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 commutator of kinetic energy and position operator 8!> \par History 9!> JGH: from qs_kinetic 10!> \author Juerg Hutter 11! ************************************************************************************************** 12MODULE commutator_rkinetic 13 USE ai_contraction, ONLY: block_add,& 14 contraction 15 USE ai_kinetic, ONLY: kinetic 16 USE basis_set_types, ONLY: gto_basis_set_p_type,& 17 gto_basis_set_type 18 USE dbcsr_api, ONLY: dbcsr_get_block_p,& 19 dbcsr_p_type 20 USE kinds, ONLY: dp 21 USE orbital_pointers, ONLY: coset,& 22 ncoset 23 USE qs_integral_utils, ONLY: basis_set_list_setup,& 24 get_memory_usage 25 USE qs_kind_types, ONLY: qs_kind_type 26 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 27 get_neighbor_list_set_p,& 28 neighbor_list_iterate,& 29 neighbor_list_iterator_create,& 30 neighbor_list_iterator_p_type,& 31 neighbor_list_iterator_release,& 32 neighbor_list_set_p_type 33 34!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num 35#include "./base/base_uses.f90" 36 37 IMPLICIT NONE 38 39 PRIVATE 40 41! *** Global parameters *** 42 43 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'commutator_rkinetic' 44 45! *** Public subroutines *** 46 47 PUBLIC :: build_com_tr_matrix 48 49CONTAINS 50 51! ************************************************************************************************** 52!> \brief Calculation of commutator [T,r] over Cartesian Gaussian functions. 53!> \param matrix_tr ... 54!> \param qs_kind_set ... 55!> \param basis_type basis set to be used 56!> \param sab_nl pair list (must be consistent with basis sets!) 57!> \date 11.10.2010 58!> \par History 59!> Ported from qs_overlap, replaces code in build_core_hamiltonian 60!> Refactoring [07.2014] JGH 61!> Simplify options and use new kinetic energy integral routine 62!> Adapted from qs_kinetic [07.2016] 63!> \author JGH 64!> \version 1.0 65! ************************************************************************************************** 66 SUBROUTINE build_com_tr_matrix(matrix_tr, qs_kind_set, basis_type, sab_nl) 67 68 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_tr 69 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 70 CHARACTER(LEN=*), INTENT(IN) :: basis_type 71 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 72 POINTER :: sab_nl 73 74 CHARACTER(len=*), PARAMETER :: routineN = 'build_com_tr_matrix', & 75 routineP = moduleN//':'//routineN 76 77 INTEGER :: handle, iatom, icol, ikind, ir, irow, & 78 iset, jatom, jkind, jset, ldsab, ltab, & 79 mepos, ncoa, ncob, nkind, nseta, & 80 nsetb, nthread, sgfa, sgfb 81 INTEGER, DIMENSION(3) :: cell 82 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 83 npgfb, nsgfa, nsgfb 84 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 85 LOGICAL :: do_symmetric, found, trans 86 REAL(KIND=dp) :: rab2, tab 87 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: qab, tkab 88 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: kab 89 REAL(KIND=dp), DIMENSION(3) :: rab 90 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 91 REAL(KIND=dp), DIMENSION(:, :), POINTER :: kx_block, ky_block, kz_block, rpgfa, & 92 rpgfb, scon_a, scon_b, zeta, zetb 93 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list 94 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b 95 TYPE(neighbor_list_iterator_p_type), & 96 DIMENSION(:), POINTER :: nl_iterator 97 98 CALL timeset(routineN, handle) 99 100 nkind = SIZE(qs_kind_set) 101 102 ! check for symmetry 103 CPASSERT(SIZE(sab_nl) > 0) 104 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric) 105 106 ! prepare basis set 107 ALLOCATE (basis_set_list(nkind)) 108 CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set) 109 110 ! *** Allocate work storage *** 111 ldsab = get_memory_usage(qs_kind_set, basis_type) 112 113 nthread = 1 114!$ nthread = omp_get_max_threads() 115 ! Iterate of neighbor list 116 CALL neighbor_list_iterator_create(nl_iterator, sab_nl, nthread=nthread) 117 118!$OMP PARALLEL DEFAULT(NONE) & 119!$OMP SHARED (nthread,ldsab,nl_iterator, do_symmetric,& 120!$OMP ncoset,matrix_tr,basis_set_list) & 121!$OMP PRIVATE (kx_block,ky_block,kz_block,mepos,kab,qab,tab,ikind,jkind,iatom,jatom,rab,cell,& 122!$OMP basis_set_a,basis_set_b,& 123!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, & 124!$OMP zeta, first_sgfb, lb_max, lb_min, ltab, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, tkab, & 125!$OMP zetb, scon_a, scon_b, irow, icol, found, trans, rab2, sgfa, sgfb, iset, jset) 126 127 mepos = 0 128!$ mepos = omp_get_thread_num() 129 130 ALLOCATE (kab(ldsab, ldsab, 3), qab(ldsab, ldsab)) 131 132 DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0) 133 CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, & 134 iatom=iatom, jatom=jatom, r=rab, cell=cell) 135 basis_set_a => basis_set_list(ikind)%gto_basis_set 136 IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE 137 basis_set_b => basis_set_list(jkind)%gto_basis_set 138 IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE 139 ! basis ikind 140 first_sgfa => basis_set_a%first_sgf 141 la_max => basis_set_a%lmax 142 la_min => basis_set_a%lmin 143 npgfa => basis_set_a%npgf 144 nseta = basis_set_a%nset 145 nsgfa => basis_set_a%nsgf_set 146 rpgfa => basis_set_a%pgf_radius 147 set_radius_a => basis_set_a%set_radius 148 scon_a => basis_set_a%scon 149 zeta => basis_set_a%zet 150 ! basis jkind 151 first_sgfb => basis_set_b%first_sgf 152 lb_max => basis_set_b%lmax 153 lb_min => basis_set_b%lmin 154 npgfb => basis_set_b%npgf 155 nsetb = basis_set_b%nset 156 nsgfb => basis_set_b%nsgf_set 157 rpgfb => basis_set_b%pgf_radius 158 set_radius_b => basis_set_b%set_radius 159 scon_b => basis_set_b%scon 160 zetb => basis_set_b%zet 161 162 IF (do_symmetric) THEN 163 IF (iatom <= jatom) THEN 164 irow = iatom 165 icol = jatom 166 ELSE 167 irow = jatom 168 icol = iatom 169 END IF 170 ELSE 171 irow = iatom 172 icol = jatom 173 END IF 174 NULLIFY (kx_block) 175 CALL dbcsr_get_block_p(matrix=matrix_tr(1)%matrix, & 176 row=irow, col=icol, BLOCK=kx_block, found=found) 177 CPASSERT(found) 178 NULLIFY (ky_block) 179 CALL dbcsr_get_block_p(matrix=matrix_tr(2)%matrix, & 180 row=irow, col=icol, BLOCK=ky_block, found=found) 181 CPASSERT(found) 182 NULLIFY (kz_block) 183 CALL dbcsr_get_block_p(matrix=matrix_tr(3)%matrix, & 184 row=irow, col=icol, BLOCK=kz_block, found=found) 185 CPASSERT(found) 186 187 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 188 tab = SQRT(rab2) 189 trans = do_symmetric .AND. (iatom > jatom) 190 191 DO iset = 1, nseta 192 193 ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1)) 194 sgfa = first_sgfa(1, iset) 195 196 DO jset = 1, nsetb 197 198 IF (set_radius_a(iset) + set_radius_b(jset) < tab) CYCLE 199 200 ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1)) 201 sgfb = first_sgfb(1, jset) 202 203 ! calclulate integrals 204 ltab = MAX(npgfa(iset)*ncoset(la_max(iset) + 1), npgfb(jset)*ncoset(lb_max(jset) + 1)) 205 ALLOCATE (tkab(ltab, ltab)) 206 CALL kinetic(la_max(iset) + 1, la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 207 lb_max(jset) + 1, lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), & 208 rab, tkab) 209 ! commutator 210 CALL comab_opr(la_max(iset), npgfa(iset), rpgfa(:, iset), la_min(iset), & 211 lb_max(jset), npgfb(jset), rpgfb(:, jset), lb_min(jset), & 212 tab, tkab, kab) 213 DEALLOCATE (tkab) 214 ! Contraction step 215 DO ir = 1, 3 216 CALL contraction(kab(:, :, ir), qab, ca=scon_a(:, sgfa:), na=ncoa, ma=nsgfa(iset), & 217 cb=scon_b(:, sgfb:), nb=ncob, mb=nsgfb(jset), trans=trans) 218!$OMP CRITICAL(blockadd) 219 SELECT CASE (ir) 220 CASE (1) 221 CALL block_add("IN", qab, nsgfa(iset), nsgfb(jset), kx_block, sgfa, sgfb, trans=trans) 222 CASE (2) 223 CALL block_add("IN", qab, nsgfa(iset), nsgfb(jset), ky_block, sgfa, sgfb, trans=trans) 224 CASE (3) 225 CALL block_add("IN", qab, nsgfa(iset), nsgfb(jset), kz_block, sgfa, sgfb, trans=trans) 226 END SELECT 227!$OMP END CRITICAL(blockadd) 228 END DO 229 230 END DO 231 END DO 232 233 END DO 234 DEALLOCATE (kab, qab) 235!$OMP END PARALLEL 236 CALL neighbor_list_iterator_release(nl_iterator) 237 238 ! Release work storage 239 DEALLOCATE (basis_set_list) 240 241 CALL timestop(handle) 242 243 END SUBROUTINE build_com_tr_matrix 244 245! ************************************************************************************************** 246!> \brief Calculate the commutator [O,r] from the integrals [a|O|b]. 247!> We assume that on input all integrals [a+1|O|b+1] are available. 248!> [a|[O,ri]|b] = [a|O|b+1i] - [a+1i|O|b] 249!> \param la_max ... 250!> \param npgfa ... 251!> \param rpgfa ... 252!> \param la_min ... 253!> \param lb_max ... 254!> \param npgfb ... 255!> \param rpgfb ... 256!> \param lb_min ... 257!> \param dab ... 258!> \param ab ... 259!> \param comabr ... 260!> 261!> \date 25.07.2016 262!> \par Literature 263!> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986) 264!> \par Parameters 265!> - ax,ay,az : Angular momentum index numbers of orbital a. 266!> - bx,by,bz : Angular momentum index numbers of orbital b. 267!> - coset : Cartesian orbital set pointer. 268!> - l{a,b} : Angular momentum quantum number of shell a or b. 269!> - l{a,b}_max: Maximum angular momentum quantum number of shell a or b. 270!> - l{a,b}_min: Minimum angular momentum quantum number of shell a or b. 271!> - ncoset : Number of orbitals in a Cartesian orbital set. 272!> - npgf{a,b} : Degree of contraction of shell a or b. 273!> - rab : Distance vector between the atomic centers a and b. 274!> - rab2 : Square of the distance between the atomic centers a and b. 275!> - rac : Distance vector between the atomic centers a and c. 276!> - rac2 : Square of the distance between the atomic centers a and c. 277!> - rbc : Distance vector between the atomic centers b and c. 278!> - rbc2 : Square of the distance between the atomic centers b and c. 279!> - rpgf{a,b} : Radius of the primitive Gaussian-type function a or b. 280!> - zet{a,b} : Exponents of the Gaussian-type functions a or b. 281!> - zetp : Reciprocal of the sum of the exponents of orbital a and b. 282!> 283!> \author JGH 284!> \version 1.0 285! ************************************************************************************************** 286 SUBROUTINE comab_opr(la_max, npgfa, rpgfa, la_min, lb_max, npgfb, rpgfb, lb_min, & 287 dab, ab, comabr) 288 INTEGER, INTENT(IN) :: la_max, npgfa 289 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa 290 INTEGER, INTENT(IN) :: la_min, lb_max, npgfb 291 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb 292 INTEGER, INTENT(IN) :: lb_min 293 REAL(KIND=dp), INTENT(IN) :: dab 294 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ab 295 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: comabr 296 297 INTEGER :: ax, ay, az, bx, by, bz, coa, coap, & 298 coapx, coapy, coapz, cob, cobp, cobpx, & 299 cobpy, cobpz, ipgf, jpgf, la, lb, na, & 300 nap, nb, nbp, ofa, ofb 301 302 comabr = 0.0_dp 303 304 ofa = ncoset(la_min - 1) 305 ofb = ncoset(lb_min - 1) 306 307 na = 0 308 nap = 0 309 DO ipgf = 1, npgfa 310 nb = 0 311 nbp = 0 312 DO jpgf = 1, npgfb 313 IF (rpgfa(ipgf) + rpgfb(jpgf) > dab) THEN 314 DO la = la_min, la_max 315 DO ax = 0, la 316 DO ay = 0, la - ax 317 az = la - ax - ay 318 coa = na + coset(ax, ay, az) - ofa 319 coap = nap + coset(ax, ay, az) - ofa 320 coapx = nap + coset(ax + 1, ay, az) - ofa 321 coapy = nap + coset(ax, ay + 1, az) - ofa 322 coapz = nap + coset(ax, ay, az + 1) - ofa 323 DO lb = lb_min, lb_max 324 DO bx = 0, lb 325 DO by = 0, lb - bx 326 bz = lb - bx - by 327 cob = nb + coset(bx, by, bz) - ofb 328 cobp = nbp + coset(bx, by, bz) - ofb 329 cobpx = nbp + coset(bx + 1, by, bz) - ofb 330 cobpy = nbp + coset(bx, by + 1, bz) - ofb 331 cobpz = nbp + coset(bx, by, bz + 1) - ofb 332 ! [a|[O,ri]|b] = [a|O|b+1i] - [a+1i|O|b] 333 comabr(coa, cob, 1) = ab(coap, cobpx) - ab(coapx, cobp) 334 comabr(coa, cob, 2) = ab(coap, cobpy) - ab(coapy, cobp) 335 comabr(coa, cob, 3) = ab(coap, cobpz) - ab(coapz, cobp) 336 END DO 337 END DO 338 END DO 339 END DO 340 END DO 341 END DO 342 END IF 343 nb = nb + ncoset(lb_max) - ofb 344 nbp = nbp + ncoset(lb_max + 1) - ofb 345 END DO 346 na = na + ncoset(la_max) - ofa 347 nap = nap + ncoset(la_max + 1) - ofa 348 END DO 349 350 END SUBROUTINE comab_opr 351 352END MODULE commutator_rkinetic 353 354