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 Coordinate and basis-center derivative ladders for Cartesian Gaussian integrals.
10 : !> Derivatives act on a basis function with the integral operator held fixed.
11 : !> Electronic-coordinate derivatives have the opposite sign.
12 : ! **************************************************************************************************
13 : MODULE ai_cartesian_operators
14 : USE kinds, ONLY: dp
15 : USE mathlib, ONLY: binomial_gen
16 : USE orbital_pointers, ONLY: coset, &
17 : indco, &
18 : ncoset
19 : #include "../base/base_uses.f90"
20 :
21 : IMPLICIT NONE
22 : PRIVATE
23 :
24 : PUBLIC :: cartesian_ladder_map, cartesian_ladder, apply_cartesian_operator_pair, apply_cartesian_operator, apply_cartesian_1d, &
25 : cartesian_polynomial_1d
26 :
27 : INTERFACE cartesian_ladder
28 : MODULE PROCEDURE cartesian_ladder_real, cartesian_ladder_complex
29 : END INTERFACE cartesian_ladder
30 :
31 : INTERFACE apply_cartesian_operator_pair
32 : MODULE PROCEDURE apply_cartesian_operator_pair_real, apply_cartesian_operator_pair_complex
33 : END INTERFACE apply_cartesian_operator_pair
34 :
35 : INTERFACE apply_cartesian_1d
36 : MODULE PROCEDURE apply_cartesian_1d_real, apply_cartesian_1d_complex
37 : END INTERFACE apply_cartesian_1d
38 :
39 : CONTAINS
40 :
41 : ! **************************************************************************************************
42 : !> \brief Construct base, raised, lowered and degree indices for a Cartesian ladder.
43 : !> \param lmax highest output shell or one-dimensional degree
44 : !> \param lmin lowest output shell (zero for one-dimensional tables)
45 : !> \param axis Cartesian axis, or zero for a one-dimensional table
46 : !> \param map rows contain base, raised, lowered and degree indices
47 : !> \param source_lmin optional lowest input shell; defaults to lmin for Cartesian blocks
48 : ! **************************************************************************************************
49 152590623 : SUBROUTINE cartesian_ladder_map(lmax, lmin, axis, map, source_lmin)
50 : INTEGER, INTENT(IN) :: lmax, lmin, axis
51 : INTEGER, DIMENSION(:, :), INTENT(OUT) :: map
52 : INTEGER, INTENT(IN), OPTIONAL :: source_lmin
53 :
54 : INTEGER :: i, ico, lowered(3), offset, power(3), &
55 : raised(3), source_offset
56 :
57 152590623 : IF (axis == 0) THEN
58 122582140 : DO i = 0, lmax
59 427842752 : map(:, i + 1) = [i + 1, i + 2, MAX(1, i), i]
60 : END DO
61 : ELSE
62 106323636 : offset = ncoset(lmin - 1)
63 106323636 : source_offset = offset
64 106323636 : IF (PRESENT(source_lmin)) source_offset = ncoset(source_lmin - 1)
65 567330225 : DO ico = offset + 1, ncoset(lmax)
66 1844026356 : power = indco(:, ico)
67 461006589 : raised = power
68 461006589 : lowered = power
69 461006589 : raised(axis) = raised(axis) + 1
70 461006589 : lowered(axis) = MAX(0, lowered(axis) - 1)
71 : map(:, ico - offset) = [ico - source_offset, &
72 : coset(raised(1), raised(2), raised(3)) - source_offset, &
73 2411356581 : MAX(1, coset(lowered(1), lowered(2), lowered(3)) - source_offset), power(axis)]
74 : END DO
75 : END IF
76 152590623 : END SUBROUTINE cartesian_ladder_map
77 :
78 : #:for suffix, integral_type in [('real', 'REAL'), ('complex', 'COMPLEX')]
79 : ! **************************************************************************************************
80 : !> \brief Apply one coordinate or basis-center derivative ladder to indexed integrals.
81 : !> \param source integral block including the raised functions
82 : !> \param values transformed integral block
83 : !> \param map base, raised, lowered and degree indices
84 : !> \param derivative select a basis-center derivative instead of coordinate multiplication
85 : !> \param exponent exponent of the differentiated Gaussian
86 : !> \param origin displacement of the Gaussian center from the coordinate reference point
87 : !> \param right apply the map to columns rather than rows
88 : ! **************************************************************************************************
89 440903700 : PURE SUBROUTINE cartesian_ladder_${suffix}$ (source, values, map, derivative, exponent, origin, right)
90 : ${integral_type}$ (KIND=dp), DIMENSION(:, :), INTENT(IN) :: source
91 : ${integral_type}$ (KIND=dp), DIMENSION(:, :), INTENT(OUT) :: values
92 : INTEGER, DIMENSION(:, :), INTENT(IN) :: map
93 : LOGICAL, INTENT(IN) :: derivative
94 : REAL(KIND=dp), INTENT(IN) :: exponent, origin
95 : LOGICAL, INTENT(IN) :: right
96 :
97 : INTEGER :: i, j, second_index
98 881807400 : REAL(KIND=dp) :: scale, second_scale(SIZE(map, 2))
99 :
100 440903700 : IF (derivative) THEN
101 439150056 : scale = 2.0_dp*exponent
102 1759218687 : second_scale = -REAL(map(4, :), dp)
103 : second_index = 3
104 : ELSE
105 1753644 : scale = 1.0_dp
106 4786998 : second_scale = origin
107 : second_index = 1
108 : END IF
109 :
110 : ! Keep the first (contiguous) dimension innermost on both sides.
111 440903700 : IF (right) THEN
112 906334554 : DO j = 1, SIZE(map, 2)
113 3904870353 : DO i = 1, SIZE(values, 1)
114 : values(i, j) = scale*source(i, map(2, j)) + &
115 3696870393 : second_scale(j)*source(i, map(second_index, j))
116 : END DO
117 : END DO
118 : ELSE
119 893316069 : DO j = 1, SIZE(values, 2)
120 3384782178 : DO i = 1, SIZE(map, 2)
121 : values(i, j) = scale*source(map(2, i), j) + &
122 3151878438 : second_scale(i)*source(map(second_index, i), j)
123 : END DO
124 : END DO
125 : END IF
126 440903700 : END SUBROUTINE cartesian_ladder_${suffix}$
127 :
128 : ! **************************************************************************************************
129 : !> \brief Apply an ordered pair of coordinate or derivative ladders to integral blocks.
130 : !> Maps describe source-to-work and work-to-output transformations. The first result
131 : !> is reused for all requested second axes. On the same basis function the resulting
132 : !> operator product is O_first O_second: the second operator acts on the Gaussian first.
133 : !> A derivative represents minus the electronic gradient, equal to a basis-center
134 : !> derivative on an unmodified Gaussian. Exponents and supplied origins are fixed
135 : !> coefficients; total center derivatives of coordinate prefactors require their
136 : !> parameter dependence to be accounted for by the caller.
137 : !> \param source integral block including all required raised and lowered functions
138 : !> \param values transformed components in the order specified by components
139 : !> \param map_first first maps, indexed by map entry, output function and axis
140 : !> \param map_second second maps into work, with the same axis convention
141 : !> \param components first and second map-axis indices; use a single map axis for a 1D table
142 : !> \param derivative select a derivative instead of coordinate multiplication for each step
143 : !> \param exponents Gaussian exponents for the two steps, unused for coordinate multiplication
144 : !> \param right select columns instead of rows for each step
145 : !> \param work intermediate block, retained after the last requested first axis
146 : !> \param origin optional center-to-reference displacements indexed by map axis and step; defaults to zero
147 : ! **************************************************************************************************
148 80405740 : PURE SUBROUTINE apply_cartesian_operator_pair_${suffix}$ (source, values, map_first, map_second, &
149 40202870 : components, derivative, exponents, right, work, origin)
150 : ${integral_type}$ (KIND=dp), DIMENSION(:, :), INTENT(IN) :: source
151 : ${integral_type}$ (KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: values
152 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: map_first, map_second
153 : INTEGER, DIMENSION(:, :), INTENT(IN) :: components
154 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: exponents
155 : LOGICAL, DIMENSION(2), INTENT(IN) :: derivative, right
156 : ${integral_type}$ (KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: work
157 : REAL(KIND=dp), DIMENSION(3, 2), INTENT(IN), OPTIONAL :: origin
158 :
159 : INTEGER :: axis, icomp
160 : REAL(KIND=dp) :: center(3, 2)
161 :
162 40202870 : center = 0.0_dp
163 40202870 : IF (PRESENT(origin)) center = origin
164 80642792 : DO axis = 1, SIZE(map_first, 3)
165 40973976 : IF (.NOT. ANY(components(1, :) == axis)) CYCLE
166 : CALL cartesian_ladder(source, work, map_first(:, :, axis), &
167 40439922 : derivative(1), exponents(1), center(axis, 1), right(1))
168 123213566 : DO icomp = 1, SIZE(components, 2)
169 42570774 : IF (components(1, icomp) /= axis) CYCLE
170 : CALL cartesian_ladder(work, values(:, :, icomp), map_second(:, :, components(2, icomp)), &
171 83010696 : derivative(2), exponents(2), center(components(2, icomp), 2), right(2))
172 : END DO
173 : END DO
174 40202870 : END SUBROUTINE apply_cartesian_operator_pair_${suffix}$
175 : #:endfor
176 :
177 : ! **************************************************************************************************
178 : !> \brief Apply a Cartesian operator on either side of a primitive integral block.
179 : !> The input includes one extra shell on both sides. Coordinate multiplication
180 : !> uses the supplied shell-range layout; derivative inputs extend the lowest
181 : !> shell by one on the selected side for lowering. Screened primitive blocks
182 : !> must be zero in the input, including their raised and lowered functions.
183 : !> \param la_max highest output shell on A
184 : !> \param la_min lowest output shell on A
185 : !> \param npgfa number of primitives on A
186 : !> \param lb_max highest output shell on B
187 : !> \param lb_min lowest output shell on B
188 : !> \param npgfb number of primitives on B
189 : !> \param ab input integral block with raised shells
190 : !> \param values three Cartesian components in the output shell-range layout
191 : !> \param right act on the ket rather than the bra
192 : !> \param derivative select a basis-center derivative
193 : !> \param exponents Gaussian exponents on the selected side, required for derivatives
194 : !> \param origin displacement of the selected center from the coordinate reference point
195 : ! **************************************************************************************************
196 5213474 : SUBROUTINE apply_cartesian_operator(la_max, la_min, npgfa, lb_max, lb_min, npgfb, &
197 5213474 : ab, values, right, derivative, exponents, origin)
198 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa, lb_max, lb_min, &
199 : npgfb
200 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ab
201 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: values
202 : LOGICAL, INTENT(IN) :: right, derivative
203 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: exponents
204 : REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: origin
205 :
206 : INTEGER :: axis, ia, iap, ib, ibp, ipgf, jpgf, n, &
207 : na, nap, nb, nbp, source_min_a, &
208 : source_min_b
209 5213474 : INTEGER, ALLOCATABLE :: map(:, :, :)
210 : LOGICAL :: valid_shape
211 : REAL(KIND=dp) :: center(3), exponent
212 :
213 5213474 : na = ncoset(la_max) - ncoset(la_min - 1)
214 5213474 : nb = ncoset(lb_max) - ncoset(lb_min - 1)
215 5213474 : source_min_a = la_min
216 5213474 : source_min_b = lb_min
217 5213474 : IF (derivative) THEN
218 5212277 : IF (right) THEN
219 5169079 : source_min_b = MAX(0, lb_min - 1)
220 : ELSE
221 43198 : source_min_a = MAX(0, la_min - 1)
222 : END IF
223 : END IF
224 5213474 : nap = ncoset(la_max + 1) - ncoset(source_min_a - 1)
225 5213474 : nbp = ncoset(lb_max + 1) - ncoset(source_min_b - 1)
226 5213474 : IF (derivative) THEN
227 5212277 : CPASSERT(PRESENT(exponents))
228 : END IF
229 15640422 : valid_shape = ALL(SHAPE(ab) >= [nap*npgfa, nbp*npgfb])
230 0 : CPASSERT(valid_shape)
231 20853896 : valid_shape = ALL(SHAPE(values) >= [na*npgfa, nb*npgfb, 3])
232 0 : CPASSERT(valid_shape)
233 5213474 : center = 0.0_dp
234 5213474 : IF (PRESENT(origin)) center = origin
235 5213474 : n = MERGE(nb, na, right)
236 20853896 : ALLOCATE (map(4, n, 3))
237 20853896 : DO axis = 1, 3
238 20853896 : IF (right) THEN
239 15508092 : CALL cartesian_ladder_map(lb_max, lb_min, axis, map(:, :, axis), source_min_b)
240 : ELSE
241 132330 : CALL cartesian_ladder_map(la_max, la_min, axis, map(:, :, axis), source_min_a)
242 : END IF
243 : END DO
244 :
245 1202207294 : values = 0.0_dp
246 16385728 : DO ipgf = 1, npgfa
247 11172254 : ia = (ipgf - 1)*na
248 11172254 : iap = (ipgf - 1)*nap
249 45843562 : DO jpgf = 1, npgfb
250 29457834 : ib = (jpgf - 1)*nb
251 29457834 : ibp = (jpgf - 1)*nbp
252 29457834 : exponent = 0.0_dp
253 29711406 : IF (derivative) exponent = exponents(MERGE(jpgf, ipgf, right))
254 129003590 : DO axis = 1, 3
255 : CALL cartesian_ladder(ab(iap + 1:iap + nap, ibp + 1:ibp + nbp), &
256 : values(ia + 1:ia + na, ib + 1:ib + nb, axis), map(:, :, axis), &
257 117831336 : derivative, exponent, center(axis), right)
258 : END DO
259 : END DO
260 : END DO
261 5213474 : END SUBROUTINE apply_cartesian_operator
262 :
263 : #:for suffix, integral_type in [('real', 'REAL'), ('complex', 'COMPLEX')]
264 : ! **************************************************************************************************
265 : !> \brief Apply the common ladder to separable one-dimensional integral tables.
266 : !> \param source input tables indexed by the degrees on A and B and the Cartesian axis
267 : !> \param values transformed tables; entries outside the requested range are untouched
268 : !> \param la_max highest output degree on A
269 : !> \param lb_max highest output degree on B
270 : !> \param right act on B rather than A
271 : !> \param derivative select a basis-center derivative
272 : !> \param exponent Gaussian exponent on the selected side
273 : !> \param origin displacement of that center from the coordinate reference point
274 : ! **************************************************************************************************
275 36687303 : SUBROUTINE apply_cartesian_1d_${suffix}$ (source, values, la_max, lb_max, right, derivative, exponent, origin)
276 : ${integral_type}$ (KIND=dp), DIMENSION(0:, 0:, :), INTENT(IN) :: source
277 : ${integral_type}$ (KIND=dp), DIMENSION(0:, 0:, :), INTENT(INOUT) :: values
278 : INTEGER, INTENT(IN) :: la_max, lb_max
279 : LOGICAL, INTENT(IN) :: right, derivative
280 : REAL(KIND=dp), INTENT(IN) :: exponent, origin(3)
281 :
282 36687303 : INTEGER :: axis, map(4, MERGE(lb_max, la_max, right) + 1)
283 :
284 72810752 : CALL cartesian_ladder_map(MERGE(lb_max, la_max, right), 0, 0, map)
285 146749212 : DO axis = 1, 3
286 : CALL cartesian_ladder(source(:, :, axis), values(0:la_max, 0:lb_max, axis), map, &
287 146749212 : derivative, exponent, origin(axis), right)
288 : END DO
289 36687303 : END SUBROUTINE apply_cartesian_1d_${suffix}$
290 : #:endfor
291 :
292 : ! **************************************************************************************************
293 : !> \brief Form requested coordinate powers from separable Cartesian integral tables.
294 : !> Only requested powers are stored. With the origin on A, each power is a
295 : !> direct shift of the bra degree; other origins use the binomial expansion.
296 : !> \param source input tables with enough raised bra degrees for the largest power
297 : !> \param values tables indexed by bra degree, ket degree, requested power and axis
298 : !> \param la_max highest output degree on A
299 : !> \param lb_max highest output degree on B
300 : !> \param powers nonnegative coordinate powers to evaluate
301 : !> \param origin displacement of A from the coordinate reference point
302 : ! **************************************************************************************************
303 1462329 : SUBROUTINE cartesian_polynomial_1d(source, values, la_max, lb_max, powers, origin)
304 : REAL(KIND=dp), DIMENSION(0:, 0:, :), INTENT(IN) :: source
305 : REAL(KIND=dp), DIMENSION(0:, 0:, 0:, :), &
306 : INTENT(OUT) :: values
307 : INTEGER, INTENT(IN) :: la_max, lb_max, powers(0:)
308 : REAL(KIND=dp), INTENT(IN) :: origin(3)
309 :
310 : INTEGER :: axis, i, k, p
311 : REAL(KIND=dp) :: coefficient
312 :
313 5849316 : DO axis = 1, 3
314 17905656 : DO i = 0, SIZE(powers) - 1
315 12056340 : p = powers(i)
316 16443327 : IF (ABS(origin(axis)) <= 0.0_dp) THEN
317 7194888 : values(:, :, i, axis) = source(p:p + la_max, 0:lb_max, axis)
318 : ELSE
319 113336541 : values(:, :, i, axis) = 0.0_dp
320 32960719 : DO k = 0, p
321 21654423 : coefficient = binomial_gen(REAL(p, dp), k)*origin(axis)**(p - k)
322 : values(:, :, i, axis) = values(:, :, i, axis) + &
323 232379941 : coefficient*source(k:k + la_max, 0:lb_max, axis)
324 : END DO
325 : END IF
326 : END DO
327 : END DO
328 1462329 : END SUBROUTINE cartesian_polynomial_1d
329 :
330 : END MODULE ai_cartesian_operators
|
