1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculation of contracted, spherical Gaussian integrals using the (OS) integral 8!> scheme. Routines for the following two-center integrals: 9!> i) (a|O(r12)|b) where O(r12) is the overlap, coulomb operator etc. 10!> ii) (aba) and (abb) s-overlaps 11!> \par History 12!> created [06.2015] 13!> 05.2019: Added truncated coulomb operator (A. Bussy) 14!> \author Dorothea Golze 15! ************************************************************************************************** 16MODULE generic_os_integrals 17 USE ai_contraction_sphi, ONLY: ab_contract,& 18 abc_contract,& 19 abcd_contract 20 USE ai_derivatives, ONLY: dabdr_noscreen 21 USE ai_operator_ra2m, ONLY: operator_ra2m 22 USE ai_operators_r12, ONLY: ab_sint_os,& 23 cps_coulomb2,& 24 cps_gauss2,& 25 cps_truncated2,& 26 cps_verf2,& 27 cps_verfc2,& 28 cps_vgauss2,& 29 operator2 30 USE ai_overlap, ONLY: overlap_aab,& 31 overlap_ab,& 32 overlap_abb 33 USE ai_overlap_aabb, ONLY: overlap_aabb 34 USE basis_set_types, ONLY: get_gto_basis_set,& 35 gto_basis_set_type 36 USE constants_operator, ONLY: operator_coulomb,& 37 operator_gauss,& 38 operator_truncated,& 39 operator_verf,& 40 operator_verfc,& 41 operator_vgauss 42 USE debug_os_integrals, ONLY: overlap_aabb_test,& 43 overlap_ab_test,& 44 overlap_abc_test 45 USE kinds, ONLY: dp 46 USE orbital_pointers, ONLY: ncoset 47#include "./base/base_uses.f90" 48 49 IMPLICIT NONE 50 51 PRIVATE 52 53! ************************************************************************************************** 54 55 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'generic_os_integrals' 56 57 PUBLIC :: int_operators_r12_ab_os, int_overlap_ab_os, int_ra2m_ab_os, int_overlap_aba_os, & 58 int_overlap_abb_os, int_overlap_aabb_os 59 60CONTAINS 61 62! ************************************************************************************************** 63!> \brief Calcululates the two-center integrals of the type (a|O(r12)|b) using the OS scheme 64!> \param r12_operator the integral operator, which depends on r12=|r1-r2| 65!> \param vab integral matrix of spherical contracted Gaussian functions 66!> \param dvab derivative of the integrals 67!> \param rab distance vector between center A and B 68!> \param fba basis at center A 69!> \param fbb basis at center B 70!> \param omega parameter in the operator 71!> \param r_cutoff the cutoff in case of truncated coulomb operator 72!> \param calculate_forces ... 73! ************************************************************************************************** 74 SUBROUTINE int_operators_r12_ab_os(r12_operator, vab, dvab, rab, fba, fbb, omega, r_cutoff, & 75 calculate_forces) 76 77 INTEGER, INTENT(IN) :: r12_operator 78 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab 79 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), & 80 OPTIONAL :: dvab 81 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 82 TYPE(gto_basis_set_type), POINTER :: fba, fbb 83 REAL(KIND=dp), INTENT(IN), OPTIONAL :: omega, r_cutoff 84 LOGICAL, INTENT(IN) :: calculate_forces 85 86 CHARACTER(LEN=*), PARAMETER :: routineN = 'int_operators_r12_ab_os', & 87 routineP = moduleN//':'//routineN 88 89 INTEGER :: handle 90 REAL(KIND=dp) :: my_omega, my_r_cutoff 91 92 PROCEDURE(ab_sint_os), POINTER :: cps_operator2 93 94 NULLIFY (cps_operator2) 95 CALL timeset(routineN, handle) 96 97 my_omega = 1.0_dp 98 my_r_cutoff = 1.0_dp 99 100 SELECT CASE (r12_operator) 101 CASE (operator_coulomb) 102 cps_operator2 => cps_coulomb2 103 CASE (operator_verf) 104 cps_operator2 => cps_verf2 105 IF (PRESENT(omega)) my_omega = omega 106 CASE (operator_verfc) 107 cps_operator2 => cps_verfc2 108 IF (PRESENT(omega)) my_omega = omega 109 CASE (operator_vgauss) 110 cps_operator2 => cps_vgauss2 111 IF (PRESENT(omega)) my_omega = omega 112 CASE (operator_gauss) 113 cps_operator2 => cps_gauss2 114 IF (PRESENT(omega)) my_omega = omega 115 CASE (operator_truncated) 116 cps_operator2 => cps_truncated2 117 IF (PRESENT(r_cutoff)) my_r_cutoff = r_cutoff 118 CASE DEFAULT 119 CPABORT("Operator not available") 120 END SELECT 121 122 CALL int_operator_ab_os_low(cps_operator2, vab, dvab, rab, fba, fbb, my_omega, my_r_cutoff, & 123 calculate_forces) 124 125 CALL timestop(handle) 126 127 END SUBROUTINE int_operators_r12_ab_os 128 129! ************************************************************************************************** 130!> \brief calculate integrals (a|O(r12)|b) 131!> \param cps_operator2 procedure pointer for the respective operator. 132!> \param vab integral matrix of spherical contracted Gaussian functions 133!> \param dvab derivative of the integrals 134!> \param rab distance vector between center A and B 135!> \param fba basis at center A 136!> \param fbb basis at center B 137!> \param omega parameter in the operator 138!> \param r_cutoff ... 139!> \param calculate_forces ... 140! ************************************************************************************************** 141 SUBROUTINE int_operator_ab_os_low(cps_operator2, vab, dvab, rab, fba, fbb, omega, r_cutoff, & 142 calculate_forces) 143 144 PROCEDURE(ab_sint_os), POINTER :: cps_operator2 145 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab 146 REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, & 147 INTENT(INOUT) :: dvab 148 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 149 TYPE(gto_basis_set_type), POINTER :: fba, fbb 150 REAL(KIND=dp), INTENT(IN) :: omega, r_cutoff 151 LOGICAL, INTENT(IN) :: calculate_forces 152 153 CHARACTER(LEN=*), PARAMETER :: routineN = 'int_operator_ab_os_low', routineP = moduleN//':'//routineN 154 155 INTEGER :: handle, i, iset, jset, lds, m1, m2, & 156 maxco, maxcoa, maxcob, maxl, maxla, & 157 maxlb, ncoa, ncoap, ncob, ncobp, & 158 nseta, nsetb, sgfa, sgfb 159 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 160 npgfb, nsgfa, nsgfb 161 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 162 REAL(KIND=dp) :: dab, rab2 163 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vac, vac_plus 164 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: devab, vwork 165 REAL(KIND=dp), DIMENSION(:, :), POINTER :: sphi_a, sphi_b, zeta, zetb 166 167 CALL timeset(routineN, handle) 168 NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, & 169 first_sgfa, first_sgfb, sphi_a, sphi_b, zeta, zetb) 170 171 ! basis ikind 172 first_sgfa => fba%first_sgf 173 la_max => fba%lmax 174 la_min => fba%lmin 175 npgfa => fba%npgf 176 nseta = fba%nset 177 nsgfa => fba%nsgf_set 178 sphi_a => fba%sphi 179 zeta => fba%zet 180 ! basis jkind 181 first_sgfb => fbb%first_sgf 182 lb_max => fbb%lmax 183 lb_min => fbb%lmin 184 npgfb => fbb%npgf 185 nsetb = fbb%nset 186 nsgfb => fbb%nsgf_set 187 sphi_b => fbb%sphi 188 zetb => fbb%zet 189 190 CALL get_gto_basis_set(fba, maxco=maxcoa, maxl=maxla) 191 CALL get_gto_basis_set(fbb, maxco=maxcob, maxl=maxlb) 192 maxco = MAX(maxcoa, maxcob) 193 IF (calculate_forces) THEN 194 maxl = MAX(maxla + 1, maxlb) 195 ELSE 196 maxl = MAX(maxla, maxlb) 197 ENDIF 198 lds = ncoset(maxl) 199 200 rab2 = SUM(rab*rab) 201 dab = SQRT(rab2) 202 203 vab = 0.0_dp 204 IF (calculate_forces) dvab = 0.0_dp 205 206 DO iset = 1, nseta 207 208 ncoa = npgfa(iset)*ncoset(la_max(iset)) 209 ncoap = npgfa(iset)*ncoset(la_max(iset) + 1) 210 sgfa = first_sgfa(1, iset) 211 212 DO jset = 1, nsetb 213 214 ncob = npgfb(jset)*ncoset(lb_max(jset)) 215 ncobp = npgfb(jset)*ncoset(lb_max(jset) + 1) 216 sgfb = first_sgfb(1, jset) 217 m1 = sgfa + nsgfa(iset) - 1 218 m2 = sgfb + nsgfb(jset) - 1 219 220 ! calculate integrals 221 IF (calculate_forces) THEN 222 ALLOCATE (vwork(ncoap, ncobp, la_max(iset) + lb_max(jset) + 3), & 223 vac(ncoa, ncob), vac_plus(ncoap, ncobp), devab(ncoa, ncob, 3)) 224 devab = 0._dp 225 vwork = 0.0_dp 226 vac = 0.0_dp 227 CALL operator2(cps_operator2, la_max(iset) + 1, npgfa(iset), zeta(:, iset), la_min(iset), & 228 lb_max(jset) + 1, npgfb(jset), zetb(:, jset), lb_min(jset), & 229 omega, r_cutoff, rab, rab2, vac, vwork, maxder=1, vac_plus=vac_plus) 230 CALL dabdr_noscreen(la_max(iset), npgfa(iset), zeta(:, iset), lb_max(jset), npgfb(jset), & 231 vac_plus, devab(:, :, 1), devab(:, :, 2), devab(:, :, 3)) 232 DO i = 1, 3 233 CALL ab_contract(dvab(sgfa:m1, sgfb:m2, i), devab(:, :, i), sphi_a(:, sgfa:), & 234 sphi_b(:, sgfb:), ncoa, ncob, nsgfa(iset), nsgfb(jset)) 235 ENDDO 236 237 ELSE 238 ALLOCATE (vwork(ncoa, ncob, la_max(iset) + lb_max(jset) + 1), & 239 vac(ncoa, ncob), vac_plus(ncoap, ncobp), devab(ncoa, ncob, 3)) 240 vwork = 0.0_dp 241 vac = 0.0_dp 242 CALL operator2(cps_operator2, la_max(iset), npgfa(iset), zeta(:, iset), la_min(iset), & 243 lb_max(jset), npgfb(jset), zetb(:, jset), lb_min(jset), & 244 omega, r_cutoff, rab, rab2, vac, vwork) 245 ENDIF 246 247 CALL ab_contract(vab(sgfa:m1, sgfb:m2), vac(1:ncoa, 1:ncob), sphi_a(:, sgfa:), sphi_b(:, sgfb:), & 248 ncoa, ncob, nsgfa(iset), nsgfb(jset)) 249 DEALLOCATE (vwork, vac, vac_plus, devab) 250 END DO 251 END DO 252 253 CALL timestop(handle) 254 255 END SUBROUTINE int_operator_ab_os_low 256 257! ************************************************************************************************** 258!> \brief calculate overlap integrals (a,b) 259!> \param sab integral (a,b) 260!> \param dsab derivative of sab with respect to A 261!> \param rab distance vector between center A and B 262!> \param fba basis at center A 263!> \param fbb basis at center B 264!> \param calculate_forces ... 265!> \param debug integrals are debugged by recursive routines if requested 266!> \param dmax maximal deviation between integrals when debugging 267! ************************************************************************************************** 268 SUBROUTINE int_overlap_ab_os(sab, dsab, rab, fba, fbb, calculate_forces, debug, dmax) 269 270 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: sab 271 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), & 272 OPTIONAL :: dsab 273 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 274 TYPE(gto_basis_set_type), POINTER :: fba, fbb 275 LOGICAL, INTENT(IN) :: calculate_forces, debug 276 REAL(KIND=dp), INTENT(INOUT) :: dmax 277 278 CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_ab_os', & 279 routineP = moduleN//':'//routineN 280 281 INTEGER :: handle, i, iset, jset, m1, m2, maxco, & 282 maxcoa, maxcob, maxl, maxla, maxlb, & 283 ncoa, ncob, nseta, nsetb, sgfa, sgfb 284 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 285 npgfb, nsgfa, nsgfb 286 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 287 LOGICAL :: failure 288 REAL(KIND=dp) :: dab, ra(3), rb(3) 289 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sint 290 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dsint 291 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 292 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, scon_a, scon_b, zeta, zetb 293 294 failure = .FALSE. 295 CALL timeset(routineN, handle) 296 NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, & 297 first_sgfa, first_sgfb, scon_a, scon_b, zeta, zetb) 298 299 ! basis ikind 300 first_sgfa => fba%first_sgf 301 la_max => fba%lmax 302 la_min => fba%lmin 303 npgfa => fba%npgf 304 nseta = fba%nset 305 nsgfa => fba%nsgf_set 306 rpgfa => fba%pgf_radius 307 set_radius_a => fba%set_radius 308 scon_a => fba%scon 309 zeta => fba%zet 310 ! basis jkind 311 first_sgfb => fbb%first_sgf 312 lb_max => fbb%lmax 313 lb_min => fbb%lmin 314 npgfb => fbb%npgf 315 nsetb = fbb%nset 316 nsgfb => fbb%nsgf_set 317 rpgfb => fbb%pgf_radius 318 set_radius_b => fbb%set_radius 319 scon_b => fbb%scon 320 zetb => fbb%zet 321 322 CALL get_gto_basis_set(fba, maxco=maxcoa, maxl=maxla) 323 CALL get_gto_basis_set(fbb, maxco=maxcob, maxl=maxlb) 324 maxco = MAX(maxcoa, maxcob) 325 maxl = MAX(maxla, maxlb) 326 ALLOCATE (sint(maxco, maxco)) 327 ALLOCATE (dsint(maxco, maxco, 3)) 328 329 dab = SQRT(SUM(rab**2)) 330 331 sab = 0.0_dp 332 IF (calculate_forces) THEN 333 IF (PRESENT(dsab)) dsab = 0.0_dp 334 END IF 335 336 DO iset = 1, nseta 337 338 ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1)) 339 sgfa = first_sgfa(1, iset) 340 341 DO jset = 1, nsetb 342 343 IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE 344 345 ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1)) 346 sgfb = first_sgfb(1, jset) 347 m1 = sgfa + nsgfa(iset) - 1 348 m2 = sgfb + nsgfb(jset) - 1 349 350 IF (calculate_forces) THEN 351 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 352 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), & 353 rab, sint, dsint) 354 ELSE 355 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 356 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), & 357 rab, sint) 358 359 ENDIF 360 361 ! debug if requested 362 IF (debug) THEN 363 ra = 0.0_dp 364 rb = rab 365 CALL overlap_ab_test(la_max(iset), la_min(iset), npgfa(iset), zeta(:, iset), & 366 lb_max(jset), lb_min(jset), npgfb(jset), zetb(:, jset), & 367 ra, rb, sint, dmax) 368 ENDIF 369 370 CALL ab_contract(sab(sgfa:m1, sgfb:m2), sint(1:ncoa, 1:ncob), scon_a(:, sgfa:), scon_b(:, sgfb:), & 371 ncoa, ncob, nsgfa(iset), nsgfb(jset)) 372 IF (calculate_forces) THEN 373 DO i = 1, 3 374 CALL ab_contract(dsab(sgfa:m1, sgfb:m2, i), dsint(1:ncoa, 1:ncob, i), scon_a(:, sgfa:), & 375 scon_b(:, sgfb:), ncoa, ncob, nsgfa(iset), nsgfb(jset)) 376 ENDDO 377 ENDIF 378 END DO 379 END DO 380 381 DEALLOCATE (sint, dsint) 382 383 CALL timestop(handle) 384 385 END SUBROUTINE int_overlap_ab_os 386 387! ************************************************************************************************** 388!> \brief calculate integrals (a|(r-Ra)^(2m)|b) 389!> \param sab integral (a|(r-Ra)^(2m)|b) 390!> \param dsab derivative of sab with respect to A 391!> \param rab distance vector between center A and B 392!> \param fba fba basis at center A 393!> \param fbb fbb basis at center B 394!> \param m exponent m in operator (r-Ra)^(2m) 395!> \param calculate_forces ... 396! ************************************************************************************************** 397 SUBROUTINE int_ra2m_ab_os(sab, dsab, rab, fba, fbb, m, calculate_forces) 398 399 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: sab 400 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), & 401 OPTIONAL :: dsab 402 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 403 TYPE(gto_basis_set_type), POINTER :: fba, fbb 404 INTEGER, INTENT(IN) :: m 405 LOGICAL, INTENT(IN) :: calculate_forces 406 407 CHARACTER(LEN=*), PARAMETER :: routineN = 'int_ra2m_ab_os', routineP = moduleN//':'//routineN 408 409 INTEGER :: handle, i, iset, jset, m1, m2, maxco, & 410 maxcoa, maxcob, maxl, maxla, maxlb, & 411 ncoa, ncob, nseta, nsetb, sgfa, sgfb 412 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 413 npgfb, nsgfa, nsgfb 414 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 415 LOGICAL :: failure 416 REAL(KIND=dp) :: dab 417 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sint 418 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dsint 419 REAL(KIND=dp), DIMENSION(:, :), POINTER :: scon_a, scon_b, zeta, zetb 420 421 failure = .FALSE. 422 CALL timeset(routineN, handle) 423 NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, & 424 first_sgfa, first_sgfb, scon_a, scon_b, zeta, zetb) 425 426 ! basis ikind 427 first_sgfa => fba%first_sgf 428 la_max => fba%lmax 429 la_min => fba%lmin 430 npgfa => fba%npgf 431 nseta = fba%nset 432 nsgfa => fba%nsgf_set 433 scon_a => fba%scon 434 zeta => fba%zet 435 ! basis jkind 436 first_sgfb => fbb%first_sgf 437 lb_max => fbb%lmax 438 lb_min => fbb%lmin 439 npgfb => fbb%npgf 440 nsetb = fbb%nset 441 nsgfb => fbb%nsgf_set 442 scon_b => fbb%scon 443 zetb => fbb%zet 444 445 CALL get_gto_basis_set(fba, maxco=maxcoa, maxl=maxla) 446 CALL get_gto_basis_set(fbb, maxco=maxcob, maxl=maxlb) 447 maxco = MAX(maxcoa, maxcob) 448 maxl = MAX(maxla, maxlb) 449 ALLOCATE (sint(maxco, maxco)) 450 ALLOCATE (dsint(maxco, maxco, 3)) 451 452 dab = SQRT(SUM(rab**2)) 453 454 sab = 0.0_dp 455 IF (calculate_forces) THEN 456 IF (PRESENT(dsab)) dsab = 0.0_dp 457 END IF 458 459 DO iset = 1, nseta 460 461 ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1)) 462 sgfa = first_sgfa(1, iset) 463 464 DO jset = 1, nsetb 465 466 ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1)) 467 sgfb = first_sgfb(1, jset) 468 m1 = sgfa + nsgfa(iset) - 1 469 m2 = sgfb + nsgfb(jset) - 1 470 471 CALL operator_ra2m(la_max(iset), la_min(iset), npgfa(iset), zeta(:, iset), & 472 lb_max(jset), lb_min(jset), npgfb(jset), zetb(:, jset), & 473 m, rab, sint, dsint, calculate_forces) 474 475 CALL ab_contract(sab(sgfa:m1, sgfb:m2), sint, scon_a(:, sgfa:), scon_b(:, sgfb:), & 476 ncoa, ncob, nsgfa(iset), nsgfb(jset)) 477 IF (calculate_forces) THEN 478 DO i = 1, 3 479 CALL ab_contract(dsab(sgfa:m1, sgfb:m2, i), dsint(:, :, i), scon_a(:, sgfa:), & 480 scon_b(:, sgfb:), ncoa, ncob, nsgfa(iset), nsgfb(jset)) 481 ENDDO 482 ENDIF 483 END DO 484 END DO 485 486 DEALLOCATE (sint, dsint) 487 488 CALL timestop(handle) 489 490 END SUBROUTINE int_ra2m_ab_os 491 492! ************************************************************************************************** 493!> \brief calculate integrals (a,b,fa) 494!> \param abaint integral (a,b,fa) 495!> \param dabdaint derivative of abaint with respect to A 496!> \param rab distance vector between center A and B 497!> \param oba orbital basis at center A 498!> \param obb orbital basis at center B 499!> \param fba auxiliary basis set at center A 500!> \param calculate_forces ... 501!> \param debug integrals are debugged by recursive routines if requested 502!> \param dmax maximal deviation between integrals when debugging 503! ************************************************************************************************** 504 SUBROUTINE int_overlap_aba_os(abaint, dabdaint, rab, oba, obb, fba, & 505 calculate_forces, debug, dmax) 506 507 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: abaint 508 REAL(KIND=dp), ALLOCATABLE, & 509 DIMENSION(:, :, :, :), OPTIONAL :: dabdaint 510 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 511 TYPE(gto_basis_set_type), POINTER :: oba, obb, fba 512 LOGICAL, INTENT(IN) :: calculate_forces, debug 513 REAL(KIND=dp), INTENT(INOUT) :: dmax 514 515 CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_aba_os', & 516 routineP = moduleN//':'//routineN 517 518 INTEGER :: handle, i, iset, jset, kaset, m1, m2, & 519 m3, ncoa, ncob, ncoc, nseta, nsetb, & 520 nsetca, sgfa, sgfb, sgfc 521 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, lca_max, & 522 lca_min, npgfa, npgfb, npgfca, nsgfa, & 523 nsgfb, nsgfca 524 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, first_sgfca 525 REAL(KIND=dp) :: dab, dac, dbc 526 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: saba 527 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: sdaba 528 REAL(KIND=dp), DIMENSION(3) :: ra, rac, rb, rbc 529 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b, set_radius_ca 530 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, rpgfca, scon_a, scon_b, & 531 scon_ca, zeta, zetb, zetca 532 533 CALL timeset(routineN, handle) 534 NULLIFY (la_max, la_min, lb_max, lb_min, lca_max, lca_min, npgfa, npgfb, & 535 npgfca, nsgfa, nsgfb, nsgfca) 536 NULLIFY (first_sgfa, first_sgfb, first_sgfca, set_radius_a, set_radius_b, & 537 set_radius_ca, rpgfa, rpgfb, rpgfca, scon_a, scon_b, scon_ca, & 538 zeta, zetb, zetca) 539 540 ! basis ikind 541 first_sgfa => oba%first_sgf 542 la_max => oba%lmax 543 la_min => oba%lmin 544 npgfa => oba%npgf 545 nseta = oba%nset 546 nsgfa => oba%nsgf_set 547 rpgfa => oba%pgf_radius 548 set_radius_a => oba%set_radius 549 scon_a => oba%scon 550 zeta => oba%zet 551 ! basis jkind 552 first_sgfb => obb%first_sgf 553 lb_max => obb%lmax 554 lb_min => obb%lmin 555 npgfb => obb%npgf 556 nsetb = obb%nset 557 nsgfb => obb%nsgf_set 558 rpgfb => obb%pgf_radius 559 set_radius_b => obb%set_radius 560 scon_b => obb%scon 561 zetb => obb%zet 562 563 ! basis RI A 564 first_sgfca => fba%first_sgf 565 lca_max => fba%lmax 566 lca_min => fba%lmin 567 npgfca => fba%npgf 568 nsetca = fba%nset 569 nsgfca => fba%nsgf_set 570 rpgfca => fba%pgf_radius 571 set_radius_ca => fba%set_radius 572 scon_ca => fba%scon 573 zetca => fba%zet 574 575 dab = SQRT(SUM(rab**2)) 576 577 abaint = 0.0_dp 578 IF (calculate_forces) THEN 579 IF (PRESENT(dabdaint)) dabdaint = 0.0_dp 580 END IF 581 582 DO iset = 1, nseta 583 584 ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1)) 585 sgfa = first_sgfa(1, iset) 586 587 DO jset = 1, nsetb 588 589 IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE 590 591 ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1)) 592 sgfb = first_sgfb(1, jset) 593 m1 = sgfa + nsgfa(iset) - 1 594 m2 = sgfb + nsgfb(jset) - 1 595 596 ! calculate integrals abaint and derivative [d(a,b,a)/dA] dabdaint if requested 597 rac = 0._dp 598 dac = 0._dp 599 rbc = -rab 600 dbc = dab 601 602 DO kaset = 1, nsetca 603 604 IF (set_radius_b(jset) + set_radius_ca(kaset) < dab) CYCLE 605 606 ncoc = npgfca(kaset)*(ncoset(lca_max(kaset)) - ncoset(lca_min(kaset) - 1)) 607 sgfc = first_sgfca(1, kaset) 608 m3 = sgfc + nsgfca(kaset) - 1 609 IF (ncoa*ncob*ncoc > 0) THEN 610 ALLOCATE (saba(ncoa, ncob, ncoc)) 611 ! integrals 612 IF (calculate_forces) THEN 613 ALLOCATE (sdaba(ncoa, ncob, ncoc, 3)) 614 CALL overlap_aab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 615 lca_max(kaset), lca_min(kaset), npgfca(kaset), rpgfca(:, kaset), zetca(:, kaset), & 616 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), & 617 rab, saba=saba, daba=sdaba) 618 619 DO i = 1, 3 620 CALL abc_contract(dabdaint(sgfa:m1, sgfb:m2, sgfc:m3, i), sdaba(1:ncoa, 1:ncob, 1:ncoc, i), & 621 scon_a(:, sgfa:), scon_b(:, sgfb:), scon_ca(:, sgfc:), & 622 ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfca(kaset)) 623 ENDDO 624 625 DEALLOCATE (sdaba) 626 ELSE 627 CALL overlap_aab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 628 lca_max(kaset), lca_min(kaset), npgfca(kaset), rpgfca(:, kaset), zetca(:, kaset), & 629 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), & 630 rab, saba=saba) 631 632 ENDIF 633 ! debug if requested 634 IF (debug) THEN 635 ra = 0.0_dp 636 rb = rab 637 CALL overlap_abc_test(la_max(iset), npgfa(iset), zeta(:, iset), la_min(iset), & 638 lb_max(jset), npgfb(jset), zetb(:, jset), lb_min(jset), & 639 lca_max(kaset), npgfca(kaset), zetca(:, kaset), lca_min(kaset), & 640 ra, rb, ra, saba, dmax) 641 ENDIF 642 CALL abc_contract(abaint(sgfa:m1, sgfb:m2, sgfc:m3), saba(1:ncoa, 1:ncob, 1:ncoc), & 643 scon_a(:, sgfa:), scon_b(:, sgfb:), scon_ca(:, sgfc:), & 644 ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfca(kaset)) 645 DEALLOCATE (saba) 646 END IF 647 END DO 648 END DO 649 END DO 650 651 CALL timestop(handle) 652 653 END SUBROUTINE int_overlap_aba_os 654 655! ************************************************************************************************** 656!> \brief calculate integrals (a,b,fb) 657!> \param abbint integral (a,b,fb) 658!> \param dabbint derivative of abbint with respect to A 659!> \param rab distance vector between center A and B 660!> \param oba orbital basis at center A 661!> \param obb orbital basis at center B 662!> \param fbb auxiliary basis set at center B 663!> \param calculate_forces ... 664!> \param debug integrals are debugged by recursive routines if requested 665!> \param dmax maximal deviation between integrals when debugging 666! ************************************************************************************************** 667 SUBROUTINE int_overlap_abb_os(abbint, dabbint, rab, oba, obb, fbb, calculate_forces, debug, dmax) 668 669 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: abbint 670 REAL(KIND=dp), ALLOCATABLE, & 671 DIMENSION(:, :, :, :), OPTIONAL :: dabbint 672 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 673 TYPE(gto_basis_set_type), POINTER :: oba, obb, fbb 674 LOGICAL, INTENT(IN) :: calculate_forces, debug 675 REAL(KIND=dp), INTENT(INOUT) :: dmax 676 677 CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_abb_os', & 678 routineP = moduleN//':'//routineN 679 680 INTEGER :: handle, i, iset, jset, kbset, m1, m2, & 681 m3, ncoa, ncob, ncoc, nseta, nsetb, & 682 nsetcb, sgfa, sgfb, sgfc 683 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, lcb_max, & 684 lcb_min, npgfa, npgfb, npgfcb, nsgfa, & 685 nsgfb, nsgfcb 686 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, first_sgfcb 687 REAL(KIND=dp) :: dab, dac, dbc 688 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: sabb 689 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: sdabb 690 REAL(KIND=dp), DIMENSION(3) :: ra, rac, rb, rbc 691 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b, set_radius_cb 692 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, rpgfcb, scon_a, scon_b, & 693 scon_cb, zeta, zetb, zetcb 694 695 CALL timeset(routineN, handle) 696 NULLIFY (la_max, la_min, lb_max, lb_min, lcb_max, lcb_min, npgfa, npgfb, & 697 npgfcb, nsgfa, nsgfb, nsgfcb) 698 NULLIFY (first_sgfa, first_sgfb, first_sgfcb, set_radius_a, set_radius_b, & 699 set_radius_cb, rpgfa, rpgfb, rpgfcb, scon_a, scon_b, scon_cb, & 700 zeta, zetb, zetcb) 701 702 ! basis ikind 703 first_sgfa => oba%first_sgf 704 la_max => oba%lmax 705 la_min => oba%lmin 706 npgfa => oba%npgf 707 nseta = oba%nset 708 nsgfa => oba%nsgf_set 709 rpgfa => oba%pgf_radius 710 set_radius_a => oba%set_radius 711 scon_a => oba%scon 712 zeta => oba%zet 713 ! basis jkind 714 first_sgfb => obb%first_sgf 715 lb_max => obb%lmax 716 lb_min => obb%lmin 717 npgfb => obb%npgf 718 nsetb = obb%nset 719 nsgfb => obb%nsgf_set 720 rpgfb => obb%pgf_radius 721 set_radius_b => obb%set_radius 722 scon_b => obb%scon 723 zetb => obb%zet 724 725 ! basis RI B 726 first_sgfcb => fbb%first_sgf 727 lcb_max => fbb%lmax 728 lcb_min => fbb%lmin 729 npgfcb => fbb%npgf 730 nsetcb = fbb%nset 731 nsgfcb => fbb%nsgf_set 732 rpgfcb => fbb%pgf_radius 733 set_radius_cb => fbb%set_radius 734 scon_cb => fbb%scon 735 zetcb => fbb%zet 736 737 dab = SQRT(SUM(rab**2)) 738 739 abbint = 0.0_dp 740 IF (calculate_forces) THEN 741 IF (PRESENT(dabbint)) dabbint = 0.0_dp 742 END IF 743 744 DO iset = 1, nseta 745 746 ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1)) 747 sgfa = first_sgfa(1, iset) 748 749 DO jset = 1, nsetb 750 751 IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE 752 753 ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1)) 754 sgfb = first_sgfb(1, jset) 755 m1 = sgfa + nsgfa(iset) - 1 756 m2 = sgfb + nsgfb(jset) - 1 757 758 ! calculate integrals abbint and derivative [d(a,b,b)/dA] dabbint if requested 759 rac = rab 760 dac = dab 761 rbc = 0._dp 762 dbc = 0._dp 763 764 DO kbset = 1, nsetcb 765 766 IF (set_radius_a(iset) + set_radius_cb(kbset) < dab) CYCLE 767 768 ncoc = npgfcb(kbset)*(ncoset(lcb_max(kbset)) - ncoset(lcb_min(kbset) - 1)) 769 sgfc = first_sgfcb(1, kbset) 770 m3 = sgfc + nsgfcb(kbset) - 1 771 IF (ncoa*ncob*ncoc > 0) THEN 772 ALLOCATE (sabb(ncoa, ncob, ncoc)) 773 IF (calculate_forces) THEN 774 ALLOCATE (sdabb(ncoa, ncob, ncoc, 3)) 775 CALL overlap_abb(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 776 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), & 777 lcb_max(kbset), lcb_min(kbset), npgfcb(kbset), rpgfcb(:, kbset), zetcb(:, kbset), & 778 rab, sabb, sdabb) 779 780 DO i = 1, 3 781 CALL abc_contract(dabbint(sgfa:m1, sgfb:m2, sgfc:m3, i), sdabb(1:ncoa, 1:ncob, 1:ncoc, i), & 782 scon_a(:, sgfa:), scon_b(:, sgfb:), scon_cb(:, sgfc:), & 783 ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfcb(kbset)) 784 ENDDO 785 DEALLOCATE (sdabb) 786 ELSE 787 CALL overlap_abb(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 788 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), & 789 lcb_max(kbset), lcb_min(kbset), npgfcb(kbset), rpgfcb(:, kbset), zetcb(:, kbset), & 790 rab, sabb) 791 ENDIF 792 ! debug if requested 793 IF (debug) THEN 794 ra = 0.0_dp 795 rb = rab 796 CALL overlap_abc_test(la_max(iset), npgfa(iset), zeta(:, iset), la_min(iset), & 797 lb_max(jset), npgfb(jset), zetb(:, jset), lb_min(jset), & 798 lcb_max(kbset), npgfcb(kbset), zetcb(:, kbset), lcb_min(kbset), & 799 ra, rb, rb, sabb, dmax) 800 ENDIF 801 802 CALL abc_contract(abbint(sgfa:m1, sgfb:m2, sgfc:m3), sabb(1:ncoa, 1:ncob, 1:ncoc), & 803 scon_a(:, sgfa:), scon_b(:, sgfb:), scon_cb(:, sgfc:), & 804 ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfcb(kbset)) 805 DEALLOCATE (sabb) 806 ENDIF 807 END DO 808 809 END DO 810 END DO 811 812 CALL timestop(handle) 813 814 END SUBROUTINE int_overlap_abb_os 815 816! ************************************************************************************************** 817!> \brief calculate overlap integrals (aa,bb) 818!> \param saabb integral (aa,bb) 819!> \param oba orbital basis at center A 820!> \param obb orbital basis at center B 821!> \param rab ... 822!> \param debug integrals are debugged by recursive routines if requested 823!> \param dmax maximal deviation between integrals when debugging 824! ************************************************************************************************** 825 SUBROUTINE int_overlap_aabb_os(saabb, oba, obb, rab, debug, dmax) 826 827 REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: saabb 828 TYPE(gto_basis_set_type), POINTER :: oba, obb 829 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 830 LOGICAL, INTENT(IN) :: debug 831 REAL(KIND=dp), INTENT(INOUT) :: dmax 832 833 CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_aabb_os', & 834 routineP = moduleN//':'//routineN 835 836 INTEGER :: handle, iset, isgfa1, jset, jsgfa2, kset, ksgfb1, lds, lset, lsgfb2, m1, m2, m3, & 837 m4, maxco, maxcoa, maxcob, maxl, maxla, maxlb, ncoa1, ncoa2, ncob1, ncob2, nseta, nsetb, & 838 sgfa1, sgfa2, sgfb1, sgfb2 839 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 840 npgfb, nsgfa, nsgfb 841 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 842 LOGICAL :: asets_equal, bsets_equal 843 REAL(KIND=dp) :: dab, ra(3), rb(3) 844 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: swork 845 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: sint 846 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 847 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb 848 849 CALL timeset(routineN, handle) 850 NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, & 851 first_sgfa, first_sgfb, set_radius_a, set_radius_b, rpgfa, rpgfb, & 852 sphi_a, sphi_b, zeta, zetb) 853 854 ! basis ikind 855 first_sgfa => oba%first_sgf 856 la_max => oba%lmax 857 la_min => oba%lmin 858 npgfa => oba%npgf 859 nseta = oba%nset 860 nsgfa => oba%nsgf_set 861 rpgfa => oba%pgf_radius 862 set_radius_a => oba%set_radius 863 sphi_a => oba%sphi 864 zeta => oba%zet 865 ! basis jkind 866 first_sgfb => obb%first_sgf 867 lb_max => obb%lmax 868 lb_min => obb%lmin 869 npgfb => obb%npgf 870 nsetb = obb%nset 871 nsgfb => obb%nsgf_set 872 rpgfb => obb%pgf_radius 873 set_radius_b => obb%set_radius 874 sphi_b => obb%sphi 875 zetb => obb%zet 876 877 CALL get_gto_basis_set(oba, maxco=maxcoa, maxl=maxla) 878 CALL get_gto_basis_set(obb, maxco=maxcob, maxl=maxlb) 879 maxco = MAX(maxcoa, maxcob) 880 maxla = 2*maxla 881 maxlb = 2*maxlb 882 maxl = MAX(maxla, maxlb) 883 lds = ncoset(maxl) 884 ALLOCATE (sint(maxco, maxco, maxco, maxco)) 885 ALLOCATE (swork(lds, lds)) 886 sint = 0._dp 887 swork = 0._dp 888 889 dab = SQRT(SUM(rab**2)) 890 891 DO iset = 1, nseta 892 893 ncoa1 = npgfa(iset)*ncoset(la_max(iset)) 894 sgfa1 = first_sgfa(1, iset) 895 m1 = sgfa1 + nsgfa(iset) - 1 896 897 DO jset = iset, nseta 898 899 ncoa2 = npgfa(jset)*ncoset(la_max(jset)) 900 sgfa2 = first_sgfa(1, jset) 901 m2 = sgfa2 + nsgfa(jset) - 1 902 903 DO kset = 1, nsetb 904 905 ncob1 = npgfb(kset)*ncoset(lb_max(kset)) 906 sgfb1 = first_sgfb(1, kset) 907 m3 = sgfb1 + nsgfb(kset) - 1 908 909 DO lset = kset, nsetb 910 911 ncob2 = npgfb(lset)*ncoset(lb_max(lset)) 912 sgfb2 = first_sgfb(1, lset) 913 m4 = sgfb2 + nsgfb(lset) - 1 914 915 ! check if sets are identical to spare some integral evaluation 916 asets_equal = .FALSE. 917 IF (iset == jset) asets_equal = .TRUE. 918 bsets_equal = .FALSE. 919 IF (kset == lset) bsets_equal = .TRUE. 920 ! calculate integrals 921 CALL overlap_aabb(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 922 la_max(jset), la_min(jset), npgfa(jset), rpgfa(:, jset), zeta(:, jset), & 923 lb_max(kset), lb_min(kset), npgfb(kset), rpgfb(:, kset), zetb(:, kset), & 924 lb_max(lset), lb_min(lset), npgfb(lset), rpgfb(:, lset), zetb(:, lset), & 925 asets_equal, bsets_equal, rab, dab, sint, swork, lds) 926 ! debug if requested 927 IF (debug) THEN 928 ra = 0.0_dp 929 rb = rab 930 CALL overlap_aabb_test(la_max(iset), la_min(iset), npgfa(iset), zeta(:, iset), & 931 la_max(jset), la_min(jset), npgfa(jset), zeta(:, jset), & 932 lb_max(kset), lb_min(kset), npgfb(kset), zetb(:, kset), & 933 lb_max(lset), lb_min(lset), npgfb(lset), zetb(:, lset), & 934 ra, rb, sint, dmax) 935 ENDIF 936 937 CALL abcd_contract(saabb(sgfa1:m1, sgfa2:m2, sgfb1:m3, sgfb2:m4), sint, sphi_a(:, sgfa1:), & 938 sphi_a(:, sgfa2:), sphi_b(:, sgfb1:), sphi_b(:, sgfb2:), ncoa1, ncoa2, & 939 ncob1, ncob2, nsgfa(iset), nsgfa(jset), nsgfb(kset), nsgfb(lset)) 940 941 ! account for the fact that some integrals are alike 942 DO isgfa1 = sgfa1, m1 943 DO jsgfa2 = sgfa2, m2 944 DO ksgfb1 = sgfb1, m3 945 DO lsgfb2 = sgfb2, m4 946 saabb(jsgfa2, isgfa1, ksgfb1, lsgfb2) = saabb(isgfa1, jsgfa2, ksgfb1, lsgfb2) 947 saabb(isgfa1, jsgfa2, lsgfb2, ksgfb1) = saabb(isgfa1, jsgfa2, ksgfb1, lsgfb2) 948 saabb(jsgfa2, isgfa1, lsgfb2, ksgfb1) = saabb(isgfa1, jsgfa2, ksgfb1, lsgfb2) 949 END DO 950 END DO 951 END DO 952 END DO 953 954 END DO 955 END DO 956 END DO 957 END DO 958 959 DEALLOCATE (sint, swork) 960 961 CALL timestop(handle) 962 963 END SUBROUTINE int_overlap_aabb_os 964 965END MODULE generic_os_integrals 966