1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Analytical derivatives of Integrals for semi-empirical methods 8!> \author Teodoro Laino - Zurich University 04.2007 [tlaino] 9!> \par History 10!> 23.11.2007 jhu short range version of integrals 11!> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver 12!> for computing integrals 13!> Teodoro Laino (05.2008) [tlaino] - University of Zurich : analytical 14!> derivatives for d-orbitals 15! ************************************************************************************************** 16MODULE semi_empirical_int_ana 17 18 USE input_constants, ONLY: do_method_am1,& 19 do_method_pchg,& 20 do_method_pdg,& 21 do_method_pm3,& 22 do_method_pm6,& 23 do_method_pm6fm,& 24 do_method_undef,& 25 do_se_IS_kdso_d 26 USE kinds, ONLY: dp 27 USE multipole_types, ONLY: do_multipole_none 28 USE physcon, ONLY: angstrom,& 29 evolt 30 USE semi_empirical_int_arrays, ONLY: & 31 fac_x_to_z, ijkl_ind, ijkl_sym, inddd, inddp, indexa, indexb, indpp, int2c_type, l_index, & 32 map_x_to_z, rij_threshold 33 USE semi_empirical_int_num, ONLY: nucint_d_num,& 34 nucint_sp_num,& 35 terep_d_num,& 36 terep_sp_num 37 USE semi_empirical_int_utils, ONLY: d_ijkl_d,& 38 d_ijkl_sp,& 39 rot_2el_2c_first,& 40 rotmat,& 41 store_2el_2c_diag 42 USE semi_empirical_types, ONLY: rotmat_create,& 43 rotmat_release,& 44 rotmat_type,& 45 se_int_control_type,& 46 se_int_screen_type,& 47 se_taper_type,& 48 semi_empirical_type,& 49 setup_se_int_control_type 50 USE taper_types, ONLY: dtaper_eval,& 51 taper_eval 52#include "./base/base_uses.f90" 53 54 IMPLICIT NONE 55 PRIVATE 56#include "semi_empirical_int_debug.h" 57#include "semi_empirical_int_args.h" 58 59 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_int_ana' 60 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. 61 PUBLIC :: rotnuc_ana, rotint_ana, corecore_ana, corecore_el_ana 62 63CONTAINS 64 65! ************************************************************************************************** 66!> \brief Computes analytical gradients for semiempirical integrals 67!> \param sepi Atomic parameters of first atom 68!> \param sepj Atomic parameters of second atom 69!> \param rijv Coordinate vector i -> j 70!> \param itype ... 71!> \param e1b Array of electron-nuclear attraction integrals, Electron on atom ni attracting nucleus of nj. 72!> \param e2a Array of electron-nuclear attraction integrals, Electron on atom nj attracting nucleus of ni. 73!> \param de1b derivative of e1b term 74!> \param de2a derivative of e2a term 75!> \param se_int_control input parameters that control the calculation of SE 76!> integrals (shortrange, R3 residual, screening type) 77!> \param se_taper ... 78!> \par History 79!> 04.2007 created [tlaino] 80!> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver 81!> for computing integrals 82!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : removed the core-core part 83!> \author Teodoro Laino [tlaino] - Zurich University 84!> \note 85!> Analytical version of the MOPAC rotnuc routine 86! ************************************************************************************************** 87 RECURSIVE SUBROUTINE rotnuc_ana(sepi, sepj, rijv, itype, e1b, e2a, de1b, de2a, & 88 se_int_control, se_taper) 89 TYPE(semi_empirical_type), POINTER :: sepi, sepj 90 REAL(dp), DIMENSION(3), INTENT(IN) :: rijv 91 INTEGER, INTENT(IN) :: itype 92 REAL(dp), DIMENSION(45), INTENT(OUT), OPTIONAL :: e1b, e2a 93 REAL(dp), DIMENSION(3, 45), INTENT(OUT), OPTIONAL :: de1b, de2a 94 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 95 TYPE(se_taper_type), POINTER :: se_taper 96 97 CHARACTER(len=*), PARAMETER :: routineN = 'rotnuc_ana', routineP = moduleN//':'//routineN 98 99 INTEGER :: i, idd, idp, ind1, ind2, ipp, j, & 100 last_orbital(2), m, n 101 LOGICAL :: invert, l_de1b, l_de2a, l_e1b, l_e2a, & 102 lgrad, task(2) 103 REAL(KIND=dp) :: rij, xtmp 104 REAL(KIND=dp), DIMENSION(10, 2) :: core, dcore 105 REAL(KIND=dp), DIMENSION(3) :: drij 106 REAL(KIND=dp), DIMENSION(3, 45) :: tmp_d 107 REAL(KIND=dp), DIMENSION(45) :: tmp 108 TYPE(rotmat_type), POINTER :: ij_matrix 109 110 NULLIFY (ij_matrix) 111 rij = DOT_PRODUCT(rijv, rijv) 112 ! Initialization 113 l_e1b = PRESENT(e1b) 114 l_e2a = PRESENT(e2a) 115 l_de1b = PRESENT(de1b) 116 l_de2a = PRESENT(de2a) 117 lgrad = l_de1b .OR. l_de2a 118 119 IF (rij > rij_threshold) THEN 120 ! Compute Integrals in diatomic frame opportunely inverted 121 rij = SQRT(rij) 122 ! Create the rotation matrix 123 CALL rotmat_create(ij_matrix) 124 CALL rotmat(sepi, sepj, rijv, rij, ij_matrix, do_derivatives=lgrad, do_invert=invert) 125 126 IF (lgrad) THEN 127 drij(1) = rijv(1)/rij 128 drij(2) = rijv(2)/rij 129 drij(3) = rijv(3)/rij 130 ! Possibly Invert Frame 131 IF (invert) THEN 132 xtmp = drij(3) 133 drij(3) = drij(1) 134 drij(1) = xtmp 135 END IF 136 END IF 137 138 CALL dcore_nucint_ana(sepi, sepj, rij, core=core, dcore=dcore, itype=itype, se_taper=se_taper, & 139 se_int_control=se_int_control, lgrad=lgrad) 140 141 ! Copy parameters over to arrays for do loop. 142 last_orbital(1) = sepi%natorb 143 last_orbital(2) = sepj%natorb 144 task(1) = l_e1b 145 task(2) = l_e2a 146 DO n = 1, 2 147 IF (.NOT. task(n)) CYCLE 148 DO i = 1, last_orbital(n) 149 ind1 = i - 1 150 DO j = 1, i 151 ind2 = j - 1 152 m = (i*(i - 1))/2 + j 153 ! Perform Rotations ... 154 IF (ind2 == 0) THEN 155 IF (ind1 == 0) THEN 156 ! Type of Integral (SS/) 157 tmp(m) = core(1, n) 158 ELSE IF (ind1 < 4) THEN 159 ! Type of Integral (SP/) 160 tmp(m) = ij_matrix%sp(1, ind1)*core(2, n) 161 ELSE 162 ! Type of Integral (SD/) 163 tmp(m) = ij_matrix%sd(1, ind1 - 3)*core(5, n) 164 END IF 165 ELSE IF (ind2 < 4) THEN 166 IF (ind1 < 4) THEN 167 ! Type of Integral (PP/) 168 ipp = indpp(ind1, ind2) 169 tmp(m) = core(3, n)*ij_matrix%pp(ipp, 1, 1) + & 170 core(4, n)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3)) 171 ELSE 172 ! Type of Integral (PD/) 173 idp = inddp(ind1 - 3, ind2) 174 tmp(m) = core(6, n)*ij_matrix%pd(idp, 1, 1) + & 175 core(8, n)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3)) 176 END IF 177 ELSE 178 ! Type of Integral (DD/) 179 idd = inddd(ind1 - 3, ind2 - 3) 180 tmp(m) = core(7, n)*ij_matrix%dd(idd, 1, 1) + & 181 core(9, n)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + & 182 core(10, n)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5)) 183 END IF 184 END DO 185 END DO 186 IF (n == 1) THEN 187 DO i = 1, sepi%atm_int_size 188 e1b(i) = -tmp(i) 189 END DO 190 END IF 191 IF (n == 2) THEN 192 DO i = 1, sepj%atm_int_size 193 e2a(i) = -tmp(i) 194 END DO 195 END IF 196 END DO 197 IF (invert .AND. l_e1b) CALL invert_integral(sepi, sepi, int1el=e1b) 198 IF (invert .AND. l_e2a) CALL invert_integral(sepj, sepj, int1el=e2a) 199 200 ! Possibly compute derivatives 201 task(1) = l_de1b 202 task(2) = l_de2a 203 DO n = 1, 2 204 IF (.NOT. task(n)) CYCLE 205 DO i = 1, last_orbital(n) 206 ind1 = i - 1 207 DO j = 1, i 208 ind2 = j - 1 209 m = (i*(i - 1))/2 + j 210 ! Perform Rotations ... 211 IF (ind2 == 0) THEN 212 IF (ind1 == 0) THEN 213 ! Type of Integral (SS/) 214 tmp_d(1, m) = dcore(1, n)*drij(1) 215 tmp_d(2, m) = dcore(1, n)*drij(2) 216 tmp_d(3, m) = dcore(1, n)*drij(3) 217 ELSE IF (ind1 < 4) THEN 218 ! Type of Integral (SP/) 219 tmp_d(1, m) = ij_matrix%sp_d(1, 1, ind1)*core(2, n) + & 220 ij_matrix%sp(1, ind1)*dcore(2, n)*drij(1) 221 222 tmp_d(2, m) = ij_matrix%sp_d(2, 1, ind1)*core(2, n) + & 223 ij_matrix%sp(1, ind1)*dcore(2, n)*drij(2) 224 225 tmp_d(3, m) = ij_matrix%sp_d(3, 1, ind1)*core(2, n) + & 226 ij_matrix%sp(1, ind1)*dcore(2, n)*drij(3) 227 ELSE 228 ! Type of Integral (SD/) 229 tmp_d(1, m) = ij_matrix%sd_d(1, 1, ind1 - 3)*core(5, n) + & 230 ij_matrix%sd(1, ind1 - 3)*dcore(5, n)*drij(1) 231 232 tmp_d(2, m) = ij_matrix%sd_d(2, 1, ind1 - 3)*core(5, n) + & 233 ij_matrix%sd(1, ind1 - 3)*dcore(5, n)*drij(2) 234 235 tmp_d(3, m) = ij_matrix%sd_d(3, 1, ind1 - 3)*core(5, n) + & 236 ij_matrix%sd(1, ind1 - 3)*dcore(5, n)*drij(3) 237 END IF 238 ELSE IF (ind2 < 4) THEN 239 IF (ind1 < 4) THEN 240 ! Type of Integral (PP/) 241 ipp = indpp(ind1, ind2) 242 tmp_d(1, m) = dcore(3, n)*drij(1)*ij_matrix%pp(ipp, 1, 1) + & 243 core(3, n)*ij_matrix%pp_d(1, ipp, 1, 1) + & 244 dcore(4, n)*drij(1)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3)) + & 245 core(4, n)*(ij_matrix%pp_d(1, ipp, 2, 2) + ij_matrix%pp_d(1, ipp, 3, 3)) 246 247 tmp_d(2, m) = dcore(3, n)*drij(2)*ij_matrix%pp(ipp, 1, 1) + & 248 core(3, n)*ij_matrix%pp_d(2, ipp, 1, 1) + & 249 dcore(4, n)*drij(2)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3)) + & 250 core(4, n)*(ij_matrix%pp_d(2, ipp, 2, 2) + ij_matrix%pp_d(2, ipp, 3, 3)) 251 252 tmp_d(3, m) = dcore(3, n)*drij(3)*ij_matrix%pp(ipp, 1, 1) + & 253 core(3, n)*ij_matrix%pp_d(3, ipp, 1, 1) + & 254 dcore(4, n)*drij(3)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3)) + & 255 core(4, n)*(ij_matrix%pp_d(3, ipp, 2, 2) + ij_matrix%pp_d(3, ipp, 3, 3)) 256 ELSE 257 ! Type of Integral (PD/) 258 idp = inddp(ind1 - 3, ind2) 259 tmp_d(1, m) = dcore(6, n)*drij(1)*ij_matrix%pd(idp, 1, 1) + & 260 core(6, n)*ij_matrix%pd_d(1, idp, 1, 1) + & 261 dcore(8, n)*drij(1)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3)) + & 262 core(8, n)*(ij_matrix%pd_d(1, idp, 2, 2) + ij_matrix%pd_d(1, idp, 3, 3)) 263 264 tmp_d(2, m) = dcore(6, n)*drij(2)*ij_matrix%pd(idp, 1, 1) + & 265 core(6, n)*ij_matrix%pd_d(2, idp, 1, 1) + & 266 dcore(8, n)*drij(2)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3)) + & 267 core(8, n)*(ij_matrix%pd_d(2, idp, 2, 2) + ij_matrix%pd_d(2, idp, 3, 3)) 268 269 tmp_d(3, m) = dcore(6, n)*drij(3)*ij_matrix%pd(idp, 1, 1) + & 270 core(6, n)*ij_matrix%pd_d(3, idp, 1, 1) + & 271 dcore(8, n)*drij(3)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3)) + & 272 core(8, n)*(ij_matrix%pd_d(3, idp, 2, 2) + ij_matrix%pd_d(3, idp, 3, 3)) 273 END IF 274 ELSE 275 ! Type of Integral (DD/) 276 idd = inddd(ind1 - 3, ind2 - 3) 277 tmp_d(1, m) = dcore(7, n)*drij(1)*ij_matrix%dd(idd, 1, 1) + & 278 core(7, n)*ij_matrix%dd_d(1, idd, 1, 1) + & 279 dcore(9, n)*drij(1)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + & 280 core(9, n)*(ij_matrix%dd_d(1, idd, 2, 2) + ij_matrix%dd_d(1, idd, 3, 3)) + & 281 dcore(10, n)*drij(1)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5)) + & 282 core(10, n)*(ij_matrix%dd_d(1, idd, 4, 4) + ij_matrix%dd_d(1, idd, 5, 5)) 283 284 tmp_d(2, m) = dcore(7, n)*drij(2)*ij_matrix%dd(idd, 1, 1) + & 285 core(7, n)*ij_matrix%dd_d(2, idd, 1, 1) + & 286 dcore(9, n)*drij(2)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + & 287 core(9, n)*(ij_matrix%dd_d(2, idd, 2, 2) + ij_matrix%dd_d(2, idd, 3, 3)) + & 288 dcore(10, n)*drij(2)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5)) + & 289 core(10, n)*(ij_matrix%dd_d(2, idd, 4, 4) + ij_matrix%dd_d(2, idd, 5, 5)) 290 291 tmp_d(3, m) = dcore(7, n)*drij(3)*ij_matrix%dd(idd, 1, 1) + & 292 core(7, n)*ij_matrix%dd_d(3, idd, 1, 1) + & 293 dcore(9, n)*drij(3)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + & 294 core(9, n)*(ij_matrix%dd_d(3, idd, 2, 2) + ij_matrix%dd_d(3, idd, 3, 3)) + & 295 dcore(10, n)*drij(3)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5)) + & 296 core(10, n)*(ij_matrix%dd_d(3, idd, 4, 4) + ij_matrix%dd_d(3, idd, 5, 5)) 297 END IF 298 END DO 299 END DO 300 IF (n == 1) THEN 301 DO i = 1, sepi%atm_int_size 302 de1b(1, i) = -tmp_d(1, i) 303 de1b(2, i) = -tmp_d(2, i) 304 de1b(3, i) = -tmp_d(3, i) 305 END DO 306 END IF 307 IF (n == 2) THEN 308 DO i = 1, sepj%atm_int_size 309 de2a(1, i) = -tmp_d(1, i) 310 de2a(2, i) = -tmp_d(2, i) 311 de2a(3, i) = -tmp_d(3, i) 312 END DO 313 END IF 314 END DO 315 IF (invert .AND. l_de1b) CALL invert_derivative(sepi, sepi, dint1el=de1b) 316 IF (invert .AND. l_de2a) CALL invert_derivative(sepj, sepj, dint1el=de2a) 317 CALL rotmat_release(ij_matrix) 318 319 ! Possibly debug the analytical values versus the numerical ones 320 IF (debug_this_module) THEN 321 CALL check_drotnuc_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, e1b, e2a, de1b, de2a) 322 END IF 323 END IF 324 END SUBROUTINE rotnuc_ana 325 326! ************************************************************************************************** 327!> \brief Computes analytical gradients for semiempirical core-core interaction. 328!> \param sepi Atomic parameters of first atom 329!> \param sepj Atomic parameters of second atom 330!> \param rijv Coordinate vector i -> j 331!> \param itype ... 332!> \param enuc nuclear-nuclear repulsion term. 333!> \param denuc derivative of nuclear-nuclear repulsion term. 334!> \param se_int_control input parameters that control the calculation of SE 335!> integrals (shortrange, R3 residual, screening type) 336!> \param se_taper ... 337!> \par History 338!> 04.2007 created [tlaino] 339!> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver 340!> for computing integrals 341!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : removed the 342!> core-core part 343!> \author Teodoro Laino [tlaino] - Zurich University 344!> \note 345!> Analytical version of the MOPAC rotnuc routine 346! ************************************************************************************************** 347 RECURSIVE SUBROUTINE corecore_ana(sepi, sepj, rijv, itype, enuc, denuc, se_int_control, se_taper) 348 TYPE(semi_empirical_type), POINTER :: sepi, sepj 349 REAL(dp), DIMENSION(3), INTENT(IN) :: rijv 350 INTEGER, INTENT(IN) :: itype 351 REAL(dp), INTENT(OUT), OPTIONAL :: enuc 352 REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: denuc 353 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 354 TYPE(se_taper_type), POINTER :: se_taper 355 356 CHARACTER(len=*), PARAMETER :: routineN = 'corecore_ana', routineP = moduleN//':'//routineN 357 358 INTEGER :: ig, nt 359 LOGICAL :: l_denuc, l_enuc 360 REAL(dp) :: aab, alpi, alpj, apdg, ax, dai, daj, dax, dbi, dbj, denuc_loc, dqcorr, drija, & 361 dscale, dssss, dssss_sr, dtmp, dzz, enuc_loc, pai, paj, pbi, pbj, qcorr, rij, rija, & 362 scale, ssss, ssss_sr, tmp, xab, xtmp, zaf, zbf, zz 363 REAL(dp), DIMENSION(3) :: drij 364 REAL(dp), DIMENSION(4) :: fni1, fni2, fni3, fnj1, fnj2, fnj3 365 TYPE(se_int_control_type) :: se_int_control_off 366 367 rij = DOT_PRODUCT(rijv, rijv) 368 ! Initialization 369 l_enuc = PRESENT(enuc) 370 l_denuc = PRESENT(denuc) 371 IF ((rij > rij_threshold) .AND. (l_enuc .OR. l_denuc)) THEN 372 ! Compute Integrals in diatomic frame 373 rij = SQRT(rij) 374 CALL setup_se_int_control_type(se_int_control_off, shortrange=.FALSE., do_ewald_r3=.FALSE., & 375 do_ewald_gks=.FALSE., integral_screening=se_int_control%integral_screening, & 376 max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.) 377 CALL dssss_nucint_ana(sepi, sepj, rij, ssss=ssss, dssss=dssss, itype=itype, se_taper=se_taper, & 378 se_int_control=se_int_control_off, lgrad=l_denuc) 379 ! In case let's compute the short-range part of the (ss|ss) integral 380 IF (se_int_control%shortrange) THEN 381 CALL dssss_nucint_ana(sepi, sepj, rij, ssss=ssss_sr, dssss=dssss_sr, itype=itype, & 382 se_taper=se_taper, se_int_control=se_int_control, lgrad=l_denuc) 383 ELSE 384 ssss_sr = ssss 385 dssss_sr = dssss 386 END IF 387 ! Zeroing local method dependent core-core corrections 388 enuc_loc = 0.0_dp 389 denuc_loc = 0.0_dp 390 qcorr = 0.0_dp 391 scale = 0.0_dp 392 dscale = 0.0_dp 393 dqcorr = 0.0_dp 394 zz = sepi%zeff*sepj%zeff 395 ! Core Core electrostatic contribution 396 IF (l_enuc) enuc_loc = zz*ssss_sr 397 IF (l_denuc) denuc_loc = zz*dssss_sr 398 ! Method dependent code 399 tmp = zz*ssss 400 IF (l_denuc) dtmp = zz*dssss 401 IF (itype /= do_method_pm6 .AND. itype /= do_method_pm6fm) THEN 402 alpi = sepi%alp 403 alpj = sepj%alp 404 scale = EXP(-alpi*rij) + EXP(-alpj*rij) 405 IF (l_denuc) THEN 406 dscale = -alpi*EXP(-alpi*rij) - alpj*EXP(-alpj*rij) 407 END IF 408 nt = sepi%z + sepj%z 409 IF (nt == 8 .OR. nt == 9) THEN 410 IF (sepi%z == 7 .OR. sepi%z == 8) THEN 411 scale = scale + (angstrom*rij - 1._dp)*EXP(-alpi*rij) 412 IF (l_denuc) THEN 413 dscale = dscale + angstrom*EXP(-alpi*rij) - (angstrom*rij - 1._dp)*alpi*EXP(-alpi*rij) 414 END IF 415 END IF 416 IF (sepj%z == 7 .OR. sepj%z == 8) THEN 417 scale = scale + (angstrom*rij - 1._dp)*EXP(-alpj*rij) 418 IF (l_denuc) THEN 419 dscale = dscale + angstrom*EXP(-alpj*rij) - (angstrom*rij - 1._dp)*alpj*EXP(-alpj*rij) 420 END IF 421 END IF 422 ENDIF 423 IF (l_denuc) THEN 424 dscale = SIGN(1.0_dp, scale*tmp)*(dscale*tmp + scale*dtmp) 425 dzz = -zz/rij**2 426 END IF 427 scale = ABS(scale*tmp) 428 zz = zz/rij 429 IF (itype == do_method_am1 .OR. itype == do_method_pm3 .OR. itype == do_method_pdg) THEN 430 IF (itype == do_method_am1 .AND. sepi%z == 5) THEN 431 !special case AM1 Boron 432 SELECT CASE (sepj%z) 433 CASE DEFAULT 434 nt = 1 435 CASE (1) 436 nt = 2 437 CASE (6) 438 nt = 3 439 CASE (9, 17, 35, 53) 440 nt = 4 441 END SELECT 442 fni1(1) = sepi%bfn1(1, nt) 443 fni1(2) = sepi%bfn1(2, nt) 444 fni1(3) = sepi%bfn1(3, nt) 445 fni1(4) = sepi%bfn1(4, nt) 446 fni2(1) = sepi%bfn2(1, nt) 447 fni2(2) = sepi%bfn2(2, nt) 448 fni2(3) = sepi%bfn2(3, nt) 449 fni2(4) = sepi%bfn2(4, nt) 450 fni3(1) = sepi%bfn3(1, nt) 451 fni3(2) = sepi%bfn3(2, nt) 452 fni3(3) = sepi%bfn3(3, nt) 453 fni3(4) = sepi%bfn3(4, nt) 454 ELSE 455 fni1(1) = sepi%fn1(1) 456 fni1(2) = sepi%fn1(2) 457 fni1(3) = sepi%fn1(3) 458 fni1(4) = sepi%fn1(4) 459 fni2(1) = sepi%fn2(1) 460 fni2(2) = sepi%fn2(2) 461 fni2(3) = sepi%fn2(3) 462 fni2(4) = sepi%fn2(4) 463 fni3(1) = sepi%fn3(1) 464 fni3(2) = sepi%fn3(2) 465 fni3(3) = sepi%fn3(3) 466 fni3(4) = sepi%fn3(4) 467 END IF 468 IF (itype == do_method_am1 .AND. sepj%z == 5) THEN 469 !special case AM1 Boron 470 SELECT CASE (sepi%z) 471 CASE DEFAULT 472 nt = 1 473 CASE (1) 474 nt = 2 475 CASE (6) 476 nt = 3 477 CASE (9, 17, 35, 53) 478 nt = 4 479 END SELECT 480 fnj1(1) = sepj%bfn1(1, nt) 481 fnj1(2) = sepj%bfn1(2, nt) 482 fnj1(3) = sepj%bfn1(3, nt) 483 fnj1(4) = sepj%bfn1(4, nt) 484 fnj2(1) = sepj%bfn2(1, nt) 485 fnj2(2) = sepj%bfn2(2, nt) 486 fnj2(3) = sepj%bfn2(3, nt) 487 fnj2(4) = sepj%bfn2(4, nt) 488 fnj3(1) = sepj%bfn3(1, nt) 489 fnj3(2) = sepj%bfn3(2, nt) 490 fnj3(3) = sepj%bfn3(3, nt) 491 fnj3(4) = sepj%bfn3(4, nt) 492 ELSE 493 fnj1(1) = sepj%fn1(1) 494 fnj1(2) = sepj%fn1(2) 495 fnj1(3) = sepj%fn1(3) 496 fnj1(4) = sepj%fn1(4) 497 fnj2(1) = sepj%fn2(1) 498 fnj2(2) = sepj%fn2(2) 499 fnj2(3) = sepj%fn2(3) 500 fnj2(4) = sepj%fn2(4) 501 fnj3(1) = sepj%fn3(1) 502 fnj3(2) = sepj%fn3(2) 503 fnj3(3) = sepj%fn3(3) 504 fnj3(4) = sepj%fn3(4) 505 END IF 506 ! AM1/PM3/PDG correction to nuclear repulsion 507 DO ig = 1, SIZE(fni1) 508 IF (ABS(fni1(ig)) > 0._dp) THEN 509 ax = fni2(ig)*(rij - fni3(ig))**2 510 IF (ax <= 25._dp) THEN 511 scale = scale + zz*fni1(ig)*EXP(-ax) 512 IF (l_denuc) THEN 513 dax = fni2(ig)*2.0_dp*(rij - fni3(ig)) 514 dscale = dscale + dzz*fni1(ig)*EXP(-ax) - dax*zz*fni1(ig)*EXP(-ax) 515 END IF 516 ENDIF 517 ENDIF 518 IF (ABS(fnj1(ig)) > 0._dp) THEN 519 ax = fnj2(ig)*(rij - fnj3(ig))**2 520 IF (ax <= 25._dp) THEN 521 scale = scale + zz*fnj1(ig)*EXP(-ax) 522 IF (l_denuc) THEN 523 dax = fnj2(ig)*2.0_dp*(rij - fnj3(ig)) 524 dscale = dscale + dzz*fnj1(ig)*EXP(-ax) - dax*zz*fnj1(ig)*EXP(-ax) 525 END IF 526 ENDIF 527 ENDIF 528 END DO 529 ENDIF 530 IF (itype == do_method_pdg) THEN 531 ! PDDG function 532 zaf = sepi%zeff/nt 533 zbf = sepj%zeff/nt 534 pai = sepi%pre(1) 535 pbi = sepi%pre(2) 536 paj = sepj%pre(1) 537 pbj = sepj%pre(2) 538 dai = sepi%d(1) 539 dbi = sepi%d(2) 540 daj = sepj%d(1) 541 dbj = sepj%d(2) 542 apdg = 10._dp*angstrom**2 543 qcorr = (zaf*pai + zbf*paj)*EXP(-apdg*(rij - dai - daj)**2) + & 544 (zaf*pai + zbf*pbj)*EXP(-apdg*(rij - dai - dbj)**2) + & 545 (zaf*pbi + zbf*paj)*EXP(-apdg*(rij - dbi - daj)**2) + & 546 (zaf*pbi + zbf*pbj)*EXP(-apdg*(rij - dbi - dbj)**2) 547 IF (l_denuc) THEN 548 dqcorr = (zaf*pai + zbf*paj)*EXP(-apdg*(rij - dai - daj)**2)*(-2.0_dp*apdg*(rij - dai - daj)) + & 549 (zaf*pai + zbf*pbj)*EXP(-apdg*(rij - dai - dbj)**2)*(-2.0_dp*apdg*(rij - dai - dbj)) + & 550 (zaf*pbi + zbf*paj)*EXP(-apdg*(rij - dbi - daj)**2)*(-2.0_dp*apdg*(rij - dbi - daj)) + & 551 (zaf*pbi + zbf*pbj)*EXP(-apdg*(rij - dbi - dbj)**2)*(-2.0_dp*apdg*(rij - dbi - dbj)) 552 END IF 553 ELSEIF (itype == do_method_pchg) THEN 554 qcorr = 0.0_dp 555 scale = 0.0_dp 556 dscale = 0.0_dp 557 dqcorr = 0.0_dp 558 ELSE 559 qcorr = 0.0_dp 560 dqcorr = 0.0_dp 561 END IF 562 ELSE 563 ! PM6 core-core terms 564 scale = tmp 565 IF (l_denuc) dscale = dtmp 566 drija = angstrom 567 rija = rij*drija 568 xab = sepi%xab(sepj%z) 569 aab = sepi%aab(sepj%z) 570 IF ((sepi%z == 1 .AND. (sepj%z == 6 .OR. sepj%z == 7 .OR. sepj%z == 8)) .OR. & 571 (sepj%z == 1 .AND. (sepi%z == 6 .OR. sepi%z == 7 .OR. sepi%z == 8))) THEN 572 ! Special Case O-H or N-H or C-H 573 IF (l_denuc) dscale = dscale*(2._dp*xab*EXP(-aab*rija*rija)) - & 574 scale*2._dp*xab*EXP(-aab*rija*rija)*(2.0_dp*aab*rija)*drija 575 IF (l_enuc) scale = scale*(2._dp*xab*EXP(-aab*rija*rija)) 576 ELSEIF (sepi%z == 6 .AND. sepj%z == 6) THEN 577 ! Special Case C-C 578 IF (l_denuc) dscale = & 579 dscale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6)) + 9.28_dp*EXP(-5.98_dp*rija)) & 580 - scale*2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6))*aab*(1.0_dp + 6.0_dp*0.0003_dp*rija**5)*drija & 581 - scale*9.28_dp*EXP(-5.98_dp*rija)*5.98_dp*drija 582 IF (l_enuc) scale = scale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6)) + 9.28_dp*EXP(-5.98_dp*rija)) 583 ELSEIF ((sepi%z == 8 .AND. sepj%z == 14) .OR. & 584 (sepj%z == 8 .AND. sepi%z == 14)) THEN 585 ! Special Case Si-O 586 IF (l_denuc) dscale = & 587 dscale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6)) - 0.0007_dp*EXP(-(rija - 2.9_dp)**2)) & 588 - scale*2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6))*aab*(1.0_dp + 6.0_dp*0.0003_dp*rija**5)*drija + & 589 scale*0.0007_dp*EXP(-(rija - 2.9_dp)**2)*(2.0_dp*(rija - 2.9_dp)*drija) 590 IF (l_enuc) scale = scale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6)) - 0.0007_dp*EXP(-(rija - 2.9_dp)**2)) 591 ELSE 592 ! General Case 593 ! Factor of 2 found by experiment 594 IF (l_denuc) dscale = dscale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6))) & 595 - scale*2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6))*aab*(1.0_dp + 6.0_dp*0.0003_dp*rija**5)*drija 596 IF (l_enuc) scale = scale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6))) 597 END IF 598 ! General correction term a*exp(-b*(rij-c)^2) 599 xtmp = 1.e-8_dp/evolt*((REAL(sepi%z, dp)**(1._dp/3._dp) + REAL(sepj%z, dp)**(1._dp/3._dp))/rija)**12 600 IF (l_enuc) THEN 601 qcorr = (sepi%a*EXP(-sepi%b*(rij - sepi%c)**2))*zz/rij + & 602 (sepj%a*EXP(-sepj%b*(rij - sepj%c)**2))*zz/rij + & 603 ! Hard core repulsion 604 xtmp 605 END IF 606 IF (l_denuc) THEN 607 dqcorr = (sepi%a*EXP(-sepi%b*(rij - sepi%c)**2)*(-2.0_dp*sepi%b*(rij - sepi%c)))*zz/rij - & 608 (sepi%a*EXP(-sepi%b*(rij - sepi%c)**2))*zz/rij**2 + & 609 (sepj%a*EXP(-sepj%b*(rij - sepj%c)**2)*(-2.0_dp*sepj%b*(rij - sepj%c)))*zz/rij - & 610 (sepj%a*EXP(-sepj%b*(rij - sepj%c)**2))*zz/rij**2 + & 611 ! Hard core repulsion 612 (-12.0_dp*xtmp/rija*drija) 613 END IF 614 ENDIF 615 616 ! Only at the very end let's sum-up the several contributions energy/derivatives 617 ! This assignment should be method independent 618 IF (l_enuc) THEN 619 enuc = enuc_loc + scale + qcorr 620 END IF 621 IF (l_denuc) THEN 622 drij(1) = rijv(1)/rij 623 drij(2) = rijv(2)/rij 624 drij(3) = rijv(3)/rij 625 denuc = (denuc_loc + dscale + dqcorr)*drij 626 END IF 627 ! Debug statement 628 IF (debug_this_module) THEN 629 CALL check_dcorecore_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, enuc, denuc) 630 END IF 631 ENDIF 632 END SUBROUTINE corecore_ana 633 634! ************************************************************************************************** 635!> \brief Computes analytical gradients for semiempirical core-core electrostatic 636!> interaction only. 637!> \param sepi Atomic parameters of first atom 638!> \param sepj Atomic parameters of second atom 639!> \param rijv Coordinate vector i -> j 640!> \param itype ... 641!> \param enuc nuclear-nuclear electrostatic repulsion term. 642!> \param denuc derivative of nuclear-nuclear electrostatic 643!> repulsion term. 644!> \param se_int_control input parameters that control the calculation of SE 645!> integrals (shortrange, R3 residual, screening type) 646!> \param se_taper ... 647!> \par History 648!> 04.2007 created [tlaino] 649!> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver 650!> for computing integrals 651!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : removed the 652!> core-core part 653!> \author Teodoro Laino [tlaino] - Zurich University 654!> \note 655!> Analytical version of the MOPAC rotnuc routine 656! ************************************************************************************************** 657 RECURSIVE SUBROUTINE corecore_el_ana(sepi, sepj, rijv, itype, enuc, denuc, & 658 se_int_control, se_taper) 659 TYPE(semi_empirical_type), POINTER :: sepi, sepj 660 REAL(dp), DIMENSION(3), INTENT(IN) :: rijv 661 INTEGER, INTENT(IN) :: itype 662 REAL(dp), INTENT(OUT), OPTIONAL :: enuc 663 REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: denuc 664 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 665 TYPE(se_taper_type), POINTER :: se_taper 666 667 CHARACTER(len=*), PARAMETER :: routineN = 'corecore_el_ana', & 668 routineP = moduleN//':'//routineN 669 670 LOGICAL :: l_denuc, l_enuc 671 REAL(dp) :: drij(3), dssss, dssss_sr, rij, ssss, & 672 ssss_sr, tmp, zz 673 TYPE(se_int_control_type) :: se_int_control_off 674 675 rij = DOT_PRODUCT(rijv, rijv) 676 ! Initialization 677 l_enuc = PRESENT(enuc) 678 l_denuc = PRESENT(denuc) 679 IF ((rij > rij_threshold) .AND. (l_enuc .OR. l_denuc)) THEN 680 ! Compute Integrals in diatomic frame 681 rij = SQRT(rij) 682 CALL setup_se_int_control_type(se_int_control_off, shortrange=.FALSE., do_ewald_r3=.FALSE., & 683 do_ewald_gks=.FALSE., integral_screening=se_int_control%integral_screening, & 684 max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.) 685 CALL dssss_nucint_ana(sepi, sepj, rij, ssss=ssss, dssss=dssss, itype=itype, se_taper=se_taper, & 686 se_int_control=se_int_control_off, lgrad=l_denuc) 687 ! In case let's compute the short-range part of the (ss|ss) integral 688 IF (se_int_control%shortrange .OR. se_int_control%pc_coulomb_int) THEN 689 CALL dssss_nucint_ana(sepi, sepj, rij, ssss=ssss_sr, dssss=dssss_sr, itype=itype, & 690 se_taper=se_taper, se_int_control=se_int_control, lgrad=l_denuc) 691 ELSE 692 ssss_sr = ssss 693 dssss_sr = dssss 694 END IF 695 zz = sepi%zeff*sepj%zeff 696 ! Core Core electrostatic contribution 697 IF (l_enuc) enuc = zz*ssss_sr 698 IF (l_denuc) THEN 699 drij(1) = rijv(1)/rij 700 drij(2) = rijv(2)/rij 701 drij(3) = rijv(3)/rij 702 tmp = zz*dssss_sr 703 denuc = tmp*drij 704 END IF 705 END IF 706 END SUBROUTINE corecore_el_ana 707 708! ************************************************************************************************** 709!> \brief Exploits inversion symmetry to avoid divergence 710!> \param sepi ... 711!> \param sepj ... 712!> \param int1el ... 713!> \param int2el ... 714!> \par History 715!> 04.2007 created [tlaino] 716!> 05.2008 New driver for integral invertion (supports d-orbitals) 717!> \author Teodoro Laino - Zurich University 718! ************************************************************************************************** 719 SUBROUTINE invert_integral(sepi, sepj, int1el, int2el) 720 TYPE(semi_empirical_type), POINTER :: sepi, sepj 721 REAL(dp), DIMENSION(:), INTENT(INOUT), OPTIONAL :: int1el, int2el 722 723 CHARACTER(len=*), PARAMETER :: routineN = 'invert_integral', & 724 routineP = moduleN//':'//routineN 725 726 INTEGER :: fdim, gind, gknd, i, imap, ind, j, jmap, & 727 jnd, k, kmap, knd, l, lmap, lnd, ndim, & 728 sdim, tdim, tind 729 REAL(KIND=dp) :: ifac, jfac, kfac, lfac 730 REAL(KIND=dp), DIMENSION(2025) :: tmp2el 731 REAL(KIND=dp), DIMENSION(45) :: tmp1el 732 733! One-electron integral 734 735 IF (PRESENT(int1el)) THEN 736 fdim = sepi%atm_int_size 737 ndim = 0 738 DO i = 1, fdim 739 tmp1el(i) = 0.0_dp 740 END DO 741 DO i = 1, sepi%natorb 742 DO j = 1, i 743 ndim = ndim + 1 744 745 ! Get the integral in the original frame (along z) 746 DO ind = 1, 2 747 imap = map_x_to_z(ind, i) 748 IF (imap == 0) CYCLE 749 ifac = fac_x_to_z(ind, i) 750 DO jnd = 1, 2 751 jmap = map_x_to_z(jnd, j) 752 IF (jmap == 0) CYCLE 753 jfac = fac_x_to_z(jnd, j) 754 gind = indexb(imap, jmap) 755 756 tmp1el(ndim) = tmp1el(ndim) + ifac*jfac*int1el(gind) 757 END DO 758 END DO 759 END DO 760 END DO 761 DO i = 1, fdim 762 int1el(i) = tmp1el(i) 763 END DO 764 END IF 765 766 ! Two electron integrals 767 IF (PRESENT(int2el)) THEN 768 sdim = sepi%atm_int_size 769 tdim = sepj%atm_int_size 770 fdim = sdim*tdim 771 ndim = 0 772 DO i = 1, fdim 773 tmp2el(i) = 0.0_dp 774 END DO 775 DO i = 1, sepi%natorb 776 DO j = 1, i 777 DO k = 1, sepj%natorb 778 DO l = 1, k 779 ndim = ndim + 1 780 781 ! Get the integral in the original frame (along z) 782 DO ind = 1, 2 783 imap = map_x_to_z(ind, i) 784 IF (imap == 0) CYCLE 785 ifac = fac_x_to_z(ind, i) 786 DO jnd = 1, 2 787 jmap = map_x_to_z(jnd, j) 788 IF (jmap == 0) CYCLE 789 jfac = fac_x_to_z(jnd, j) 790 gind = indexb(imap, jmap) 791 792 ! Get the integral in the original frame (along z) 793 DO knd = 1, 2 794 kmap = map_x_to_z(knd, k) 795 IF (kmap == 0) CYCLE 796 kfac = fac_x_to_z(knd, k) 797 DO lnd = 1, 2 798 lmap = map_x_to_z(lnd, l) 799 IF (lmap == 0) CYCLE 800 lfac = fac_x_to_z(lnd, l) 801 gknd = indexb(kmap, lmap) 802 803 tind = (gind - 1)*tdim + gknd 804 tmp2el(ndim) = tmp2el(ndim) + ifac*jfac*lfac*kfac*int2el(tind) 805 END DO 806 END DO 807 808 END DO 809 END DO 810 811 END DO 812 END DO 813 END DO 814 END DO 815 DO i = 1, fdim 816 int2el(i) = tmp2el(i) 817 END DO 818 END IF 819 END SUBROUTINE invert_integral 820 821! ************************************************************************************************** 822!> \brief Exploits inversion symmetry to avoid divergence 823!> \param sepi ... 824!> \param sepj ... 825!> \param dint1el ... 826!> \param dint2el ... 827!> \par History 828!> 04.2007 created [tlaino] 829!> \author Teodoro Laino - Zurich University 830! ************************************************************************************************** 831 SUBROUTINE invert_derivative(sepi, sepj, dint1el, dint2el) 832 TYPE(semi_empirical_type), POINTER :: sepi, sepj 833 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), & 834 OPTIONAL :: dint1el, dint2el 835 836 CHARACTER(len=*), PARAMETER :: routineN = 'invert_derivative', & 837 routineP = moduleN//':'//routineN 838 839 INTEGER :: i, m 840 REAL(KIND=dp) :: tmp 841 842! Integral part 843 844 DO i = 1, 3 845 IF (PRESENT(dint1el)) THEN 846 CALL invert_integral(sepi, sepj, int1el=dint1el(i, :)) 847 END IF 848 IF (PRESENT(dint2el)) THEN 849 CALL invert_integral(sepi, sepj, int2el=dint2el(i, :)) 850 END IF 851 END DO 852 853 ! Derivatives part 854 IF (PRESENT(dint1el)) THEN 855 DO m = 1, SIZE(dint1el, 2) 856 tmp = dint1el(3, m) 857 dint1el(3, m) = dint1el(1, m) 858 dint1el(1, m) = tmp 859 END DO 860 END IF 861 IF (PRESENT(dint2el)) THEN 862 DO m = 1, SIZE(dint2el, 2) 863 tmp = dint2el(3, m) 864 dint2el(3, m) = dint2el(1, m) 865 dint2el(1, m) = tmp 866 END DO 867 END IF 868 END SUBROUTINE invert_derivative 869 870! ************************************************************************************************** 871!> \brief Calculates the ssss integral and analytical derivatives (main driver) 872!> \param sepi parameters of atom i 873!> \param sepj parameters of atom j 874!> \param rij interatomic distance 875!> \param ssss ... 876!> \param dssss derivative of (ssss) integral 877!> derivatives are intended w.r.t. rij 878!> \param itype ... 879!> \param se_taper ... 880!> \param se_int_control input parameters that control the calculation of SE 881!> integrals (shortrange, R3 residual, screening type) 882!> \param lgrad ... 883!> \par History 884!> 03.2007 created [tlaino] 885!> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver 886!> for computing integrals 887!> \author Teodoro Laino - Zurich University 888!> \note 889!> Analytical version - Analytical evaluation of gradients 890!> Teodoro Laino - Zurich University 04.2007 891!> 892! ************************************************************************************************** 893 SUBROUTINE dssss_nucint_ana(sepi, sepj, rij, ssss, dssss, itype, se_taper, se_int_control, & 894 lgrad) 895 TYPE(semi_empirical_type), POINTER :: sepi, sepj 896 REAL(dp), INTENT(IN) :: rij 897 REAL(dp), INTENT(OUT) :: ssss, dssss 898 INTEGER, INTENT(IN) :: itype 899 TYPE(se_taper_type), POINTER :: se_taper 900 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 901 LOGICAL, INTENT(IN) :: lgrad 902 903 CHARACTER(len=*), PARAMETER :: routineN = 'dssss_nucint_ana', & 904 routineP = moduleN//':'//routineN 905 906 REAL(KIND=dp) :: dft, ft 907 TYPE(se_int_screen_type) :: se_int_screen 908 909! Compute the Tapering function 910 911 ft = 1.0_dp 912 dft = 0.0_dp 913 IF (itype /= do_method_pchg) THEN 914 ft = taper_eval(se_taper%taper, rij) 915 dft = dtaper_eval(se_taper%taper, rij) 916 END IF 917 ! Evaluate additional taper function for dumped integrals 918 IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN 919 se_int_screen%ft = 1.0_dp 920 se_int_screen%dft = 0.0_dp 921 IF (itype /= do_method_pchg) THEN 922 se_int_screen%ft = taper_eval(se_taper%taper_add, rij) 923 se_int_screen%dft = dtaper_eval(se_taper%taper_add, rij) 924 END IF 925 END IF 926 927 ! Value of the integrals for sp shell 928 CALL nucint_sp_num(sepi, sepj, rij, ssss=ssss, itype=itype, se_int_control=se_int_control, & 929 se_int_screen=se_int_screen) 930 931 IF (lgrad) THEN 932 ! Integrals derivatives for sp shell 933 CALL dnucint_sp_ana(sepi, sepj, rij, dssss=dssss, itype=itype, se_int_control=se_int_control, & 934 se_int_screen=se_int_screen) 935 END IF 936 937 ! Tapering the value of the integrals 938 IF (lgrad) THEN 939 dssss = ft*dssss + dft*ssss 940 END IF 941 ssss = ft*ssss 942 943 ! Debug Procedure.. Check valifity of analytical gradients of nucint 944 IF (debug_this_module .AND. lgrad) THEN 945 CALL check_dssss_nucint_ana(sepi, sepj, rij, dssss, itype, se_int_control, se_taper=se_taper) 946 END IF 947 END SUBROUTINE dssss_nucint_ana 948 949! ************************************************************************************************** 950!> \brief Calculates the nuclear attraction integrals and analytical integrals (main driver) 951!> \param sepi parameters of atom i 952!> \param sepj parameters of atom j 953!> \param rij interatomic distance 954!> \param core ... 955!> \param dcore derivative of 4 X 2 array of electron-core attraction integrals 956!> derivatives are intended w.r.t. rij 957!> The storage of the nuclear attraction integrals core(kl/ij) iS 958!> (SS/)=1, (SO/)=2, (OO/)=3, (PP/)=4 959!> where ij=1 if the orbitals centred on atom i, =2 if on atom j. 960!> \param itype type of semi_empirical model 961!> extension to the original routine to compute qm/mm 962!> integrals 963!> \param se_taper ... 964!> \param se_int_control input parameters that control the calculation of SE 965!> integrals (shortrange, R3 residual, screening type) 966!> \param lgrad ... 967!> \par History 968!> 03.2007 created [tlaino] 969!> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver 970!> for computing integrals 971!> \author Teodoro Laino - Zurich University 972!> \note 973!> Analytical version - Analytical evaluation of gradients 974!> Teodoro Laino - Zurich University 04.2007 975!> 976! ************************************************************************************************** 977 SUBROUTINE dcore_nucint_ana(sepi, sepj, rij, core, dcore, itype, se_taper, & 978 se_int_control, lgrad) 979 TYPE(semi_empirical_type), POINTER :: sepi, sepj 980 REAL(dp), INTENT(IN) :: rij 981 REAL(dp), DIMENSION(10, 2), INTENT(OUT) :: core, dcore 982 INTEGER, INTENT(IN) :: itype 983 TYPE(se_taper_type), POINTER :: se_taper 984 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 985 LOGICAL, INTENT(IN) :: lgrad 986 987 CHARACTER(len=*), PARAMETER :: routineN = 'dcore_nucint_ana', & 988 routineP = moduleN//':'//routineN 989 990 INTEGER :: i 991 REAL(KIND=dp) :: dft, ft 992 TYPE(se_int_screen_type) :: se_int_screen 993 994! Compute the Tapering function 995 996 ft = 1.0_dp 997 dft = 0.0_dp 998 IF (itype /= do_method_pchg) THEN 999 ft = taper_eval(se_taper%taper, rij) 1000 dft = dtaper_eval(se_taper%taper, rij) 1001 END IF 1002 ! Evaluate additional taper function for dumped integrals 1003 IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN 1004 se_int_screen%ft = 1.0_dp 1005 se_int_screen%dft = 0.0_dp 1006 IF (itype /= do_method_pchg) THEN 1007 se_int_screen%ft = taper_eval(se_taper%taper_add, rij) 1008 se_int_screen%dft = dtaper_eval(se_taper%taper_add, rij) 1009 END IF 1010 END IF 1011 1012 ! Value of the integrals for sp shell 1013 CALL nucint_sp_num(sepi, sepj, rij, core=core, itype=itype, & 1014 se_int_control=se_int_control, se_int_screen=se_int_screen) 1015 1016 IF (sepi%dorb .OR. sepj%dorb) THEN 1017 ! Compute the contribution from d-orbitals 1018 CALL nucint_d_num(sepi, sepj, rij, core, itype, & 1019 se_int_control=se_int_control, se_int_screen=se_int_screen) 1020 END IF 1021 1022 IF (lgrad) THEN 1023 ! Integrals derivatives for sp shell 1024 CALL dnucint_sp_ana(sepi, sepj, rij, dcore=dcore, itype=itype, & 1025 se_int_control=se_int_control, se_int_screen=se_int_screen) 1026 1027 IF (sepi%dorb .OR. sepj%dorb) THEN 1028 ! Integral derivatives involving d-orbitals 1029 CALL dnucint_d_ana(sepi, sepj, rij, dcore=dcore, itype=itype, & 1030 se_int_control=se_int_control, se_int_screen=se_int_screen) 1031 END IF 1032 END IF 1033 1034 ! Tapering the value of the integrals 1035 IF (lgrad) THEN 1036 DO i = 1, sepi%core_size 1037 dcore(i, 1) = ft*dcore(i, 1) + dft*core(i, 1) 1038 END DO 1039 DO i = 1, sepj%core_size 1040 dcore(i, 2) = ft*dcore(i, 2) + dft*core(i, 2) 1041 END DO 1042 END IF 1043 DO i = 1, sepi%core_size 1044 core(i, 1) = ft*core(i, 1) 1045 END DO 1046 DO i = 1, sepj%core_size 1047 core(i, 2) = ft*core(i, 2) 1048 END DO 1049 1050 ! Debug Procedure.. Check valifity of analytical gradients of nucint 1051 IF (debug_this_module .AND. lgrad) THEN 1052 CALL check_dcore_nucint_ana(sepi, sepj, rij, dcore, itype, se_int_control, se_taper=se_taper) 1053 END IF 1054 END SUBROUTINE dcore_nucint_ana 1055 1056! ************************************************************************************************** 1057!> \brief Calculates the nuclear attraction integrals and derivatives for sp basis 1058!> \param sepi parameters of atom i 1059!> \param sepj parameters of atom j 1060!> \param rij interatomic distance 1061!> \param dssss derivative of (ssss) integral 1062!> derivatives are intended w.r.t. rij 1063!> where ij=1 if the orbitals centred on atom i, =2 if on atom j. 1064!> \param dcore derivative of 4 X 2 array of electron-core attraction integrals 1065!> The storage of the nuclear attraction integrals core(kl/ij) iS 1066!> (SS/)=1, (SP/)=2, (PP/)=3, (P+P+/)=4, (SD/)=5, 1067!> (DP/)=6, (DD/)=7, (D+P+)=8, (D+D+/)=9, (D#D#)=10 1068!> \param itype type of semi_empirical model 1069!> extension to the original routine to compute qm/mm 1070!> integrals 1071!> \param se_int_control input parameters that control the calculation of SE 1072!> integrals (shortrange, R3 residual, screening type) 1073!> \param se_int_screen ... 1074!> \par History 1075!> 04.2007 created [tlaino] 1076!> Teodoro Laino (03.2008) [tlaino] - University of Zurich : new driver 1077!> for computing integrals 1078!> 05.2008 Teodoro Laino [tlaino] - University of Zurich: major rewriting 1079!> \author Teodoro Laino - Zurich University 1080!> \note 1081!> Analytical version - Analytical evaluation of gradients 1082!> Teodoro Laino - Zurich University 04.2007 1083!> routine adapted from mopac7 (repp) 1084!> vector version written by Ernest R. Davidson, Indiana University 1085! ************************************************************************************************** 1086 SUBROUTINE dnucint_sp_ana(sepi, sepj, rij, dssss, dcore, itype, se_int_control, & 1087 se_int_screen) 1088 TYPE(semi_empirical_type), POINTER :: sepi, sepj 1089 REAL(dp), INTENT(IN) :: rij 1090 REAL(dp), INTENT(INOUT), OPTIONAL :: dssss 1091 REAL(dp), DIMENSION(10, 2), INTENT(INOUT), & 1092 OPTIONAL :: dcore 1093 INTEGER, INTENT(IN) :: itype 1094 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 1095 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 1096 1097 CHARACTER(len=*), PARAMETER :: routineN = 'dnucint_sp_ana', routineP = moduleN//':'//routineN 1098 1099 INTEGER :: ij, kl 1100 LOGICAL :: l_core, l_ssss, si, sj 1101 1102 l_core = PRESENT(dcore) 1103 l_ssss = PRESENT(dssss) 1104 IF (.NOT. (l_core .OR. l_ssss)) RETURN 1105 1106 si = (sepi%natorb > 1) 1107 sj = (sepj%natorb > 1) 1108 1109 ij = indexa(1, 1) 1110 IF (l_ssss) THEN 1111 ! Store the value for the derivative of <S S | S S > (Used for computing the core-core interactions) 1112 dssss = d_ijkl_sp(sepi, sepj, ij, ij, 0, 0, 0, 0, -1, rij, CPPint_args) 1113 END IF 1114 1115 IF (l_core) THEN 1116 ! <S S | S S > 1117 kl = indexa(1, 1) 1118 dcore(1, 1) = d_ijkl_sp(sepi, sepj, kl, ij, 0, 0, 0, 0, 2, rij, CPPint_args)*sepj%zeff 1119 IF (si) THEN 1120 ! <S P | S S > 1121 kl = indexa(2, 1) 1122 dcore(2, 1) = d_ijkl_sp(sepi, sepj, kl, ij, 0, 1, 0, 0, 2, rij, CPPint_args)*sepj%zeff 1123 ! <P P | S S > 1124 kl = indexa(2, 2) 1125 dcore(3, 1) = d_ijkl_sp(sepi, sepj, kl, ij, 1, 1, 0, 0, 2, rij, CPPint_args)*sepj%zeff 1126 ! <P+ P+ | S S > 1127 kl = indexa(3, 3) 1128 dcore(4, 1) = d_ijkl_sp(sepi, sepj, kl, ij, 1, 1, 0, 0, 2, rij, CPPint_args)*sepj%zeff 1129 END IF 1130 1131 ! <S S | S S > 1132 kl = indexa(1, 1) 1133 dcore(1, 2) = d_ijkl_sp(sepi, sepj, ij, kl, 0, 0, 0, 0, 1, rij, CPPint_args)*sepi%zeff 1134 IF (sj) THEN 1135 ! <S S | S P > 1136 kl = indexa(2, 1) 1137 dcore(2, 2) = d_ijkl_sp(sepi, sepj, ij, kl, 0, 0, 0, 1, 1, rij, CPPint_args)*sepi%zeff 1138 ! <S S | P P > 1139 kl = indexa(2, 2) 1140 dcore(3, 2) = d_ijkl_sp(sepi, sepj, ij, kl, 0, 0, 1, 1, 1, rij, CPPint_args)*sepi%zeff 1141 ! <S S | P+ P+ > 1142 kl = indexa(3, 3) 1143 dcore(4, 2) = d_ijkl_sp(sepi, sepj, ij, kl, 0, 0, 1, 1, 1, rij, CPPint_args)*sepi%zeff 1144 END IF 1145 END IF 1146 END SUBROUTINE dnucint_sp_ana 1147 1148! ************************************************************************************************** 1149!> \brief Calculates the analytical derivative of the nuclear attraction 1150!> integrals involving d orbitals 1151!> \param sepi parameters of atom i 1152!> \param sepj parameters of atom j 1153!> \param rij interatomic distance 1154!> \param dcore 4 X 2 array of electron-core attraction integrals 1155!> The storage of the nuclear attraction integrals core(kl/ij) iS 1156!> (SS/)=1, (SP/)=2, (PP/)=3, (P+P+/)=4, (SD/)=5, 1157!> (DP/)=6, (DD/)=7, (D+P+)=8, (D+D+/)=9, (D#D#)=10 1158!> 1159!> where ij=1 if the orbitals centred on atom i, =2 if on atom j. 1160!> \param itype type of semi_empirical model 1161!> extension to the original routine to compute qm/mm 1162!> integrals 1163!> \param se_int_control input parameters that control the calculation of SE 1164!> integrals (shortrange, R3 residual, screening type) 1165!> \param se_int_screen ... 1166!> \author 1167!> Teodoro Laino (05.2008) [tlaino] - University of Zurich: created 1168! ************************************************************************************************** 1169 SUBROUTINE dnucint_d_ana(sepi, sepj, rij, dcore, itype, se_int_control, & 1170 se_int_screen) 1171 TYPE(semi_empirical_type), POINTER :: sepi, sepj 1172 REAL(dp), INTENT(IN) :: rij 1173 REAL(dp), DIMENSION(10, 2), INTENT(INOUT) :: dcore 1174 INTEGER, INTENT(IN) :: itype 1175 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 1176 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 1177 1178 CHARACTER(len=*), PARAMETER :: routineN = 'dnucint_d_ana', routineP = moduleN//':'//routineN 1179 1180 INTEGER :: ij, kl 1181 1182! Check if d-orbitals are present 1183 1184 IF (sepi%dorb .OR. sepj%dorb) THEN 1185 ij = indexa(1, 1) 1186 IF (sepj%dorb) THEN 1187 ! <S S | D S> 1188 kl = indexa(5, 1) 1189 dcore(5, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 0, 1, rij, CPPint_args)*sepi%zeff 1190 ! <S S | D P > 1191 kl = indexa(5, 2) 1192 dcore(6, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 1, 1, rij, CPPint_args)*sepi%zeff 1193 ! <S S | D D > 1194 kl = indexa(5, 5) 1195 dcore(7, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, CPPint_args)*sepi%zeff 1196 ! <S S | D+P+> 1197 kl = indexa(6, 3) 1198 dcore(8, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 1, 1, rij, CPPint_args)*sepi%zeff 1199 ! <S S | D+D+> 1200 kl = indexa(6, 6) 1201 dcore(9, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, CPPint_args)*sepi%zeff 1202 ! <S S | D#D#> 1203 kl = indexa(8, 8) 1204 dcore(10, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, CPPint_args)*sepi%zeff 1205 END IF 1206 IF (sepi%dorb) THEN 1207 ! <D S | S S> 1208 kl = indexa(5, 1) 1209 dcore(5, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 0, 0, 0, 2, rij, CPPint_args)*sepj%zeff 1210 ! <D P | S S > 1211 kl = indexa(5, 2) 1212 dcore(6, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 1, 0, 0, 2, rij, CPPint_args)*sepj%zeff 1213 ! <D D | S S > 1214 kl = indexa(5, 5) 1215 dcore(7, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, CPPint_args)*sepj%zeff 1216 ! <D+P+| S S > 1217 kl = indexa(6, 3) 1218 dcore(8, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 1, 0, 0, 2, rij, CPPint_args)*sepj%zeff 1219 ! <D+D+| S S > 1220 kl = indexa(6, 6) 1221 dcore(9, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, CPPint_args)*sepj%zeff 1222 ! <D#D#| S S > 1223 kl = indexa(8, 8) 1224 dcore(10, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, CPPint_args)*sepj%zeff 1225 END IF 1226 END IF 1227 END SUBROUTINE dnucint_d_ana 1228 1229! ************************************************************************************************** 1230!> \brief calculates the derivative of the two-particle interactions 1231!> \param sepi Atomic parameters of first atom 1232!> \param sepj Atomic parameters of second atom 1233!> \param rijv Coordinate vector i -> j 1234!> \param w Array of two-electron repulsion integrals. 1235!> \param dw ... 1236!> \param se_int_control ... 1237!> \param se_taper ... 1238!> \par History 1239!> 04.2007 created [tlaino] 1240!> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver 1241!> for computing integrals 1242!> \author Teodoro Laino - Zurich University 1243!> \note 1244!> Analytical version - Analytical evaluation of gradients 1245!> Teodoro Laino - Zurich University 04.2007 1246!> routine adapted from mopac7 (repp) 1247!> vector version written by Ernest R. Davidson, Indiana University 1248! ************************************************************************************************** 1249 RECURSIVE SUBROUTINE rotint_ana(sepi, sepj, rijv, w, dw, se_int_control, se_taper) 1250 TYPE(semi_empirical_type), POINTER :: sepi, sepj 1251 REAL(dp), DIMENSION(3), INTENT(IN) :: rijv 1252 REAL(dp), DIMENSION(2025), INTENT(OUT), OPTIONAL :: w 1253 REAL(dp), DIMENSION(3, 2025), INTENT(OUT), & 1254 OPTIONAL :: dw 1255 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 1256 TYPE(se_taper_type), POINTER :: se_taper 1257 1258 CHARACTER(len=*), PARAMETER :: routineN = 'rotint_ana', routineP = moduleN//':'//routineN 1259 1260 INTEGER :: i, i1, ii, ij, ij1, iminus, istep, & 1261 iw_loc, j, j1, jj, k, kk, kl, l, & 1262 limij, limkl, mm 1263 LOGICAL :: invert, l_w, lgrad 1264 LOGICAL, DIMENSION(45, 45) :: logv, logv_d 1265 REAL(dp) :: rij, xtmp 1266 REAL(dp), DIMENSION(3) :: drij 1267 REAL(KIND=dp) :: cc, cc_d(3), wrepp, wrepp_d(3) 1268 REAL(KIND=dp), DIMENSION(2025) :: ww 1269 REAL(KIND=dp), DIMENSION(3, 2025) :: ww_d 1270 REAL(KIND=dp), DIMENSION(3, 45, 45) :: v_d 1271 REAL(KIND=dp), DIMENSION(45, 45) :: v 1272 REAL(KIND=dp), DIMENSION(491) :: rep, rep_d 1273 TYPE(rotmat_type), POINTER :: ij_matrix 1274 1275 NULLIFY (ij_matrix) 1276 l_w = PRESENT(w) 1277 lgrad = PRESENT(dw) 1278 IF (.NOT. (l_w .OR. lgrad)) RETURN 1279 1280 rij = DOT_PRODUCT(rijv, rijv) 1281 IF (rij > rij_threshold) THEN 1282 ! The repulsion integrals over molecular frame (w) are stored in the 1283 ! order in which they will later be used. ie. (i,j/k,l) where 1284 ! j.le.i and l.le.k and l varies most rapidly and i least 1285 ! rapidly. (anti-normal computer storage) 1286 rij = SQRT(rij) 1287 1288 ! Create the rotation matrix 1289 CALL rotmat_create(ij_matrix) 1290 CALL rotmat(sepi, sepj, rijv, rij, ij_matrix, do_derivatives=lgrad, do_invert=invert) 1291 1292 ! Compute integrals in diatomic frame as well their derivatives (if requested) 1293 CALL dterep_ana(sepi, sepj, rij, rep, rep_d, se_taper, se_int_control, lgrad=lgrad) 1294 1295 IF (lgrad) THEN 1296 drij(1) = rijv(1)/rij 1297 drij(2) = rijv(2)/rij 1298 drij(3) = rijv(3)/rij 1299 ! Possibly Invert Frame 1300 IF (invert) THEN 1301 xtmp = drij(3) 1302 drij(3) = drij(1) 1303 drij(1) = xtmp 1304 END IF 1305 END IF 1306 1307 ii = sepi%natorb 1308 kk = sepj%natorb 1309 ! First step in rotation of integrals 1310 CALL rot_2el_2c_first(sepi, sepj, rijv, se_int_control, se_taper, invert, ii, kk, rep, logv, ij_matrix, & 1311 v, lgrad, rep_d, v_d, logv_d, drij) 1312 1313 ! Integrals if requested 1314 IF (l_w) THEN 1315 ! Rotate Integrals 1316 IF (ii*kk > 0) THEN 1317 limij = sepi%atm_int_size 1318 limkl = sepj%atm_int_size 1319 istep = limkl*limij 1320 DO i1 = 1, istep 1321 ww(i1) = 0.0_dp 1322 END DO 1323 ! Second step in rotation of integrals 1324 DO i1 = 1, ii 1325 DO j1 = 1, i1 1326 ij = indexa(i1, j1) 1327 jj = indexb(i1, j1) 1328 mm = int2c_type(jj) 1329 DO k = 1, kk 1330 DO l = 1, k 1331 kl = indexb(k, l) 1332 IF (logv(ij, kl)) THEN 1333 wrepp = v(ij, kl) 1334 SELECT CASE (mm) 1335 CASE (1) 1336 ! (SS/) 1337 i = 1 1338 j = 1 1339 iw_loc = (indexb(i, j) - 1)*limkl + kl 1340 ww(iw_loc) = wrepp 1341 CASE (2) 1342 ! (SP/) 1343 j = 1 1344 DO i = 1, 3 1345 iw_loc = (indexb(i + 1, j) - 1)*limkl + kl 1346 ww(iw_loc) = ww(iw_loc) + ij_matrix%sp(i1 - 1, i)*wrepp 1347 END DO 1348 CASE (3) 1349 ! (PP/) 1350 DO i = 1, 3 1351 cc = ij_matrix%pp(i, i1 - 1, j1 - 1) 1352 iw_loc = (indexb(i + 1, i + 1) - 1)*limkl + kl 1353 ww(iw_loc) = ww(iw_loc) + cc*wrepp 1354 iminus = i - 1 1355 IF (iminus /= 0) THEN 1356 DO j = 1, iminus 1357 cc = ij_matrix%pp(1 + i + j, i1 - 1, j1 - 1) 1358 iw_loc = (indexb(i + 1, j + 1) - 1)*limkl + kl 1359 ww(iw_loc) = ww(iw_loc) + cc*wrepp 1360 END DO 1361 END IF 1362 END DO 1363 CASE (4) 1364 ! (SD/) 1365 j = 1 1366 DO i = 1, 5 1367 iw_loc = (indexb(i + 4, j) - 1)*limkl + kl 1368 ww(iw_loc) = ww(iw_loc) + ij_matrix%sd(i1 - 4, i)*wrepp 1369 END DO 1370 CASE (5) 1371 ! (DP/) 1372 DO i = 1, 5 1373 DO j = 1, 3 1374 iw_loc = (indexb(i + 4, j + 1) - 1)*limkl + kl 1375 ij1 = 3*(i - 1) + j 1376 ww(iw_loc) = ww(iw_loc) + ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp 1377 END DO 1378 END DO 1379 CASE (6) 1380 ! (DD/) 1381 DO i = 1, 5 1382 cc = ij_matrix%dd(i, i1 - 4, j1 - 4) 1383 iw_loc = (indexb(i + 4, i + 4) - 1)*limkl + kl 1384 ww(iw_loc) = ww(iw_loc) + cc*wrepp 1385 iminus = i - 1 1386 IF (iminus /= 0) THEN 1387 DO j = 1, iminus 1388 ij1 = inddd(i, j) 1389 cc = ij_matrix%dd(ij1, i1 - 4, j1 - 4) 1390 iw_loc = (indexb(i + 4, j + 4) - 1)*limkl + kl 1391 ww(iw_loc) = ww(iw_loc) + cc*wrepp 1392 END DO 1393 END IF 1394 END DO 1395 END SELECT 1396 END IF 1397 END DO 1398 END DO 1399 END DO 1400 END DO 1401 ! Store two electron integrals in the triangular format 1402 CALL store_2el_2c_diag(limij, limkl, ww(1:istep), w) 1403 IF (invert) CALL invert_integral(sepi, sepj, int2el=w) 1404 END IF 1405 1406 IF (debug_this_module) THEN 1407 ! Check value of integrals 1408 CALL check_rotint_ana(sepi, sepj, rijv, w, se_int_control=se_int_control, se_taper=se_taper) 1409 END IF 1410 END IF 1411 1412 ! Gradients if requested 1413 IF (lgrad) THEN 1414 ! Rotate Integrals derivatives 1415 IF (ii*kk > 0) THEN 1416 limij = sepi%atm_int_size 1417 limkl = sepj%atm_int_size 1418 istep = limkl*limij 1419 DO i1 = 1, istep 1420 ww_d(1, i1) = 0.0_dp 1421 ww_d(2, i1) = 0.0_dp 1422 ww_d(3, i1) = 0.0_dp 1423 END DO 1424 1425 ! Second step in rotation of integrals 1426 DO i1 = 1, ii 1427 DO j1 = 1, i1 1428 ij = indexa(i1, j1) 1429 jj = indexb(i1, j1) 1430 mm = int2c_type(jj) 1431 DO k = 1, kk 1432 DO l = 1, k 1433 kl = indexb(k, l) 1434 IF (logv_d(ij, kl)) THEN 1435 wrepp_d(1) = v_d(1, ij, kl) 1436 wrepp_d(2) = v_d(2, ij, kl) 1437 wrepp_d(3) = v_d(3, ij, kl) 1438 wrepp = v(ij, kl) 1439 SELECT CASE (mm) 1440 CASE (1) 1441 ! (SS/) 1442 i = 1 1443 j = 1 1444 iw_loc = (indexb(i, j) - 1)*limkl + kl 1445 ww_d(1, iw_loc) = wrepp_d(1) 1446 ww_d(2, iw_loc) = wrepp_d(2) 1447 ww_d(3, iw_loc) = wrepp_d(3) 1448 CASE (2) 1449 ! (SP/) 1450 j = 1 1451 DO i = 1, 3 1452 iw_loc = (indexb(i + 1, j) - 1)*limkl + kl 1453 ww_d(1, iw_loc) = ww_d(1, iw_loc) + ij_matrix%sp_d(1, i1 - 1, i)*wrepp + & 1454 ij_matrix%sp(i1 - 1, i)*wrepp_d(1) 1455 1456 ww_d(2, iw_loc) = ww_d(2, iw_loc) + ij_matrix%sp_d(2, i1 - 1, i)*wrepp + & 1457 ij_matrix%sp(i1 - 1, i)*wrepp_d(2) 1458 1459 ww_d(3, iw_loc) = ww_d(3, iw_loc) + ij_matrix%sp_d(3, i1 - 1, i)*wrepp + & 1460 ij_matrix%sp(i1 - 1, i)*wrepp_d(3) 1461 END DO 1462 CASE (3) 1463 ! (PP/) 1464 DO i = 1, 3 1465 cc = ij_matrix%pp(i, i1 - 1, j1 - 1) 1466 cc_d(1) = ij_matrix%pp_d(1, i, i1 - 1, j1 - 1) 1467 cc_d(2) = ij_matrix%pp_d(2, i, i1 - 1, j1 - 1) 1468 cc_d(3) = ij_matrix%pp_d(3, i, i1 - 1, j1 - 1) 1469 iw_loc = (indexb(i + 1, i + 1) - 1)*limkl + kl 1470 ww_d(1, iw_loc) = ww_d(1, iw_loc) + cc_d(1)*wrepp + cc*wrepp_d(1) 1471 ww_d(2, iw_loc) = ww_d(2, iw_loc) + cc_d(2)*wrepp + cc*wrepp_d(2) 1472 ww_d(3, iw_loc) = ww_d(3, iw_loc) + cc_d(3)*wrepp + cc*wrepp_d(3) 1473 iminus = i - 1 1474 IF (iminus /= 0) THEN 1475 DO j = 1, iminus 1476 cc = ij_matrix%pp(1 + i + j, i1 - 1, j1 - 1) 1477 cc_d(1) = ij_matrix%pp_d(1, 1 + i + j, i1 - 1, j1 - 1) 1478 cc_d(2) = ij_matrix%pp_d(2, 1 + i + j, i1 - 1, j1 - 1) 1479 cc_d(3) = ij_matrix%pp_d(3, 1 + i + j, i1 - 1, j1 - 1) 1480 iw_loc = (indexb(i + 1, j + 1) - 1)*limkl + kl 1481 ww_d(1, iw_loc) = ww_d(1, iw_loc) + cc_d(1)*wrepp + cc*wrepp_d(1) 1482 ww_d(2, iw_loc) = ww_d(2, iw_loc) + cc_d(2)*wrepp + cc*wrepp_d(2) 1483 ww_d(3, iw_loc) = ww_d(3, iw_loc) + cc_d(3)*wrepp + cc*wrepp_d(3) 1484 END DO 1485 END IF 1486 END DO 1487 CASE (4) 1488 ! (SD/) 1489 j = 1 1490 DO i = 1, 5 1491 iw_loc = (indexb(i + 4, j) - 1)*limkl + kl 1492 ww_d(1, iw_loc) = ww_d(1, iw_loc) + ij_matrix%sd_d(1, i1 - 4, i)*wrepp + & 1493 ij_matrix%sd(i1 - 4, i)*wrepp_d(1) 1494 1495 ww_d(2, iw_loc) = ww_d(2, iw_loc) + ij_matrix%sd_d(2, i1 - 4, i)*wrepp + & 1496 ij_matrix%sd(i1 - 4, i)*wrepp_d(2) 1497 1498 ww_d(3, iw_loc) = ww_d(3, iw_loc) + ij_matrix%sd_d(3, i1 - 4, i)*wrepp + & 1499 ij_matrix%sd(i1 - 4, i)*wrepp_d(3) 1500 END DO 1501 CASE (5) 1502 ! (DP/) 1503 DO i = 1, 5 1504 DO j = 1, 3 1505 iw_loc = (indexb(i + 4, j + 1) - 1)*limkl + kl 1506 ij1 = 3*(i - 1) + j 1507 ww_d(1, iw_loc) = ww_d(1, iw_loc) + ij_matrix%pd_d(1, ij1, i1 - 4, j1 - 1)*wrepp + & 1508 ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp_d(1) 1509 1510 ww_d(2, iw_loc) = ww_d(2, iw_loc) + ij_matrix%pd_d(2, ij1, i1 - 4, j1 - 1)*wrepp + & 1511 ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp_d(2) 1512 1513 ww_d(3, iw_loc) = ww_d(3, iw_loc) + ij_matrix%pd_d(3, ij1, i1 - 4, j1 - 1)*wrepp + & 1514 ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp_d(3) 1515 END DO 1516 END DO 1517 CASE (6) 1518 ! (DD/) 1519 DO i = 1, 5 1520 cc = ij_matrix%dd(i, i1 - 4, j1 - 4) 1521 cc_d = ij_matrix%dd_d(:, i, i1 - 4, j1 - 4) 1522 iw_loc = (indexb(i + 4, i + 4) - 1)*limkl + kl 1523 ww_d(1, iw_loc) = ww_d(1, iw_loc) + cc_d(1)*wrepp + cc*wrepp_d(1) 1524 ww_d(2, iw_loc) = ww_d(2, iw_loc) + cc_d(2)*wrepp + cc*wrepp_d(2) 1525 ww_d(3, iw_loc) = ww_d(3, iw_loc) + cc_d(3)*wrepp + cc*wrepp_d(3) 1526 iminus = i - 1 1527 IF (iminus /= 0) THEN 1528 DO j = 1, iminus 1529 ij1 = inddd(i, j) 1530 cc = ij_matrix%dd(ij1, i1 - 4, j1 - 4) 1531 cc_d(1) = ij_matrix%dd_d(1, ij1, i1 - 4, j1 - 4) 1532 cc_d(2) = ij_matrix%dd_d(2, ij1, i1 - 4, j1 - 4) 1533 cc_d(3) = ij_matrix%dd_d(3, ij1, i1 - 4, j1 - 4) 1534 iw_loc = (indexb(i + 4, j + 4) - 1)*limkl + kl 1535 ww_d(1, iw_loc) = ww_d(1, iw_loc) + cc_d(1)*wrepp + cc*wrepp_d(1) 1536 ww_d(2, iw_loc) = ww_d(2, iw_loc) + cc_d(2)*wrepp + cc*wrepp_d(2) 1537 ww_d(3, iw_loc) = ww_d(3, iw_loc) + cc_d(3)*wrepp + cc*wrepp_d(3) 1538 END DO 1539 END IF 1540 END DO 1541 END SELECT 1542 END IF 1543 END DO 1544 END DO 1545 END DO 1546 END DO 1547 ! Store two electron integrals in the triangular format 1548 CALL store_2el_2c_diag(limij, limkl, ww_dx=ww_d(1, 1:istep), ww_dy=ww_d(2, 1:istep), ww_dz=ww_d(3, 1:istep), & 1549 dw=dw) 1550 IF (invert) CALL invert_derivative(sepi, sepj, dint2el=dw) 1551 END IF 1552 1553 IF (debug_this_module) THEN 1554 ! Check derivatives 1555 CALL check_rotint_ana(sepi, sepj, rijv, dw=dw, se_int_control=se_int_control, se_taper=se_taper) 1556 END IF 1557 END IF 1558 CALL rotmat_release(ij_matrix) 1559 ENDIF 1560 END SUBROUTINE rotint_ana 1561 1562! ************************************************************************************************** 1563!> \brief Calculates the derivative and the value of two-electron repulsion 1564!> integrals and the nuclear attraction integrals w.r.t. |r| 1565!> \param sepi parameters of atom i 1566!> \param sepj parameters of atom j 1567!> \param rij interatomic distance 1568!> \param rep rray of two-electron repulsion integrals 1569!> \param rep_d array of two-electron repulsion integrals derivatives 1570!> \param se_taper ... 1571!> \param se_int_control input parameters that control the calculation of SE 1572!> integrals (shortrange, R3 residual, screening type) 1573!> \param lgrad ... 1574!> \par History 1575!> 03.2008 created [tlaino] 1576!> \author Teodoro Laino [tlaino] - Zurich University 1577! ************************************************************************************************** 1578 RECURSIVE SUBROUTINE dterep_ana(sepi, sepj, rij, rep, rep_d, se_taper, & 1579 se_int_control, lgrad) 1580 TYPE(semi_empirical_type), POINTER :: sepi, sepj 1581 REAL(KIND=dp), INTENT(IN) :: rij 1582 REAL(KIND=dp), DIMENSION(491), INTENT(OUT) :: rep, rep_d 1583 TYPE(se_taper_type), POINTER :: se_taper 1584 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 1585 LOGICAL, INTENT(IN) :: lgrad 1586 1587 CHARACTER(len=*), PARAMETER :: routineN = 'dterep_ana', routineP = moduleN//':'//routineN 1588 1589 INTEGER :: i, ij, j, k, kl, l, lasti, lastj, li, & 1590 lj, lk, ll, numb 1591 REAL(KIND=dp) :: dft, ft, ft1 1592 TYPE(se_int_screen_type) :: se_int_screen 1593 1594! Compute the tapering function and its derivatives 1595 1596 ft = taper_eval(se_taper%taper, rij) 1597 dft = 0.0_dp 1598 ft1 = ft 1599 IF (lgrad) THEN 1600 ft1 = 1.0_dp 1601 dft = dtaper_eval(se_taper%taper, rij) 1602 END IF 1603 ! Evaluate additional taper function for dumped integrals 1604 IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN 1605 se_int_screen%ft = taper_eval(se_taper%taper_add, rij) 1606 IF (lgrad) & 1607 se_int_screen%dft = dtaper_eval(se_taper%taper_add, rij) 1608 END IF 1609 1610 ! Integral Values for sp shells only 1611 CALL terep_sp_num(sepi, sepj, rij, rep, se_int_control=se_int_control, & 1612 se_int_screen=se_int_screen, ft=ft1) 1613 1614 IF (sepi%dorb .OR. sepj%dorb) THEN 1615 ! Compute the contribution from d-orbitals 1616 CALL terep_d_num(sepi, sepj, rij, rep, se_int_control=se_int_control, & 1617 se_int_screen=se_int_screen, ft=ft1) 1618 END IF 1619 1620 IF (lgrad) THEN 1621 ! Integral Derivatives 1622 CALL dterep_sp_ana(sepi, sepj, rij, rep_d, rep, se_int_control, & 1623 se_int_screen, ft, dft) 1624 1625 IF (sepi%dorb .OR. sepj%dorb) THEN 1626 ! Compute the derivatives from d-orbitals 1627 CALL dterep_d_ana(sepi, sepj, rij, rep_d, rep, se_int_control, & 1628 se_int_screen, ft, dft) 1629 END IF 1630 1631 ! Tapering Integral values 1632 lasti = sepi%natorb 1633 lastj = sepj%natorb 1634 DO i = 1, lasti 1635 li = l_index(i) 1636 DO j = 1, i 1637 lj = l_index(j) 1638 ij = indexa(i, j) 1639 DO k = 1, lastj 1640 lk = l_index(k) 1641 DO l = 1, k 1642 ll = l_index(l) 1643 kl = indexa(k, l) 1644 numb = ijkl_ind(ij, kl) 1645 IF (numb > 0) rep(numb) = rep(numb)*ft 1646 END DO 1647 END DO 1648 END DO 1649 END DO 1650 END IF 1651 1652 ! Possibly debug 2el 2cent integrals and derivatives 1653 IF (debug_this_module) THEN 1654 CALL check_dterep_ana(sepi, sepj, rij, rep, rep_d, se_int_control, se_taper=se_taper, & 1655 lgrad=lgrad) 1656 END IF 1657 END SUBROUTINE dterep_ana 1658 1659! ************************************************************************************************** 1660!> \brief Calculates the derivative and the value of two-electron repulsion 1661!> integrals and the nuclear attraction integrals w.r.t. |r| - sp shells only 1662!> \param sepi parameters of atom i 1663!> \param sepj parameters of atom j 1664!> \param rij interatomic distance 1665!> \param drep array of derivatives of two-electron repulsion integrals 1666!> \param rep array of two-electron repulsion integrals 1667!> \param se_int_control input parameters that control the calculation of SE 1668!> integrals (shortrange, R3 residual, screening type) 1669!> \param se_int_screen ... 1670!> \param ft ... 1671!> \param dft ... 1672!> \par History 1673!> 04.2007 created [tlaino] 1674!> Teodoro Laino (03.2008) [tlaino] - University of Zurich : new driver 1675!> for computing integrals 1676!> 05.2008 Teodoro Laino [tlaino] - University of Zurich: major rewriting 1677!> \author Teodoro Laino - Zurich University 1678!> \note 1679!> Analytical version - Analytical evaluation of gradients 1680!> Teodoro Laino - Zurich University 04.2007 1681!> routine adapted from mopac7 (repp) 1682!> vector version written by Ernest R. Davidson, Indiana University 1683! ************************************************************************************************** 1684 SUBROUTINE dterep_sp_ana(sepi, sepj, rij, drep, rep, se_int_control, & 1685 se_int_screen, ft, dft) 1686 TYPE(semi_empirical_type), POINTER :: sepi, sepj 1687 REAL(dp), INTENT(IN) :: rij 1688 REAL(dp), DIMENSION(491), INTENT(OUT) :: drep 1689 REAL(dp), DIMENSION(491), INTENT(IN) :: rep 1690 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 1691 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 1692 REAL(dp), INTENT(IN) :: ft, dft 1693 1694 CHARACTER(len=*), PARAMETER :: routineN = 'dterep_sp_ana', routineP = moduleN//':'//routineN 1695 1696 INTEGER :: i, ij, j, k, kl, l, lasti, lastj, li, & 1697 lj, lk, ll, nold, numb 1698 REAL(KIND=dp) :: tmp 1699 1700 lasti = sepi%natorb 1701 lastj = sepj%natorb 1702 DO i = 1, MIN(lasti, 4) 1703 li = l_index(i) 1704 DO j = 1, i 1705 lj = l_index(j) 1706 ij = indexa(i, j) 1707 DO k = 1, MIN(lastj, 4) 1708 lk = l_index(k) 1709 DO l = 1, k 1710 ll = l_index(l) 1711 kl = indexa(k, l) 1712 1713 numb = ijkl_ind(ij, kl) 1714 IF (numb > 0) THEN 1715 nold = ijkl_sym(numb) 1716 IF (nold > 0) THEN 1717 drep(numb) = drep(nold) 1718 ELSE IF (nold < 0) THEN 1719 drep(numb) = -drep(-nold) 1720 ELSE IF (nold == 0) THEN 1721 tmp = d_ijkl_sp(sepi, sepj, ij, kl, li, lj, lk, ll, 0, rij, & 1722 se_int_control, se_int_screen, do_method_undef) 1723 drep(numb) = dft*rep(numb) + ft*tmp 1724 END IF 1725 END IF 1726 END DO 1727 END DO 1728 END DO 1729 END DO 1730 END SUBROUTINE dterep_sp_ana 1731 1732! ************************************************************************************************** 1733!> \brief Calculates the derivatives of the two-electron repulsion integrals - d shell only 1734!> \param sepi parameters of atom i 1735!> \param sepj parameters of atom j 1736!> \param rij interatomic distance 1737!> \param drep ... 1738!> \param rep array of two-electron repulsion integrals 1739!> \param se_int_control input parameters that control the calculation of 1740!> integrals (shortrange, R3 residual, screening type) 1741!> \param se_int_screen ... 1742!> \param ft ... 1743!> \param dft ... 1744!> \par History 1745!> Teodoro Laino (05.2008) [tlaino] - University of Zurich : new driver 1746!> for computing integral derivatives for d-orbitals 1747! ************************************************************************************************** 1748 SUBROUTINE dterep_d_ana(sepi, sepj, rij, drep, rep, se_int_control, & 1749 se_int_screen, ft, dft) 1750 TYPE(semi_empirical_type), POINTER :: sepi, sepj 1751 REAL(dp), INTENT(IN) :: rij 1752 REAL(dp), DIMENSION(491), INTENT(INOUT) :: drep 1753 REAL(dp), DIMENSION(491), INTENT(IN) :: rep 1754 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 1755 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 1756 REAL(dp), INTENT(IN) :: ft, dft 1757 1758 CHARACTER(len=*), PARAMETER :: routineN = 'dterep_d_ana', routineP = moduleN//':'//routineN 1759 1760 INTEGER :: i, ij, j, k, kl, l, lasti, lastj, li, & 1761 lj, lk, ll, nold, numb 1762 REAL(KIND=dp) :: tmp 1763 1764 lasti = sepi%natorb 1765 lastj = sepj%natorb 1766 DO i = 1, lasti 1767 li = l_index(i) 1768 DO j = 1, i 1769 lj = l_index(j) 1770 ij = indexa(i, j) 1771 DO k = 1, lastj 1772 lk = l_index(k) 1773 DO l = 1, k 1774 ll = l_index(l) 1775 kl = indexa(k, l) 1776 1777 numb = ijkl_ind(ij, kl) 1778 ! From 1 to 34 we store integrals involving sp shells 1779 IF (numb > 34) THEN 1780 nold = ijkl_sym(numb) 1781 IF (nold > 34) THEN 1782 drep(numb) = drep(nold) 1783 ELSE IF (nold < -34) THEN 1784 drep(numb) = -drep(-nold) 1785 ELSE IF (nold == 0) THEN 1786 tmp = d_ijkl_d(sepi, sepj, ij, kl, li, lj, lk, ll, 0, rij, & 1787 se_int_control, se_int_screen, do_method_undef) 1788 drep(numb) = dft*rep(numb) + ft*tmp 1789 END IF 1790 END IF 1791 END DO 1792 END DO 1793 END DO 1794 END DO 1795 END SUBROUTINE dterep_d_ana 1796 1797END MODULE semi_empirical_int_ana 1798