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 Local UPF matrix elements from the original radial potential.
10 : !> The Gaussian core electrostatics supplies -Z erf(sqrt(alpha)*r)/r.
11 : !> Its complement is integrated with analytic angles and radial quadrature.
12 : !> Gaussian products, center derivatives and RI integrals share this operator.
13 : ! **************************************************************************************************
14 : MODULE upf_local_integrals
15 : USE kinds, ONLY: dp
16 : USE mathconstants, ONLY: fourpi,&
17 : pi
18 : USE mathlib, ONLY: binomial
19 : USE orbital_pointers, ONLY: coset,&
20 : indco,&
21 : ncoset
22 : USE upf_projector_integrals, ONLY: upf_projector_init,&
23 : upf_projector_overlap,&
24 : upf_projector_type
25 : USE upf_types, ONLY: atom_upfpot_type
26 : #include "../base/base_uses.f90"
27 :
28 : IMPLICIT NONE
29 : PRIVATE
30 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'upf_local_integrals'
31 :
32 : TYPE upf_local_type
33 : TYPE(upf_projector_type) :: radial
34 : REAL(KIND=dp) :: eps = 0.0_dp
35 : END TYPE upf_local_type
36 :
37 : PUBLIC :: upf_local_type, upf_local_init, upf_local_matrix, upf_local_ri
38 :
39 : CONTAINS
40 :
41 : ! **************************************************************************************************
42 : !> \brief Prepare the short-range local potential outside parallel regions.
43 : !> \param pot Original UPF data, with PP_LOCAL already converted to Hartree.
44 : !> \param alpha Gaussian core-charge exponent, not its screening radius.
45 : !> \param degree Largest Gaussian product degree, including derivatives or moments.
46 : !> \param local Immutable integration data. Reinitialization retains its threshold.
47 : !> \param eps Optional threshold for estimating the potential's neighbor-list radius.
48 : !> This is NOT an integration tolerance: EPS_PPL is normally 1e-2.
49 : !> All original radial samples, including the tail, enter the quadrature.
50 : ! **************************************************************************************************
51 56 : SUBROUTINE upf_local_init(pot, alpha, degree, local, eps)
52 : TYPE(atom_upfpot_type), INTENT(IN) :: pot
53 : REAL(KIND=dp), INTENT(IN) :: alpha
54 : INTEGER, INTENT(IN) :: degree
55 : TYPE(upf_local_type), INTENT(INOUT) :: local
56 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps
57 :
58 : INTEGER :: i, last, n
59 56 : REAL(KIND=dp), ALLOCATABLE :: short_range(:)
60 5880 : TYPE(atom_upfpot_type) :: scalar
61 :
62 56 : CPASSERT(alpha > 0.0_dp)
63 56 : CPASSERT(SIZE(pot%r) == SIZE(pot%vlocal))
64 56 : CPASSERT(SIZE(pot%r) >= 3)
65 56 : IF (PRESENT(eps)) local%eps = eps
66 56 : CPASSERT(local%eps >= 0.0_dp)
67 168 : ALLOCATE (short_range(SIZE(pot%vlocal)))
68 68562 : short_range(:) = pot%vlocal
69 68562 : DO i = 1, SIZE(pot%r)
70 68562 : IF (pot%r(i) > 0.0_dp) THEN
71 68486 : short_range(i) = short_range(i) + pot%zion*ERF(SQRT(alpha)*pot%r(i))/pot%r(i)
72 : ELSE
73 20 : short_range(i) = short_range(i) + 2.0_dp*pot%zion*SQRT(alpha/pi)
74 : END IF
75 : END DO
76 56 : last = 0
77 28366 : DO i = SIZE(short_range), 1, -1
78 28366 : IF (ABS(short_range(i)) > local%eps) THEN
79 : last = i
80 : EXIT
81 : END IF
82 : END DO
83 56 : n = SIZE(short_range)
84 56 : scalar%number_of_proj = 1
85 224 : ALLOCATE (scalar%r(n), scalar%rab(n), scalar%lbeta(1))
86 68562 : scalar%r(:) = pot%r(:n)
87 68562 : scalar%rab(:) = pot%rab(:n)
88 112 : scalar%lbeta = 0
89 112 : ALLOCATE (scalar%beta(n, 1), scalar%dion(1, 1))
90 : ! Y00 = 1/sqrt(4*pi); absorb it to integrate a scalar radial potential.
91 68562 : scalar%beta(:, 1) = SQRT(fourpi)*short_range(:n)
92 168 : scalar%dion = 0.0_dp
93 56 : CALL upf_projector_init(scalar, degree, local%radial)
94 : ! Screening estimates the support in a pair-list bound. It must never
95 : ! truncate the function actually integrated within a retained pair.
96 56 : local%radial%radius = 0.0_dp
97 56 : IF (last > 0) local%radial%radius = pot%r(last)
98 56 : END SUBROUTINE upf_local_init
99 :
100 : ! **************************************************************************************************
101 : !> \brief Evaluate primitive three-center local integrals by the Gaussian product theorem.
102 : !> \param local Prepared radial potential.
103 : !> \param la Largest Cartesian degree at A; all degrees 0..la are returned.
104 : !> \param alpha Primitive exponents at A.
105 : !> \param lb Largest Cartesian degree at B; all degrees 0..lb are returned.
106 : !> \param beta Primitive exponents at B.
107 : !> \param rab B minus A.
108 : !> \param rac Potential center C minus A.
109 : !> \param v Matrix elements in primitive-major Cartesian order.
110 : !> \param da Optional derivatives with respect to A, at fixed B and C.
111 : !> \param db Optional derivatives with respect to B, at fixed A and C.
112 : !> The derivative with respect to C is -(da+db).
113 : ! **************************************************************************************************
114 10306 : SUBROUTINE upf_local_matrix(local, la, alpha, lb, beta, rab, rac, v, da, db)
115 : TYPE(upf_local_type), INTENT(IN) :: local
116 : INTEGER, INTENT(IN) :: la
117 : REAL(KIND=dp), INTENT(IN) :: alpha(:)
118 : INTEGER, INTENT(IN) :: lb
119 : REAL(KIND=dp), INTENT(IN) :: beta(:), rab(3), rac(3)
120 : REAL(KIND=dp), INTENT(OUT) :: v(:, :)
121 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: da(:, :, :), db(:, :, :)
122 :
123 : INTEGER :: a, aa, axis, b, bb, da_max, db_max, &
124 : degree, i, ia, ib, ix, iy, iz, j, k, &
125 : pa(3), pb(3), shifted(3)
126 : REAL(KIND=dp) :: center(3), p, pref, shift_a(3), &
127 : shift_b(3), value
128 : REAL(KIND=dp), ALLOCATABLE :: moments(:, :, :), &
129 10306 : polynomial(:, :, :, :), work(:, :)
130 :
131 10306 : CPASSERT(la >= 0 .AND. lb >= 0)
132 54376 : CPASSERT(ALL(alpha > 0.0_dp) .AND. ALL(beta > 0.0_dp))
133 10306 : CPASSERT(SIZE(v, 1) == SIZE(alpha)*ncoset(la))
134 10306 : CPASSERT(SIZE(v, 2) == SIZE(beta)*ncoset(lb))
135 10306 : degree = la + lb
136 10306 : IF (PRESENT(da) .OR. PRESENT(db)) degree = degree + 1
137 10306 : da_max = la
138 10306 : db_max = lb
139 10306 : IF (PRESENT(da)) THEN
140 1544 : CPASSERT(ALL(SHAPE(da) == [SIZE(v, 1), SIZE(v, 2), 3]))
141 1698524 : da = 0.0_dp
142 386 : da_max = la + 1
143 : END IF
144 10306 : IF (PRESENT(db)) THEN
145 1544 : CPASSERT(ALL(SHAPE(db) == [SIZE(v, 1), SIZE(v, 2), 3]))
146 1698524 : db = 0.0_dp
147 386 : db_max = lb + 1
148 : END IF
149 82448 : ALLOCATE (moments(ncoset(degree), 1, 1), polynomial(0:da_max, 0:db_max, 0:degree, 3))
150 41224 : ALLOCATE (work(ncoset(da_max), ncoset(db_max)))
151 1318201 : v = 0.0_dp
152 27337 : DO ib = 1, SIZE(beta)
153 60701 : DO ia = 1, SIZE(alpha)
154 33364 : p = alpha(ia) + beta(ib)
155 133456 : shift_a = beta(ib)*rab/p
156 133456 : shift_b = -alpha(ia)*rab/p
157 133456 : center = rac - shift_a
158 133456 : pref = EXP(-alpha(ia)*beta(ib)/p*SUM(rab**2))
159 33364 : IF (pref == 0.0_dp) CYCLE
160 66600 : CALL upf_projector_overlap(local%radial, degree, [p], center, 0, moments)
161 33300 : polynomial = 0.0_dp
162 133200 : DO axis = 1, 3
163 354732 : DO a = 0, da_max
164 851283 : DO b = 0, MIN(db_max, degree - a)
165 1816509 : DO i = 0, a
166 3712851 : DO j = 0, b
167 : polynomial(a, b, i + j, axis) = polynomial(a, b, i + j, axis) &
168 3183000 : + binomial(a, i)*binomial(b, j)*shift_a(axis)**(a - i)*shift_b(axis)**(b - j)
169 : END DO
170 : END DO
171 : END DO
172 : END DO
173 : END DO
174 235558 : DO b = 1, ncoset(db_max)
175 809032 : pb = indco(:, b)
176 3059791 : DO a = 1, ncoset(da_max)
177 11296932 : pa = indco(:, a)
178 11296932 : IF (SUM(pa + pb) > degree) CYCLE
179 2294865 : value = 0.0_dp
180 7689722 : DO iz = 0, pa(3) + pb(3)
181 19878329 : DO iy = 0, pa(2) + pb(2)
182 44190549 : DO ix = 0, pa(1) + pb(1)
183 26607085 : k = coset(ix, iy, iz)
184 : value = value + polynomial(pa(1), pb(1), ix, 1)*polynomial(pa(2), pb(2), iy, 2) &
185 38795692 : *polynomial(pa(3), pb(3), iz, 3)*moments(k, 1, 1)
186 : END DO
187 : END DO
188 : END DO
189 3026491 : work(a, b) = pref*value
190 : END DO
191 : END DO
192 198361 : DO b = 1, ncoset(lb)
193 592120 : pb = indco(:, b)
194 148030 : bb = (ib - 1)*ncoset(lb) + b
195 1426739 : DO a = 1, ncoset(la)
196 4981380 : pa = indco(:, a)
197 1245345 : aa = (ia - 1)*ncoset(la) + a
198 1245345 : v(aa, bb) = work(a, b)
199 5129410 : DO axis = 1, 3
200 3736035 : IF (PRESENT(da)) THEN
201 1656960 : shifted = pa
202 1656960 : shifted(axis) = shifted(axis) + 1
203 1656960 : da(aa, bb, axis) = 2*alpha(ia)*work(coset(shifted(1), shifted(2), shifted(3)), b)
204 1656960 : IF (pa(axis) > 0) THEN
205 730680 : shifted = pa
206 730680 : shifted(axis) = shifted(axis) - 1
207 : da(aa, bb, axis) = da(aa, bb, axis) &
208 730680 : - pa(axis)*work(coset(shifted(1), shifted(2), shifted(3)), b)
209 : END IF
210 : END IF
211 4981380 : IF (PRESENT(db)) THEN
212 1656960 : shifted = pb
213 1656960 : shifted(axis) = shifted(axis) + 1
214 1656960 : db(aa, bb, axis) = 2*beta(ib)*work(a, coset(shifted(1), shifted(2), shifted(3)))
215 1656960 : IF (pb(axis) > 0) THEN
216 692880 : shifted = pb
217 692880 : shifted(axis) = shifted(axis) - 1
218 : db(aa, bb, axis) = db(aa, bb, axis) &
219 692880 : - pb(axis)*work(a, coset(shifted(1), shifted(2), shifted(3)))
220 : END IF
221 : END IF
222 : END DO
223 : END DO
224 : END DO
225 : END DO
226 : END DO
227 10306 : END SUBROUTINE upf_local_matrix
228 :
229 : ! **************************************************************************************************
230 : !> \brief Local potential against a single Cartesian Gaussian for RI.
231 : !> \param local Prepared radial potential.
232 : !> \param la Largest Cartesian degree.
233 : !> \param zeta Primitive exponents.
234 : !> \param rac Potential center minus Gaussian center.
235 : !> \param v Primitive-major Cartesian integrals.
236 : !> \param dv Optional Gaussian-center derivatives, in x,y,z order.
237 : ! **************************************************************************************************
238 64 : SUBROUTINE upf_local_ri(local, la, zeta, rac, v, dv)
239 : TYPE(upf_local_type), INTENT(IN) :: local
240 : INTEGER, INTENT(IN) :: la
241 : REAL(KIND=dp), INTENT(IN) :: zeta(:), rac(3)
242 : REAL(KIND=dp), INTENT(OUT) :: v(:)
243 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: dv(:, :)
244 :
245 : INTEGER :: n, nder
246 : REAL(KIND=dp), ALLOCATABLE :: integrals(:, :, :)
247 :
248 64 : nder = 0
249 64 : IF (PRESENT(dv)) nder = 1
250 64 : n = SIZE(zeta)*ncoset(la)
251 64 : CPASSERT(SIZE(v) == n)
252 256 : ALLOCATE (integrals(n, 1, ncoset(nder)))
253 64 : CALL upf_projector_overlap(local%radial, la, zeta, rac, nder, integrals)
254 4544 : v = integrals(:, 1, 1)
255 64 : IF (PRESENT(dv)) THEN
256 96 : CPASSERT(ALL(SHAPE(dv) == [n, 3]))
257 6848 : dv = integrals(:, 1, 2:4)
258 : END IF
259 64 : END SUBROUTINE upf_local_ri
260 0 : END MODULE upf_local_integrals
|
