1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5! ************************************************************************************************** 6!> \brief Calculation of the non-local pseudopotential contribution to the core Hamiltonian 7!> <a|V(non-local)|b> = <a|p(l,i)>*h(i,j)*<p(l,j)|b> 8!> \par History 9!> - refactered from qs_core_hamiltian [Joost VandeVondele, 2008-11-01] 10!> - full rewrite [jhu, 2009-01-23] 11! ************************************************************************************************** 12MODULE commutator_rpnl 13 USE ai_moments, ONLY: moment 14 USE ai_overlap, ONLY: overlap 15 USE basis_set_types, ONLY: gto_basis_set_p_type,& 16 gto_basis_set_type 17 USE dbcsr_api, ONLY: dbcsr_get_block_p,& 18 dbcsr_p_type 19 USE external_potential_types, ONLY: gth_potential_p_type,& 20 gth_potential_type,& 21 sgp_potential_p_type,& 22 sgp_potential_type 23 USE kinds, ONLY: dp 24 USE orbital_pointers, ONLY: init_orbital_pointers,& 25 nco,& 26 ncoset 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: get_iterator_info,& 31 neighbor_list_iterate,& 32 neighbor_list_iterator_create,& 33 neighbor_list_iterator_p_type,& 34 neighbor_list_iterator_release,& 35 neighbor_list_set_p_type 36 USE sap_kind_types, ONLY: alist_type,& 37 clist_type,& 38 get_alist,& 39 release_sap_int,& 40 sap_int_type,& 41 sap_sort 42 43!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads 44 45#include "./base/base_uses.f90" 46 47 IMPLICIT NONE 48 49 PRIVATE 50 51 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'commutator_rpnl' 52 53 PUBLIC :: build_com_rpnl 54 55CONTAINS 56 57! ************************************************************************************************** 58!> \brief ... 59!> \param matrix_rv ... 60!> \param qs_kind_set ... 61!> \param sab_orb ... 62!> \param sap_ppnl ... 63!> \param eps_ppnl ... 64! ************************************************************************************************** 65 SUBROUTINE build_com_rpnl(matrix_rv, qs_kind_set, sab_orb, sap_ppnl, eps_ppnl) 66 67 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_rv 68 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 69 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 70 POINTER :: sab_orb, sap_ppnl 71 REAL(KIND=dp), INTENT(IN) :: eps_ppnl 72 73 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_com_rpnl', routineP = moduleN//':'//routineN 74 75 INTEGER :: handle, i, iab, iac, iatom, ibc, icol, ikind, ilist, inode, irow, iset, jatom, & 76 jkind, jneighbor, kac, katom, kbc, kkind, l, lc_max, lc_min, ldai, ldsab, lppnl, maxco, & 77 maxder, maxl, maxlgto, maxlppnl, maxppnl, maxsgf, mepos, na, nb, ncoa, ncoc, nkind, & 78 nlist, nneighbor, nnode, np, nppnl, nprjc, nseta, nsgfa, nthread, prjc, sgfa 79 INTEGER, DIMENSION(3) :: cell_b, cell_c 80 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nprj_ppnl, & 81 nsgf_seta 82 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa 83 LOGICAL :: found, gpot, ppnl_present, spot 84 REAL(KIND=dp) :: dac, ppnl_radius 85 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work, sab, work 86 REAL(KIND=dp), DIMENSION(1) :: rprjc, zetc 87 REAL(KIND=dp), DIMENSION(3) :: rab, rac 88 REAL(KIND=dp), DIMENSION(:), POINTER :: alpha_ppnl, set_radius_a 89 REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, rpgfa, sphi_a, vprj_ppnl, x_block, & 90 y_block, z_block, zeta 91 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint 92 TYPE(alist_type), POINTER :: alist_ac, alist_bc 93 TYPE(clist_type), POINTER :: clist 94 TYPE(gth_potential_p_type), DIMENSION(:), POINTER :: gpotential 95 TYPE(gth_potential_type), POINTER :: gth_potential 96 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set 97 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 98 TYPE(neighbor_list_iterator_p_type), & 99 DIMENSION(:), POINTER :: nl_iterator 100 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int 101 TYPE(sgp_potential_p_type), DIMENSION(:), POINTER :: spotential 102 TYPE(sgp_potential_type), POINTER :: sgp_potential 103 104 CALL timeset(routineN, handle) 105 106 ppnl_present = ASSOCIATED(sap_ppnl) 107 108 IF (ppnl_present) THEN 109 110 nkind = SIZE(qs_kind_set) 111 112 CALL get_qs_kind_set(qs_kind_set, & 113 maxco=maxco, & 114 maxlgto=maxlgto, & 115 maxsgf=maxsgf, & 116 maxlppnl=maxlppnl, & 117 maxppnl=maxppnl) 118 119 maxl = MAX(maxlgto, maxlppnl) 120 CALL init_orbital_pointers(maxl + 1) 121 122 ldsab = MAX(maxco, ncoset(maxlppnl), maxsgf, maxppnl) 123 ldai = ncoset(maxl + 1) 124 125 !sap_int needs to be shared as multiple threads need to access this 126 ALLOCATE (sap_int(nkind*nkind)) 127 DO i = 1, nkind*nkind 128 NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex) 129 sap_int(i)%nalist = 0 130 END DO 131 132 !set up direct access to basis and potential 133 ALLOCATE (basis_set(nkind), gpotential(nkind), spotential(nkind)) 134 DO ikind = 1, nkind 135 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) 136 IF (ASSOCIATED(orb_basis_set)) THEN 137 basis_set(ikind)%gto_basis_set => orb_basis_set 138 ELSE 139 NULLIFY (basis_set(ikind)%gto_basis_set) 140 END IF 141 CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, & 142 sgp_potential=sgp_potential) 143 IF (ASSOCIATED(gth_potential)) THEN 144 gpotential(ikind)%gth_potential => gth_potential 145 NULLIFY (spotential(ikind)%sgp_potential) 146 ELSE IF (ASSOCIATED(sgp_potential)) THEN 147 spotential(ikind)%sgp_potential => sgp_potential 148 NULLIFY (gpotential(ikind)%gth_potential) 149 ELSE 150 NULLIFY (gpotential(ikind)%gth_potential) 151 NULLIFY (spotential(ikind)%sgp_potential) 152 END IF 153 END DO 154 155 maxder = 4 156 nthread = 1 157!$ nthread = omp_get_max_threads() 158 159 !calculate the overlap integrals <a|p> 160 CALL neighbor_list_iterator_create(nl_iterator, sap_ppnl, nthread=nthread) 161!$OMP PARALLEL & 162!$OMP DEFAULT (NONE) & 163!$OMP SHARED (nl_iterator, basis_set, spotential, gpotential, maxder, ncoset, & 164!$OMP sap_int, nkind, ldsab, ldai, nco ) & 165!$OMP PRIVATE (mepos, ikind, kkind, iatom, katom, nlist, ilist, nneighbor, jneighbor, & 166!$OMP cell_c, rac, iac, first_sgfa, la_max, la_min, npgfa, nseta, nsgfa, nsgf_seta, & 167!$OMP sphi_a, zeta, cprj, lppnl, nppnl, nprj_ppnl, & 168!$OMP clist, iset, ncoa, sgfa, prjc, work, sab, ai_work, nprjc, ppnl_radius, & 169!$OMP ncoc, rpgfa, vprj_ppnl, i, l, gpot, spot, & 170!$OMP set_radius_a, rprjc, dac, lc_max, lc_min, zetc, alpha_ppnl) 171 mepos = 0 172!$ mepos = omp_get_thread_num() 173 174 ALLOCATE (sab(ldsab, ldsab, maxder), work(ldsab, ldsab, maxder)) 175 sab = 0.0_dp 176 ALLOCATE (ai_work(ldai, ldai, 1)) 177 ai_work = 0.0_dp 178 179 DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0) 180 CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=kkind, iatom=iatom, & 181 jatom=katom, nlist=nlist, ilist=ilist, nnode=nneighbor, inode=jneighbor, cell=cell_c, r=rac) 182 iac = ikind + nkind*(kkind - 1) 183 IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE 184 gpot = ASSOCIATED(gpotential(kkind)%gth_potential) 185 spot = ASSOCIATED(spotential(kkind)%sgp_potential) 186 IF ((.NOT. gpot) .AND. (.NOT. spot)) CYCLE 187 ! get definition of basis set 188 first_sgfa => basis_set(ikind)%gto_basis_set%first_sgf 189 la_max => basis_set(ikind)%gto_basis_set%lmax 190 la_min => basis_set(ikind)%gto_basis_set%lmin 191 npgfa => basis_set(ikind)%gto_basis_set%npgf 192 nseta = basis_set(ikind)%gto_basis_set%nset 193 nsgfa = basis_set(ikind)%gto_basis_set%nsgf 194 nsgf_seta => basis_set(ikind)%gto_basis_set%nsgf_set 195 rpgfa => basis_set(ikind)%gto_basis_set%pgf_radius 196 set_radius_a => basis_set(ikind)%gto_basis_set%set_radius 197 sphi_a => basis_set(ikind)%gto_basis_set%sphi 198 zeta => basis_set(ikind)%gto_basis_set%zet 199 ! get definition of PP projectors 200 IF (gpot) THEN 201 alpha_ppnl => gpotential(kkind)%gth_potential%alpha_ppnl 202 cprj => gpotential(kkind)%gth_potential%cprj 203 lppnl = gpotential(kkind)%gth_potential%lppnl 204 nppnl = gpotential(kkind)%gth_potential%nppnl 205 nprj_ppnl => gpotential(kkind)%gth_potential%nprj_ppnl 206 ppnl_radius = gpotential(kkind)%gth_potential%ppnl_radius 207 vprj_ppnl => gpotential(kkind)%gth_potential%vprj_ppnl 208 ELSEIF (spot) THEN 209 CPABORT('SGP not implemented') 210 ELSE 211 CPABORT('PPNL unknown') 212 END IF 213!$OMP CRITICAL(sap_int_critical) 214 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN 215 sap_int(iac)%a_kind = ikind 216 sap_int(iac)%p_kind = kkind 217 sap_int(iac)%nalist = nlist 218 ALLOCATE (sap_int(iac)%alist(nlist)) 219 DO i = 1, nlist 220 NULLIFY (sap_int(iac)%alist(i)%clist) 221 sap_int(iac)%alist(i)%aatom = 0 222 sap_int(iac)%alist(i)%nclist = 0 223 END DO 224 END IF 225 IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN 226 sap_int(iac)%alist(ilist)%aatom = iatom 227 sap_int(iac)%alist(ilist)%nclist = nneighbor 228 ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor)) 229 DO i = 1, nneighbor 230 sap_int(iac)%alist(ilist)%clist(i)%catom = 0 231 END DO 232 END IF 233!$OMP END CRITICAL(sap_int_critical) 234 dac = SQRT(SUM(rac*rac)) 235 clist => sap_int(iac)%alist(ilist)%clist(jneighbor) 236 clist%catom = katom 237 clist%cell = cell_c 238 clist%rac = rac 239 ALLOCATE (clist%acint(nsgfa, nppnl, maxder), & 240 clist%achint(nsgfa, nppnl, maxder)) 241 clist%acint = 0._dp 242 clist%achint = 0._dp 243 clist%nsgf_cnt = 0 244 NULLIFY (clist%sgf_list) 245 DO iset = 1, nseta 246 ncoa = npgfa(iset)*ncoset(la_max(iset)) 247 sgfa = first_sgfa(1, iset) 248 work = 0._dp 249 prjc = 1 250 DO l = 0, lppnl 251 nprjc = nprj_ppnl(l)*nco(l) 252 IF (nprjc == 0) CYCLE 253 rprjc(1) = ppnl_radius 254 IF (set_radius_a(iset) + rprjc(1) < dac) CYCLE 255 lc_max = l + 2*(nprj_ppnl(l) - 1) 256 lc_min = l 257 zetc(1) = alpha_ppnl(l) 258 ncoc = ncoset(lc_max) 259 ! Calculate the primitive overlap and dipole moment integrals 260 CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 261 lc_max, lc_min, 1, rprjc, zetc, rac, dac, sab(:, :, 1), 0, .FALSE., ai_work, ldai) 262 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), & 263 lc_max, 1, zetc, rprjc, 1, rac, (/0._dp, 0._dp, 0._dp/), sab(:, :, 2:4)) 264 ! *** Transformation step projector functions (cartesian->spherical) *** 265 DO i = 1, maxder 266 CALL dgemm("N", "N", ncoa, nprjc, ncoc, 1.0_dp, sab(1, 1, i), ldsab, & 267 cprj(1, prjc), SIZE(cprj, 1), 0.0_dp, work(1, 1, i), ldsab) 268 END DO 269 prjc = prjc + nprjc 270 END DO 271 DO i = 1, maxder 272 ! Contraction step (basis functions) 273 CALL dgemm("T", "N", nsgf_seta(iset), nppnl, ncoa, 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), & 274 work(1, 1, i), ldsab, 0.0_dp, clist%acint(sgfa, 1, i), nsgfa) 275 ! Multiply with interaction matrix(h) 276 CALL dgemm("N", "N", nsgf_seta(iset), nppnl, nppnl, 1.0_dp, clist%acint(sgfa, 1, i), nsgfa, & 277 vprj_ppnl(1, 1), SIZE(vprj_ppnl, 1), 0.0_dp, clist%achint(sgfa, 1, i), nsgfa) 278 END DO 279 END DO 280 clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1))) 281 clist%maxach = MAXVAL(ABS(clist%achint(:, :, 1))) 282 END DO 283 284 DEALLOCATE (sab, ai_work, work) 285!$OMP END PARALLEL 286 CALL neighbor_list_iterator_release(nl_iterator) 287 288 ! *** Set up a sorting index 289 CALL sap_sort(sap_int) 290 ! *** All integrals needed have been calculated and stored in sap_int 291 ! *** We now calculate the Hamiltonian matrix elements 292 CALL neighbor_list_iterator_create(nl_iterator, sab_orb, nthread=nthread) 293 294!$OMP PARALLEL & 295!$OMP DEFAULT (NONE) & 296!$OMP SHARED (nl_iterator, basis_set, matrix_rv, & 297!$OMP sap_int, nkind, eps_ppnl ) & 298!$OMP PRIVATE (mepos, ikind, jkind, iatom, jatom, nlist, ilist, nnode, inode, cell_b, rab, & 299!$OMP iab, irow, icol, x_block, y_block, z_block, & 300!$OMP found, iac, ibc, alist_ac, alist_bc, acint, bcint, & 301!$OMP achint, bchint, na, np, nb, katom, rac, kkind, kac, kbc, i) 302 303 mepos = 0 304!$ mepos = omp_get_thread_num() 305 306 DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0) 307 CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, & 308 jatom=jatom, nlist=nlist, ilist=ilist, nnode=nnode, inode=inode, cell=cell_b, r=rab) 309 IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE 310 IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE 311 iab = ikind + nkind*(jkind - 1) 312 313 ! *** Create matrix blocks for a new matrix block column *** 314 IF (iatom <= jatom) THEN 315 irow = iatom 316 icol = jatom 317 ELSE 318 irow = jatom 319 icol = iatom 320 END IF 321 CALL dbcsr_get_block_p(matrix_rv(1)%matrix, irow, icol, x_block, found) 322 CALL dbcsr_get_block_p(matrix_rv(2)%matrix, irow, icol, y_block, found) 323 CALL dbcsr_get_block_p(matrix_rv(3)%matrix, irow, icol, z_block, found) 324 325 ! loop over all kinds for projector atom 326 IF (ASSOCIATED(x_block) .AND. ASSOCIATED(y_block) .AND. ASSOCIATED(z_block)) THEN 327 DO kkind = 1, nkind 328 iac = ikind + nkind*(kkind - 1) 329 ibc = jkind + nkind*(kkind - 1) 330 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE 331 IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) CYCLE 332 CALL get_alist(sap_int(iac), alist_ac, iatom) 333 CALL get_alist(sap_int(ibc), alist_bc, jatom) 334 IF (.NOT. ASSOCIATED(alist_ac)) CYCLE 335 IF (.NOT. ASSOCIATED(alist_bc)) CYCLE 336 DO kac = 1, alist_ac%nclist 337 DO kbc = 1, alist_bc%nclist 338 IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE 339 IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN 340 IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) CYCLE 341 acint => alist_ac%clist(kac)%acint 342 bcint => alist_bc%clist(kbc)%acint 343 achint => alist_ac%clist(kac)%achint 344 bchint => alist_bc%clist(kbc)%achint 345 na = SIZE(acint, 1) 346 np = SIZE(acint, 2) 347 nb = SIZE(bcint, 1) 348!$OMP CRITICAL(h_block_critical) 349 IF (iatom <= jatom) THEN 350 ! Vnl*r 351 CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 1), na, & 352 bcint(1, 1, 2), nb, 1.0_dp, x_block, SIZE(x_block, 1)) 353 CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 1), na, & 354 bcint(1, 1, 3), nb, 1.0_dp, y_block, SIZE(y_block, 1)) 355 CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 1), na, & 356 bcint(1, 1, 4), nb, 1.0_dp, z_block, SIZE(z_block, 1)) 357 ! -r*Vnl 358 CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 2), na, & 359 bcint(1, 1, 1), nb, 1.0_dp, x_block, SIZE(x_block, 1)) 360 CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 3), na, & 361 bcint(1, 1, 1), nb, 1.0_dp, y_block, SIZE(y_block, 1)) 362 CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 4), na, & 363 bcint(1, 1, 1), nb, 1.0_dp, z_block, SIZE(z_block, 1)) 364 ELSE 365 ! Vnl*r 366 CALL dgemm("N", "T", nb, na, np, 1.0_dp, bchint(1, 1, 2), nb, & 367 acint(1, 1, 1), na, 1.0_dp, x_block, SIZE(x_block, 1)) 368 CALL dgemm("N", "T", nb, na, np, 1.0_dp, bchint(1, 1, 3), nb, & 369 acint(1, 1, 1), na, 1.0_dp, y_block, SIZE(y_block, 1)) 370 CALL dgemm("N", "T", nb, na, np, 1.0_dp, bchint(1, 1, 4), nb, & 371 acint(1, 1, 1), na, 1.0_dp, z_block, SIZE(z_block, 1)) 372 ! -r*Vnl 373 CALL dgemm("N", "T", nb, na, np, -1.0_dp, bchint(1, 1, 1), nb, & 374 acint(1, 1, 2), na, 1.0_dp, x_block, SIZE(x_block, 1)) 375 CALL dgemm("N", "T", nb, na, np, -1.0_dp, bchint(1, 1, 1), nb, & 376 acint(1, 1, 3), na, 1.0_dp, y_block, SIZE(y_block, 1)) 377 CALL dgemm("N", "T", nb, na, np, -1.0_dp, bchint(1, 1, 1), nb, & 378 acint(1, 1, 4), na, 1.0_dp, z_block, SIZE(z_block, 1)) 379 END IF 380!$OMP END CRITICAL(h_block_critical) 381 EXIT ! We have found a match and there can be only one single match 382 END IF 383 END DO 384 END DO 385 END DO 386 ENDIF 387 END DO 388!$OMP END PARALLEL 389 CALL neighbor_list_iterator_release(nl_iterator) 390 391 CALL release_sap_int(sap_int) 392 393 DEALLOCATE (basis_set, gpotential, spotential) 394 395 END IF !ppnl_present 396 397 CALL timestop(handle) 398 399 END SUBROUTINE build_com_rpnl 400 401END MODULE commutator_rpnl 402