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