1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Integrals for semi-empiric methods 8!> \par History 9!> JGH (27.10.2004) : separate routine for nuclear attraction integrals 10!> JGH (13.03.2006) : tapering function 11!> Teodoro Laino (03.2008) [tlaino] - University of Zurich : new driver 12!> for computing integrals 13!> \author JGH (11.10.2004) 14! ************************************************************************************************** 15MODULE semi_empirical_int_num 16 17 USE input_constants, ONLY: do_method_am1,& 18 do_method_pchg,& 19 do_method_pdg,& 20 do_method_pm3,& 21 do_method_pm6,& 22 do_method_pm6fm,& 23 do_method_undef,& 24 do_se_IS_kdso_d 25 USE kinds, ONLY: dp 26 USE multipole_types, ONLY: do_multipole_none 27 USE physcon, ONLY: angstrom,& 28 evolt 29 USE semi_empirical_int_arrays, ONLY: & 30 ijkl_ind, ijkl_sym, inddd, inddp, indexa, indexb, indpp, int2c_type, l_index, rij_threshold 31 USE semi_empirical_int_utils, ONLY: ijkl_d,& 32 ijkl_sp,& 33 rot_2el_2c_first,& 34 rotmat,& 35 store_2el_2c_diag 36 USE semi_empirical_types, ONLY: rotmat_create,& 37 rotmat_release,& 38 rotmat_type,& 39 se_int_control_type,& 40 se_int_screen_type,& 41 se_taper_type,& 42 semi_empirical_type,& 43 setup_se_int_control_type 44 USE taper_types, ONLY: taper_eval 45#include "./base/base_uses.f90" 46 47 IMPLICIT NONE 48 49 PRIVATE 50#include "semi_empirical_int_args.h" 51 52 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. 53 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_int_num' 54 PUBLIC :: rotint_num, rotnuc_num, corecore_num, & 55 drotint_num, drotnuc_num, dcorecore_num, & 56 ssss_nucint_num, core_nucint_num, terep_num, & 57 nucint_sp_num, terep_sp_num, terep_d_num, & 58 nucint_d_num, corecore_el_num, dcorecore_el_num 59 60CONTAINS 61 62! ************************************************************************************************** 63!> \brief Computes the two particle interactions in the lab frame 64!> \param sepi Atomic parameters of first atom 65!> \param sepj Atomic parameters of second atom 66!> \param rijv Coordinate vector i -> j 67!> \param w Array of two-electron repulsion integrals. 68!> \param se_int_control input parameters that control the calculation of SE 69!> integrals (shortrange, R3 residual, screening type) 70!> \param se_taper ... 71!> \note routine adapted from mopac7 (rotate) 72!> written by Ernest R. Davidson, Indiana University. 73!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : major rewriting 74! ************************************************************************************************** 75 SUBROUTINE rotint_num(sepi, sepj, rijv, w, se_int_control, se_taper) 76 TYPE(semi_empirical_type), POINTER :: sepi, sepj 77 REAL(dp), DIMENSION(3), INTENT(IN) :: rijv 78 REAL(dp), DIMENSION(2025), INTENT(OUT) :: w 79 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 80 TYPE(se_taper_type), POINTER :: se_taper 81 82 CHARACTER(len=*), PARAMETER :: routineN = 'rotint_num', routineP = moduleN//':'//routineN 83 84 INTEGER :: i, i1, ii, ij, ij1, iminus, istep, & 85 iw_loc, j, j1, jj, k, kk, kl, l, & 86 limij, limkl, mm 87 LOGICAL, DIMENSION(45, 45) :: logv 88 REAL(dp) :: rij 89 REAL(KIND=dp) :: cc, wrepp 90 REAL(KIND=dp), DIMENSION(2025) :: ww 91 REAL(KIND=dp), DIMENSION(45, 45) :: v 92 REAL(KIND=dp), DIMENSION(491) :: rep 93 TYPE(rotmat_type), POINTER :: ij_matrix 94 95 NULLIFY (ij_matrix) 96 rij = DOT_PRODUCT(rijv, rijv) 97 IF (rij > rij_threshold) THEN 98 ! The repulsion integrals over molecular frame (w) are stored in the 99 ! order in which they will later be used. ie. (i,j/k,l) where 100 ! j.le.i and l.le.k and l varies most rapidly and i least 101 ! rapidly. (anti-normal computer storage) 102 rij = SQRT(rij) 103 104 ! Create the rotation matrix 105 CALL rotmat_create(ij_matrix) 106 CALL rotmat(sepi, sepj, rijv, rij, ij_matrix, do_derivatives=.FALSE.) 107 108 ! Compute Integrals in Diatomic Frame 109 CALL terep_num(sepi, sepj, rij, rep, se_taper=se_taper, se_int_control=se_int_control) 110 111 ! Rotate Integrals 112 ii = sepi%natorb 113 kk = sepj%natorb 114 IF (ii*kk > 0) THEN 115 limij = sepi%atm_int_size 116 limkl = sepj%atm_int_size 117 istep = limkl*limij 118 DO i1 = 1, istep 119 ww(i1) = 0.0_dp 120 END DO 121 122 ! First step in rotation of integrals 123 CALL rot_2el_2c_first(sepi, sepj, rijv, se_int_control, se_taper, .FALSE., ii, kk, rep, & 124 logv, ij_matrix, v, lgrad=.FALSE.) 125 126 ! Second step in rotation of integrals 127 DO i1 = 1, ii 128 DO j1 = 1, i1 129 ij = indexa(i1, j1) 130 jj = indexb(i1, j1) 131 mm = int2c_type(jj) 132 DO k = 1, kk 133 DO l = 1, k 134 kl = indexb(k, l) 135 IF (logv(ij, kl)) THEN 136 wrepp = v(ij, kl) 137 SELECT CASE (mm) 138 CASE (1) 139 ! (SS/) 140 i = 1 141 j = 1 142 iw_loc = (indexb(i, j) - 1)*limkl + kl 143 ww(iw_loc) = wrepp 144 CASE (2) 145 ! (SP/) 146 j = 1 147 DO i = 1, 3 148 iw_loc = (indexb(i + 1, j) - 1)*limkl + kl 149 ww(iw_loc) = ww(iw_loc) + ij_matrix%sp(i1 - 1, i)*wrepp 150 END DO 151 CASE (3) 152 ! (PP/) 153 DO i = 1, 3 154 cc = ij_matrix%pp(i, i1 - 1, j1 - 1) 155 iw_loc = (indexb(i + 1, i + 1) - 1)*limkl + kl 156 ww(iw_loc) = ww(iw_loc) + cc*wrepp 157 iminus = i - 1 158 IF (iminus /= 0) THEN 159 DO j = 1, iminus 160 cc = ij_matrix%pp(1 + i + j, i1 - 1, j1 - 1) 161 iw_loc = (indexb(i + 1, j + 1) - 1)*limkl + kl 162 ww(iw_loc) = ww(iw_loc) + cc*wrepp 163 END DO 164 END IF 165 END DO 166 CASE (4) 167 ! (SD/) 168 j = 1 169 DO i = 1, 5 170 iw_loc = (indexb(i + 4, j) - 1)*limkl + kl 171 ww(iw_loc) = ww(iw_loc) + ij_matrix%sd(i1 - 4, i)*wrepp 172 END DO 173 CASE (5) 174 ! (DP/) 175 DO i = 1, 5 176 DO j = 1, 3 177 iw_loc = (indexb(i + 4, j + 1) - 1)*limkl + kl 178 ij1 = 3*(i - 1) + j 179 ww(iw_loc) = ww(iw_loc) + ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp 180 END DO 181 END DO 182 CASE (6) 183 ! (DD/) 184 DO i = 1, 5 185 cc = ij_matrix%dd(i, i1 - 4, j1 - 4) 186 iw_loc = (indexb(i + 4, i + 4) - 1)*limkl + kl 187 ww(iw_loc) = ww(iw_loc) + cc*wrepp 188 iminus = i - 1 189 IF (iminus /= 0) THEN 190 DO j = 1, iminus 191 ij1 = inddd(i, j) 192 cc = ij_matrix%dd(ij1, i1 - 4, j1 - 4) 193 iw_loc = (indexb(i + 4, j + 4) - 1)*limkl + kl 194 ww(iw_loc) = ww(iw_loc) + cc*wrepp 195 END DO 196 END IF 197 END DO 198 END SELECT 199 END IF 200 END DO 201 END DO 202 END DO 203 END DO 204 END IF 205 CALL rotmat_release(ij_matrix) 206 207 ! Store two electron integrals in the triangular format 208 CALL store_2el_2c_diag(limij, limkl, ww, w) 209 ENDIF 210 END SUBROUTINE rotint_num 211 212! ************************************************************************************************** 213!> \brief Calculates the derivative pf two-electron repulsion integrals and the 214!> nuclear attraction integrals w.r.t. |r| 215!> \param sepi paramters of atom i 216!> \param sepj paramters of atom j 217!> \param rij interatomic distance 218!> \param rep array of two-electron repulsion integrals 219!> \param se_taper ... 220!> \param se_int_control input parameters that control the calculation of SE 221!> integrals (shortrange, R3 residual, screening type) 222!> \par History 223!> 03.2008 created [tlaino] 224!> \author Teodoro Laino [tlaino] - Zurich University 225! ************************************************************************************************** 226 SUBROUTINE terep_num(sepi, sepj, rij, rep, se_taper, se_int_control) 227 TYPE(semi_empirical_type), POINTER :: sepi, sepj 228 REAL(dp), INTENT(IN) :: rij 229 REAL(dp), DIMENSION(491), INTENT(OUT) :: rep 230 TYPE(se_taper_type), POINTER :: se_taper 231 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 232 233 CHARACTER(len=*), PARAMETER :: routineN = 'terep_num', routineP = moduleN//':'//routineN 234 235 REAL(KIND=dp) :: ft 236 TYPE(se_int_screen_type) :: se_int_screen 237 238 ft = taper_eval(se_taper%taper, rij) 239 ! In case of dumped integrals compute an additional taper term 240 IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN 241 se_int_screen%ft = taper_eval(se_taper%taper_add, rij) 242 END IF 243 244 ! Contribution from sp shells 245 CALL terep_sp_num(sepi, sepj, rij, rep, se_int_control, se_int_screen, ft) 246 247 IF (sepi%dorb .OR. sepj%dorb) THEN 248 ! Compute the contribution from d shells 249 CALL terep_d_num(sepi, sepj, rij, rep, se_int_control, se_int_screen, & 250 ft) 251 END IF 252 253 END SUBROUTINE terep_num 254 255! ************************************************************************************************** 256!> \brief Calculates the two-electron repulsion integrals - sp shell only 257!> \param sepi paramters of atom i 258!> \param sepj paramters of atom j 259!> \param rij ... 260!> \param rep array of two-electron repulsion integrals 261!> \param se_int_control input parameters that control the calculation of SE 262!> integrals (shortrange, R3 residual, screening type) 263!> \param se_int_screen contains information for computing the screened 264!> integrals KDSO-D 265!> \param ft ... 266!> \par History 267!> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver 268!> for computing integrals 269! ************************************************************************************************** 270 SUBROUTINE terep_sp_num(sepi, sepj, rij, rep, se_int_control, se_int_screen, & 271 ft) 272 TYPE(semi_empirical_type), POINTER :: sepi, sepj 273 REAL(dp), INTENT(IN) :: rij 274 REAL(dp), DIMENSION(491), INTENT(OUT) :: rep 275 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 276 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 277 REAL(dp), INTENT(IN) :: ft 278 279 CHARACTER(len=*), PARAMETER :: routineN = 'terep_sp_num', routineP = moduleN//':'//routineN 280 281 INTEGER :: i, ij, j, k, kl, l, lasti, lastj, li, & 282 lj, lk, ll, nold, numb 283 REAL(KIND=dp) :: tmp 284 285 lasti = sepi%natorb 286 lastj = sepj%natorb 287 DO i = 1, MIN(lasti, 4) 288 li = l_index(i) 289 DO j = 1, i 290 lj = l_index(j) 291 ij = indexa(i, j) 292 DO k = 1, MIN(lastj, 4) 293 lk = l_index(k) 294 DO l = 1, k 295 ll = l_index(l) 296 kl = indexa(k, l) 297 298 numb = ijkl_ind(ij, kl) 299 IF (numb > 0) THEN 300 nold = ijkl_sym(numb) 301 IF (nold > 0) THEN 302 rep(numb) = rep(nold) 303 ELSE IF (nold < 0) THEN 304 rep(numb) = -rep(-nold) 305 ELSE IF (nold == 0) THEN 306 tmp = ijkl_sp(sepi, sepj, ij, kl, li, lj, lk, ll, 0, rij, & 307 se_int_control, se_int_screen, do_method_undef) & 308 *ft 309 rep(numb) = tmp 310 END IF 311 END IF 312 END DO 313 END DO 314 END DO 315 END DO 316 END SUBROUTINE terep_sp_num 317 318! ************************************************************************************************** 319!> \brief Calculates the two-electron repulsion integrals - d shell only 320!> \param sepi ... 321!> \param sepj ... 322!> \param rij interatomic distance 323!> \param rep array of two-electron repulsion integrals 324!> \param se_int_control input parameters that control the calculation of SE 325!> integrals (shortrange, R3 residual, screening type) 326!> \param se_int_screen contains information for computing the screened 327!> integrals KDSO-D 328!> \param ft ... 329!> \par History 330!> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver 331!> for computing integrals 332! ************************************************************************************************** 333 SUBROUTINE terep_d_num(sepi, sepj, rij, rep, se_int_control, se_int_screen, & 334 ft) 335 TYPE(semi_empirical_type), POINTER :: sepi, sepj 336 REAL(dp), INTENT(IN) :: rij 337 REAL(dp), DIMENSION(491), INTENT(INOUT) :: rep 338 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 339 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 340 REAL(dp), INTENT(IN) :: ft 341 342 CHARACTER(len=*), PARAMETER :: routineN = 'terep_d_num', routineP = moduleN//':'//routineN 343 344 INTEGER :: i, ij, j, k, kl, l, lasti, lastj, li, & 345 lj, lk, ll, nold, numb 346 REAL(KIND=dp) :: tmp 347 348 lasti = sepi%natorb 349 lastj = sepj%natorb 350 DO i = 1, lasti 351 li = l_index(i) 352 DO j = 1, i 353 lj = l_index(j) 354 ij = indexa(i, j) 355 DO k = 1, lastj 356 lk = l_index(k) 357 DO l = 1, k 358 ll = l_index(l) 359 kl = indexa(k, l) 360 361 numb = ijkl_ind(ij, kl) 362 ! From 1 to 34 we store integrals involving sp shells 363 IF (numb > 34) THEN 364 nold = ijkl_sym(numb) 365 IF (nold > 34) THEN 366 rep(numb) = rep(nold) 367 ELSE IF (nold < -34) THEN 368 rep(numb) = -rep(-nold) 369 ELSE IF (nold == 0) THEN 370 tmp = ijkl_d(sepi, sepj, ij, kl, li, lj, lk, ll, 0, rij, & 371 se_int_control, se_int_screen, do_method_undef) & 372 *ft 373 rep(numb) = tmp 374 END IF 375 END IF 376 END DO 377 END DO 378 END DO 379 END DO 380 381 END SUBROUTINE terep_d_num 382 383! ************************************************************************************************** 384!> \brief Computes the two-particle interactions. 385!> \param sepi Atomic parameters of first atom 386!> \param sepj Atomic parameters of second atom 387!> \param rijv Coordinate vector i -> j 388!> \param e1b Array of electron-nuclear attraction integrals, 389!> e1b = Electron on atom ni attracting nucleus of nj. 390!> \param e2a Array of electron-nuclear attraction integrals, 391!> e2a = Electron on atom nj attracting nucleus of ni. 392!> \param itype ... 393!> \param se_int_control ... 394!> \param se_taper ... 395!> \note routine adapted from mopac7 (rotate) 396!> written by Ernest R. Davidson, Indiana University. 397!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : major rewriting 398!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : removed the core-core part 399! ************************************************************************************************** 400 SUBROUTINE rotnuc_num(sepi, sepj, rijv, e1b, e2a, itype, se_int_control, se_taper) 401 TYPE(semi_empirical_type), POINTER :: sepi, sepj 402 REAL(dp), DIMENSION(3), INTENT(IN) :: rijv 403 REAL(dp), DIMENSION(45), INTENT(OUT), OPTIONAL :: e1b, e2a 404 INTEGER, INTENT(IN) :: itype 405 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 406 TYPE(se_taper_type), POINTER :: se_taper 407 408 CHARACTER(len=*), PARAMETER :: routineN = 'rotnuc_num', routineP = moduleN//':'//routineN 409 410 INTEGER :: i, idd, idp, ind1, ind2, ipp, j, & 411 last_orbital(2), m, n 412 LOGICAL :: l_e1b, l_e2a, task(2) 413 REAL(dp) :: rij 414 REAL(dp), DIMENSION(10, 2) :: core 415 REAL(dp), DIMENSION(45) :: tmp 416 TYPE(rotmat_type), POINTER :: ij_matrix 417 418 NULLIFY (ij_matrix) 419 l_e1b = PRESENT(e1b) 420 l_e2a = PRESENT(e2a) 421 rij = DOT_PRODUCT(rijv, rijv) 422 423 IF (rij > rij_threshold) THEN 424 rij = SQRT(rij) 425 ! Create the rotation matrix 426 CALL rotmat_create(ij_matrix) 427 CALL rotmat(sepi, sepj, rijv, rij, ij_matrix, do_derivatives=.FALSE.) 428 429 ! Compute Integrals in Diatomic Frame 430 CALL core_nucint_num(sepi, sepj, rij, core=core, itype=itype, se_taper=se_taper, & 431 se_int_control=se_int_control) 432 433 ! Copy parameters over to arrays for do loop. 434 last_orbital(1) = sepi%natorb 435 last_orbital(2) = sepj%natorb 436 task(1) = l_e1b 437 task(2) = l_e2a 438 DO n = 1, 2 439 IF (.NOT. task(n)) CYCLE 440 DO i = 1, last_orbital(n) 441 ind1 = i - 1 442 DO j = 1, i 443 ind2 = j - 1 444 m = (i*(i - 1))/2 + j 445 ! Perform Rotations ... 446 IF (ind2 == 0) THEN 447 IF (ind1 == 0) THEN 448 ! Type of Integral (SS/) 449 tmp(m) = core(1, n) 450 ELSE IF (ind1 < 4) THEN 451 ! Type of Integral (SP/) 452 tmp(m) = ij_matrix%sp(1, ind1)*core(2, n) 453 ELSE 454 ! Type of Integral (SD/) 455 tmp(m) = ij_matrix%sd(1, ind1 - 3)*core(5, n) 456 END IF 457 ELSE IF (ind2 < 4) THEN 458 IF (ind1 < 4) THEN 459 ! Type of Integral (PP/) 460 ipp = indpp(ind1, ind2) 461 tmp(m) = core(3, n)*ij_matrix%pp(ipp, 1, 1) + & 462 core(4, n)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3)) 463 ELSE 464 ! Type of Integral (PD/) 465 idp = inddp(ind1 - 3, ind2) 466 tmp(m) = core(6, n)*ij_matrix%pd(idp, 1, 1) + & 467 core(8, n)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3)) 468 END IF 469 ELSE 470 ! Type of Integral (DD/) 471 idd = inddd(ind1 - 3, ind2 - 3) 472 tmp(m) = core(7, n)*ij_matrix%dd(idd, 1, 1) + & 473 core(9, n)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + & 474 core(10, n)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5)) 475 END IF 476 END DO 477 END DO 478 IF (n == 1) THEN 479 DO i = 1, sepi%atm_int_size 480 e1b(i) = -tmp(i) 481 END DO 482 END IF 483 IF (n == 2) THEN 484 DO i = 1, sepj%atm_int_size 485 e2a(i) = -tmp(i) 486 END DO 487 END IF 488 END DO 489 CALL rotmat_release(ij_matrix) 490 END IF 491 END SUBROUTINE rotnuc_num 492 493! ************************************************************************************************** 494!> \brief Computes the core-core interactions. 495!> \param sepi Atomic parameters of first atom 496!> \param sepj Atomic parameters of second atom 497!> \param rijv Coordinate vector i -> j 498!> \param enuc nuclear-nuclear repulsion term. 499!> \param itype ... 500!> \param se_int_control input parameters that control the calculation of SE 501!> integrals (shortrange, R3 residual, screening type) 502!> \param se_taper ... 503!> \note routine adapted from mopac7 (rotate) 504!> written by Ernest R. Davidson, Indiana University. 505!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : major rewriting 506!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : splitted from rotnuc 507! ************************************************************************************************** 508 SUBROUTINE corecore_num(sepi, sepj, rijv, enuc, itype, se_int_control, se_taper) 509 TYPE(semi_empirical_type), POINTER :: sepi, sepj 510 REAL(dp), DIMENSION(3), INTENT(IN) :: rijv 511 REAL(dp), INTENT(OUT) :: enuc 512 INTEGER, INTENT(IN) :: itype 513 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 514 TYPE(se_taper_type), POINTER :: se_taper 515 516 CHARACTER(len=*), PARAMETER :: routineN = 'corecore_num', routineP = moduleN//':'//routineN 517 518 INTEGER :: ig, nt 519 REAL(dp) :: aab, alpi, alpj, apdg, ax, dai, daj, & 520 dbi, dbj, pai, paj, pbi, pbj, qcorr, & 521 rij, rija, scale, ssss, ssss_sr, tmp, & 522 xab, zaf, zbf, zz 523 REAL(dp), DIMENSION(4) :: fni1, fni2, fni3, fnj1, fnj2, fnj3 524 TYPE(se_int_control_type) :: se_int_control_off 525 526 rij = DOT_PRODUCT(rijv, rijv) 527 enuc = 0.0_dp 528 IF (rij > rij_threshold) THEN 529 rij = SQRT(rij) 530 CALL setup_se_int_control_type(se_int_control_off, shortrange=.FALSE., do_ewald_r3=.FALSE., & 531 do_ewald_gks=.FALSE., integral_screening=se_int_control%integral_screening, & 532 max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.) 533 CALL ssss_nucint_num(sepi, sepj, rij, ssss=ssss, itype=itype, se_taper=se_taper, & 534 se_int_control=se_int_control_off) 535 ! In case let's compute the short-range part of the (ss|ss) integral 536 IF (se_int_control%shortrange) THEN 537 CALL ssss_nucint_num(sepi, sepj, rij, ssss=ssss_sr, itype=itype, se_taper=se_taper, & 538 se_int_control=se_int_control) 539 ELSE 540 ssss_sr = ssss 541 END IF 542 zz = sepi%zeff*sepj%zeff 543 ! Nuclear-Nuclear electrostatic contribution 544 enuc = zz*ssss_sr 545 ! Method dependent correction terms 546 IF (itype /= do_method_pm6 .AND. itype /= do_method_pm6fm) THEN 547 alpi = sepi%alp 548 alpj = sepj%alp 549 scale = EXP(-alpi*rij) + EXP(-alpj*rij) 550 551 nt = sepi%z + sepj%z 552 IF (nt == 8 .OR. nt == 9) THEN 553 IF (sepi%z == 7 .OR. sepi%z == 8) scale = scale + (angstrom*rij - 1._dp)*EXP(-alpi*rij) 554 IF (sepj%z == 7 .OR. sepj%z == 8) scale = scale + (angstrom*rij - 1._dp)*EXP(-alpj*rij) 555 ENDIF 556 scale = ABS(scale*zz*ssss) 557 zz = zz/rij 558 IF (itype == do_method_am1 .OR. itype == do_method_pm3 .OR. itype == do_method_pdg) THEN 559 IF (itype == do_method_am1 .AND. sepi%z == 5) THEN 560 !special case AM1 Boron 561 SELECT CASE (sepj%z) 562 CASE DEFAULT 563 nt = 1 564 CASE (1) 565 nt = 2 566 CASE (6) 567 nt = 3 568 CASE (9, 17, 35, 53) 569 nt = 4 570 END SELECT 571 fni1(:) = sepi%bfn1(:, nt) 572 fni2(:) = sepi%bfn2(:, nt) 573 fni3(:) = sepi%bfn3(:, nt) 574 ELSE 575 fni1(:) = sepi%fn1(:) 576 fni2(:) = sepi%fn2(:) 577 fni3(:) = sepi%fn3(:) 578 END IF 579 IF (itype == do_method_am1 .AND. sepj%z == 5) THEN 580 !special case AM1 Boron 581 SELECT CASE (sepi%z) 582 CASE DEFAULT 583 nt = 1 584 CASE (1) 585 nt = 2 586 CASE (6) 587 nt = 3 588 CASE (9, 17, 35, 53) 589 nt = 4 590 END SELECT 591 fnj1(:) = sepj%bfn1(:, nt) 592 fnj2(:) = sepj%bfn2(:, nt) 593 fnj3(:) = sepj%bfn3(:, nt) 594 ELSE 595 fnj1(:) = sepj%fn1(:) 596 fnj2(:) = sepj%fn2(:) 597 fnj3(:) = sepj%fn3(:) 598 END IF 599 ! AM1/PM3/PDG correction to nuclear repulsion 600 DO ig = 1, SIZE(fni1) 601 IF (ABS(fni1(ig)) > 0._dp) THEN 602 ax = fni2(ig)*(rij - fni3(ig))**2 603 IF (ax <= 25._dp) THEN 604 scale = scale + zz*fni1(ig)*EXP(-ax) 605 ENDIF 606 ENDIF 607 IF (ABS(fnj1(ig)) > 0._dp) THEN 608 ax = fnj2(ig)*(rij - fnj3(ig))**2 609 IF (ax <= 25._dp) THEN 610 scale = scale + zz*fnj1(ig)*EXP(-ax) 611 ENDIF 612 ENDIF 613 END DO 614 ENDIF 615 IF (itype == do_method_pdg) THEN 616 ! PDDG function 617 zaf = sepi%zeff/nt 618 zbf = sepj%zeff/nt 619 pai = sepi%pre(1) 620 pbi = sepi%pre(2) 621 paj = sepj%pre(1) 622 pbj = sepj%pre(2) 623 dai = sepi%d(1) 624 dbi = sepi%d(2) 625 daj = sepj%d(1) 626 dbj = sepj%d(2) 627 apdg = 10._dp*angstrom**2 628 qcorr = & 629 (zaf*pai + zbf*paj)*EXP(-apdg*(rij - dai - daj)**2) + & 630 (zaf*pai + zbf*pbj)*EXP(-apdg*(rij - dai - dbj)**2) + & 631 (zaf*pbi + zbf*paj)*EXP(-apdg*(rij - dbi - daj)**2) + & 632 (zaf*pbi + zbf*pbj)*EXP(-apdg*(rij - dbi - dbj)**2) 633 ELSEIF (itype == do_method_pchg) THEN 634 qcorr = 0.0_dp 635 scale = 0.0_dp 636 ELSE 637 qcorr = 0.0_dp 638 END IF 639 ELSE 640 rija = rij*angstrom 641 scale = zz*ssss 642 ! PM6 core-core terms 643 xab = sepi%xab(sepj%z) 644 aab = sepi%aab(sepj%z) 645 IF ((sepi%z == 1 .AND. (sepj%z == 6 .OR. sepj%z == 7 .OR. sepj%z == 8)) .OR. & 646 (sepj%z == 1 .AND. (sepi%z == 6 .OR. sepi%z == 7 .OR. sepi%z == 8))) THEN 647 ! Special Case O-H or N-H or C-H 648 scale = scale*(2._dp*xab*EXP(-aab*rija*rija)) 649 ELSEIF (sepi%z == 6 .AND. sepj%z == 6) THEN 650 ! Special Case C-C 651 scale = scale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6)) + 9.28_dp*EXP(-5.98_dp*rija)) 652 ELSEIF ((sepi%z == 8 .AND. sepj%z == 14) .OR. & 653 (sepj%z == 8 .AND. sepi%z == 14)) THEN 654 ! Special Case Si-O 655 scale = scale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6)) - 0.0007_dp*EXP(-(rija - 2.9_dp)**2)) 656 ELSE 657 ! General Case 658 ! Factor of 2 found by experiment 659 scale = scale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6))) 660 END IF 661 ! General correction term a*exp(-b*(rij-c)^2) 662 qcorr = (sepi%a*EXP(-sepi%b*(rij - sepi%c)**2))*sepi%zeff*sepj%zeff/rij + & 663 (sepj%a*EXP(-sepj%b*(rij - sepj%c)**2))*sepi%zeff*sepj%zeff/rij 664 ! Hard core repulsion 665 tmp = (REAL(sepi%z, dp)**(1._dp/3._dp) + REAL(sepj%z, dp)**(1._dp/3._dp)) 666 qcorr = qcorr + 1.e-8_dp/evolt*(tmp/rija)**12 667 END IF 668 enuc = enuc + scale + qcorr 669 END IF 670 END SUBROUTINE corecore_num 671 672! ************************************************************************************************** 673!> \brief Computes the electrostatic core-core interactions only. 674!> \param sepi Atomic parameters of first atom 675!> \param sepj Atomic parameters of second atom 676!> \param rijv Coordinate vector i -> j 677!> \param enuc nuclear-nuclear repulsion term. 678!> \param itype ... 679!> \param se_int_control input parameters that control the calculation of SE 680!> integrals (shortrange, R3 residual, screening type) 681!> \param se_taper ... 682!> \author Teodoro Laino [tlaino] - 05.2009 683! ************************************************************************************************** 684 SUBROUTINE corecore_el_num(sepi, sepj, rijv, enuc, itype, se_int_control, se_taper) 685 TYPE(semi_empirical_type), POINTER :: sepi, sepj 686 REAL(dp), DIMENSION(3), INTENT(IN) :: rijv 687 REAL(dp), INTENT(OUT) :: enuc 688 INTEGER, INTENT(IN) :: itype 689 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 690 TYPE(se_taper_type), POINTER :: se_taper 691 692 CHARACTER(len=*), PARAMETER :: routineN = 'corecore_el_num', & 693 routineP = moduleN//':'//routineN 694 695 REAL(dp) :: rij, ssss, ssss_sr, zz 696 TYPE(se_int_control_type) :: se_int_control_off 697 698 rij = DOT_PRODUCT(rijv, rijv) 699 enuc = 0.0_dp 700 IF (rij > rij_threshold) THEN 701 rij = SQRT(rij) 702 CALL setup_se_int_control_type(se_int_control_off, shortrange=.FALSE., do_ewald_r3=.FALSE., & 703 do_ewald_gks=.FALSE., integral_screening=se_int_control%integral_screening, & 704 max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.) 705 CALL ssss_nucint_num(sepi, sepj, rij, ssss=ssss, itype=itype, se_taper=se_taper, & 706 se_int_control=se_int_control_off) 707 ! In case let's compute the short-range part of the (ss|ss) integral 708 IF (se_int_control%shortrange .OR. se_int_control%pc_coulomb_int) THEN 709 CALL ssss_nucint_num(sepi, sepj, rij, ssss=ssss_sr, itype=itype, se_taper=se_taper, & 710 se_int_control=se_int_control) 711 ELSE 712 ssss_sr = ssss 713 END IF 714 zz = sepi%zeff*sepj%zeff 715 ! Nuclear-Nuclear electrostatic contribution 716 enuc = zz*ssss_sr 717 END IF 718 END SUBROUTINE corecore_el_num 719 720! ************************************************************************************************** 721!> \brief Calculates the SSSS integrals (main driver) 722!> \param sepi paramters of atom i 723!> \param sepj paramters of atom j 724!> \param rij interatomic distance 725!> \param ssss derivative of (ssss) integral 726!> derivatives are intended w.r.t. rij 727!> \param itype type of semi_empirical model 728!> extension to the original routine to compute qm/mm integrals 729!> \param se_taper ... 730!> \param se_int_control input parameters that control the calculation of SE 731!> integrals (shortrange, R3 residual, screening type) 732!> \par History 733!> 03.2008 created [tlaino] 734!> \author Teodoro Laino - Zurich University 735! ************************************************************************************************** 736 SUBROUTINE ssss_nucint_num(sepi, sepj, rij, ssss, itype, se_taper, se_int_control) 737 TYPE(semi_empirical_type), POINTER :: sepi, sepj 738 REAL(dp), INTENT(IN) :: rij 739 REAL(dp), INTENT(OUT) :: ssss 740 INTEGER, INTENT(IN) :: itype 741 TYPE(se_taper_type), POINTER :: se_taper 742 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 743 744 CHARACTER(len=*), PARAMETER :: routineN = 'ssss_nucint_num', & 745 routineP = moduleN//':'//routineN 746 747 REAL(KIND=dp) :: ft 748 TYPE(se_int_screen_type) :: se_int_screen 749 750! Computing Tapering function 751 752 ft = 1.0_dp 753 IF (itype /= do_method_pchg) THEN 754 ft = taper_eval(se_taper%taper, rij) 755 END IF 756 757 ! In case of dumped integrals compute an additional taper term 758 IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN 759 se_int_screen%ft = 1.0_dp 760 IF (itype /= do_method_pchg) THEN 761 se_int_screen%ft = taper_eval(se_taper%taper_add, rij) 762 END IF 763 END IF 764 765 ! Contribution from the sp shells 766 CALL nucint_sp_num(sepi, sepj, rij, ssss=ssss, itype=itype, & 767 se_int_control=se_int_control, se_int_screen=se_int_screen) 768 769 ! Tapering the integrals 770 ssss = ft*ssss 771 772 END SUBROUTINE ssss_nucint_num 773 774! ************************************************************************************************** 775!> \brief Calculates the nuclear attraction integrals CORE (main driver) 776!> \param sepi paramters of atom i 777!> \param sepj paramters of atom j 778!> \param rij interatomic distance 779!> \param core derivative of 4 X 2 array of electron-core attraction integrals 780!> derivatives are intended w.r.t. rij 781!> \param itype type of semi_empirical model 782!> extension to the original routine to compute qm/mm integrals 783!> \param se_taper ... 784!> \param se_int_control input parameters that control the calculation of SE 785!> integrals (shortrange, R3 residual, screening type) 786!> \par History 787!> 03.2008 created [tlaino] 788!> \author Teodoro Laino - Zurich University 789! ************************************************************************************************** 790 SUBROUTINE core_nucint_num(sepi, sepj, rij, core, itype, se_taper, se_int_control) 791 TYPE(semi_empirical_type), POINTER :: sepi, sepj 792 REAL(dp), INTENT(IN) :: rij 793 REAL(dp), DIMENSION(10, 2), INTENT(OUT) :: core 794 INTEGER, INTENT(IN) :: itype 795 TYPE(se_taper_type), POINTER :: se_taper 796 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 797 798 CHARACTER(len=*), PARAMETER :: routineN = 'core_nucint_num', & 799 routineP = moduleN//':'//routineN 800 801 INTEGER :: i 802 REAL(KIND=dp) :: ft 803 TYPE(se_int_screen_type) :: se_int_screen 804 805! Computing the Tapering function 806 807 ft = 1.0_dp 808 IF (itype /= do_method_pchg) THEN 809 ft = taper_eval(se_taper%taper, rij) 810 END IF 811 812 ! In case of dumped integrals compute an additional taper term 813 IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN 814 se_int_screen%ft = 1.0_dp 815 IF (itype /= do_method_pchg) THEN 816 se_int_screen%ft = taper_eval(se_taper%taper_add, rij) 817 END IF 818 END IF 819 820 ! Contribution from the sp shells 821 CALL nucint_sp_num(sepi, sepj, rij, core=core, itype=itype, & 822 se_int_control=se_int_control, se_int_screen=se_int_screen) 823 824 IF (sepi%dorb .OR. sepj%dorb) THEN 825 ! Compute the contribution from d shells 826 CALL nucint_d_num(sepi, sepj, rij, core, itype, se_int_control, & 827 se_int_screen) 828 END IF 829 830 ! Tapering the integrals 831 DO i = 1, sepi%core_size 832 core(i, 1) = ft*core(i, 1) 833 END DO 834 DO i = 1, sepj%core_size 835 core(i, 2) = ft*core(i, 2) 836 END DO 837 838 END SUBROUTINE core_nucint_num 839 840! ************************************************************************************************** 841!> \brief ... 842!> \param sepi ... 843!> \param sepj ... 844!> \param rij ... 845!> \param ssss ... 846!> \param core ... 847!> \param itype ... 848!> \param se_int_control ... 849!> \param se_int_screen ... 850!> \par History 851!> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver 852!> for computing integrals 853!> Teodoro Laino (04.2008) [tlaino] - Totally rewritten: nothing to do with 854!> the old version 855! ************************************************************************************************** 856 857 SUBROUTINE nucint_sp_num(sepi, sepj, rij, ssss, core, itype, se_int_control, & 858 se_int_screen) 859 TYPE(semi_empirical_type), POINTER :: sepi, sepj 860 REAL(dp), INTENT(IN) :: rij 861 REAL(dp), INTENT(INOUT), OPTIONAL :: ssss 862 REAL(dp), DIMENSION(10, 2), INTENT(INOUT), & 863 OPTIONAL :: core 864 INTEGER, INTENT(IN) :: itype 865 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 866 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 867 868 CHARACTER(len=*), PARAMETER :: routineN = 'nucint_sp_num', routineP = moduleN//':'//routineN 869 870 INTEGER :: ij, kl 871 LOGICAL :: l_core, l_ssss, si, sj 872 873 l_core = PRESENT(core) 874 l_ssss = PRESENT(ssss) 875 IF (.NOT. (l_core .OR. l_ssss)) RETURN 876 si = (sepi%natorb > 1) 877 sj = (sepj%natorb > 1) 878 879 ij = indexa(1, 1) 880 IF (l_ssss) THEN 881 ! Store the value for <S S | S S > (Used for computing the core-core interactions) 882 ssss = ijkl_sp(sepi, sepj, ij, ij, 0, 0, 0, 0, -1, rij, CPPint_args) 883 END IF 884 885 IF (l_core) THEN 886 ! <S S | S S > 887 kl = indexa(1, 1) 888 core(1, 1) = ijkl_sp(sepi, sepj, kl, ij, 0, 0, 0, 0, 2, rij, CPPint_args)*sepj%zeff 889 IF (si) THEN 890 ! <S P | S S > 891 kl = indexa(2, 1) 892 core(2, 1) = ijkl_sp(sepi, sepj, kl, ij, 0, 1, 0, 0, 2, rij, CPPint_args)*sepj%zeff 893 ! <P P | S S > 894 kl = indexa(2, 2) 895 core(3, 1) = ijkl_sp(sepi, sepj, kl, ij, 1, 1, 0, 0, 2, rij, CPPint_args)*sepj%zeff 896 ! <P+ P+ | S S > 897 kl = indexa(3, 3) 898 core(4, 1) = ijkl_sp(sepi, sepj, kl, ij, 1, 1, 0, 0, 2, rij, CPPint_args)*sepj%zeff 899 END IF 900 901 ! <S S | S S > 902 kl = indexa(1, 1) 903 core(1, 2) = ijkl_sp(sepi, sepj, ij, kl, 0, 0, 0, 0, 1, rij, CPPint_args)*sepi%zeff 904 IF (sj) THEN 905 ! <S S | S P > 906 kl = indexa(2, 1) 907 core(2, 2) = ijkl_sp(sepi, sepj, ij, kl, 0, 0, 0, 1, 1, rij, CPPint_args)*sepi%zeff 908 ! <S S | P P > 909 kl = indexa(2, 2) 910 core(3, 2) = ijkl_sp(sepi, sepj, ij, kl, 0, 0, 1, 1, 1, rij, CPPint_args)*sepi%zeff 911 ! <S S | P+ P+ > 912 kl = indexa(3, 3) 913 core(4, 2) = ijkl_sp(sepi, sepj, ij, kl, 0, 0, 1, 1, 1, rij, CPPint_args)*sepi%zeff 914 END IF 915 END IF 916 917 END SUBROUTINE nucint_sp_num 918 919! ************************************************************************************************** 920!> \brief Calculates the nuclear attraction integrals involving d orbitals 921!> \param sepi paramters of atom i 922!> \param sepj paramters of atom j 923!> \param rij interatomic distance 924!> \param core 4 X 2 array of electron-core attraction integrals 925!> 926!> The storage of the nuclear attraction integrals core(kl/ij) iS 927!> (SS/)=1, (SP/)=2, (PP/)=3, (P+P+/)=4, (SD/)=5, 928!> (DP/)=6, (DD/)=7, (D+P+)=8, (D+D+/)=9, (D#D#)=10 929!> 930!> where ij=1 if the orbitals centred on atom i, =2 if on atom j. 931!> \param itype type of semi_empirical model 932!> extension to the original routine to compute qm/mm integrals 933!> \param se_int_control input parameters that control the calculation of SE 934!> integrals (shortrange, R3 residual, screening type) 935!> \param se_int_screen contains information for computing the screened 936!> integrals KDSO-D 937!> \author 938!> Teodoro Laino (03.2008) [tlaino] - University of Zurich 939! ************************************************************************************************** 940 SUBROUTINE nucint_d_num(sepi, sepj, rij, core, itype, se_int_control, & 941 se_int_screen) 942 TYPE(semi_empirical_type), POINTER :: sepi, sepj 943 REAL(dp), INTENT(IN) :: rij 944 REAL(dp), DIMENSION(10, 2), INTENT(INOUT) :: core 945 INTEGER, INTENT(IN) :: itype 946 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 947 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 948 949 CHARACTER(len=*), PARAMETER :: routineN = 'nucint_d_num', routineP = moduleN//':'//routineN 950 951 INTEGER :: ij, kl 952 953! Check if d-orbitals are present 954 955 IF (sepi%dorb .OR. sepj%dorb) THEN 956 ij = indexa(1, 1) 957 IF (sepj%dorb) THEN 958 ! <S S | D S> 959 kl = indexa(5, 1) 960 core(5, 2) = ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 0, 1, rij, CPPint_args)*sepi%zeff 961 ! <S S | D P > 962 kl = indexa(5, 2) 963 core(6, 2) = ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 1, 1, rij, CPPint_args)*sepi%zeff 964 ! <S S | D D > 965 kl = indexa(5, 5) 966 core(7, 2) = ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, CPPint_args)*sepi%zeff 967 ! <S S | D+P+> 968 kl = indexa(6, 3) 969 core(8, 2) = ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 1, 1, rij, CPPint_args)*sepi%zeff 970 ! <S S | D+D+> 971 kl = indexa(6, 6) 972 core(9, 2) = ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, CPPint_args)*sepi%zeff 973 ! <S S | D#D#> 974 kl = indexa(8, 8) 975 core(10, 2) = ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, CPPint_args)*sepi%zeff 976 END IF 977 IF (sepi%dorb) THEN 978 ! <D S | S S> 979 kl = indexa(5, 1) 980 core(5, 1) = ijkl_d(sepi, sepj, kl, ij, 2, 0, 0, 0, 2, rij, CPPint_args)*sepj%zeff 981 ! <D P | S S > 982 kl = indexa(5, 2) 983 core(6, 1) = ijkl_d(sepi, sepj, kl, ij, 2, 1, 0, 0, 2, rij, CPPint_args)*sepj%zeff 984 ! <D D | S S > 985 kl = indexa(5, 5) 986 core(7, 1) = ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, CPPint_args)*sepj%zeff 987 ! <D+P+| S S > 988 kl = indexa(6, 3) 989 core(8, 1) = ijkl_d(sepi, sepj, kl, ij, 2, 1, 0, 0, 2, rij, CPPint_args)*sepj%zeff 990 ! <D+D+| S S > 991 kl = indexa(6, 6) 992 core(9, 1) = ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, CPPint_args)*sepj%zeff 993 ! <D#D#| S S > 994 kl = indexa(8, 8) 995 core(10, 1) = ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, CPPint_args)*sepj%zeff 996 END IF 997 998 END IF 999 END SUBROUTINE nucint_d_num 1000 1001! ************************************************************************************************** 1002!> \brief Numerical Derivatives for rotint 1003!> \param sepi ... 1004!> \param sepj ... 1005!> \param r ... 1006!> \param dw ... 1007!> \param delta ... 1008!> \param se_int_control ... 1009!> \param se_taper ... 1010! ************************************************************************************************** 1011 SUBROUTINE drotint_num(sepi, sepj, r, dw, delta, se_int_control, se_taper) 1012 TYPE(semi_empirical_type), POINTER :: sepi, sepj 1013 REAL(dp), DIMENSION(3), INTENT(IN) :: r 1014 REAL(dp), DIMENSION(3, 2025), INTENT(OUT) :: dw 1015 REAL(dp), INTENT(IN) :: delta 1016 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 1017 TYPE(se_taper_type), POINTER :: se_taper 1018 1019 CHARACTER(len=*), PARAMETER :: routineN = 'drotint_num', routineP = moduleN//':'//routineN 1020 1021 INTEGER :: i, j, nsize 1022 REAL(dp) :: od 1023 REAL(dp), DIMENSION(2025) :: wm, wp 1024 REAL(dp), DIMENSION(3) :: rr 1025 1026 od = 0.5_dp/delta 1027 nsize = sepi%atm_int_size*sepj%atm_int_size 1028 DO i = 1, 3 1029 rr = r 1030 rr(i) = rr(i) + delta 1031 CALL rotint_num(sepi, sepj, rr, wp, se_int_control, se_taper=se_taper) 1032 rr(i) = rr(i) - 2._dp*delta 1033 CALL rotint_num(sepi, sepj, rr, wm, se_int_control, se_taper=se_taper) 1034 DO j = 1, nsize 1035 dw(i, j) = od*(wp(j) - wm(j)) 1036 END DO 1037 END DO 1038 1039 END SUBROUTINE drotint_num 1040 1041! ************************************************************************************************** 1042!> \brief Numerical Derivatives for rotnuc 1043!> \param sepi ... 1044!> \param sepj ... 1045!> \param r ... 1046!> \param de1b ... 1047!> \param de2a ... 1048!> \param itype ... 1049!> \param delta ... 1050!> \param se_int_control ... 1051!> \param se_taper ... 1052! ************************************************************************************************** 1053 SUBROUTINE drotnuc_num(sepi, sepj, r, de1b, de2a, itype, delta, se_int_control, se_taper) 1054 TYPE(semi_empirical_type), POINTER :: sepi, sepj 1055 REAL(dp), DIMENSION(3), INTENT(IN) :: r 1056 REAL(dp), DIMENSION(3, 45), INTENT(OUT), OPTIONAL :: de1b, de2a 1057 INTEGER, INTENT(IN) :: itype 1058 REAL(dp), INTENT(IN) :: delta 1059 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 1060 TYPE(se_taper_type), POINTER :: se_taper 1061 1062 CHARACTER(len=*), PARAMETER :: routineN = 'drotnuc_num', routineP = moduleN//':'//routineN 1063 1064 INTEGER :: i, j 1065 LOGICAL :: l_de1b, l_de2a 1066 REAL(dp) :: od 1067 REAL(dp), DIMENSION(3) :: rr 1068 REAL(dp), DIMENSION(45) :: e1m, e1p, e2m, e2p 1069 1070 l_de1b = PRESENT(de1b) 1071 l_de2a = PRESENT(de2a) 1072 IF (.NOT. (l_de1b .OR. l_de2a)) RETURN 1073 od = 0.5_dp/delta 1074 DO i = 1, 3 1075 rr = r 1076 rr(i) = rr(i) + delta 1077 CALL rotnuc_num(sepi, sepj, rr, e1p, e2p, itype, se_int_control, se_taper=se_taper) 1078 rr(i) = rr(i) - 2._dp*delta 1079 CALL rotnuc_num(sepi, sepj, rr, e1m, e2m, itype, se_int_control, se_taper=se_taper) 1080 IF (l_de1b) THEN 1081 DO j = 1, sepi%atm_int_size 1082 de1b(i, j) = od*(e1p(j) - e1m(j)) 1083 END DO 1084 END IF 1085 IF (l_de2a) THEN 1086 DO j = 1, sepj%atm_int_size 1087 de2a(i, j) = od*(e2p(j) - e2m(j)) 1088 END DO 1089 END IF 1090 END DO 1091 END SUBROUTINE drotnuc_num 1092 1093! ************************************************************************************************** 1094!> \brief Numerical Derivatives for corecore 1095!> \param sepi ... 1096!> \param sepj ... 1097!> \param r ... 1098!> \param denuc ... 1099!> \param itype ... 1100!> \param delta ... 1101!> \param se_int_control ... 1102!> \param se_taper ... 1103! ************************************************************************************************** 1104 SUBROUTINE dcorecore_num(sepi, sepj, r, denuc, itype, delta, se_int_control, se_taper) 1105 TYPE(semi_empirical_type), POINTER :: sepi, sepj 1106 REAL(dp), DIMENSION(3), INTENT(IN) :: r 1107 REAL(dp), DIMENSION(3), INTENT(OUT) :: denuc 1108 INTEGER, INTENT(IN) :: itype 1109 REAL(dp), INTENT(IN) :: delta 1110 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 1111 TYPE(se_taper_type), POINTER :: se_taper 1112 1113 CHARACTER(len=*), PARAMETER :: routineN = 'dcorecore_num', routineP = moduleN//':'//routineN 1114 1115 INTEGER :: i 1116 REAL(dp) :: enucm, enucp, od 1117 REAL(dp), DIMENSION(3) :: rr 1118 1119 od = 0.5_dp/delta 1120 DO i = 1, 3 1121 rr = r 1122 rr(i) = rr(i) + delta 1123 CALL corecore_num(sepi, sepj, rr, enucp, itype, se_int_control, se_taper=se_taper) 1124 rr(i) = rr(i) - 2._dp*delta 1125 CALL corecore_num(sepi, sepj, rr, enucm, itype, se_int_control, se_taper=se_taper) 1126 denuc(i) = od*(enucp - enucm) 1127 END DO 1128 END SUBROUTINE dcorecore_num 1129 1130! ************************************************************************************************** 1131!> \brief Numerical Derivatives for corecore 1132!> \param sepi ... 1133!> \param sepj ... 1134!> \param r ... 1135!> \param denuc ... 1136!> \param itype ... 1137!> \param delta ... 1138!> \param se_int_control ... 1139!> \param se_taper ... 1140! ************************************************************************************************** 1141 SUBROUTINE dcorecore_el_num(sepi, sepj, r, denuc, itype, delta, se_int_control, se_taper) 1142 TYPE(semi_empirical_type), POINTER :: sepi, sepj 1143 REAL(dp), DIMENSION(3), INTENT(IN) :: r 1144 REAL(dp), DIMENSION(3), INTENT(OUT) :: denuc 1145 INTEGER, INTENT(IN) :: itype 1146 REAL(dp), INTENT(IN) :: delta 1147 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 1148 TYPE(se_taper_type), POINTER :: se_taper 1149 1150 CHARACTER(len=*), PARAMETER :: routineN = 'dcorecore_el_num', & 1151 routineP = moduleN//':'//routineN 1152 1153 INTEGER :: i 1154 REAL(dp) :: enucm, enucp, od 1155 REAL(dp), DIMENSION(3) :: rr 1156 1157 od = 0.5_dp/delta 1158 DO i = 1, 3 1159 rr = r 1160 rr(i) = rr(i) + delta 1161 CALL corecore_el_num(sepi, sepj, rr, enucp, itype, se_int_control, se_taper=se_taper) 1162 rr(i) = rr(i) - 2._dp*delta 1163 CALL corecore_el_num(sepi, sepj, rr, enucm, itype, se_int_control, se_taper=se_taper) 1164 denuc(i) = od*(enucp - enucm) 1165 END DO 1166 END SUBROUTINE dcorecore_el_num 1167 1168END MODULE semi_empirical_int_num 1169