1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5! ************************************************************************************************** 6!> \brief General overlap type integrals containers 7!> \par History 8!> - rewrite of PPNL and OCE integrals 9! ************************************************************************************************** 10MODULE sap_kind_types 11 12 USE ai_moments, ONLY: moment 13 USE ai_overlap, ONLY: overlap 14 USE basis_set_types, ONLY: gto_basis_set_p_type,& 15 gto_basis_set_type 16 USE cell_types, ONLY: cell_type,& 17 pbc 18 USE external_potential_types, ONLY: gth_potential_p_type,& 19 gth_potential_type,& 20 sgp_potential_p_type,& 21 sgp_potential_type 22 USE kinds, ONLY: dp 23 USE orbital_pointers, ONLY: init_orbital_pointers,& 24 nco,& 25 ncoset 26 USE particle_types, ONLY: particle_type 27 USE qs_kind_types, ONLY: get_qs_kind,& 28 get_qs_kind_set,& 29 qs_kind_type 30 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type 31 USE util, ONLY: locate,& 32 sort 33#include "./base/base_uses.f90" 34 35 IMPLICIT NONE 36 37 PRIVATE 38 39 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'sap_kind_types' 40 41 TYPE clist_type 42 INTEGER :: catom, nsgf_cnt 43 INTEGER, DIMENSION(:), POINTER :: sgf_list 44 INTEGER, DIMENSION(3) :: cell 45 LOGICAL :: sgf_soft_only 46 REAL(KIND=dp) :: maxac, maxach 47 REAL(KIND=dp), DIMENSION(3) :: rac 48 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: acint 49 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint 50 END TYPE clist_type 51 52 TYPE alist_type 53 INTEGER :: aatom 54 INTEGER :: nclist 55 TYPE(clist_type), DIMENSION(:), POINTER :: clist 56 END TYPE alist_type 57 58 TYPE sap_int_type 59 INTEGER :: a_kind, p_kind 60 INTEGER :: nalist 61 TYPE(alist_type), DIMENSION(:), POINTER :: alist 62 INTEGER, DIMENSION(:), POINTER :: asort, aindex 63 END TYPE sap_int_type 64 65 PUBLIC :: sap_int_type, clist_type, alist_type, & 66 release_sap_int, get_alist, alist_pre_align_blk, & 67 alist_post_align_blk, sap_sort, build_sap_ints 68 69CONTAINS 70 71!========================================================================================================== 72 73! ************************************************************************************************** 74!> \brief ... 75!> \param sap_int ... 76! ************************************************************************************************** 77 SUBROUTINE release_sap_int(sap_int) 78 79 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int 80 81 CHARACTER(LEN=*), PARAMETER :: routineN = 'release_sap_int', & 82 routineP = moduleN//':'//routineN 83 84 INTEGER :: i, j, k 85 TYPE(clist_type), POINTER :: clist 86 87 CPASSERT(ASSOCIATED(sap_int)) 88 89 DO i = 1, SIZE(sap_int) 90 IF (ASSOCIATED(sap_int(i)%alist)) THEN 91 DO j = 1, SIZE(sap_int(i)%alist) 92 IF (ASSOCIATED(sap_int(i)%alist(j)%clist)) THEN 93 DO k = 1, SIZE(sap_int(i)%alist(j)%clist) 94 clist => sap_int(i)%alist(j)%clist(k) 95 IF (ASSOCIATED(clist%acint)) THEN 96 DEALLOCATE (clist%acint) 97 END IF 98 IF (ASSOCIATED(clist%sgf_list)) THEN 99 DEALLOCATE (clist%sgf_list) 100 END IF 101 IF (ASSOCIATED(clist%achint)) THEN 102 DEALLOCATE (clist%achint) 103 END IF 104 END DO 105 DEALLOCATE (sap_int(i)%alist(j)%clist) 106 END IF 107 END DO 108 DEALLOCATE (sap_int(i)%alist) 109 END IF 110 IF (ASSOCIATED(sap_int(i)%asort)) THEN 111 DEALLOCATE (sap_int(i)%asort) 112 END IF 113 IF (ASSOCIATED(sap_int(i)%aindex)) THEN 114 DEALLOCATE (sap_int(i)%aindex) 115 END IF 116 END DO 117 118 DEALLOCATE (sap_int) 119 120 END SUBROUTINE release_sap_int 121 122! ************************************************************************************************** 123!> \brief ... 124!> \param sap_int ... 125!> \param alist ... 126!> \param atom ... 127! ************************************************************************************************** 128 SUBROUTINE get_alist(sap_int, alist, atom) 129 130 TYPE(sap_int_type), INTENT(IN) :: sap_int 131 TYPE(alist_type), INTENT(OUT), POINTER :: alist 132 INTEGER, INTENT(IN) :: atom 133 134 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_alist', routineP = moduleN//':'//routineN 135 136 INTEGER :: i 137 138 NULLIFY (alist) 139 i = locate(sap_int%asort, atom) 140 IF (i > 0 .AND. i <= SIZE(sap_int%alist)) THEN 141 i = sap_int%aindex(i) 142 alist => sap_int%alist(i) 143 ELSE IF (i == 0) THEN 144 NULLIFY (alist) 145 ELSE 146 CPABORT("") 147 END IF 148 149 END SUBROUTINE get_alist 150 151! ************************************************************************************************** 152!> \brief ... 153!> \param blk_in ... 154!> \param ldin ... 155!> \param blk_out ... 156!> \param ldout ... 157!> \param ilist ... 158!> \param in ... 159!> \param jlist ... 160!> \param jn ... 161! ************************************************************************************************** 162 SUBROUTINE alist_pre_align_blk(blk_in, ldin, blk_out, ldout, ilist, in, jlist, jn) 163 INTEGER, INTENT(IN) :: in, ilist(*), ldout 164 REAL(dp), INTENT(INOUT) :: blk_out(ldout, *) 165 INTEGER, INTENT(IN) :: ldin 166 REAL(dp), INTENT(IN) :: blk_in(ldin, *) 167 INTEGER, INTENT(IN) :: jlist(*), jn 168 169 INTEGER :: i, i0, i1, i2, i3, inn, inn1, j, j0 170 171 inn = MOD(in, 4) 172 inn1 = inn + 1 173 DO j = 1, jn 174 j0 = jlist(j) 175 DO i = 1, inn 176 i0 = ilist(i) 177 blk_out(i, j) = blk_in(i0, j0) 178 ENDDO 179 DO i = inn1, in, 4 180 i0 = ilist(i) 181 i1 = ilist(i + 1) 182 i2 = ilist(i + 2) 183 i3 = ilist(i + 3) 184 blk_out(i, j) = blk_in(i0, j0) 185 blk_out(i + 1, j) = blk_in(i1, j0) 186 blk_out(i + 2, j) = blk_in(i2, j0) 187 blk_out(i + 3, j) = blk_in(i3, j0) 188 ENDDO 189 ENDDO 190 END SUBROUTINE alist_pre_align_blk 191 192! ************************************************************************************************** 193!> \brief ... 194!> \param blk_in ... 195!> \param ldin ... 196!> \param blk_out ... 197!> \param ldout ... 198!> \param ilist ... 199!> \param in ... 200!> \param jlist ... 201!> \param jn ... 202! ************************************************************************************************** 203 SUBROUTINE alist_post_align_blk(blk_in, ldin, blk_out, ldout, ilist, in, jlist, jn) 204 INTEGER, INTENT(IN) :: in, ilist(*), ldout 205 REAL(dp), INTENT(INOUT) :: blk_out(ldout, *) 206 INTEGER, INTENT(IN) :: ldin 207 REAL(dp), INTENT(IN) :: blk_in(ldin, *) 208 INTEGER, INTENT(IN) :: jlist(*), jn 209 210 INTEGER :: i, i0, i1, i2, i3, inn, inn1, j, j0 211 212 inn = MOD(in, 4) 213 inn1 = inn + 1 214 DO j = 1, jn 215 j0 = jlist(j) 216 DO i = 1, inn 217 i0 = ilist(i) 218 blk_out(i0, j0) = blk_out(i0, j0) + blk_in(i, j) 219 ENDDO 220 DO i = inn1, in, 4 221 i0 = ilist(i) 222 i1 = ilist(i + 1) 223 i2 = ilist(i + 2) 224 i3 = ilist(i + 3) 225 blk_out(i0, j0) = blk_out(i0, j0) + blk_in(i, j) 226 blk_out(i1, j0) = blk_out(i1, j0) + blk_in(i + 1, j) 227 blk_out(i2, j0) = blk_out(i2, j0) + blk_in(i + 2, j) 228 blk_out(i3, j0) = blk_out(i3, j0) + blk_in(i + 3, j) 229 ENDDO 230 ENDDO 231 END SUBROUTINE alist_post_align_blk 232 233! ************************************************************************************************** 234!> \brief ... 235!> \param sap_int ... 236! ************************************************************************************************** 237 SUBROUTINE sap_sort(sap_int) 238 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int 239 240 CHARACTER(LEN=*), PARAMETER :: routineN = 'sap_sort', routineP = moduleN//':'//routineN 241 242 INTEGER :: iac, na 243 244! *** Set up a sorting index 245 246!$OMP DO 247 DO iac = 1, SIZE(sap_int) 248 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE 249 na = SIZE(sap_int(iac)%alist) 250 ALLOCATE (sap_int(iac)%asort(na), sap_int(iac)%aindex(na)) 251 sap_int(iac)%asort(1:na) = sap_int(iac)%alist(1:na)%aatom 252 CALL sort(sap_int(iac)%asort, na, sap_int(iac)%aindex) 253 END DO 254!$OMP END DO 255 256 END SUBROUTINE sap_sort 257 258!========================================================================================================== 259 260! ************************************************************************************************** 261!> \brief Calculate overlap and optionally momenta <a|x^n|p> between GTOs and nl. pseudo potential projectors 262!> adapted from core_ppnl.F::build_core_ppnl 263!> \param sap_int allocated in parent routine (nkind*nkind) 264!> \param sap_ppnl ... 265!> \param qs_kind_set ... 266!> \param nder Either number of derivatives or order of moments 267!> \param moment_mode if present and true, moments are calculated instead of derivatives 268!> \param refpoint optionally the reference point for moment calculation 269!> \param particle_set needed if refpoint is present 270!> \param cell needed if refpoint is present 271! ************************************************************************************************** 272 SUBROUTINE build_sap_ints(sap_int, sap_ppnl, qs_kind_set, nder, moment_mode, refpoint, particle_set, cell) 273 TYPE(sap_int_type), DIMENSION(:), INTENT(INOUT), & 274 POINTER :: sap_int 275 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 276 INTENT(IN), POINTER :: sap_ppnl 277 TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), & 278 POINTER :: qs_kind_set 279 INTEGER, INTENT(IN) :: nder 280 LOGICAL, INTENT(IN), OPTIONAL :: moment_mode 281 REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: refpoint 282 TYPE(particle_type), DIMENSION(:), INTENT(IN), & 283 OPTIONAL, POINTER :: particle_set 284 TYPE(cell_type), INTENT(IN), OPTIONAL, POINTER :: cell 285 286 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_sap_ints', routineP = moduleN//':'//routineN 287 288 INTEGER :: first_col, handle, i, iac, iatom, ikind, ilist, iset, j, jneighbor, katom, kkind, & 289 l, lc_max, lc_min, ldai, ldsab, lppnl, maxco, maxder, maxl, maxlgto, maxlppnl, maxppnl, & 290 maxsgf, na, nb, ncoa, ncoc, nkind, nlist, nneighbor, nnl, np, nppnl, nprjc, nseta, nsgfa, & 291 prjc, sgfa, slot 292 INTEGER, DIMENSION(3) :: cell_c 293 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nprj_ppnl, & 294 nsgf_seta 295 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa 296 LOGICAL :: dogth, my_momentmode, my_ref 297 LOGICAL, DIMENSION(0:9) :: is_nonlocal 298 REAL(KIND=dp) :: dac, ppnl_radius 299 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: radp 300 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sab, work 301 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work 302 REAL(KIND=dp), DIMENSION(1) :: rprjc, zetc 303 REAL(KIND=dp), DIMENSION(3) :: ra, rac, raf, rc, rcf, rf 304 REAL(KIND=dp), DIMENSION(:), POINTER :: a_nl, alpha_ppnl, hprj, set_radius_a 305 REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, h_nl, rpgfa, sphi_a, vprj_ppnl, & 306 zeta 307 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: c_nl 308 TYPE(clist_type), POINTER :: clist 309 TYPE(gth_potential_p_type), DIMENSION(:), POINTER :: gpotential 310 TYPE(gth_potential_type), POINTER :: gth_potential 311 TYPE(gto_basis_set_p_type), ALLOCATABLE, & 312 DIMENSION(:) :: basis_set 313 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 314 TYPE(sgp_potential_p_type), DIMENSION(:), POINTER :: spotential 315 TYPE(sgp_potential_type), POINTER :: sgp_potential 316 317 CALL timeset(routineN, handle) 318 319 nkind = SIZE(qs_kind_set) 320 maxder = ncoset(nder) 321 322 ! determine whether moments or derivatives should be calculated 323 my_momentmode = .FALSE. 324 IF (PRESENT(moment_mode)) THEN 325 my_momentmode = moment_mode 326 ENDIF 327 328 my_ref = .FALSE. 329 IF (PRESENT(refpoint)) THEN 330 CPASSERT((PRESENT(cell) .AND. PRESENT(particle_set))) ! need these as well if refpoint is provided 331 rf = refpoint 332 my_ref = .TRUE. 333 ENDIF 334 335 CALL get_qs_kind_set(qs_kind_set, & 336 maxco=maxco, & 337 maxlgto=maxlgto, & 338 maxsgf=maxsgf, & 339 maxlppnl=maxlppnl, & 340 maxppnl=maxppnl) 341 342 maxl = MAX(maxlgto, maxlppnl) 343 CALL init_orbital_pointers(maxl + nder + 1) 344 345 ldsab = MAX(maxco, ncoset(maxlppnl), maxsgf, maxppnl) 346 ldai = ncoset(maxl + nder + 1) 347 348 !set up direct access to basis and potential 349 NULLIFY (gpotential, spotential) 350 ALLOCATE (basis_set(nkind), gpotential(nkind), spotential(nkind)) 351 DO ikind = 1, nkind 352 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) 353 IF (ASSOCIATED(orb_basis_set)) THEN 354 basis_set(ikind)%gto_basis_set => orb_basis_set 355 ELSE 356 NULLIFY (basis_set(ikind)%gto_basis_set) 357 END IF 358 CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, sgp_potential=sgp_potential) 359 NULLIFY (gpotential(ikind)%gth_potential) 360 NULLIFY (spotential(ikind)%sgp_potential) 361 IF (ASSOCIATED(gth_potential)) THEN 362 gpotential(ikind)%gth_potential => gth_potential 363 ELSE IF (ASSOCIATED(sgp_potential)) THEN 364 spotential(ikind)%sgp_potential => sgp_potential 365 END IF 366 END DO 367 368 !allocate sap int 369 DO slot = 1, sap_ppnl(1)%nl_size 370 371 ikind = sap_ppnl(1)%nlist_task(slot)%ikind 372 kkind = sap_ppnl(1)%nlist_task(slot)%jkind 373 iatom = sap_ppnl(1)%nlist_task(slot)%iatom 374 katom = sap_ppnl(1)%nlist_task(slot)%jatom 375 nlist = sap_ppnl(1)%nlist_task(slot)%nlist 376 ilist = sap_ppnl(1)%nlist_task(slot)%ilist 377 nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode 378 379 iac = ikind + nkind*(kkind - 1) 380 IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE 381 IF (.NOT. ASSOCIATED(gpotential(kkind)%gth_potential) .AND. & 382 .NOT. ASSOCIATED(spotential(kkind)%sgp_potential)) CYCLE 383 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN 384 sap_int(iac)%a_kind = ikind 385 sap_int(iac)%p_kind = kkind 386 sap_int(iac)%nalist = nlist 387 ALLOCATE (sap_int(iac)%alist(nlist)) 388 DO i = 1, nlist 389 NULLIFY (sap_int(iac)%alist(i)%clist) 390 sap_int(iac)%alist(i)%aatom = 0 391 sap_int(iac)%alist(i)%nclist = 0 392 END DO 393 END IF 394 IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN 395 sap_int(iac)%alist(ilist)%aatom = iatom 396 sap_int(iac)%alist(ilist)%nclist = nneighbor 397 ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor)) 398 DO i = 1, nneighbor 399 sap_int(iac)%alist(ilist)%clist(i)%catom = 0 400 END DO 401 END IF 402 END DO 403 404 !calculate the overlap integrals <a|pp> 405!$OMP PARALLEL & 406!$OMP DEFAULT (NONE) & 407!$OMP SHARED (basis_set, gpotential, spotential, maxder, ncoset, my_momentmode, & 408!$OMP sap_ppnl, sap_int, nkind, ldsab, ldai, nder, nco, my_ref, cell, particle_set, rf) & 409!$OMP PRIVATE (ikind, kkind, iatom, katom, nlist, ilist, nneighbor, jneighbor, & 410!$OMP cell_c, rac, iac, first_sgfa, la_max, la_min, npgfa, nseta, nsgfa, nsgf_seta, & 411!$OMP slot, sphi_a, zeta, cprj, hprj, lppnl, nppnl, nprj_ppnl, & 412!$OMP clist, iset, ncoa, sgfa, prjc, work, sab, ai_work, nprjc, ppnl_radius, & 413!$OMP ncoc, rpgfa, first_col, vprj_ppnl, i, j, l, dogth, & 414!$OMP set_radius_a, rprjc, dac, lc_max, lc_min, zetc, alpha_ppnl, & 415!$OMP na, nb, np, nnl, is_nonlocal, a_nl, h_nl, c_nl, radp, raf, rcf, ra, rc) 416 417 ! allocate work storage 418 ALLOCATE (sab(ldsab, ldsab*maxder), work(ldsab, ldsab*maxder)) 419 sab = 0.0_dp 420 ALLOCATE (ai_work(ldai, ldai, ncoset(nder + 1))) 421 ai_work = 0.0_dp 422 423!$OMP DO SCHEDULE(GUIDED) 424 ! loop over neighbourlist 425 DO slot = 1, sap_ppnl(1)%nl_size 426 ikind = sap_ppnl(1)%nlist_task(slot)%ikind 427 kkind = sap_ppnl(1)%nlist_task(slot)%jkind 428 iatom = sap_ppnl(1)%nlist_task(slot)%iatom 429 katom = sap_ppnl(1)%nlist_task(slot)%jatom 430 nlist = sap_ppnl(1)%nlist_task(slot)%nlist 431 ilist = sap_ppnl(1)%nlist_task(slot)%ilist 432 nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode 433 jneighbor = sap_ppnl(1)%nlist_task(slot)%inode 434 cell_c(:) = sap_ppnl(1)%nlist_task(slot)%cell(:) 435 rac(1:3) = sap_ppnl(1)%nlist_task(slot)%r(1:3) 436 437 iac = ikind + nkind*(kkind - 1) 438 IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE 439 ! get definition of basis set 440 first_sgfa => basis_set(ikind)%gto_basis_set%first_sgf 441 la_max => basis_set(ikind)%gto_basis_set%lmax 442 la_min => basis_set(ikind)%gto_basis_set%lmin 443 npgfa => basis_set(ikind)%gto_basis_set%npgf 444 nseta = basis_set(ikind)%gto_basis_set%nset 445 nsgfa = basis_set(ikind)%gto_basis_set%nsgf 446 nsgf_seta => basis_set(ikind)%gto_basis_set%nsgf_set 447 rpgfa => basis_set(ikind)%gto_basis_set%pgf_radius 448 set_radius_a => basis_set(ikind)%gto_basis_set%set_radius 449 sphi_a => basis_set(ikind)%gto_basis_set%sphi 450 zeta => basis_set(ikind)%gto_basis_set%zet 451 ! get definition of PP projectors 452 IF (ASSOCIATED(gpotential(kkind)%gth_potential)) THEN 453 ! GTH potential 454 dogth = .TRUE. 455 alpha_ppnl => gpotential(kkind)%gth_potential%alpha_ppnl 456 cprj => gpotential(kkind)%gth_potential%cprj 457 lppnl = gpotential(kkind)%gth_potential%lppnl 458 nppnl = gpotential(kkind)%gth_potential%nppnl 459 nprj_ppnl => gpotential(kkind)%gth_potential%nprj_ppnl 460 ppnl_radius = gpotential(kkind)%gth_potential%ppnl_radius 461 vprj_ppnl => gpotential(kkind)%gth_potential%vprj_ppnl 462 ELSE IF (ASSOCIATED(spotential(kkind)%sgp_potential)) THEN 463 ! SGP potential 464 dogth = .FALSE. 465 nprjc = spotential(kkind)%sgp_potential%nppnl 466 IF (nprjc == 0) CYCLE 467 nnl = spotential(kkind)%sgp_potential%n_nonlocal 468 lppnl = spotential(kkind)%sgp_potential%lmax 469 is_nonlocal = .FALSE. 470 is_nonlocal(0:lppnl) = spotential(kkind)%sgp_potential%is_nonlocal(0:lppnl) 471 a_nl => spotential(kkind)%sgp_potential%a_nonlocal 472 h_nl => spotential(kkind)%sgp_potential%h_nonlocal 473 c_nl => spotential(kkind)%sgp_potential%c_nonlocal 474 ppnl_radius = spotential(kkind)%sgp_potential%ppnl_radius 475 ALLOCATE (radp(nnl)) 476 radp(:) = ppnl_radius 477 cprj => spotential(kkind)%sgp_potential%cprj_ppnl 478 hprj => spotential(kkind)%sgp_potential%vprj_ppnl 479 nppnl = SIZE(cprj, 2) 480 ELSE 481 CYCLE 482 END IF 483 484 IF (my_ref) THEN 485 ra(:) = pbc(particle_set(iatom)%r(:) - rf, cell) + rf 486 rc(:) = pbc(particle_set(katom)%r(:) - rf, cell) + rf 487 raf(:) = ra(:) - rf(:) 488 rcf(:) = rc(:) - rf(:) 489 ELSE 490 raf(:) = -rac 491 rcf(:) = (/0._dp, 0._dp, 0._dp/) 492 ENDIF 493 494 dac = NORM2(rac) 495 clist => sap_int(iac)%alist(ilist)%clist(jneighbor) 496 clist%catom = katom 497 clist%cell = cell_c 498 clist%rac = rac 499 ALLOCATE (clist%acint(nsgfa, nppnl, maxder), & 500 clist%achint(nsgfa, nppnl, maxder)) 501 clist%acint = 0._dp 502 clist%achint = 0._dp 503 clist%nsgf_cnt = 0 504 NULLIFY (clist%sgf_list) 505 506 DO iset = 1, nseta 507 ncoa = npgfa(iset)*ncoset(la_max(iset)) 508 sgfa = first_sgfa(1, iset) 509 IF (dogth) THEN 510 ! GTH potential 511 prjc = 1 512 work = 0._dp 513 DO l = 0, lppnl 514 nprjc = nprj_ppnl(l)*nco(l) 515 IF (nprjc == 0) CYCLE 516 rprjc(1) = ppnl_radius 517 IF (set_radius_a(iset) + rprjc(1) < dac) CYCLE 518 lc_max = l + 2*(nprj_ppnl(l) - 1) 519 lc_min = l 520 zetc(1) = alpha_ppnl(l) 521 ncoc = ncoset(lc_max) 522 523 ! *** Calculate the primitive overlap integrals *** 524 IF (my_momentmode) THEN 525 CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 526 lc_max, lc_min, 1, rprjc, zetc, rac, dac, sab, 0, .FALSE., ai_work, ldai) 527 ELSE 528 CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 529 lc_max, lc_min, 1, rprjc, zetc, rac, dac, sab, nder, .TRUE., ai_work, ldai) 530 ENDIF 531 IF (my_momentmode .AND. nder >= 1) THEN 532 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), & 533 lc_max, 1, zetc, rprjc, nder, raf, rcf, ai_work) 534 ! reduce ai_work to sab 535 na = ncoa 536 np = ncoc 537 DO i = 1, maxder - 1 538 first_col = i*ldsab 539 sab(1:na, first_col + 1:first_col + np) = ai_work(1:na, 1:np, i) 540 ENDDO 541 ENDIF 542 543 ! *** Transformation step projector functions (cartesian->spherical) *** 544 na = ncoa 545 nb = nprjc 546 np = ncoc 547 DO i = 1, maxder 548 first_col = (i - 1)*ldsab 549 work(1:na, first_col + prjc:first_col + prjc + nb - 1) = & 550 MATMUL(sab(1:na, first_col + 1:first_col + np), cprj(1:np, prjc:prjc + nb - 1)) 551 END DO 552 prjc = prjc + nprjc 553 END DO 554 555 na = nsgf_seta(iset) 556 nb = nppnl 557 np = ncoa 558 DO i = 1, maxder 559 first_col = (i - 1)*ldsab + 1 560 561 ! *** Contraction step (basis functions) *** 562 clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = & 563 MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work(1:np, first_col:first_col + nb - 1)) 564 565 ! *** Multiply with interaction matrix(h) *** 566 clist%achint(sgfa:sgfa + na - 1, 1:nb, i) = & 567 MATMUL(clist%acint(sgfa:sgfa + na - 1, 1:nb, i), vprj_ppnl(1:nb, 1:nb)) 568 END DO 569 ELSE 570 ! SGP potential 571 ! *** Calculate the primitive overlap integrals *** 572 IF (my_momentmode) THEN 573 CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 574 lppnl, 0, nnl, radp, a_nl, rac, dac, sab, 0, .FALSE., ai_work, ldai) 575 ELSE 576 CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 577 lppnl, 0, nnl, radp, a_nl, rac, dac, sab, nder, .TRUE., ai_work, ldai) 578 ENDIF 579 IF (my_momentmode .AND. nder >= 1) THEN 580 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), & 581 lppnl, nnl, a_nl, radp, nder, raf, rcf, ai_work) 582 ! reduce ai_work to sab 583 na = ncoa 584 DO i = 1, maxder - 1 585 first_col = i*ldsab 586 sab(1:na, first_col:first_col + nprjc - 1) = ai_work(1:na, 1:nprjc, i) 587 ENDDO 588 ENDIF 589 590 na = nsgf_seta(iset) 591 nb = nppnl 592 np = ncoa 593 DO i = 1, maxder 594 first_col = (i - 1)*ldsab + 1 595 ! *** Transformation step projector functions (cartesian->spherical) *** 596 work(1:np, 1:nb) = MATMUL(sab(1:np, first_col:first_col + nprjc - 1), cprj(1:nprjc, 1:nb)) 597 598 ! *** Contraction step (basis functions) *** 599 clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = & 600 MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work(1:np, 1:nb)) 601 602 ! *** Multiply with interaction matrix(h) *** 603 ncoc = sgfa + nsgf_seta(iset) - 1 604 DO j = 1, nppnl 605 clist%achint(sgfa:ncoc, j, i) = clist%acint(sgfa:ncoc, j, i)*hprj(j) 606 END DO 607 END DO 608 END IF 609 END DO 610 clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1))) 611 clist%maxach = MAXVAL(ABS(clist%achint(:, :, 1))) 612 IF (.NOT. dogth) DEALLOCATE (radp) 613 614 ENDDO 615 616 DEALLOCATE (sab, ai_work, work) 617 618!$OMP END PARALLEL 619 620 DEALLOCATE (basis_set, gpotential, spotential) 621 622 CALL timestop(handle) 623 624 END SUBROUTINE build_sap_ints 625 626END MODULE sap_kind_types 627