1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Utilities to post-process semi-empirical parameters 8!> \par History 9!> [tlaino] 03.2008 - Splitting from semi_empirical_parameters and 10!> keeping there only the setting of the SE params 11!> \author Teodoro Laino [tlaino] - University of Zurich 12!> \date 03.2008 [tlaino] 13! ************************************************************************************************** 14MODULE semi_empirical_par_utils 15 16 USE kinds, ONLY: dp 17 USE mathconstants, ONLY: fac 18 USE mathlib, ONLY: binomial 19 USE physcon, ONLY: bohr,& 20 evolt 21 USE semi_empirical_int_arrays, ONLY: int_ij,& 22 int_kl,& 23 int_onec2el 24 USE semi_empirical_types, ONLY: get_se_param,& 25 semi_empirical_type 26#include "./base/base_uses.f90" 27 28 IMPLICIT NONE 29 30 PRIVATE 31 32 INTEGER, PARAMETER, PRIVATE :: nelem = 106 33 34 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_par_utils' 35 36! STANDARD MOPAC PARAMETERS USED FOR AM1, RM1, MNDO, PM3, PM6, 37! PM6-FM 38! 39! H He 40! Li Be B C N O F Ne 41! Na Mg Al Si P S Cl Ar 42! K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr 43! Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe 44! Cs Ba La Ce-Lu Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn 45! Fr Ra Ac Th Pa U Np Pu Am Cm Bk Cf Es Fm Md No Lr Rf Ha 106 46 47! "s" shell 48 INTEGER, DIMENSION(0:nelem), PRIVATE :: Nos = (/-1, & ! 0 49 1, 2, & ! 2 50 1, 2, 2, 2, 2, 2, 2, 0, & ! 10 51 1, 2, 2, 2, 2, 2, 2, 0, & ! 18 52 1, 2, 2, 2, 2, 1, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 0, & ! 36 53 1, 2, 2, 2, 1, 1, 2, 1, 1, 0, 1, 2, 2, 2, 2, 2, 2, 0, & ! 54 54 1, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, & 55 2, 2, 1, 2, 2, 2, 1, 1, 2, 2, 2, 2, 2, 2, 0, & ! 86 56 1, 1, 2, 4, 2, 2, 2, 2, 2, 2, 2, 2, 0, 3, -3, 1, 2, 1, -2, -1/) 57 58! "p" shell 59 INTEGER, DIMENSION(0:nelem), PRIVATE :: Nop = (/-1, & ! 0 60 0, 0, & ! 2 61 0, 0, 1, 2, 3, 4, 5, 6, & ! 10 62 0, 0, 1, 2, 3, 4, 5, 6, & ! 18 63 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 4, 5, 6, & ! 36 64 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 4, 5, 6, & ! 54 65 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 66 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 4, 5, 6, & ! 86 67 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/) 68 69! "d" shell 70 INTEGER, DIMENSION(0:nelem), PRIVATE :: Nod = (/-1, & ! 0 71 0, 0, & ! 2 72 0, 0, 0, 0, 0, 0, 0, 0, & ! 10 73 0, 0, 0, 0, 0, 0, 0, 0, & ! 18 74 0, 0, 1, 2, 3, 5, 5, 6, 7, 8, 10, 0, 0, 0, 0, 0, 0, 0, & ! 36 75 0, 0, 1, 2, 4, 5, 5, 7, 8, 10, 10, 0, 0, 0, 0, 0, 0, 0, & ! 54 76 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, & 77 2, 3, 5, 5, 6, 7, 9, 10, 0, 0, 0, 0, 0, 0, 0, & ! 86 78 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/) 79 80! H <Quantum Numbers for s, p, d and f orbitals> He 81! Li Be B C N O F Ne 82! Na Mg Al Si P S Cl Ar 83! K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr 84! Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe 85! Cs Ba La Ce-Lu Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn 86! Fr Ra Ac Th Pa U Np Pu Am Cm Bk Cf Es Fm Md No Lr Rf Ha 106 87 88 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: nqs = (/-1, & ! 0 89 1, 1, & ! 2 90 2, 2, 2, 2, 2, 2, 2, 2, & ! 10 91 3, 3, 3, 3, 3, 3, 3, 3, & ! 18 92 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, & ! 36 93 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & ! 54 94 6, 6, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 6, & 95 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, -1, -1, -1, & ! 86 96 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/) 97 98 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: nqp = (/-1, & ! 0 99 -1, -1, & ! 2 100 2, 2, 2, 2, 2, 2, 2, 2, & ! 10 101 3, 3, 3, 3, 3, 3, 3, 3, & ! 18 102 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, & ! 36 103 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & ! 54 104 6, 6, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 6, & 105 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, -1, -1, -1, & ! 86 106 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/) 107 108 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: nqd = (/-1, & ! 0 109 -1, -1, & ! 2 110 -1, -1, -1, -1, -1, -1, -1, -1, & ! 10 111 -1, -1, 3, 3, 3, 3, 3, -1, & ! 18 112 -1, -1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, -1, & ! 36 113 -1, -1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, -1, & ! 54 114 -1, -1, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 5, & 115 5, 5, 5, 5, 5, 5, 5, 5, 5, -1, -1, -1, -1, -1, -1, & ! 86 116 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/) 117 118 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: nqf = (/-1, & ! 0 119 -1, -1, & ! 2 120 -1, -1, -1, -1, -1, -1, -1, -1, & ! 10 121 -1, -1, -1, -1, -1, -1, -1, -1, & ! 18 122 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, & ! 36 123 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, & ! 54 124 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, & 125 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, & ! 86 126 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/) 127 128 ! Element Valence 129 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: zval = (/-1, & ! 0 130 1, 2, & ! 2 131 1, 2, 3, 4, 5, 6, 7, 8, & ! 10 132 1, 2, 3, 4, 5, 6, 7, 8, & ! 18 133 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 2, 3, 4, 5, 6, 7, 8, & ! 36 134 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 2, 3, 4, 5, 6, 7, 8, & ! 54 135 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 3, & 136 4, 5, 6, 7, 8, 9, 10, 11, 2, 3, 4, 5, 6, 7, -1, & ! 86 137 -1, -1, -1, 4, -1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/) 138 139 ! Number of 1 center 2 electron integrals involving partially filled d shells 140 ! r016: <SS|DD> 141 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: ir016 = (/0, & ! 0 142 0, 0, & ! 2 143 0, 0, 0, 0, 0, 0, 0, 0, & ! 10 144 0, 0, 0, 0, 0, 0, 0, 0, & ! 18 145 0, 0, 2, 4, 6, 5, 10, 12, 14, 16, 10, 0, 0, 0, 0, 0, 0, 0, & ! 36 146 0, 0, 4, 4, 4, 5, 10, 7, 8, 0, 10, 0, 0, 0, 0, 0, 0, 0, & ! 54 147 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 148 4, 6, 8, 10, 12, 14, 9, 10, 0, 0, 0, 0, 0, 0, 0, & ! 86 149 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/) 150 151 ! r066: <DD|DD> "0" term 152 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: ir066 = (/0, & ! 0 153 0, 0, & ! 2 154 0, 0, 0, 0, 0, 0, 0, 0, & ! 10 155 0, 0, 0, 0, 0, 0, 0, 0, & ! 18 156 0, 0, 0, 1, 3, 10, 10, 15, 21, 28, 45, 0, 0, 0, 0, 0, 0, 0, & ! 36 157 0, 0, 0, 1, 6, 10, 10, 21, 28, 45, 45, 0, 0, 0, 0, 0, 0, 0, & ! 54 158 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 159 1, 3, 6, 10, 15, 21, 36, 45, 0, 0, 0, 0, 0, 0, 0, & ! 86 160 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/) 161 162 ! r244: <SD|SD> 163 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: ir244 = (/0, & ! 0 164 0, 0, & ! 2 165 0, 0, 0, 0, 0, 0, 0, 0, & ! 10 166 0, 0, 0, 0, 0, 0, 0, 0, & ! 18 167 0, 0, 1, 2, 3, 5, 5, 6, 7, 8, 5, 0, 0, 0, 0, 0, 0, 0, & ! 36 168 0, 0, 1, 2, 4, 5, 5, 5, 5, 0, 5, 0, 0, 0, 0, 0, 0, 0, & ! 54 169 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 170 2, 3, 4, 5, 6, 7, 5, 5, 0, 0, 0, 0, 0, 0, 0, & ! 86 171 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/) 172 173 ! r266: <DD|DD> "2" term 174 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: ir266 = (/0, & ! 0 175 0, 0, & ! 2 176 0, 0, 0, 0, 0, 0, 0, 0, & ! 10 177 0, 0, 0, 0, 0, 0, 0, 0, & ! 18 178 0, 0, 0, 8, 15, 35, 35, 35, 43, 50, 70, 0, 0, 0, 0, 0, 0, 0, & ! 36 179 0, 0, 0, 8, 21, 35, 35, 43, 50, 70, 70, 0, 0, 0, 0, 0, 0, 0, & ! 54 180 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 181 8, 15, 21, 35, 35, 43, 56, 70, 0, 0, 0, 0, 0, 0, 0, & ! 86 182 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/) 183 184 ! r466: <DD|DD> "4" term 185 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: ir466 = (/0, & ! 0 186 0, 0, & ! 2 187 0, 0, 0, 0, 0, 0, 0, 0, & ! 10 188 0, 0, 0, 0, 0, 0, 0, 0, & ! 18 189 0, 0, 0, 1, 8, 35, 35, 35, 36, 43, 70, 0, 0, 0, 0, 0, 0, 0, & ! 36 190 0, 0, 0, 1, 21, 35, 35, 36, 43, 70, 70, 0, 0, 0, 0, 0, 0, 0, & ! 54 191 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 192 1, 8, 21, 35, 35, 36, 56, 70, 0, 0, 0, 0, 0, 0, 0, & ! 86 193 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/) 194 195 INTERFACE amn_l 196 MODULE PROCEDURE amn_l1, amn_l2 197 END INTERFACE 198 199 PUBLIC :: convert_param_to_cp2k, calpar, valence_electrons, get_se_basis, & 200 setup_1c_2el_int, amn_l 201 202CONTAINS 203 204! ************************************************************************************************** 205!> \brief Gives back the number of valence electrons for element z and also the 206!> number of atomic orbitals for that specific element 207!> \param sep ... 208!> \param extended_basis_set ... 209! ************************************************************************************************** 210 SUBROUTINE valence_electrons(sep, extended_basis_set) 211 TYPE(semi_empirical_type), POINTER :: sep 212 LOGICAL, INTENT(IN) :: extended_basis_set 213 214 CHARACTER(len=*), PARAMETER :: routineN = 'valence_electrons', & 215 routineP = moduleN//':'//routineN 216 217 INTEGER :: natorb, z 218 LOGICAL :: check, use_p_orbitals 219 REAL(KIND=dp) :: zeff 220 221 use_p_orbitals = .TRUE. 222 z = sep%z 223 CPASSERT(z >= 0) 224 ! Special case for Hydrogen.. If requested allow p-orbitals on it.. 225 SELECT CASE (z) 226 CASE (0, 2) 227 use_p_orbitals = .FALSE. 228 CASE (1) 229 use_p_orbitals = sep%p_orbitals_on_h 230 CASE DEFAULT 231 ! Nothing to do.. 232 END SELECT 233 ! Determine the number of atomic orbitals 234 natorb = 0 235 IF (nqs(z) > 0) natorb = natorb + 1 236 IF ((nqp(z) > 0) .OR. use_p_orbitals) natorb = natorb + 3 237 IF (extended_basis_set .AND. element_has_d(sep)) natorb = natorb + 5 238 IF (extended_basis_set .AND. element_has_f(sep)) natorb = natorb + 7 239 ! Check and assignemnt 240 check = (natorb <= 4) .OR. (extended_basis_set) 241 CPASSERT(check) 242 sep%natorb = natorb 243 sep%extended_basis_set = extended_basis_set 244 ! Determine the Z eff 245 zeff = REAL(zval(z), KIND=dp) 246 sep%zeff = zeff 247 END SUBROUTINE valence_electrons 248 249! ************************************************************************************************** 250!> \brief Gives back the number of basis function for each l 251!> \param sep ... 252!> \param l ... 253!> \return ... 254! ************************************************************************************************** 255 FUNCTION get_se_basis(sep, l) RESULT(n) 256 TYPE(semi_empirical_type), POINTER :: sep 257 INTEGER, INTENT(IN) :: l 258 INTEGER :: n 259 260 CHARACTER(len=*), PARAMETER :: routineN = 'get_se_basis', routineP = moduleN//':'//routineN 261 262 IF (sep%z < 0 .OR. sep%z > nelem) THEN 263 CPABORT("Invalid atomic number !") 264 ELSE 265 IF (l == 0) THEN 266 n = nqs(sep%z) 267 ELSEIF (l == 1) THEN 268 ! Special case for Hydrogen.. If requested allow p-orbitals on it.. 269 IF ((sep%z == 1) .AND. sep%p_orbitals_on_h) THEN 270 n = 1 271 ELSE 272 n = nqp(sep%z) 273 END IF 274 ELSEIF (l == 2) THEN 275 n = nqd(sep%z) 276 ELSEIF (l == 3) THEN 277 n = nqf(sep%z) 278 ELSE 279 CPABORT("Invalid l quantum number !") 280 END IF 281 IF (n < 0) THEN 282 CPABORT("Invalid n quantum number !") 283 END IF 284 END IF 285 END FUNCTION get_se_basis 286 287! ************************************************************************************************** 288!> \brief Converts parameter units to internal 289!> \param sep ... 290!> \date 03.2008 [tlaino] 291!> \author Teodoro Laino [tlaino] - University of Zurich 292! ************************************************************************************************** 293 SUBROUTINE convert_param_to_cp2k(sep) 294 TYPE(semi_empirical_type), POINTER :: sep 295 296 CHARACTER(len=*), PARAMETER :: routineN = 'convert_param_to_cp2k', & 297 routineP = moduleN//':'//routineN 298 299 sep%beta = sep%beta/evolt 300 sep%uss = sep%uss/evolt 301 sep%upp = sep%upp/evolt 302 sep%udd = sep%udd/evolt 303 sep%alp = sep%alp/bohr 304 sep%eisol = sep%eisol/evolt 305 sep%gss = sep%gss/evolt 306 sep%gsp = sep%gsp/evolt 307 sep%gpp = sep%gpp/evolt 308 sep%gp2 = sep%gp2/evolt 309 sep%gsd = sep%gsd/evolt 310 sep%gpd = sep%gpd/evolt 311 sep%gdd = sep%gdd/evolt 312 sep%hsp = sep%hsp/evolt 313 sep%fn1 = sep%fn1*bohr/evolt 314 sep%fn2 = sep%fn2/bohr/bohr 315 sep%fn3 = sep%fn3*bohr 316 sep%bfn1 = sep%bfn1*bohr/evolt 317 sep%bfn2 = sep%bfn2/bohr/bohr 318 sep%bfn3 = sep%bfn3*bohr 319 sep%f0sd = sep%f0sd 320 sep%g2sd = sep%g2sd 321 sep%a = sep%a*bohr/evolt 322 sep%b = sep%b/bohr/bohr 323 sep%c = sep%c*bohr 324 sep%pre = sep%pre/evolt 325 sep%d = sep%d/bohr 326 327 END SUBROUTINE convert_param_to_cp2k 328 329! ************************************************************************************************** 330!> \brief Calculates missing parameters 331!> \param z ... 332!> \param sep ... 333!> \date 03.2008 [tlaino] 334!> \author Teodoro Laino [tlaino] - University of Zurich 335! ************************************************************************************************** 336 SUBROUTINE calpar(z, sep) 337 INTEGER :: z 338 TYPE(semi_empirical_type), POINTER :: sep 339 340 CHARACTER(len=*), PARAMETER :: routineN = 'calpar', routineP = moduleN//':'//routineN 341 342 INTEGER :: iod, iop, ios, j, jmax, k, l 343 REAL(KIND=dp) :: ad, am, aq, d1, d2, d3, dd, df, eisol, gdd1, gp2, gp2c, gpp, gppc, gqq, & 344 gsp, gspc, gss, gssc, hpp, hpp1, hpp2, hsp, hsp1, hsp2, hspc, p, p4, q1, q2, q3, qf, qn, & 345 qq, udd, upp, uss, zp, zs 346 347 IF (.NOT. sep%defined) RETURN 348 uss = sep%uss 349 upp = sep%upp 350 udd = sep%udd 351 gss = sep%gss 352 gpp = sep%gpp 353 gsp = sep%gsp 354 gp2 = sep%gp2 355 hsp = sep%hsp 356 zs = sep%sto_exponents(0) 357 zp = sep%sto_exponents(1) 358 ios = Nos(z) 359 iop = Nop(z) 360 iod = Nod(z) 361 362 p = 2.0_dp 363 p4 = p**4 364 ! GSSC is the number of two-electron terms of type <SS|SS> 365 gssc = REAL(MAX(ios - 1, 0), KIND=dp) 366 k = iop 367 ! GSPC is the number of two-electron terms of type <SS|PP> 368 gspc = REAL(ios*k, KIND=dp) 369 l = MIN(k, 6 - k) 370 ! GP2C is the number of two-electron terms of type <PP|PP> 371 ! plus 0.5 of the number of HPP integrals. 372 ! (HPP is not used; instead it is replaced by 0.5(GPP-GP2)) 373 gp2c = REAL((k*(k - 1))/2, KIND=dp) + 0.5_dp*REAL((l*(l - 1))/2, KIND=dp) 374 ! GPPC is minus 0.5 times the number of HPP integrals. 375 gppc = -0.5_dp*REAL((l*(l - 1))/2, KIND=dp) 376 ! HSPC is the number of two-electron terms of type <SP|SP>. 377 ! (S and P must have the same spin. In all cases, if 378 ! P is non-zero, there are two S electrons) 379 hspc = REAL(-k, KIND=dp) 380 ! Constraint the value of the STO exponent 381 zp = MAX(0.3_dp, zp) 382 ! Take into account constraints on the values of the integrals 383 hpp = 0.5_dp*(gpp - gp2) 384 hpp = MAX(0.1_dp, hpp) 385 hsp = MAX(1.E-7_dp, hsp) 386 387 ! Evaluation of EISOL 388 eisol = uss*ios + upp*iop + udd*iod + gss*gssc + gpp*gppc + gsp*gspc + gp2*gp2c + hsp*hspc 389 390 ! Principal quantum number 391 qn = REAL(nqs(z), KIND=dp) 392 CPASSERT(qn > 0) 393 394 ! Charge separation evaluation 395 dd = (2.0_dp*qn + 1)*(4.0_dp*zs*zp)**(qn + 0.5_dp)/(zs + zp)**(2.0_dp*qn + 2)/SQRT(3.0_dp) 396 qq = SQRT((4.0_dp*qn*qn + 6.0_dp*qn + 2.0_dp)/20.0_dp)/zp 397 398 ! Calculation of the additive terms in atomic units 399 jmax = 5 400 gdd1 = (hsp/(evolt*dd**2))**(1.0_dp/3.0_dp) 401 d1 = gdd1 402 d2 = gdd1 + 0.04_dp 403 DO j = 1, jmax 404 df = d2 - d1 405 hsp1 = 0.5_dp*d1 - 0.5_dp/SQRT(4.0_dp*dd**2 + 1.0_dp/d1**2) 406 hsp2 = 0.5_dp*d2 - 0.5_dp/SQRT(4.0_dp*dd**2 + 1.0_dp/d2**2) 407 IF (ABS(hsp2 - hsp1) < EPSILON(0.0_dp)) EXIT 408 d3 = d1 + df*(hsp/evolt - hsp1)/(hsp2 - hsp1) 409 d1 = d2 410 d2 = d3 411 END DO 412 gqq = (p4*hpp/(evolt*48.0_dp*qq**4))**0.2_dp 413 q1 = gqq 414 q2 = gqq + 0.04_dp 415 DO j = 1, jmax 416 qf = q2 - q1 417 hpp1 = 0.25_dp*q1 - 0.5_dp/SQRT(4.0_dp*qq**2 + 1.0_dp/q1**2) + 0.25_dp/SQRT(8.0_dp*qq**2 + 1.0_dp/q1**2) 418 hpp2 = 0.25_dp*q2 - 0.5_dp/SQRT(4.0_dp*qq**2 + 1.0_dp/q2**2) + 0.25_dp/SQRT(8.0_dp*qq**2 + 1.0_dp/q2**2) 419 IF (ABS(hpp2 - hpp1) < EPSILON(0.0_dp)) EXIT 420 q3 = q1 + qf*(hpp/evolt - hpp1)/(hpp2 - hpp1) 421 q1 = q2 422 q2 = q3 423 END DO 424 am = gss/evolt 425 ad = d2 426 aq = q2 427 IF (z == 1) THEN 428 ad = am 429 aq = am 430 dd = 0.0_dp 431 qq = 0.0_dp 432 END IF 433 ! Overwrite these parameters if they were undefined.. otherwise keep the defined 434 ! value 435 IF (ABS(sep%eisol) < EPSILON(0.0_dp)) sep%eisol = eisol 436 IF (ABS(sep%dd) < EPSILON(0.0_dp)) sep%dd = dd 437 IF (ABS(sep%qq) < EPSILON(0.0_dp)) sep%qq = qq 438 IF (ABS(sep%am) < EPSILON(0.0_dp)) sep%am = am 439 IF (ABS(sep%ad) < EPSILON(0.0_dp)) sep%ad = ad 440 IF (ABS(sep%aq) < EPSILON(0.0_dp)) sep%aq = aq 441 ! Proceed with d-orbitals and fill the Kolpman-Ohno and Charge Separation 442 ! arrays 443 CALL calpar_d(sep) 444 END SUBROUTINE calpar 445 446! ************************************************************************************************** 447!> \brief Finalize the initialization of parameters, defining additional 448!> parameters for d-orbitals 449!> 450!> \param sep ... 451!> \date 03.2008 [tlaino] 452!> \author Teodoro Laino [tlaino] - University of Zurich 453! ************************************************************************************************** 454 SUBROUTINE calpar_d(sep) 455 TYPE(semi_empirical_type), POINTER :: sep 456 457 CHARACTER(len=*), PARAMETER :: routineN = 'calpar_d', routineP = moduleN//':'//routineN 458 459 REAL(KIND=dp), DIMENSION(6) :: amn 460 461! Determine if this element owns d-orbitals (only if the parametrization 462! supports the d-orbitals) 463 464 IF (sep%extended_basis_set) sep%dorb = element_has_d(sep) 465 IF (sep%dorb) THEN 466 CALL amn_l(sep, amn) 467 CALL eval_1c_2el_spd(sep) 468 CALL eval_cs_ko(sep, amn) 469 END IF 470 IF (.NOT. sep%dorb) THEN 471 ! Use the old integral module 472 IF (ABS(sep%am) > EPSILON(0.0_dp)) THEN 473 sep%ko(1) = 0.5_dp/sep%am 474 END IF 475 IF (ABS(sep%ad) > EPSILON(0.0_dp) .AND. (sep%z /= 1)) THEN 476 sep%ko(2) = 0.5_dp/sep%ad 477 END IF 478 IF (ABS(sep%aq) > EPSILON(0.0_dp) .AND. (sep%z /= 1)) THEN 479 sep%ko(3) = 0.5_dp/sep%aq 480 END IF 481 sep%ko(7) = sep%ko(1) 482 sep%ko(9) = sep%ko(1) 483 sep%cs(2) = sep%dd 484 sep%cs(3) = sep%qq*SQRT(2.0_dp) 485 ELSE 486 ! Use the new integral module 487 sep%ko(9) = sep%ko(1) 488 sep%aq = 0.5_dp/sep%ko(3) 489 END IF 490 ! In case the Klopman-Ohno CORE therm is provided let's overwrite the 491 ! computed one 492 IF (ABS(sep%rho) > EPSILON(0.0_dp)) THEN 493 sep%ko(9) = sep%rho 494 END IF 495 END SUBROUTINE calpar_d 496 497! ************************************************************************************************** 498!> \brief Determines if the elements has d-orbitals 499!> 500!> \param sep ... 501!> \return ... 502!> \date 05.2008 [tlaino] 503!> \author Teodoro Laino [tlaino] - University of Zurich 504! ************************************************************************************************** 505 FUNCTION element_has_d(sep) RESULT(res) 506 TYPE(semi_empirical_type), POINTER :: sep 507 LOGICAL :: res 508 509 CHARACTER(len=*), PARAMETER :: routineN = 'element_has_d', routineP = moduleN//':'//routineN 510 511 res = (nqd(sep%z) > 0 .AND. sep%sto_exponents(2) > EPSILON(0.0_dp)) 512 END FUNCTION element_has_d 513 514! ************************************************************************************************** 515!> \brief Determines if the elements has f-orbitals 516!> 517!> \param sep ... 518!> \return ... 519!> \date 05.2008 [tlaino] 520!> \author Teodoro Laino [tlaino] - University of Zurich 521! ************************************************************************************************** 522 FUNCTION element_has_f(sep) RESULT(res) 523 TYPE(semi_empirical_type), POINTER :: sep 524 LOGICAL :: res 525 526 CHARACTER(len=*), PARAMETER :: routineN = 'element_has_f', routineP = moduleN//':'//routineN 527 528 res = (nqf(sep%z) > 0 .AND. sep%sto_exponents(3) > EPSILON(0.0_dp)) 529 END FUNCTION element_has_f 530 531! ************************************************************************************************** 532!> \brief Computes the A^{\mu \nu}_l values for the evaluation of the two-center 533!> two-electron integrals. The term is the one reported in Eq.(7) of TCA 534!> 535!> \param sep ... 536!> \param amn ... 537!> \date 03.2008 [tlaino] 538!> \par Notation Index: 1 (SS), 2 (SP), 3 (SD), 4 (PP), 5 (PD), 6 (DD) 539!> 540!> \author Teodoro Laino [tlaino] - University of Zurich 541! ************************************************************************************************** 542 SUBROUTINE amn_l1(sep, amn) 543 TYPE(semi_empirical_type), POINTER :: sep 544 REAL(KIND=dp), DIMENSION(6), INTENT(OUT) :: amn 545 546 CHARACTER(len=*), PARAMETER :: routineN = 'amn_l1', routineP = moduleN//':'//routineN 547 548 INTEGER :: nd, nsp 549 REAL(KIND=dp) :: z1, z2, z3 550 551 z1 = sep%sto_exponents(0) 552 z2 = sep%sto_exponents(1) 553 z3 = sep%sto_exponents(2) 554 IF (z1 <= 0.0_dp) & 555 CALL cp_abort(__LOCATION__, & 556 "Trying to use s-orbitals, but the STO exponents is set to 0. "// & 557 "Please check if your parameterization supports the usage of s orbitals! ") 558 amn = 0.0_dp 559 nsp = nqs(sep%z) 560 IF (sep%natorb >= 4) THEN 561 IF (z2 <= 0.0_dp) & 562 CALL cp_abort(__LOCATION__, & 563 "Trying to use p-orbitals, but the STO exponents is set to 0. "// & 564 "Please check if your parameterization supports the usage of p orbitals! ") 565 amn(2) = amn_l_low(z1, z2, nsp, nsp, 1) 566 amn(3) = amn_l_low(z2, z2, nsp, nsp, 2) 567 IF (sep%dorb) THEN 568 IF (z3 <= 0.0_dp) & 569 CALL cp_abort(__LOCATION__, & 570 "Trying to use d-orbitals, but the STO exponents is set to 0. "// & 571 "Please check if your parameterization supports the usage of d orbitals! ") 572 nd = nqd(sep%z) 573 amn(4) = amn_l_low(z1, z3, nsp, nd, 2) 574 amn(5) = amn_l_low(z2, z3, nsp, nd, 1) 575 amn(6) = amn_l_low(z3, z3, nd, nd, 2) 576 END IF 577 END IF 578 END SUBROUTINE amn_l1 579 580! ************************************************************************************************** 581!> \brief Computes the A^{\mu \nu}_l values for the evaluation of the two-center 582!> two-electron integrals. The term is the one reported in Eq.(7) of TCA 583!> 584!> \param sep ... 585!> \param amn ... 586!> \date 09.2008 [tlaino] 587!> \par 588!> 589!> \author Teodoro Laino [tlaino] - University of Zurich 590! ************************************************************************************************** 591 SUBROUTINE amn_l2(sep, amn) 592 TYPE(semi_empirical_type), POINTER :: sep 593 REAL(KIND=dp), DIMENSION(6, 0:2), INTENT(OUT) :: amn 594 595 CHARACTER(len=*), PARAMETER :: routineN = 'amn_l2', routineP = moduleN//':'//routineN 596 597 INTEGER :: nd, nsp 598 REAL(KIND=dp) :: z1, z2, z3 599 600 z1 = sep%sto_exponents(0) 601 z2 = sep%sto_exponents(1) 602 z3 = sep%sto_exponents(2) 603 CPASSERT(z1 > 0.0_dp) 604 amn = 0.0_dp 605 nsp = nqs(sep%z) 606 amn(1, 0) = amn_l_low(z1, z1, nsp, nsp, 0) 607 IF (sep%natorb >= 4) THEN 608 CPASSERT(z2 > 0.0_dp) 609 amn(2, 1) = amn_l_low(z1, z2, nsp, nsp, 1) 610 amn(3, 0) = amn_l_low(z2, z2, nsp, nsp, 0) 611 amn(3, 2) = amn_l_low(z2, z2, nsp, nsp, 2) 612 IF (sep%dorb) THEN 613 CPASSERT(z3 > 0.0_dp) 614 nd = nqd(sep%z) 615 amn(4, 2) = amn_l_low(z1, z3, nsp, nd, 2) 616 amn(5, 1) = amn_l_low(z2, z3, nsp, nd, 1) 617 amn(6, 0) = amn_l_low(z3, z3, nd, nd, 0) 618 amn(6, 2) = amn_l_low(z3, z3, nd, nd, 2) 619 END IF 620 END IF 621 END SUBROUTINE amn_l2 622 623! ************************************************************************************************** 624!> \brief Low level for computing Eq.(7) of TCA 625!> \param z1 ... 626!> \param z2 ... 627!> \param n1 ... 628!> \param n2 ... 629!> \param l ... 630!> \return ... 631!> \date 03.2008 [tlaino] 632!> \author Teodoro Laino [tlaino] - University of Zurich 633! ************************************************************************************************** 634 FUNCTION amn_l_low(z1, z2, n1, n2, l) RESULT(amnl) 635 REAL(KIND=dp), INTENT(IN) :: z1, z2 636 INTEGER, INTENT(IN) :: n1, n2, l 637 REAL(KIND=dp) :: amnl 638 639 amnl = fac(n1 + n2 + l)/SQRT(fac(2*n1)*fac(2*n2))*(2.0_dp*z1/(z1 + z2))**n1* & 640 (2.0_dp*z2/(z1 + z2))**n2*2.0_dp*SQRT(z1*z2)/(z1 + z2)**(l + 1) 641 642 END FUNCTION amn_l_low 643 644! ************************************************************************************************** 645!> \brief Calculation of chare separations and additive terms used for computing 646!> the two-center two-electron integrals with d-orbitals 647!> \param sep ... 648!> \param amn ... 649!> \date 03.2008 [tlaino] 650!> \par Notation 651!> -) Charge separations [sep%cs(1:6)] [see equations (12)-(16) of TCA] 652!> -) Additive terms of Klopman-Ohno terms [sep%ko(1:9)] [see equations 653!> (19)-(26) of TCA] 654!> -) Atomic core additive term stored in sep%ko(9): used in the calculation 655!> of the core-electron attractions and core-core repulsions 656!> \author Teodoro Laino [tlaino] - University of Zurich 657! ************************************************************************************************** 658 SUBROUTINE eval_cs_ko(sep, amn) 659 TYPE(semi_empirical_type), POINTER :: sep 660 REAL(KIND=dp), DIMENSION(6), INTENT(IN) :: amn 661 662 CHARACTER(len=*), PARAMETER :: routineN = 'eval_cs_ko', routineP = moduleN//':'//routineN 663 664 REAL(KIND=dp) :: d, fg 665 666! SS term 667 668 fg = sep%gss 669 sep%ko(1) = ko_ij(0, 1.0_dp, fg) 670 IF (sep%natorb >= 4) THEN 671 ! Other terms for SP basis 672 ! SP 673 d = amn(2)/SQRT(3.0_dp) 674 fg = sep%hsp 675 sep%cs(2) = d 676 sep%ko(2) = ko_ij(1, d, fg) 677 ! PP 678 sep%ko(7) = sep%ko(1) 679 d = SQRT(amn(3)*2.0_dp/5.0_dp) 680 fg = 0.5_dp*(sep%gpp - sep%gp2) 681 sep%cs(3) = d 682 sep%ko(3) = ko_ij(2, d, fg) 683 ! Terms involving d-orbitals 684 IF (sep%dorb) THEN 685 ! SD 686 d = SQRT(amn(4)*2.0_dp/SQRT(15.0_dp)) 687 fg = sep%onec2el(19) 688 sep%cs(4) = d 689 sep%ko(4) = ko_ij(2, d, fg) 690 ! PD 691 d = amn(5)/SQRT(5.0_dp) 692 fg = sep%onec2el(23) - 1.8_dp*sep%onec2el(35) 693 sep%cs(5) = d 694 sep%ko(5) = ko_ij(1, d, fg) 695 ! DD 696 fg = 0.2_dp*(sep%onec2el(29) + 2.0_dp*sep%onec2el(30) + 2.0_dp*sep%onec2el(31)) 697 sep%ko(8) = ko_ij(0, 1.0_dp, fg) 698 d = SQRT(amn(6)*2.0_dp/7.0_dp) 699 fg = sep%onec2el(44) - (20.0_dp/35.0_dp)*sep%onec2el(52) 700 sep%cs(6) = d 701 sep%ko(6) = ko_ij(2, d, fg) 702 END IF 703 END IF 704 END SUBROUTINE eval_cs_ko 705 706! ************************************************************************************************** 707!> \brief Computes the 1 center two-electrons integrals for a SPD basis 708!> \param sep ... 709!> \date 03.2008 [tlaino] 710!> \author Teodoro Laino [tlaino] - University of Zurich 711! ************************************************************************************************** 712 SUBROUTINE eval_1c_2el_spd(sep) 713 TYPE(semi_empirical_type), POINTER :: sep 714 715 CHARACTER(len=*), PARAMETER :: routineN = 'eval_1c_2el_spd', & 716 routineP = moduleN//':'//routineN 717 718 REAL(KIND=dp) :: r016, r036, r066, r125, r155, r234, & 719 r236, r244, r246, r266, r355, r466, & 720 s15, s3, s5 721 722 IF (sep%dorb) THEN 723 s3 = SQRT(3.0_dp) 724 s5 = SQRT(5.0_dp) 725 s15 = SQRT(15.0_dp) 726 727 ! We evaluate now the Slater-Condon parameters (Rlij) 728 CALL sc_param(sep, r066, r266, r466, r016, r244, r036, r236, r155, r355, r125, & 729 r234, r246) 730 731 IF (ABS(sep%f0sd) > EPSILON(0.0_dp)) THEN 732 r016 = sep%f0sd 733 END IF 734 IF (ABS(sep%g2sd) > EPSILON(0.0_dp)) THEN 735 r244 = sep%g2sd 736 END IF 737 CALL eisol_corr(sep, r016, r066, r244, r266, r466) 738 sep%onec2el(1) = r016 739 sep%onec2el(2) = 2.0_dp/(3.0_dp*s5)*r125 740 sep%onec2el(3) = 1.0_dp/s15*r125 741 sep%onec2el(4) = 2.0_dp/(5.0_dp*s5)*r234 742 sep%onec2el(5) = r036 + 4.0_dp/35.0_dp*r236 743 sep%onec2el(6) = r036 + 2.0_dp/35.0_dp*r236 744 sep%onec2el(7) = r036 - 4.0_dp/35.0_dp*r236 745 sep%onec2el(8) = -1.0_dp/(3.0_dp*s5)*r125 746 sep%onec2el(9) = SQRT(3.0_dp/125.0_dp)*r234 747 sep%onec2el(10) = s3/35.0_dp*r236 748 sep%onec2el(11) = 3.0_dp/35.0_dp*r236 749 sep%onec2el(12) = -1.0_dp/(5.0_dp*s5)*r234 750 sep%onec2el(13) = r036 - 2.0_dp/35.0_dp*r236 751 sep%onec2el(14) = -2.0_dp*s3/35.0_dp*r236 752 sep%onec2el(15) = -sep%onec2el(3) 753 sep%onec2el(16) = -sep%onec2el(11) 754 sep%onec2el(17) = -sep%onec2el(9) 755 sep%onec2el(18) = -sep%onec2el(14) 756 sep%onec2el(19) = 1.0_dp/5.0_dp*r244 757 sep%onec2el(20) = 2.0_dp/(7.0_dp*s5)*r246 758 sep%onec2el(21) = sep%onec2el(20)/2.0_dp 759 sep%onec2el(22) = -sep%onec2el(20) 760 sep%onec2el(23) = 4.0_dp/15.0_dp*r155 + 27.0_dp/245.0_dp*r355 761 sep%onec2el(24) = 2.0_dp*s3/15.0_dp*r155 - 9.0_dp*s3/245.0_dp*r355 762 sep%onec2el(25) = 1.0_dp/15.0_dp*r155 + 18.0_dp/245.0_dp*r355 763 sep%onec2el(26) = -s3/15.0_dp*r155 + 12.0_dp*s3/245.0_dp*r355 764 sep%onec2el(27) = -s3/15.0_dp*r155 - 3.0_dp*s3/245.0_dp*r355 765 sep%onec2el(28) = -sep%onec2el(27) 766 sep%onec2el(29) = r066 + 4.0_dp/49.0_dp*r266 + 4.0_dp/49.0_dp*r466 767 sep%onec2el(30) = r066 + 2.0_dp/49.0_dp*r266 - 24.0_dp/441.0_dp*r466 768 sep%onec2el(31) = r066 - 4.0_dp/49.0_dp*r266 + 6.0_dp/441.0_dp*r466 769 sep%onec2el(32) = SQRT(3.0_dp/245.0_dp)*r246 770 sep%onec2el(33) = 1.0_dp/5.0_dp*r155 + 24.0_dp/245.0_dp*r355 771 sep%onec2el(34) = 1.0_dp/5.0_dp*r155 - 6.0_dp/245.0_dp*r355 772 sep%onec2el(35) = 3.0_dp/49.0_dp*r355 773 sep%onec2el(36) = 1.0_dp/49.0_dp*r266 + 30.0_dp/441.0_dp*r466 774 sep%onec2el(37) = s3/49.0_dp*r266 - 5.0_dp*s3/441.0_dp*r466 775 sep%onec2el(38) = r066 - 2.0_dp/49.0_dp*r266 - 4.0_dp/441.0_dp*r466 776 sep%onec2el(39) = -2.0_dp*s3/49.0_dp*r266 + 10.0_dp*s3/441.0_dp*r466 777 sep%onec2el(40) = -sep%onec2el(32) 778 sep%onec2el(41) = -sep%onec2el(34) 779 sep%onec2el(42) = -sep%onec2el(35) 780 sep%onec2el(43) = -sep%onec2el(37) 781 sep%onec2el(44) = 3.0_dp/49.0_dp*r266 + 20.0_dp/441.0_dp*r466 782 sep%onec2el(45) = -sep%onec2el(39) 783 sep%onec2el(46) = 1.0_dp/5.0_dp*r155 - 3.0_dp/35.0_dp*r355 784 sep%onec2el(47) = -sep%onec2el(46) 785 sep%onec2el(48) = 4.0_dp/49.0_dp*r266 + 15.0_dp/441.0_dp*r466 786 sep%onec2el(49) = 3.0_dp/49.0_dp*r266 - 5.0_dp/147.0_dp*r466 787 sep%onec2el(50) = -sep%onec2el(49) 788 sep%onec2el(51) = r066 + 4.0_dp/49.0_dp*r266 - 34.0_dp/441.0_dp*r466 789 sep%onec2el(52) = 35.0_dp/441.0_dp*r466 790 sep%f0dd = r066 791 sep%f2dd = r266 792 sep%f4dd = r466 793 sep%f0sd = r016 794 sep%g2sd = r244 795 sep%f0pd = r036 796 sep%f2pd = r236 797 sep%g1pd = r155 798 sep%g3pd = r355 799 END IF 800 END SUBROUTINE eval_1c_2el_spd 801 802! ************************************************************************************************** 803!> \brief Slater-Condon parameters for 1 center 2 electrons integrals 804!> \param sep ... 805!> \param r066 ... 806!> \param r266 ... 807!> \param r466 ... 808!> \param r016 ... 809!> \param r244 ... 810!> \param r036 ... 811!> \param r236 ... 812!> \param r155 ... 813!> \param r355 ... 814!> \param r125 ... 815!> \param r234 ... 816!> \param r246 ... 817!> \date 03.2008 [tlaino] 818!> \author Teodoro Laino [tlaino] - University of Zurich 819! ************************************************************************************************** 820 SUBROUTINE sc_param(sep, r066, r266, r466, r016, r244, r036, r236, r155, r355, & 821 r125, r234, r246) 822 TYPE(semi_empirical_type), POINTER :: sep 823 REAL(KIND=dp), INTENT(out) :: r066, r266, r466, r016, r244, r036, & 824 r236, r155, r355, r125, r234, r246 825 826 CHARACTER(len=*), PARAMETER :: routineN = 'sc_param', routineP = moduleN//':'//routineN 827 828 INTEGER :: nd, ns 829 REAL(KIND=dp) :: ed, ep, es 830 831 ns = nqs(sep%z) 832 nd = nqd(sep%z) 833 es = sep%zn(0) 834 ep = sep%zn(1) 835 ed = sep%zn(2) 836 r016 = sc_param_low(0, ns, es, ns, es, nd, ed, nd, ed) 837 r036 = sc_param_low(0, ns, ep, ns, ep, nd, ed, nd, ed) 838 r066 = sc_param_low(0, nd, ed, nd, ed, nd, ed, nd, ed) 839 r155 = sc_param_low(1, ns, ep, nd, ed, ns, ep, nd, ed) 840 r125 = sc_param_low(1, ns, es, ns, ep, ns, ep, nd, ed) 841 r244 = sc_param_low(2, ns, es, nd, ed, ns, es, nd, ed) 842 r236 = sc_param_low(2, ns, ep, ns, ep, nd, ed, nd, ed) 843 r266 = sc_param_low(2, nd, ed, nd, ed, nd, ed, nd, ed) 844 r234 = sc_param_low(2, ns, ep, ns, ep, ns, es, nd, ed) 845 r246 = sc_param_low(2, ns, es, nd, ed, nd, ed, nd, ed) 846 r355 = sc_param_low(3, ns, ep, nd, ed, ns, ep, nd, ed) 847 r466 = sc_param_low(4, nd, ed, nd, ed, nd, ed, nd, ed) 848 END SUBROUTINE sc_param 849 850! ************************************************************************************************** 851!> \brief Slater-Condon parameters for 1 center 2 electrons integrals - Low level 852!> \param k ... 853!> \param na ... 854!> \param ea ... 855!> \param nb ... 856!> \param eb ... 857!> \param nc ... 858!> \param ec ... 859!> \param nd ... 860!> \param ed ... 861!> \return ... 862!> \date 03.2008 [tlaino] 863!> \par Notation 864!> -) k: Type of integral 865!> -) na,na: Principle Quantum Number of AO,corresponding to electron 1 866!> -) ea,eb: Exponents of AO,corresponding to electron 1 867!> -) nb,nc: Principle Quantum Number of AO,corresponding to electron 2 868!> -) ec,ed: Exponents of AO,corresponding to electron 2 869!> \author Teodoro Laino [tlaino] - University of Zurich 870! ************************************************************************************************** 871 FUNCTION sc_param_low(k, na, ea, nb, eb, nc, ec, nd, ed) RESULT(res) 872 INTEGER, INTENT(in) :: k, na 873 REAL(KIND=dp), INTENT(in) :: ea 874 INTEGER, INTENT(in) :: nb 875 REAL(KIND=dp), INTENT(in) :: eb 876 INTEGER, INTENT(in) :: nc 877 REAL(KIND=dp), INTENT(in) :: ec 878 INTEGER, INTENT(in) :: nd 879 REAL(KIND=dp), INTENT(in) :: ed 880 REAL(KIND=dp) :: res 881 882 CHARACTER(len=*), PARAMETER :: routineN = 'sc_param_low', routineP = moduleN//':'//routineN 883 884 INTEGER :: i, m, m1, m2, n, nab, ncd 885 REAL(KIND=dp) :: a2, aab, acd, ae, aea, aeb, aec, aed, c, & 886 e, eab, ecd, ff, s0, s1, s2, s3, tmp 887 888 CPASSERT(ea > 0.0_dp) 889 CPASSERT(eb > 0.0_dp) 890 CPASSERT(ec > 0.0_dp) 891 CPASSERT(ed > 0.0_dp) 892 aea = LOG(ea) 893 aeb = LOG(eb) 894 aec = LOG(ec) 895 aed = LOG(ed) 896 nab = na + nb 897 ncd = nc + nd 898 ecd = ec + ed 899 eab = ea + eb 900 e = ecd + eab 901 n = nab + ncd 902 ae = LOG(e) 903 a2 = LOG(2.0_dp) 904 acd = LOG(ecd) 905 aab = LOG(eab) 906 ff = fac(n - 1)/SQRT(fac(2*na)*fac(2*nb)*fac(2*nc)*fac(2*nd)) 907 tmp = na*aea + nb*aeb + nc*aec + nd*aed + 0.5_dp*(aea + aeb + aec + aed) + a2*(n + 2) - ae*n 908 c = evolt*ff*EXP(tmp) 909 s0 = 1.0_dp/e 910 s1 = 0.0_dp 911 s2 = 0.0_dp 912 m = ncd - k 913 DO i = 1, m 914 s0 = s0*e/ecd 915 s1 = s1 + s0*(binomial(ncd - k - 1, i - 1) - binomial(ncd + k, i - 1))/binomial(n - 1, i - 1) 916 END DO 917 m1 = m 918 m2 = ncd + k 919 DO i = m1, m2 920 s0 = s0*e/ecd 921 s2 = s2 + s0*binomial(m2, i)/binomial(n - 1, i) 922 END DO 923 s3 = EXP(ae*n - acd*(m2 + 1) - aab*(nab - k))/binomial(n - 1, m2) 924 res = c*(s1 - s2 + s3) 925 END FUNCTION sc_param_low 926 927! ************************************************************************************************** 928!> \brief Corrects the EISOL fo the one-center terms coming from those atoms 929!> that have partially filled "d" shells 930!> \param sep ... 931!> \param r016 ... 932!> \param r066 ... 933!> \param r244 ... 934!> \param r266 ... 935!> \param r466 ... 936!> \date 03.2008 [tlaino] 937!> \par Notation 938!> r016: <SS|DD> 939!> r066: <DD|DD> "0" term 940!> r244: <SD|SD> 941!> r266: <DD|DD> "2" term 942!> r466: <DD|DD> "4" term 943!> 944!> \author Teodoro Laino [tlaino] - University of Zurich 945! ************************************************************************************************** 946 SUBROUTINE eisol_corr(sep, r016, r066, r244, r266, r466) 947 TYPE(semi_empirical_type), POINTER :: sep 948 REAL(KIND=dp), INTENT(in) :: r016, r066, r244, r266, r466 949 950 CHARACTER(len=*), PARAMETER :: routineN = 'eisol_corr', routineP = moduleN//':'//routineN 951 952 sep%eisol = sep%eisol + ir016(sep%z)*r016 + & 953 ir066(sep%z)*r066 - & 954 ir244(sep%z)*r244/5.0_dp - & 955 ir266(sep%z)*r266/49.0_dp - & 956 ir466(sep%z)*r466/49.0_dp 957 END SUBROUTINE eisol_corr 958 959! ************************************************************************************************** 960!> \brief Computes the Klopman-Ohno additive terms for 2-center 2-electron 961!> integrals requiring that the corresponding 1-center 2-electron integral 962!> is reproduced from the 2-center one for r->0 963!> \param l ... 964!> \param d ... 965!> \param fg ... 966!> \return ... 967!> \date 03.2008 [tlaino] 968!> \author Teodoro Laino [tlaino] - University of Zurich 969! ************************************************************************************************** 970 FUNCTION ko_ij(l, d, fg) RESULT(res) 971 INTEGER, INTENT(in) :: l 972 REAL(KIND=dp), INTENT(in) :: d, fg 973 REAL(KIND=dp) :: res 974 975 CHARACTER(len=*), PARAMETER :: routineN = 'ko_ij', routineP = moduleN//':'//routineN 976 INTEGER, PARAMETER :: niter = 100 977 REAL(KIND=dp), PARAMETER :: epsil = 1.0E-08_dp, g1 = 0.382_dp, & 978 g2 = 0.618_dp 979 980 INTEGER :: i 981 REAL(KIND=dp) :: a1, a2, delta, dsq, ev4, ev8, f1, f2, & 982 y1, y2 983 984 CPASSERT(fg /= 0.0_dp) 985 ! Term for SS 986 IF (l == 0) THEN 987 res = 0.5_dp*evolt/fg 988 RETURN 989 END IF 990 ! Term for Higher angular momentum 991 dsq = d*d 992 ev4 = evolt*0.25_dp 993 ev8 = evolt/8.0_dp 994 a1 = 0.1_dp 995 a2 = 5.0_dp 996 DO i = 1, niter 997 delta = a2 - a1 998 IF (delta < epsil) EXIT 999 y1 = a1 + delta*g1 1000 y2 = a1 + delta*g2 1001 IF (l == 1) THEN 1002 f1 = (ev4*(1/y1 - 1/SQRT(y1**2 + dsq)) - fg)**2 1003 f2 = (ev4*(1/y2 - 1/SQRT(y2**2 + dsq)) - fg)**2 1004 ELSE IF (l == 2) THEN 1005 f1 = (ev8*(1.0_dp/y1 - 2.0_dp/SQRT(y1**2 + dsq*0.5_dp) + 1.0_dp/SQRT(y1**2 + dsq)) - fg)**2 1006 f2 = (ev8*(1/y2 - 2.0_dp/SQRT(y2**2 + dsq*0.5_dp) + 1.0_dp/SQRT(y2**2 + dsq)) - fg)**2 1007 END IF 1008 IF (f1 < f2) THEN 1009 a2 = y2 1010 ELSE 1011 a1 = y1 1012 END IF 1013 END DO 1014 ! Convergence reached.. define additive terms 1015 IF (f1 >= f2) THEN 1016 res = a2 1017 ELSE 1018 res = a1 1019 END IF 1020 END FUNCTION ko_ij 1021 1022! ************************************************************************************************** 1023!> \brief Fills the 1 center 2 electron integrals for the construction of the 1024!> one-electron fock matrix 1025!> \param sep ... 1026!> \date 04.2008 [tlaino] 1027!> \author Teodoro Laino [tlaino] - University of Zurich 1028! ************************************************************************************************** 1029 SUBROUTINE setup_1c_2el_int(sep) 1030 TYPE(semi_empirical_type), POINTER :: sep 1031 1032 CHARACTER(len=*), PARAMETER :: routineN = 'setup_1c_2el_int', & 1033 routineP = moduleN//':'//routineN 1034 1035 INTEGER :: i, ij, ij0, ind, ip, ipx, ipy, ipz, & 1036 isize, kl, natorb 1037 LOGICAL :: defined 1038 REAL(KIND=dp) :: gp2, gpp, gsp, gss, hsp 1039 1040 CALL get_se_param(sep, defined=defined, natorb=natorb, & 1041 gss=gss, gsp=gsp, gpp=gpp, gp2=gp2, hsp=hsp) 1042 CPASSERT(defined) 1043 1044 isize = natorb*(natorb + 1)/2 1045 ALLOCATE (sep%w(isize, isize)) 1046 ! Initialize array 1047 sep%w = 0.0_dp 1048 ! Fill the array 1049 IF (natorb > 0) THEN 1050 ip = 1 1051 sep%w(ip, ip) = gss 1052 IF (natorb > 2) THEN 1053 ipx = ip + 2 1054 ipy = ip + 5 1055 ipz = ip + 9 1056 sep%w(ipx, ip) = gsp 1057 sep%w(ipy, ip) = gsp 1058 sep%w(ipz, ip) = gsp 1059 sep%w(ip, ipx) = gsp 1060 sep%w(ip, ipy) = gsp 1061 sep%w(ip, ipz) = gsp 1062 sep%w(ipx, ipx) = gpp 1063 sep%w(ipy, ipy) = gpp 1064 sep%w(ipz, ipz) = gpp 1065 sep%w(ipy, ipx) = gp2 1066 sep%w(ipz, ipx) = gp2 1067 sep%w(ipz, ipy) = gp2 1068 sep%w(ipx, ipy) = gp2 1069 sep%w(ipx, ipz) = gp2 1070 sep%w(ipy, ipz) = gp2 1071 sep%w(ip + 1, ip + 1) = hsp 1072 sep%w(ip + 3, ip + 3) = hsp 1073 sep%w(ip + 6, ip + 6) = hsp 1074 sep%w(ip + 4, ip + 4) = 0.5_dp*(gpp - gp2) 1075 sep%w(ip + 7, ip + 7) = 0.5_dp*(gpp - gp2) 1076 sep%w(ip + 8, ip + 8) = 0.5_dp*(gpp - gp2) 1077 IF (sep%dorb) THEN 1078 ij0 = ip - 1 1079 DO i = 1, 243 1080 ij = int_ij(i) 1081 kl = int_kl(i) 1082 ind = int_onec2el(i) 1083 sep%w(ij + ij0, kl + ij0) = sep%onec2el(ind)/evolt 1084 END DO 1085 END IF 1086 END IF 1087 END IF 1088 END SUBROUTINE setup_1c_2el_int 1089 1090END MODULE semi_empirical_par_utils 1091 1092