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 Build up the nuclear potential part of the core Hamiltonian matrix
10 : !> in the case of an allelectron calculation by computing the nuclear
11 : !> attraction integral matrix <a|-Z/r|b> and the integral matrix
12 : !> <a|Z*erf(r)/r|b>. The integrals <a|Z*erf(r)/r|b> can be rewritten
13 : !> as the three-center Coulomb integrals <ab||c> with a primitive
14 : !> s-type Gaussian function c multiplied by a factor N.
15 : !>
16 : !> erfc(r)
17 : !> <a|V|b> = -Z*<a|---------|b>
18 : !> r
19 : !>
20 : !> 1 - erf(r)
21 : !> = -Z*<a|------------|b>
22 : !> r
23 : !>
24 : !> 1 erf(r)
25 : !> = -Z*(<a|---|b> - <a|--------|b>)
26 : !> r r
27 : !>
28 : !> 1
29 : !> = -Z*(<a|---|b> - N*<ab||c>)
30 : !> r
31 : !>
32 : !> 1
33 : !> = -Z*<a|---|b> + Z*N*<ab||c>
34 : !> r
35 : !> \_______/ \_____/
36 : !> | |
37 : !> nuclear coulomb
38 : !> \par Literature
39 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
40 : !> \par Parameters
41 : !> - ax,ay,az : Angular momentum index numbers of orbital a.
42 : !> - bx,by,bz : Angular momentum index numbers of orbital b.
43 : !> - coset : Cartesian orbital set pointer.
44 : !> - dab : Distance between the atomic centers a and b.
45 : !> - dac : Distance between the atomic centers a and c.
46 : !> - dbc : Distance between the atomic centers b and c.
47 : !> - l{a,b} : Angular momentum quantum number of shell a or b.
48 : !> - l{a,b}_max: Maximum angular momentum quantum number of shell a or b.
49 : !> - l{a,b}_min: Minimum angular momentum quantum number of shell a or b.
50 : !> - ncoset : Number of orbitals in a Cartesian orbital set.
51 : !> - npgf{a,b} : Degree of contraction of shell a or b.
52 : !> - rab : Distance vector between the atomic centers a and b.
53 : !> - rab2 : Square of the distance between the atomic centers a and b.
54 : !> - rac : Distance vector between the atomic centers a and c.
55 : !> - rac2 : Square of the distance between the atomic centers a and c.
56 : !> - rbc2 : Square of the distance between the atomic centers b and c.
57 : !> - rpgf{a,b} : Radius of the primitive Gaussian-type function a or b.
58 : !> - zet{a,b} : Exponents of the Gaussian-type functions a or b.
59 : !> - zetp : Reciprocal of the sum of the exponents of orbital a and b.
60 : !> - zetw : Reciprocal of the sum of the exponents of orbital a, b and c.
61 : !> \par History
62 : !> 10.2008 added pVp matrix elements (jens)
63 : !> \author Matthias Krack (04.10.2000)
64 : ! **************************************************************************************************
65 :
66 : MODULE ai_verfc
67 :
68 : USE ai_cartesian_operators, ONLY: apply_cartesian_operator_pair,&
69 : cartesian_ladder_map
70 : USE ai_os_rr, ONLY: os_rr_radial
71 : USE gamma, ONLY: fgamma => fgamma_0
72 : USE kinds, ONLY: dp
73 : USE mathconstants, ONLY: pi
74 : USE orbital_pointers, ONLY: ncoset
75 : #include "../base/base_uses.f90"
76 :
77 : IMPLICIT NONE
78 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_verfc'
79 : PRIVATE
80 :
81 : ! *** Public subroutines ***
82 :
83 : PUBLIC :: verfc
84 :
85 : CONTAINS
86 :
87 : ! **************************************************************************************************
88 : !> \brief Calculation of the primitive three-center nuclear potential
89 : !> integrals <a|Z*erfc(r)/r|b> over Cartesian Gaussian-type
90 : !> functions.
91 : !> \param la_max1 ...
92 : !> \param npgfa ...
93 : !> \param zeta ...
94 : !> \param rpgfa ...
95 : !> \param la_min1 ...
96 : !> \param lb_max1 ...
97 : !> \param npgfb ...
98 : !> \param zetb ...
99 : !> \param rpgfb ...
100 : !> \param lb_min1 ...
101 : !> \param zetc ...
102 : !> \param rpgfc ...
103 : !> \param zc ...
104 : !> \param cerf ...
105 : !> \param rab ...
106 : !> \param rab2 ...
107 : !> \param rac ...
108 : !> \param rac2 ...
109 : !> \param rbc2 ...
110 : !> \param vabc ...
111 : !> \param verf ...
112 : !> \param vnuc ...
113 : !> \param f ...
114 : !> \param maxder ...
115 : !> \param vabc_plus ...
116 : !> \param vnabc ...
117 : !> \param pVp_sum ...
118 : !> \param pVp ...
119 : !> \param dkh_erfc ...
120 : !> \date 04.10.2000
121 : !> \author Matthias Krack
122 : !> \version 1.0
123 : ! **************************************************************************************************
124 4975432 : SUBROUTINE verfc(la_max1, npgfa, zeta, rpgfa, la_min1, lb_max1, npgfb, zetb, rpgfb, lb_min1, &
125 14926296 : zetc, rpgfc, zc, cerf, rab, rab2, rac, rac2, rbc2, vabc, verf, vnuc, f, &
126 14926296 : maxder, vabc_plus, vnabc, pVp_sum, pVp, dkh_erfc) !JT
127 : INTEGER, INTENT(IN) :: la_max1, npgfa
128 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
129 : INTEGER, INTENT(IN) :: la_min1, lb_max1, npgfb
130 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
131 : INTEGER, INTENT(IN) :: lb_min1
132 : REAL(KIND=dp), INTENT(IN) :: zetc, rpgfc, zc, cerf
133 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
134 : REAL(KIND=dp), INTENT(IN) :: rab2
135 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
136 : REAL(KIND=dp), INTENT(IN) :: rac2, rbc2
137 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vabc
138 : REAL(KIND=dp), DIMENSION(:, :, :) :: verf, vnuc
139 : REAL(KIND=dp), DIMENSION(0:) :: f
140 : INTEGER, INTENT(IN), OPTIONAL :: maxder
141 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL :: vabc_plus, vnabc, pVp_sum, pVp
142 : LOGICAL, OPTIONAL :: dkh_erfc
143 :
144 : INTEGER, DIMENSION(2, 3), PARAMETER :: components = RESHAPE([1, 1, 2, 2, 3, 3], [2, 3])
145 :
146 : INTEGER :: axis, i, ipgf, j, jpgf, la_max, lb_max, &
147 : maxder_local, n, na, nap, nb, &
148 : ncol_pvp, nmax, nrow_pvp
149 4975432 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: map_a, map_b
150 : LOGICAL :: do_dkh
151 : REAL(KIND=dp) :: dab, dac, dbc, f0, f1, f4, ferf, fnuc, &
152 : rcp2, t, zetp, zetq, zetw
153 4975432 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: potential, work
154 4975432 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: derivatives
155 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp, rcp, rpw
156 :
157 : ! ---------------------------------------------------------------------------
158 :
159 4975432 : IF (PRESENT(maxder)) THEN
160 78910605 : verf(1:ncoset(la_max1), 1:ncoset(lb_max1), 1) = 0._dp
161 78910605 : vnuc(1:ncoset(la_max1), 1:ncoset(lb_max1), 1) = 0._dp
162 : END IF
163 :
164 : ! The legacy output uses the full Cartesian B range.
165 : MARK_USED(lb_min1)
166 :
167 4975432 : do_dkh = .FALSE.
168 :
169 4975432 : IF (PRESENT(pVp)) THEN
170 0 : do_dkh = .TRUE.
171 0 : pVp = 0.0_dp
172 : ! pVp matrices requires angular momentum quantum numbers l+1 and l-1
173 :
174 0 : la_max = la_max1 + 1
175 0 : IF (PRESENT(maxder)) THEN
176 0 : IF (maxder > 0) THEN
177 0 : la_max = la_max1
178 : END IF
179 : END IF
180 0 : lb_max = lb_max1 + 1
181 : ELSE
182 4975432 : do_dkh = .FALSE.
183 4975432 : la_max = la_max1
184 4975432 : lb_max = lb_max1
185 : END IF
186 :
187 4975432 : nmax = la_max + lb_max + 1
188 :
189 4975432 : maxder_local = 0
190 4975432 : IF (PRESENT(maxder)) maxder_local = maxder
191 :
192 4975432 : IF (do_dkh) THEN
193 : ! With maxder > 0, la_max1 includes the extra shells for vabc_plus.
194 : ! The pVp output ends at la_max1-maxder_local; its ladder needs one more shell.
195 0 : nrow_pvp = ncoset(la_max1 - maxder_local)
196 0 : ncol_pvp = ncoset(lb_max1)
197 0 : ALLOCATE (map_a(4, nrow_pvp, 3), map_b(4, ncol_pvp, 3))
198 : ALLOCATE (potential(ncoset(la_max1 - maxder_local + 1), ncoset(lb_max1 + 1)), &
199 0 : work(nrow_pvp, ncoset(lb_max1 + 1)), derivatives(nrow_pvp, ncol_pvp, 3))
200 0 : DO axis = 1, 3
201 0 : CALL cartesian_ladder_map(la_max1 - maxder_local, 0, axis, map_a(:, :, axis), source_lmin=0)
202 0 : CALL cartesian_ladder_map(lb_max1, 0, axis, map_b(:, :, axis), source_lmin=0)
203 : END DO
204 : END IF
205 :
206 : ! *** Calculate the distances of the centers a, b and c ***
207 :
208 4975432 : dab = SQRT(rab2)
209 4975432 : dac = SQRT(rac2)
210 4975432 : dbc = SQRT(rbc2)
211 :
212 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
213 :
214 4975432 : na = 0
215 4975432 : nap = 0
216 :
217 12505470 : DO ipgf = 1, npgfa
218 :
219 : ! *** Screening ***
220 :
221 7530038 : IF (do_dkh) THEN
222 0 : IF (dkh_erfc) THEN
223 0 : IF (rpgfa(ipgf) + rpgfc < dac) THEN
224 0 : na = na + ncoset(la_max1 - maxder_local) !JT
225 0 : nap = nap + ncoset(la_max1) !JT
226 0 : CYCLE
227 : END IF
228 : END IF
229 : ELSE
230 7530038 : IF (rpgfa(ipgf) + rpgfc < dac) THEN
231 1796992 : na = na + ncoset(la_max1 - maxder_local) !JT
232 1796992 : nap = nap + ncoset(la_max1) !JT
233 1796992 : CYCLE
234 : END IF
235 : END IF
236 :
237 5733046 : nb = 0
238 :
239 14661503 : DO jpgf = 1, npgfb
240 :
241 : ! *** Screening ***
242 8928457 : IF (do_dkh) THEN
243 0 : IF (dkh_erfc) THEN
244 0 : IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
245 : (rpgfa(ipgf) + rpgfb(jpgf) < dab)) THEN
246 0 : nb = nb + ncoset(lb_max1) !JT
247 0 : CYCLE
248 : END IF
249 : END IF
250 : ELSE
251 8928457 : IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
252 : (rpgfa(ipgf) + rpgfb(jpgf) < dab)) THEN
253 2557480 : nb = nb + ncoset(lb_max1) !JT
254 2557480 : CYCLE
255 : END IF
256 : END IF
257 : ! *** Calculate some prefactors ***
258 :
259 6370977 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
260 6370977 : zetq = 1.0_dp/zetc
261 6370977 : zetw = 1.0_dp/(zeta(ipgf) + zetb(jpgf) + zetc)
262 :
263 6370977 : f1 = zetb(jpgf)*zetp
264 6370977 : f4 = -zetc*zetw
265 :
266 6370977 : f0 = EXP(-zeta(ipgf)*f1*rab2)
267 :
268 6370977 : fnuc = 2.0_dp*pi*zetp*f0
269 6370977 : ferf = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq*f0
270 :
271 25483908 : rap(:) = f1*rab(:)
272 25483908 : rcp(:) = rap(:) - rac(:)
273 25483908 : rpw(:) = f4*rcp(:)
274 :
275 : ! *** Calculate the incomplete Gamma function values ***
276 :
277 6370977 : rcp2 = rcp(1)*rcp(1) + rcp(2)*rcp(2) + rcp(3)*rcp(3)
278 :
279 6370977 : t = rcp2/zetp
280 :
281 6370977 : CALL fgamma(nmax - 1, t, f)
282 :
283 : ! *** Calculate the basic nuclear attraction integrals [s|A(0)|s]{n} ***
284 :
285 24448707 : DO n = 1, nmax
286 24448707 : vnuc(1, 1, n) = fnuc*f(n - 1)
287 : END DO
288 :
289 : ! *** Calculate the incomplete Gamma function values ***
290 :
291 6370977 : t = -f4*rcp2/zetp
292 :
293 6370977 : CALL fgamma(nmax - 1, t, f)
294 :
295 : ! *** Calculate the basic three-center Coulomb integrals [ss||s]{n} ***
296 :
297 24448707 : DO n = 1, nmax
298 24448707 : verf(1, 1, n) = ferf*f(n - 1)
299 : END DO
300 :
301 25483908 : rbp = rap - rab
302 25483908 : CALL os_rr_radial(rap, la_max, rbp, lb_max, -rcp, zetp, 1.0_dp, nmax - 1, vnuc)
303 6370977 : CALL os_rr_radial(rap, la_max, rbp, lb_max, rpw, zetp, -f4, nmax - 1, verf)
304 :
305 : ! *** Add the contribution of the current pair ***
306 : ! *** of primitive Gaussian-type functions ***
307 :
308 : !JT
309 6370977 : IF (do_dkh) THEN
310 0 : DO j = 1, ncoset(lb_max1)
311 0 : DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1 - maxder_local)
312 0 : vnabc(na + i, nb + j) = vnabc(na + i, nb + j) + cerf*verf(i, j, 1)
313 : END DO
314 : END DO
315 0 : DO j = 1, ncoset(lb_max1)
316 0 : DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1 - maxder_local)
317 0 : vabc(na + i, nb + j) = vabc(na + i, nb + j) - zc*vnuc(i, j, 1)
318 : END DO
319 : END DO
320 : ELSE
321 27743689 : DO j = 1, ncoset(lb_max1) !JT
322 83551426 : DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1 - maxder_local) !JT
323 : vabc(na + i, nb + j) = vabc(na + i, nb + j) - &
324 77180449 : zc*vnuc(i, j, 1) + cerf*verf(i, j, 1)
325 : END DO
326 : END DO
327 : END IF
328 : !JT
329 :
330 6370977 : IF (PRESENT(maxder)) THEN
331 12261433 : DO j = 1, ncoset(lb_max1) !JT
332 94776709 : DO i = 1, ncoset(la_max1) !JT
333 : vabc_plus(nap + i, nb + j) = vabc_plus(nap + i, nb + j) - &
334 91996685 : zc*vnuc(i, j, 1) + cerf*verf(i, j, 1)
335 : END DO
336 : END DO
337 : END IF
338 : !JTs
339 6370977 : IF (do_dkh) THEN
340 : ! pVp = sum_i D_Ai D_Bi V; both electronic derivatives change sign.
341 0 : potential(:, :) = -zc*vnuc(1:SIZE(potential, 1), 1:SIZE(potential, 2), 1)
342 0 : IF (dkh_erfc) THEN
343 0 : potential(:, :) = potential + cerf*verf(1:SIZE(potential, 1), 1:SIZE(potential, 2), 1)
344 : END IF
345 : CALL apply_cartesian_operator_pair(potential, derivatives, map_a, map_b, components, &
346 0 : [.TRUE., .TRUE.], [zeta(ipgf), zetb(jpgf)], [.FALSE., .TRUE.], work)
347 0 : pVp(1:nrow_pvp, 1:ncol_pvp) = SUM(derivatives, DIM=3)
348 :
349 0 : DO j = 1, ncoset(lb_max1)
350 0 : DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1 - maxder_local)
351 0 : pVp_sum(na + i, nb + j) = pVp_sum(na + i, nb + j) + pVp(i, j)
352 : END DO
353 : END DO
354 : END IF
355 : !JTe
356 12104023 : nb = nb + ncoset(lb_max1) !JT
357 :
358 : END DO
359 :
360 5733046 : na = na + ncoset(la_max1 - maxder_local) !JT
361 10708478 : nap = nap + ncoset(la_max1) !JT
362 :
363 : END DO
364 :
365 4975432 : END SUBROUTINE verfc
366 :
367 : END MODULE ai_verfc
|
