Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of the moment integrals over Cartesian Gaussian-type
10 : !> functions.
11 : !> \par Literature
12 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
13 : !> \par History
14 : !> none
15 : !> \author J. Hutter (16.02.2005)
16 : ! **************************************************************************************************
17 : MODULE ai_moments
18 :
19 : ! ax,ay,az : Angular momentum index numbers of orbital a.
20 : ! bx,by,bz : Angular momentum index numbers of orbital b.
21 : ! coset : Cartesian orbital set pointer.
22 : ! dab : Distance between the atomic centers a and b.
23 : ! l{a,b} : Angular momentum quantum number of shell a or b.
24 : ! l{a,b}_max: Maximum angular momentum quantum number of shell a or b.
25 : ! l{a,b}_min: Minimum angular momentum quantum number of shell a or b.
26 : ! rac : Distance vector between the atomic center a and reference point c.
27 : ! rbc : Distance vector between the atomic center b and reference point c.
28 : ! rpgf{a,b} : Radius of the primitive Gaussian-type function a or b.
29 : ! zet{a,b} : Exponents of the Gaussian-type functions a or b.
30 : ! zetp : Reciprocal of the sum of the exponents of orbital a and b.
31 :
32 : USE ai_cartesian_operators, ONLY: apply_cartesian_1d, &
33 : apply_cartesian_operator, &
34 : apply_cartesian_operator_pair, &
35 : cartesian_ladder_map, &
36 : cartesian_polynomial_1d
37 : USE ai_os_rr, ONLY: os_rr_ovlp
38 : USE ai_overlap, ONLY: overlap_ab
39 : USE kinds, ONLY: dp
40 : USE mathconstants, ONLY: pi
41 : USE orbital_pointers, ONLY: indco, &
42 : ncoset
43 : #include "../base/base_uses.f90"
44 :
45 : IMPLICIT NONE
46 :
47 : PRIVATE
48 :
49 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_moments'
50 :
51 : PUBLIC :: expikr, moment, diff_momop, contract_expikr, dipole_force
52 :
53 : CONTAINS
54 :
55 : ! **************************************************************************************************
56 : !> \brief Contract complex primitive integrals into real cosine and sine AO blocks.
57 : !> \param cos_block ...
58 : !> \param sin_block ...
59 : !> \param iatom ...
60 : !> \param ncoa ...
61 : !> \param nsgfa ...
62 : !> \param sgfa ...
63 : !> \param sphi_a ...
64 : !> \param ldsa ...
65 : !> \param jatom ...
66 : !> \param ncob ...
67 : !> \param nsgfb ...
68 : !> \param sgfb ...
69 : !> \param sphi_b ...
70 : !> \param ldsb ...
71 : !> \param expab complex primitive plane-wave integrals
72 : !> \param work ...
73 : !> \param ldwork ...
74 : !> \param ordered retain the ordered atom pair instead of mapping to an upper block
75 : ! **************************************************************************************************
76 1539347 : SUBROUTINE contract_expikr(cos_block, sin_block, &
77 3078694 : iatom, ncoa, nsgfa, sgfa, sphi_a, ldsa, &
78 3078694 : jatom, ncob, nsgfb, sgfb, sphi_b, ldsb, &
79 1539347 : expab, work, ldwork, ordered)
80 :
81 : REAL(dp), DIMENSION(:, :), POINTER :: cos_block, sin_block
82 : INTEGER, INTENT(IN) :: iatom, ncoa, nsgfa, sgfa
83 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: sphi_a
84 : INTEGER, INTENT(IN) :: ldsa, jatom, ncob, nsgfb, sgfb
85 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: sphi_b
86 : INTEGER, INTENT(IN) :: ldsb
87 : COMPLEX(dp), DIMENSION(:, :), INTENT(IN) :: expab
88 : REAL(dp), DIMENSION(:, :) :: work
89 : INTEGER, INTENT(IN) :: ldwork
90 : LOGICAL, INTENT(IN), OPTIONAL :: ordered
91 :
92 : LOGICAL :: direct_block
93 1539347 : REAL(dp), ALLOCATABLE :: primitive(:, :)
94 :
95 1539347 : direct_block = iatom <= jatom
96 1539347 : IF (PRESENT(ordered)) direct_block = direct_block .OR. ordered
97 :
98 6157388 : ALLOCATE (primitive(ncoa, ncob))
99 : #:for component in ['cos', 'sin']
100 : #:if component == 'cos'
101 252807510 : primitive(:, :) = REAL(expab(1:ncoa, 1:ncob), KIND=dp)
102 : #:else
103 252807510 : primitive(:, :) = AIMAG(expab(1:ncoa, 1:ncob))
104 : #:endif
105 : CALL dgemm("N", "N", ncoa, nsgfb, ncob, &
106 : 1.0_dp, primitive(1, 1), ncoa, &
107 : sphi_b(1, sgfb), ldsb, &
108 3078694 : 0.0_dp, work(1, 1), ldwork)
109 :
110 3078694 : IF (direct_block) THEN
111 : CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, &
112 : 1.0_dp, sphi_a(1, sgfa), ldsa, &
113 : work(1, 1), ldwork, &
114 : 1.0_dp, ${component}$_block(sgfa, sgfb), &
115 1914320 : SIZE(${component}$_block, 1))
116 : ELSE
117 : CALL dgemm("T", "N", nsgfb, nsgfa, ncoa, &
118 : 1.0_dp, work(1, 1), ldwork, &
119 : sphi_a(1, sgfa), ldsa, &
120 : 1.0_dp, ${component}$_block(sgfb, sgfa), &
121 1164374 : SIZE(${component}$_block, 1))
122 : END IF
123 : #:endfor
124 :
125 1539347 : END SUBROUTINE contract_expikr
126 :
127 : ! **************************************************************************************************
128 : !> \brief Evaluate primitive exp(i k.r) integrals through the complex overlap recurrence.
129 : !> \param la_max_set highest bra shell
130 : !> \param npgfa number of bra primitives
131 : !> \param zeta bra Gaussian exponents
132 : !> \param rpgfa bra primitive radii
133 : !> \param la_min_set lowest bra shell
134 : !> \param lb_max highest ket shell
135 : !> \param npgfb number of ket primitives
136 : !> \param zetb ket Gaussian exponents
137 : !> \param rpgfb ket primitive radii
138 : !> \param lb_min lowest ket shell
139 : !> \param rac bra-center coordinates relative to the phase origin
140 : !> \param rbc ket-center coordinates relative to the phase origin
141 : !> \param kvec plane-wave vector
142 : !> \param expab primitive integrals in the full Cartesian layout, zero below the requested shells
143 : !> \param dexpab optional bra-center derivatives with the plane-wave operator held fixed
144 : !> \param dexpba optional ket-center derivatives in the same bra/ket layout
145 : !> \param momentab optional coordinate-weighted integrals relative to the phase origin
146 : !> Derivative and moment outputs cover the full Cartesian range.
147 : ! **************************************************************************************************
148 7770705 : SUBROUTINE expikr(la_max_set, npgfa, zeta, rpgfa, la_min_set, &
149 1554141 : lb_max, npgfb, zetb, rpgfb, lb_min, rac, rbc, kvec, expab, dexpab, dexpba, momentab)
150 : INTEGER, INTENT(IN) :: la_max_set, npgfa
151 : REAL(KIND=dp), INTENT(IN) :: zeta(:), rpgfa(:)
152 : INTEGER, INTENT(IN) :: la_min_set, lb_max, npgfb
153 : REAL(KIND=dp), INTENT(IN) :: zetb(:), rpgfb(:)
154 : INTEGER, INTENT(IN) :: lb_min
155 : REAL(KIND=dp), INTENT(IN) :: rac(3), rbc(3), kvec(3)
156 : COMPLEX(KIND=dp), INTENT(INOUT) :: expab(:, :)
157 : COMPLEX(KIND=dp), INTENT(INOUT), OPTIONAL, TARGET :: dexpab(:, :, :), dexpba(:, :, :), &
158 : momentab(:, :, :)
159 :
160 : COMPLEX(KIND=dp) :: prefactor, zrap(3), zrbp(3)
161 1554141 : COMPLEX(KIND=dp), ALLOCATABLE :: opr_rr(:, :, :), rr(:, :, :)
162 1554141 : COMPLEX(KIND=dp), POINTER :: opr_values(:, :, :)
163 : INTEGER :: avec(3), bvec(3), coa, cob, ia, ib, &
164 : idir, ipgf, jdir, jpgf, kdir, ldrr, &
165 : lma, lmb, ma, mb, na, nb, op
166 : LOGICAL :: derivative, right
167 : REAL(KIND=dp) :: dab, exponent, f0, f1, f2, k2, kdp, &
168 : origin(3), rab(3), rab2, zet, zetp
169 :
170 6216564 : rab = rbc - rac
171 6216564 : dab = NORM2(rab)
172 6216564 : rab2 = SUM(rab**2)
173 6216564 : k2 = SUM(kvec**2)
174 1554141 : lma = la_max_set
175 1554141 : lmb = lb_max
176 1554141 : IF (PRESENT(dexpab) .OR. PRESENT(momentab)) lma = lma + 1
177 1554141 : IF (PRESENT(dexpba)) lmb = lmb + 1
178 1554141 : ldrr = MAX(lma, lmb) + 1
179 7770705 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
180 1554141 : IF (PRESENT(dexpab) .OR. PRESENT(dexpba) .OR. PRESENT(momentab)) THEN
181 71800 : ALLOCATE (opr_rr(0:la_max_set, 0:lb_max, 3))
182 : END IF
183 :
184 1554141 : na = ncoset(la_max_set)
185 1554141 : nb = ncoset(lb_max)
186 253779848 : expab(1:na*npgfa, 1:nb*npgfb) = 0.0_dp
187 4352511 : IF (PRESENT(dexpab)) dexpab(1:na*npgfa, 1:nb*npgfb, 1:3) = 0.0_dp
188 4333107 : IF (PRESENT(dexpba)) dexpba(1:na*npgfa, 1:nb*npgfb, 1:3) = 0.0_dp
189 1641681 : IF (PRESENT(momentab)) momentab(1:na*npgfa, 1:nb*npgfb, 1:3) = 0.0_dp
190 :
191 5241662 : DO ipgf = 1, npgfa
192 3687521 : ma = (ipgf - 1)*na
193 17279233 : DO jpgf = 1, npgfb
194 12037571 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) CYCLE
195 4302465 : mb = (jpgf - 1)*nb
196 4302465 : zet = zeta(ipgf) + zetb(jpgf)
197 4302465 : zetp = 1.0_dp/zet
198 4302465 : f0 = (pi*zetp)**1.5_dp
199 4302465 : f1 = zetb(jpgf)*zetp
200 4302465 : f2 = 0.5_dp*zetp
201 17209860 : kdp = zetp*DOT_PRODUCT(kvec, zeta(ipgf)*rac + zetb(jpgf)*rbc)
202 4302465 : prefactor = f0*EXP(CMPLX(-zeta(ipgf)*f1*rab2 - 0.25_dp*k2*zetp, kdp, KIND=dp))
203 :
204 : ! The normalized 1D tables carry no three-dimensional prefactor.
205 17209860 : zrap = CMPLX(f1*rab, f2*kvec, KIND=dp)
206 17209860 : zrbp = zrap - rab
207 4302465 : CALL os_rr_ovlp(zrap, lma, zrbp, lmb, zet, ldrr, rr)
208 18511242 : DO cob = ncoset(lb_min - 1) + 1, nb
209 56835108 : bvec = indco(:, cob)
210 80673102 : DO coa = ncoset(la_min_set - 1) + 1, na
211 248647440 : avec = indco(:, coa)
212 : expab(ma + coa, mb + cob) = prefactor*rr(avec(1), bvec(1), 1)* &
213 76370637 : rr(avec(2), bvec(2), 2)*rr(avec(3), bvec(3), 3)
214 : END DO
215 : END DO
216 :
217 20897381 : DO op = 1, 3
218 4302465 : SELECT CASE (op)
219 : CASE (1)
220 4302465 : IF (.NOT. PRESENT(dexpab)) CYCLE
221 31092 : opr_values => dexpab
222 31092 : exponent = zeta(ipgf)
223 : CASE (2)
224 4302465 : IF (.NOT. PRESENT(dexpba)) CYCLE
225 30152 : opr_values => dexpba
226 30152 : exponent = zetb(jpgf)
227 : CASE (3)
228 4302465 : IF (.NOT. PRESENT(momentab)) CYCLE
229 4872 : opr_values => momentab
230 12912267 : exponent = 0.0_dp
231 : END SELECT
232 66116 : right = op == 2
233 66116 : derivative = op /= 3
234 66116 : origin = 0.0_dp
235 66116 : IF (.NOT. derivative) origin = rac
236 66116 : CALL apply_cartesian_1d(rr, opr_rr, la_max_set, lb_max, right, derivative, exponent, origin)
237 12258531 : DO cob = 1, nb
238 619376 : bvec = indco(:, cob)
239 154844 : ib = mb + cob
240 13530859 : DO coa = 1, na
241 1874480 : avec = indco(:, coa)
242 468620 : ia = ma + coa
243 2029324 : DO idir = 1, 3
244 1405860 : jdir = MOD(idir, 3) + 1
245 1405860 : kdir = MOD(idir + 1, 3) + 1
246 : opr_values(ia, ib, idir) = prefactor*opr_rr(avec(idir), bvec(idir), idir)* &
247 1874480 : rr(avec(jdir), bvec(jdir), jdir)*rr(avec(kdir), bvec(kdir), kdir)
248 : END DO
249 : END DO
250 : END DO
251 : END DO
252 : END DO
253 : END DO
254 3108282 : END SUBROUTINE expikr
255 :
256 : ! **************************************************************************************************
257 : !> \brief ...
258 : !> \param la_max ...
259 : !> \param npgfa ...
260 : !> \param zeta ...
261 : !> \param rpgfa ...
262 : !> \param la_min ...
263 : !> \param lb_max ...
264 : !> \param npgfb ...
265 : !> \param zetb ...
266 : !> \param rpgfb ...
267 : !> \param lc_max ...
268 : !> \param rac ...
269 : !> \param rbc ...
270 : !> \param mab ...
271 : !> \param overlap optional zeroth moment in primitive-major cumulative Cartesian order;
272 : !> entries outside the requested row and column ranges are untouched
273 : ! **************************************************************************************************
274 623512 : SUBROUTINE moment(la_max, npgfa, zeta, rpgfa, la_min, &
275 1247024 : lb_max, npgfb, zetb, rpgfb, lc_max, rac, rbc, mab, overlap)
276 : INTEGER, INTENT(IN) :: la_max, npgfa
277 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
278 : INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
279 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
280 : INTEGER, INTENT(IN) :: lc_max
281 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
282 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: mab
283 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: overlap
284 :
285 : INTEGER :: ax, ay, az, bx, by, bz, coa, cob, i, ia, &
286 : ib, imom, ipgf, jpgf, ldrr, lx, ly, &
287 623512 : lz, ma, mb, powers(0:lc_max)
288 : REAL(KIND=dp) :: a, b, dab, f0, rab(3), rab2, rap(3), &
289 : rbp(3), zet
290 : REAL(KIND=dp), ALLOCATABLE :: polynomial(:, :, :, :), rr(:, :, :)
291 :
292 : MARK_USED(la_min)
293 2494048 : rab = rbc - rac
294 2494048 : dab = NORM2(rab)
295 623512 : rab2 = dab**2
296 623512 : ldrr = MAX(la_max + lc_max, lb_max) + 1
297 6235120 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3), polynomial(0:la_max, 0:lb_max, 0:lc_max, 3))
298 4268746 : powers = [(i, i=0, lc_max)]
299 2250181978 : mab(:npgfa*ncoset(la_max), :npgfb*ncoset(lb_max), :ncoset(lc_max) - 1) = 0.0_dp
300 720778 : IF (PRESENT(overlap)) overlap(:npgfa*ncoset(la_max), :npgfb*ncoset(lb_max)) = 0.0_dp
301 2014879 : DO ipgf = 1, npgfa
302 1391367 : ma = (ipgf - 1)*ncoset(la_max)
303 5762519 : DO jpgf = 1, npgfb
304 3747640 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) CYCLE
305 1449379 : mb = (jpgf - 1)*ncoset(lb_max)
306 1449379 : a = zeta(ipgf)
307 1449379 : b = zetb(jpgf)
308 1449379 : zet = a + b
309 5797516 : rap = b*rab/zet
310 5797516 : rbp = -a*rab/zet
311 1449379 : f0 = (pi/zet)**1.5_dp*EXP(-a*b*rab2/zet)
312 1449379 : CALL os_rr_ovlp(rap, la_max + lc_max, rbp, lb_max, zet, ldrr, rr)
313 1449379 : CALL cartesian_polynomial_1d(rr, polynomial, la_max, lb_max, powers, rac)
314 :
315 12979770 : DO cob = 1, ncoset(lb_max)
316 10139024 : bx = indco(1, cob)
317 10139024 : by = indco(2, cob)
318 10139024 : bz = indco(3, cob)
319 10139024 : ib = mb + cob
320 101068488 : DO coa = 1, ncoset(la_max)
321 87181824 : ax = indco(1, coa)
322 87181824 : ay = indco(2, coa)
323 87181824 : az = indco(3, coa)
324 87181824 : ia = ma + coa
325 87181824 : IF (PRESENT(overlap)) THEN
326 : overlap(ia, ib) = f0*polynomial(ax, bx, 0, 1)* &
327 52308 : polynomial(ay, by, 0, 2)*polynomial(az, bz, 0, 3)
328 : END IF
329 : ! Preserve the coset ordering and omission of the zeroth moment.
330 844498237 : DO imom = 1, ncoset(lc_max) - 1
331 747177389 : lx = indco(1, imom + 1)
332 747177389 : ly = indco(2, imom + 1)
333 747177389 : lz = indco(3, imom + 1)
334 : mab(ia, ib, imom) = f0*polynomial(ax, bx, lx, 1)* &
335 834359213 : polynomial(ay, by, ly, 2)*polynomial(az, bz, lz, 3)
336 : END DO
337 : END DO
338 : END DO
339 : END DO
340 : END DO
341 623512 : END SUBROUTINE moment
342 :
343 : ! **************************************************************************************************
344 : !> \brief This returns the derivative of the moment integrals [a|\mu|b].
345 : !> By default, it differentiates the primitive on the right:
346 : !> [a|\mu|d/dR_bi] = 2*zetb*[a|\mu|b+1i] - Ni(b)[a|\mu|b-1i]
347 : !> A weighted derivative combines the left and right primitive derivatives
348 : !> using deltaR for the corresponding atom centers.
349 : !> order indicates the max order of the moment operator to be calculated
350 : !> 1: dipole
351 : !> 2: quadrupole
352 : !> ...
353 : !> \param la_max ...
354 : !> \param npgfa ...
355 : !> \param zeta ...
356 : !> \param rpgfa ...
357 : !> \param la_min ...
358 : !> \param lb_max ...
359 : !> \param npgfb ...
360 : !> \param zetb ...
361 : !> \param rpgfb ...
362 : !> \param lb_min ...
363 : !> \param order ...
364 : !> \param rac ...
365 : !> \param rbc ...
366 : !> \param difmab ...
367 : !> \param mab_ext optional precomputed moments with zero screened primitive blocks, including raised shells
368 : !> \param deltaR optional weights for the left and right primitive derivatives
369 : !> \param lambda optional atom selector for the factor (iatom == lambda) - (jatom == lambda)
370 : !> \param iatom atom associated with the left basis function
371 : !> \param jatom atom associated with the right basis function
372 : !> \note
373 : ! **************************************************************************************************
374 575127 : SUBROUTINE diff_momop(la_max, npgfa, zeta, rpgfa, la_min, &
375 575127 : lb_max, npgfb, zetb, rpgfb, lb_min, &
376 575127 : order, rac, rbc, difmab, mab_ext, deltaR, lambda, iatom, jatom)
377 :
378 : INTEGER, INTENT(IN) :: la_max, npgfa
379 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
380 : INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
381 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
382 : INTEGER, INTENT(IN) :: lb_min, order
383 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
384 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(OUT) :: difmab
385 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
386 : POINTER :: mab_ext
387 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: deltaR
388 : INTEGER, INTENT(IN), OPTIONAL :: lambda, iatom, jatom
389 :
390 : INTEGER :: ider, imom, lda, lda_min, ldb, ldb_min
391 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: difmab_tmp
392 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mab
393 :
394 575127 : lda_min = MAX(0, la_min - 1)
395 575127 : ldb_min = MAX(0, lb_min - 1)
396 575127 : lda = ncoset(la_max)*npgfa
397 575127 : ldb = ncoset(lb_max)*npgfb
398 2848803 : ALLOCATE (difmab_tmp(lda, ldb, 3))
399 :
400 575127 : IF (PRESENT(mab_ext)) THEN
401 575127 : mab => mab_ext
402 : ELSE
403 : ALLOCATE (mab(npgfa*ncoset(la_max + 1), npgfb*ncoset(lb_max + 1), &
404 0 : ncoset(order) - 1))
405 0 : mab = 0.0_dp
406 : ! *** Calculate the primitive overlap integrals ***
407 : CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
408 : lb_max + 1, npgfb, zetb, rpgfb, &
409 0 : order, rac, rbc, mab)
410 :
411 : END IF
412 5731812 : DO imom = 1, ncoset(order) - 1
413 : CALL apply_cartesian_operator(la_max, 0, npgfa, lb_max, 0, npgfb, mab(:, :, imom), &
414 5156685 : difmab_tmp, right=.TRUE., derivative=.TRUE., exponents=zetb)
415 5731812 : IF (PRESENT(deltaR)) THEN
416 6804 : CPASSERT(ASSOCIATED(deltaR))
417 6804 : CPASSERT(PRESENT(iatom) .AND. PRESENT(jatom))
418 27216 : DO ider = 1, 3
419 : difmab(1:lda, 1:ldb, imom, ider) = &
420 1496880 : difmab_tmp(1:lda, 1:ldb, ider)*deltaR(ider, jatom)
421 : END DO
422 :
423 : CALL apply_cartesian_operator(la_max, 0, npgfa, lb_max, 0, npgfb, mab(:, :, imom), &
424 6804 : difmab_tmp, right=.FALSE., derivative=.TRUE., exponents=zeta)
425 27216 : DO ider = 1, 3
426 : difmab(1:lda, 1:ldb, imom, ider) = difmab(1:lda, 1:ldb, imom, ider) &
427 1496880 : + difmab_tmp(1:lda, 1:ldb, ider)*deltaR(ider, iatom)
428 : END DO
429 : ELSE
430 1124774136 : difmab(1:lda, 1:ldb, imom, :) = difmab_tmp(1:lda, 1:ldb, :)
431 : END IF
432 : END DO
433 :
434 575127 : IF (PRESENT(lambda)) THEN
435 27 : CPASSERT(.NOT. PRESENT(deltaR))
436 27 : CPASSERT(PRESENT(iatom) .AND. PRESENT(jatom))
437 27 : IF (iatom == lambda .AND. jatom == lambda) THEN
438 7383 : difmab = 0.0_dp
439 24 : ELSE IF (iatom == lambda) THEN
440 : ! The right-hand derivative is selected with a positive sign.
441 18 : ELSE IF (jatom == lambda) THEN
442 14766 : difmab = -difmab
443 : ELSE
444 29532 : difmab = 0.0_dp
445 : END IF
446 : END IF
447 :
448 575127 : IF (PRESENT(mab_ext)) THEN
449 : NULLIFY (mab)
450 : ELSE
451 0 : DEALLOCATE (mab)
452 : END IF
453 575127 : DEALLOCATE (difmab_tmp)
454 :
455 575127 : END SUBROUTINE diff_momop
456 :
457 : ! **************************************************************************************************
458 : !> \brief This returns the derivative of the dipole integrals [a|x|b], with respect
459 : !> to the position of the primitive on the left and right, i.e.
460 : !> [da/dR_ai|\mu|b] = 2*zeta*[a+1i|\mu|b] - Ni(a)[a-1i|\mu|b]
461 : !> \param la_max ...
462 : !> \param npgfa ...
463 : !> \param zeta ...
464 : !> \param rpgfa ...
465 : !> \param la_min ...
466 : !> \param lb_max ...
467 : !> \param npgfb ...
468 : !> \param zetb ...
469 : !> \param rpgfb ...
470 : !> \param lb_min ...
471 : !> \param order ...
472 : !> \param rac ...
473 : !> \param rbc ...
474 : !> \param pab ...
475 : !> \param forcea ...
476 : !> \param forceb ...
477 : !> \note
478 : ! **************************************************************************************************
479 2100 : SUBROUTINE dipole_force(la_max, npgfa, zeta, rpgfa, la_min, &
480 2100 : lb_max, npgfb, zetb, rpgfb, lb_min, &
481 2100 : order, rac, rbc, pab, forcea, forceb)
482 :
483 : INTEGER, INTENT(IN) :: la_max, npgfa
484 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
485 : INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
486 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
487 : INTEGER, INTENT(IN) :: lb_min, order
488 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
489 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: pab
490 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), TARGET :: forcea, forceb
491 :
492 : INTEGER, PARAMETER :: components(2, 9) = RESHAPE([1, 1, 2, 1, 3, 1, 1, 2, 2, 2, 3, 2, 1, 3, 2, 3, 3, 3], [2, 9])
493 : INTEGER :: axis, first_a, first_b, icomp, ipgf, jpgf, ma, mb, nrow, ncol, side
494 4200 : INTEGER :: map_a(4, ncoset(la_max) - ncoset(la_min - 1), 3), &
495 4200 : map_b(4, ncoset(lb_max) - ncoset(lb_min - 1), 3)
496 4200 : REAL(KIND=dp) :: origin(3, 2), rab(3), distance, rr_work(3*(MAX(la_max, lb_max) + 2)**2)
497 2100 : REAL(KIND=dp), ALLOCATABLE :: source(:, :), work(:, :), derivatives(:, :, :)
498 2100 : REAL(KIND=dp), POINTER :: forces(:, :)
499 :
500 0 : CPASSERT(order == 1)
501 : MARK_USED(order)
502 2100 : first_a = ncoset(la_min - 1) + 1
503 2100 : first_b = ncoset(lb_min - 1) + 1
504 2100 : nrow = ncoset(la_max) - first_a + 1
505 2100 : ncol = ncoset(lb_max) - first_b + 1
506 : ALLOCATE (source(ncoset(la_max + 1), ncoset(lb_max + 1)), work(ncoset(la_max + 1), ncoset(lb_max + 1)), &
507 21000 : derivatives(nrow, ncol, 9))
508 8400 : DO axis = 1, 3
509 6300 : CALL cartesian_ladder_map(la_max, la_min, axis, map_a(:, :, axis), source_lmin=0)
510 8400 : CALL cartesian_ladder_map(lb_max, lb_min, axis, map_b(:, :, axis), source_lmin=0)
511 : END DO
512 2100 : origin = 0.0_dp
513 8400 : rab = rbc - rac
514 8400 : distance = NORM2(rab)
515 8036 : DO ipgf = 1, npgfa
516 5936 : ma = (ipgf - 1)*ncoset(la_max)
517 25033 : DO jpgf = 1, npgfb
518 16997 : IF (rpgfa(ipgf) + rpgfb(jpgf) < distance) CYCLE
519 13411 : mb = (jpgf - 1)*ncoset(lb_max)
520 : CALL overlap_ab(la_max + 1, 0, 1, rpgfa(ipgf:ipgf), zeta(ipgf:ipgf), &
521 : lb_max + 1, 0, 1, rpgfb(jpgf:jpgf), zetb(jpgf:jpgf), rab, &
522 13411 : sab=source, rr_work=rr_work)
523 : ! Multiply by the coordinate on the side opposite the differentiated Gaussian.
524 46169 : DO side = 1, 2
525 26822 : IF (side == 1) THEN
526 13411 : forces => forcea
527 53644 : origin(:, 1) = rbc
528 : CALL apply_cartesian_operator_pair(source, derivatives, map_b, map_a, components, &
529 : [.FALSE., .TRUE.], [0.0_dp, zeta(ipgf)], &
530 40233 : [.TRUE., .FALSE.], work(:, 1:ncol), origin)
531 : ELSE
532 13411 : forces => forceb
533 53644 : origin(:, 1) = rac
534 : CALL apply_cartesian_operator_pair(source, derivatives, map_a, map_b, components, &
535 : [.FALSE., .TRUE.], [0.0_dp, zetb(jpgf)], &
536 40233 : [.FALSE., .TRUE.], work(1:nrow, :), origin)
537 : END IF
538 285217 : DO icomp = 1, 9
539 : forces(components(1, icomp), components(2, icomp)) = &
540 : forces(components(1, icomp), components(2, icomp)) + &
541 2577998 : SUM(pab(ma + first_a:ma + ncoset(la_max), mb + first_b:mb + ncoset(lb_max))*derivatives(:, :, icomp))
542 : END DO
543 : END DO
544 : END DO
545 : END DO
546 :
547 4200 : END SUBROUTINE dipole_force
548 :
549 : END MODULE ai_moments
|
