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 Half-numerical overlaps with the original scalar UPF radial projectors.
10 : !> Cartesian polynomials and angular products are integrated analytically.
11 : !> Only a radial quadrature remains; no fit or rotation of projector identities
12 : !> is performed. The UPF reader preserves the tabulated beta amplitude and
13 : !> includes the Ry-to-Hartree factor in the coupling matrix D.
14 : ! **************************************************************************************************
15 : MODULE upf_projector_integrals
16 : USE ISO_C_BINDING, ONLY: C_DOUBLE,&
17 : C_INT
18 : USE kinds, ONLY: dp
19 : USE libgrpp, ONLY: libgrpp_init
20 : USE mathconstants, ONLY: dfac,&
21 : fourpi,&
22 : maxfac
23 : USE mathlib, ONLY: binomial
24 : USE orbital_pointers, ONLY: indco,&
25 : init_orbital_pointers,&
26 : nco,&
27 : ncoset
28 : USE orbital_transformation_matrices, ONLY: init_spherical_harmonics,&
29 : orbtramat
30 : USE upf_types, ONLY: atom_upfpot_type
31 : #include "../base/base_uses.f90"
32 :
33 : IMPLICIT NONE
34 : PRIVATE
35 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'upf_projector_integrals'
36 :
37 : TYPE upf_projector_type
38 : INTEGER :: degree = -1, lmax = -1, harmonic_max = -1, nproj = 0
39 : REAL(KIND=dp) :: radius = 0.0_dp
40 : INTEGER, ALLOCATABLE :: radial_index(:), angular_index(:), lbeta(:)
41 : REAL(KIND=dp), ALLOCATABLE :: r(:), weighted_beta(:, :), d(:, :), gaunt(:, :, :)
42 : END TYPE upf_projector_type
43 :
44 : PUBLIC :: upf_projector_type, upf_projector_init, upf_projector_overlap
45 :
46 : INTERFACE
47 : FUNCTION scaled_bessel(n, x) RESULT(value) BIND(C, name="libgrpp_modified_bessel_scaled")
48 : IMPORT :: C_DOUBLE, C_INT
49 : INTEGER(C_INT), VALUE :: n
50 : REAL(C_DOUBLE), VALUE :: x
51 : REAL(C_DOUBLE) :: value
52 : END FUNCTION scaled_bessel
53 : END INTERFACE
54 :
55 : CONTAINS
56 :
57 : ! **************************************************************************************************
58 : !> \brief Prepare immutable radial and angular data outside parallel regions.
59 : !> \param pot Original UPF data, including beta(r), not r*beta(r).
60 : !> \param degree Largest Cartesian degree needed (AO degree + derivatives/moments).
61 : !> \param projectors Prepared data; reinitializing also releases the previous allocation.
62 : ! **************************************************************************************************
63 92 : SUBROUTINE upf_projector_init(pot, degree, projectors)
64 : TYPE(atom_upfpot_type), INTENT(IN) :: pot
65 : INTEGER, INTENT(IN) :: degree
66 : TYPE(upf_projector_type), INTENT(OUT) :: projectors
67 :
68 : INTEGER :: a, b, i, j, k, l, lambda, m, n, &
69 : nsimpson, p, powers(3), q, s
70 92 : REAL(KIND=dp), ALLOCATABLE :: moments(:, :), weights(:)
71 :
72 92 : CPASSERT(degree >= 0)
73 92 : IF (pot%has_so) THEN
74 0 : CPABORT("Spin-orbit UPF projectors are not supported")
75 : END IF
76 92 : CPASSERT(pot%number_of_proj > 0)
77 298 : CPASSERT(ALL(pot%lbeta >= 0))
78 92 : projectors%degree = degree
79 298 : projectors%lmax = MAXVAL(pot%lbeta)
80 92 : projectors%harmonic_max = degree + projectors%lmax
81 : ! Both the factorial and bundled Bessel tables cover this range.
82 92 : CPASSERT(projectors%harmonic_max <= MIN(20, maxfac/2))
83 92 : CALL init_orbital_pointers(projectors%harmonic_max)
84 92 : CALL init_spherical_harmonics(projectors%harmonic_max, -1)
85 92 : CALL libgrpp_init()
86 :
87 92 : n = SIZE(pot%r)
88 92 : CPASSERT(n >= 3)
89 92 : CPASSERT(SIZE(pot%rab) == n)
90 92 : CPASSERT(SIZE(pot%beta, 1) == n)
91 92 : CPASSERT(SIZE(pot%beta, 2) == pot%number_of_proj)
92 116290 : CPASSERT(ALL(pot%r(2:) > pot%r(:n - 1)))
93 116382 : CPASSERT(ALL(pot%rab > 0.0_dp))
94 460 : ALLOCATE (projectors%r(n), projectors%lbeta(pot%number_of_proj))
95 116382 : projectors%r(:) = pot%r
96 298 : projectors%lbeta(:) = pot%lbeta
97 184 : ALLOCATE (weights(n))
98 92 : weights = 0.0_dp
99 : ! Composite Simpson in the mesh index; PP_RAB is dr/d(index).
100 : ! An even mesh has one final trapezoidal interval.
101 92 : nsimpson = n - 1 + MODULO(n, 2)
102 58266 : weights(1:nsimpson:2) = 2.0_dp/3.0_dp
103 58174 : weights(2:nsimpson - 1:2) = 4.0_dp/3.0_dp
104 92 : weights(1) = 1.0_dp/3.0_dp
105 92 : weights(nsimpson) = 1.0_dp/3.0_dp
106 160 : IF (nsimpson < n) weights(n - 1:n) = weights(n - 1:n) + 0.5_dp
107 116382 : weights(:) = weights*pot%rab*pot%r**2
108 368 : ALLOCATE (projectors%weighted_beta(n, pot%number_of_proj))
109 298 : DO i = 1, pot%number_of_proj
110 275152 : projectors%weighted_beta(:, i) = weights*pot%beta(:, i)
111 97426 : DO j = n, 1, -1
112 97334 : IF (pot%beta(j, i) /= 0.0_dp) THEN
113 206 : projectors%radius = MAX(projectors%radius, pot%r(j))
114 206 : EXIT
115 : END IF
116 : END DO
117 : END DO
118 :
119 298 : projectors%nproj = SUM(2*pot%lbeta + 1)
120 368 : ALLOCATE (projectors%radial_index(projectors%nproj), projectors%angular_index(projectors%nproj))
121 92 : p = 0
122 298 : DO i = 1, pot%number_of_proj
123 206 : l = pot%lbeta(i)
124 804 : DO m = 1, 2*l + 1
125 506 : p = p + 1
126 506 : projectors%radial_index(p) = i
127 712 : projectors%angular_index(p) = l*l + m
128 : END DO
129 : END DO
130 368 : ALLOCATE (projectors%d(projectors%nproj, projectors%nproj))
131 8664 : projectors%d = 0.0_dp
132 598 : DO p = 1, projectors%nproj
133 8664 : DO q = 1, projectors%nproj
134 8066 : i = projectors%radial_index(p)
135 8066 : j = projectors%radial_index(q)
136 8572 : IF (projectors%angular_index(p) == projectors%angular_index(q)) THEN
137 866 : projectors%d(p, q) = pot%dion(i, j)
138 7200 : ELSE IF (pot%lbeta(i) /= pot%lbeta(j)) THEN
139 4812 : CPASSERT(ABS(pot%dion(i, j)) < 1.e-12_dp)
140 : END IF
141 : END DO
142 : END DO
143 :
144 : ! G(s,lm,lambda*nu) = integral Omega**s Y_lm(Omega) Y_lambda,nu(Omega).
145 0 : ALLOCATE (projectors%gaunt(ncoset(degree), (projectors%lmax + 1)**2, &
146 460 : (projectors%harmonic_max + 1)**2))
147 1199920 : projectors%gaunt = 0.0_dp
148 4432 : DO s = 1, ncoset(degree)
149 17360 : k = SUM(indco(:, s))
150 10656 : DO l = 0, projectors%lmax
151 10364 : IF (.NOT. ANY(pot%lbeta == l)) CYCLE
152 52266 : DO lambda = MAX(0, l - k), l + k
153 41702 : IF (MODULO(k + l + lambda, 2) /= 0) CYCLE
154 89928 : ALLOCATE (moments(nco(l), nco(lambda)))
155 388830 : DO b = 1, nco(lambda)
156 1384704 : DO a = 1, nco(l)
157 3983496 : powers = indco(:, s) + indco(:, ncoset(l - 1) + a) + indco(:, ncoset(lambda - 1) + b)
158 995874 : moments(a, b) = 0.0_dp
159 2378980 : IF (ALL(MODULO(powers, 2) == 0)) THEN
160 3081600 : moments(a, b) = fourpi*PRODUCT(dfac(powers - 1))/dfac(SUM(powers) + 1)
161 : END IF
162 : END DO
163 : END DO
164 : projectors%gaunt(s, l*l + 1:(l + 1)**2, lambda*lambda + 1:(lambda + 1)**2) = &
165 68671416 : MATMUL(MATMUL(orbtramat(l)%slm, moments), TRANSPOSE(orbtramat(lambda)%slm))
166 28706 : DEALLOCATE (moments)
167 : END DO
168 : END DO
169 : END DO
170 92 : END SUBROUTINE upf_projector_init
171 :
172 : ! **************************************************************************************************
173 : !> \brief Evaluate <Cartesian Gaussian|beta_i Y_lm> and AO-center derivatives.
174 : !> \param projectors Immutable prepared data.
175 : !> \param la Largest AO Cartesian degree (all degrees 0..la are returned).
176 : !> \param zeta Primitive Gaussian exponents.
177 : !> \param rac Projector center minus AO center, as in ai_overlap.
178 : !> \param nder Maximum derivative degree, or moment degree in moment mode.
179 : !> \param integrals Primitive-major AO order, original radial/m order, coset derivative order.
180 : !> \param moment_origin Optional AO center minus moment reference; requests moments
181 : !> instead of derivatives, in the same coset ordering including the overlap.
182 : ! **************************************************************************************************
183 35844 : SUBROUTINE upf_projector_overlap(projectors, la, zeta, rac, nder, integrals, moment_origin)
184 : TYPE(upf_projector_type), INTENT(IN) :: projectors
185 : INTEGER, INTENT(IN) :: la
186 : REAL(KIND=dp), INTENT(IN) :: zeta(:), rac(3)
187 : INTEGER, INTENT(IN) :: nder
188 : REAL(KIND=dp), INTENT(OUT) :: integrals(:, :, :)
189 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: moment_origin(3)
190 :
191 : INTEGER :: a(3), axis, b(3), degree, deriv(3), i, &
192 : id, ip, j, k, l, lambda, nb, nc, nh, &
193 : nr, p, power, s, shifted(3), step, t
194 : REAL(KIND=dp) :: alpha, direction(3), distance, factor
195 35844 : REAL(KIND=dp), ALLOCATABLE :: angular(:, :, :), cartesian(:, :), centered(:, :), &
196 35844 : coeff(:, :), envelope(:), harmonics(:), monomials(:), next_coeff(:), radial(:, :, :), &
197 35844 : radial_power(:), radial_work(:), shift(:, :)
198 :
199 35844 : degree = la + nder
200 35844 : CPASSERT(la >= 0 .AND. nder >= 0)
201 35844 : CPASSERT(degree <= projectors%degree)
202 75296 : CPASSERT(ALL(zeta > 0.0_dp))
203 35844 : CPASSERT(SIZE(integrals, 1) == SIZE(zeta)*ncoset(la))
204 35844 : CPASSERT(SIZE(integrals, 2) == projectors%nproj)
205 35844 : CPASSERT(SIZE(integrals, 3) == ncoset(nder))
206 35844 : nc = ncoset(degree)
207 35844 : nh = degree + projectors%lmax
208 35844 : nr = SIZE(projectors%r)
209 35844 : nb = SIZE(projectors%weighted_beta, 2)
210 179220 : ALLOCATE (harmonics((nh + 1)**2), monomials(ncoset(nh)))
211 143376 : distance = NORM2(rac)
212 35844 : direction = [0.0_dp, 0.0_dp, 1.0_dp]
213 141771 : IF (distance > 0.0_dp) direction = -rac/distance
214 637915 : DO i = 1, ncoset(nh)
215 2444128 : monomials(i) = PRODUCT(direction**indco(:, i))
216 : END DO
217 144958 : DO l = 0, nh
218 : harmonics(l*l + 1:(l + 1)**2) = &
219 9860271 : MATMUL(orbtramat(l)%slm, monomials(ncoset(l - 1) + 1:ncoset(l)))
220 : END DO
221 179220 : ALLOCATE (angular(nc, (projectors%lmax + 1)**2, 0:nh))
222 144958 : DO lambda = 0, nh
223 267128 : DO i = 1, SIZE(angular, 2)
224 : angular(:, i, lambda) = MATMUL( &
225 : projectors%gaunt(:nc, i, lambda*lambda + 1:(lambda + 1)**2), &
226 28281950 : harmonics(lambda*lambda + 1:(lambda + 1)**2))
227 : END DO
228 : END DO
229 143376 : ALLOCATE (shift(nc, nc))
230 35844 : shift = 0.0_dp
231 627595 : DO t = 1, nc
232 2367004 : a = indco(:, t)
233 28475528 : DO s = 1, nc
234 111391732 : b = indco(:, s)
235 58972504 : IF (ANY(b > a)) CYCLE
236 : factor = 1.0_dp
237 16384396 : DO axis = 1, 3
238 16384396 : factor = factor*binomial(a(axis), b(axis))*rac(axis)**(a(axis) - b(axis))
239 : END DO
240 28439684 : shift(t, s) = factor
241 : END DO
242 : END DO
243 322596 : ALLOCATE (radial(0:degree, nb, 0:nh), radial_work(nr), radial_power(nr), envelope(nr))
244 215064 : ALLOCATE (centered(nc, projectors%nproj), cartesian(nc, projectors%nproj))
245 215064 : ALLOCATE (coeff(0:degree, 3), next_coeff(0:degree))
246 2528635 : integrals = 0.0_dp
247 75296 : DO ip = 1, SIZE(zeta)
248 39452 : alpha = zeta(ip)
249 41943772 : envelope(:) = EXP(-alpha*(projectors%r - distance)**2)
250 158056 : DO lambda = 0, nh
251 132472083 : DO i = 1, nr
252 : radial_work(i) = envelope(i)*scaled_bessel(INT(lambda, C_INT), &
253 132472083 : REAL(2.0_dp*alpha*projectors%r(i)*distance, C_DOUBLE))
254 : END DO
255 132472083 : radial_power(:) = radial_work
256 623106 : DO k = 0, degree
257 1047310 : radial(k, :, lambda) = MATMUL(radial_power, projectors%weighted_beta)
258 553915001 : radial_power(:) = radial_power*projectors%r
259 : END DO
260 : END DO
261 39452 : centered = 0.0_dp
262 94036 : DO p = 1, projectors%nproj
263 54584 : i = projectors%radial_index(p)
264 54584 : j = projectors%angular_index(p)
265 1203875 : DO s = 1, nc
266 4439356 : k = SUM(indco(:, s))
267 8681970 : centered(s, p) = fourpi*DOT_PRODUCT(angular(s, j, :), radial(k, i, :))
268 : END DO
269 : END DO
270 20585712 : cartesian(:, :) = MATMUL(shift, centered)
271 118054 : DO id = 1, ncoset(nder)
272 171032 : deriv = indco(:, id)
273 781265 : DO t = 1, ncoset(la)
274 2796220 : a = indco(:, t)
275 699055 : coeff = 0.0_dp
276 2796220 : DO axis = 1, 3
277 2097165 : coeff(a(axis), axis) = 1.0_dp
278 2796220 : IF (PRESENT(moment_origin)) THEN
279 1008000 : coeff(:, axis) = 0.0_dp
280 315000 : DO power = 0, deriv(axis)
281 : coeff(a(axis) + power, axis) = binomial(deriv(axis), power)* &
282 315000 : moment_origin(axis)**(deriv(axis) - power)
283 : END DO
284 : ELSE
285 2047725 : DO step = 1, deriv(axis)
286 76560 : next_coeff = 0.0_dp
287 238410 : DO power = 0, a(axis) + step - 1
288 161850 : next_coeff(power + 1) = next_coeff(power + 1) + 2.0_dp*alpha*coeff(power, axis)
289 238410 : IF (power > 0) next_coeff(power - 1) = next_coeff(power - 1) - power*coeff(power, axis)
290 : END DO
291 2554605 : coeff(:, axis) = next_coeff
292 : END DO
293 : END IF
294 : END DO
295 35950242 : DO s = 1, nc
296 140833716 : shifted = indco(:, s)
297 35208429 : factor = coeff(shifted(1), 1)*coeff(shifted(2), 2)*coeff(shifted(3), 3)
298 35208429 : IF (factor == 0.0_dp) CYCLE
299 : integrals((ip - 1)*ncoset(la) + t, :, id) = integrals((ip - 1)*ncoset(la) + t, :, id) + &
300 6197499 : factor*cartesian(s, :)
301 : END DO
302 : END DO
303 : END DO
304 : END DO
305 35844 : END SUBROUTINE upf_projector_overlap
306 :
307 22482 : END MODULE upf_projector_integrals
|
