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, 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_1d
32 : MODULE PROCEDURE apply_cartesian_1d_real, apply_cartesian_1d_complex
33 : END INTERFACE apply_cartesian_1d
34 :
35 : CONTAINS
36 :
37 : ! **************************************************************************************************
38 : !> \brief Construct base, raised, lowered and degree indices for a Cartesian ladder.
39 : !> \param lmax highest output shell or one-dimensional degree
40 : !> \param lmin lowest output shell (zero for one-dimensional tables)
41 : !> \param axis Cartesian axis, or zero for a one-dimensional table
42 : !> \param map rows contain base, raised, lowered and degree indices
43 : !> \param source_lmin optional lowest input shell; defaults to lmin for Cartesian blocks
44 : ! **************************************************************************************************
45 80900877 : SUBROUTINE cartesian_ladder_map(lmax, lmin, axis, map, source_lmin)
46 : INTEGER, INTENT(IN) :: lmax, lmin, axis
47 : INTEGER, DIMENSION(:, :), INTENT(OUT) :: map
48 : INTEGER, INTENT(IN), OPTIONAL :: source_lmin
49 :
50 : INTEGER :: i, ico, lowered(3), offset, power(3), &
51 : raised(3), source_offset
52 :
53 80900877 : IF (axis == 0) THEN
54 172714981 : DO i = 0, lmax
55 610356725 : map(:, i + 1) = [i + 1, i + 2, MAX(1, i), i]
56 : END DO
57 : ELSE
58 17596332 : offset = ncoset(lmin - 1)
59 17596332 : source_offset = offset
60 17596332 : IF (PRESENT(source_lmin)) source_offset = ncoset(source_lmin - 1)
61 81686925 : DO ico = offset + 1, ncoset(lmax)
62 256362372 : power = indco(:, ico)
63 64090593 : raised = power
64 64090593 : lowered = power
65 64090593 : raised(axis) = raised(axis) + 1
66 64090593 : lowered(axis) = MAX(0, lowered(axis) - 1)
67 : map(:, ico - offset) = [ico - source_offset, &
68 : coset(raised(1), raised(2), raised(3)) - source_offset, &
69 338049297 : MAX(1, coset(lowered(1), lowered(2), lowered(3)) - source_offset), power(axis)]
70 : END DO
71 : END IF
72 80900877 : END SUBROUTINE cartesian_ladder_map
73 :
74 : #:for suffix, integral_type in [('real', 'REAL'), ('complex', 'COMPLEX')]
75 : ! **************************************************************************************************
76 : !> \brief Apply one coordinate or basis-center derivative ladder to indexed integrals.
77 : !> \param source integral block including the raised functions
78 : !> \param values transformed integral block
79 : !> \param map base, raised, lowered and degree indices
80 : !> \param derivative select a basis-center derivative instead of coordinate multiplication
81 : !> \param exponent exponent of the differentiated Gaussian
82 : !> \param origin displacement of the Gaussian center from the coordinate reference point
83 : !> \param right apply the map to columns rather than rows
84 : ! **************************************************************************************************
85 283823646 : PURE SUBROUTINE cartesian_ladder_${suffix}$ (source, values, map, derivative, exponent, origin, right)
86 : ${integral_type}$ (KIND=dp), DIMENSION(:, :), INTENT(IN) :: source
87 : ${integral_type}$ (KIND=dp), DIMENSION(:, :), INTENT(OUT) :: values
88 : INTEGER, DIMENSION(:, :), INTENT(IN) :: map
89 : LOGICAL, INTENT(IN) :: derivative
90 : REAL(KIND=dp), INTENT(IN) :: exponent, origin
91 : LOGICAL, INTENT(IN) :: right
92 :
93 : INTEGER :: i, j, second_index
94 567647292 : REAL(KIND=dp) :: scale, second_scale(SIZE(map, 2))
95 :
96 283823646 : IF (derivative) THEN
97 282153618 : scale = 2.0_dp*exponent
98 902893935 : second_scale = -REAL(map(4, :), dp)
99 : second_index = 3
100 : ELSE
101 1670028 : scale = 1.0_dp
102 4493268 : second_scale = origin
103 : second_index = 1
104 : END IF
105 :
106 : ! Keep the first (contiguous) dimension innermost on both sides.
107 283823646 : IF (right) THEN
108 477413865 : DO j = 1, SIZE(map, 2)
109 1651479486 : DO i = 1, SIZE(values, 1)
110 : values(i, j) = scale*source(i, map(2, j)) + &
111 1522136931 : second_scale(j)*source(i, map(second_index, j))
112 : END DO
113 : END DO
114 : ELSE
115 464096466 : DO j = 1, SIZE(values, 2)
116 1131452649 : DO i = 1, SIZE(map, 2)
117 : values(i, j) = scale*source(map(2, i), j) + &
118 976971558 : second_scale(i)*source(map(second_index, i), j)
119 : END DO
120 : END DO
121 : END IF
122 283823646 : END SUBROUTINE cartesian_ladder_${suffix}$
123 : #:endfor
124 :
125 : ! **************************************************************************************************
126 : !> \brief Apply a Cartesian operator on either side of a primitive integral block.
127 : !> The input includes one extra shell on both sides. Coordinate multiplication
128 : !> uses the supplied shell-range layout; derivative inputs extend the lowest
129 : !> shell by one on the selected side for lowering. Screened primitive blocks
130 : !> must be zero in the input, including their raised and lowered functions.
131 : !> \param la_max highest output shell on A
132 : !> \param la_min lowest output shell on A
133 : !> \param npgfa number of primitives on A
134 : !> \param lb_max highest output shell on B
135 : !> \param lb_min lowest output shell on B
136 : !> \param npgfb number of primitives on B
137 : !> \param ab input integral block with raised shells
138 : !> \param values three Cartesian components in the output shell-range layout
139 : !> \param right act on the ket rather than the bra
140 : !> \param derivative select a basis-center derivative
141 : !> \param exponents Gaussian exponents on the selected side, required for derivatives
142 : !> \param origin displacement of the selected center from the coordinate reference point
143 : ! **************************************************************************************************
144 5226074 : SUBROUTINE apply_cartesian_operator(la_max, la_min, npgfa, lb_max, lb_min, npgfb, &
145 5226074 : ab, values, right, derivative, exponents, origin)
146 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa, lb_max, lb_min, &
147 : npgfb
148 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ab
149 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: values
150 : LOGICAL, INTENT(IN) :: right, derivative
151 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: exponents
152 : REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: origin
153 :
154 : INTEGER :: axis, ia, iap, ib, ibp, ipgf, jpgf, n, &
155 : na, nap, nb, nbp, source_min_a, &
156 : source_min_b
157 5226074 : INTEGER, ALLOCATABLE :: map(:, :, :)
158 : LOGICAL :: valid_shape
159 : REAL(KIND=dp) :: center(3), exponent
160 :
161 5226074 : na = ncoset(la_max) - ncoset(la_min - 1)
162 5226074 : nb = ncoset(lb_max) - ncoset(lb_min - 1)
163 5226074 : source_min_a = la_min
164 5226074 : source_min_b = lb_min
165 5226074 : IF (derivative) THEN
166 5224877 : IF (right) THEN
167 5175379 : source_min_b = MAX(0, lb_min - 1)
168 : ELSE
169 49498 : source_min_a = MAX(0, la_min - 1)
170 : END IF
171 : END IF
172 5226074 : nap = ncoset(la_max + 1) - ncoset(source_min_a - 1)
173 5226074 : nbp = ncoset(lb_max + 1) - ncoset(source_min_b - 1)
174 5226074 : IF (derivative) THEN
175 5224877 : CPASSERT(PRESENT(exponents))
176 : END IF
177 15678222 : valid_shape = ALL(SHAPE(ab) >= [nap*npgfa, nbp*npgfb])
178 0 : CPASSERT(valid_shape)
179 20904296 : valid_shape = ALL(SHAPE(values) >= [na*npgfa, nb*npgfb, 3])
180 0 : CPASSERT(valid_shape)
181 5226074 : center = 0.0_dp
182 5226074 : IF (PRESENT(origin)) center = origin
183 5226074 : n = MERGE(nb, na, right)
184 20904296 : ALLOCATE (map(4, n, 3))
185 20904296 : DO axis = 1, 3
186 20904296 : IF (right) THEN
187 15526992 : CALL cartesian_ladder_map(lb_max, lb_min, axis, map(:, :, axis), source_min_b)
188 : ELSE
189 151230 : CALL cartesian_ladder_map(la_max, la_min, axis, map(:, :, axis), source_min_a)
190 : END IF
191 : END DO
192 :
193 1205277500 : values = 0.0_dp
194 16433944 : DO ipgf = 1, npgfa
195 11207870 : ia = (ipgf - 1)*na
196 11207870 : iap = (ipgf - 1)*nap
197 45993760 : DO jpgf = 1, npgfb
198 29559816 : ib = (jpgf - 1)*nb
199 29559816 : ibp = (jpgf - 1)*nbp
200 29559816 : exponent = 0.0_dp
201 29864379 : IF (derivative) exponent = exponents(MERGE(jpgf, ipgf, right))
202 129447134 : DO axis = 1, 3
203 : CALL cartesian_ladder(ab(iap + 1:iap + nap, ibp + 1:ibp + nbp), &
204 : values(ia + 1:ia + na, ib + 1:ib + nb, axis), map(:, :, axis), &
205 118239264 : derivative, exponent, center(axis), right)
206 : END DO
207 : END DO
208 : END DO
209 5226074 : END SUBROUTINE apply_cartesian_operator
210 :
211 : #:for suffix, integral_type in [('real', 'REAL'), ('complex', 'COMPLEX')]
212 : ! **************************************************************************************************
213 : !> \brief Apply the common ladder to separable one-dimensional integral tables.
214 : !> \param source input tables indexed by the degrees on A and B and the Cartesian axis
215 : !> \param values transformed tables; entries outside the requested range are untouched
216 : !> \param la_max highest output degree on A
217 : !> \param lb_max highest output degree on B
218 : !> \param right act on B rather than A
219 : !> \param derivative select a basis-center derivative
220 : !> \param exponent Gaussian exponent on the selected side
221 : !> \param origin displacement of that center from the coordinate reference point
222 : ! **************************************************************************************************
223 63304545 : SUBROUTINE apply_cartesian_1d_${suffix}$ (source, values, la_max, lb_max, right, derivative, exponent, origin)
224 : ${integral_type}$ (KIND=dp), DIMENSION(0:, 0:, :), INTENT(IN) :: source
225 : ${integral_type}$ (KIND=dp), DIMENSION(0:, 0:, :), INTENT(INOUT) :: values
226 : INTEGER, INTENT(IN) :: la_max, lb_max
227 : LOGICAL, INTENT(IN) :: right, derivative
228 : REAL(KIND=dp), INTENT(IN) :: exponent, origin(3)
229 :
230 63304545 : INTEGER :: axis, map(4, MERGE(lb_max, la_max, right) + 1)
231 :
232 112735566 : CALL cartesian_ladder_map(MERGE(lb_max, la_max, right), 0, 0, map)
233 253218180 : DO axis = 1, 3
234 : CALL cartesian_ladder(source(:, :, axis), values(0:la_max, 0:lb_max, axis), map, &
235 253218180 : derivative, exponent, origin(axis), right)
236 : END DO
237 63304545 : END SUBROUTINE apply_cartesian_1d_${suffix}$
238 : #:endfor
239 :
240 : ! **************************************************************************************************
241 : !> \brief Form requested coordinate powers from separable Cartesian integral tables.
242 : !> Only requested powers are stored. With the origin on A, each power is a
243 : !> direct shift of the bra degree; other origins use the binomial expansion.
244 : !> \param source input tables with enough raised bra degrees for the largest power
245 : !> \param values tables indexed by bra degree, ket degree, requested power and axis
246 : !> \param la_max highest output degree on A
247 : !> \param lb_max highest output degree on B
248 : !> \param powers nonnegative coordinate powers to evaluate
249 : !> \param origin displacement of A from the coordinate reference point
250 : ! **************************************************************************************************
251 1472390 : SUBROUTINE cartesian_polynomial_1d(source, values, la_max, lb_max, powers, origin)
252 : REAL(KIND=dp), DIMENSION(0:, 0:, :), INTENT(IN) :: source
253 : REAL(KIND=dp), DIMENSION(0:, 0:, 0:, :), &
254 : INTENT(OUT) :: values
255 : INTEGER, INTENT(IN) :: la_max, lb_max, powers(0:)
256 : REAL(KIND=dp), INTENT(IN) :: origin(3)
257 :
258 : INTEGER :: axis, i, k, p
259 : REAL(KIND=dp) :: coefficient
260 :
261 5889560 : DO axis = 1, 3
262 18006266 : DO i = 0, SIZE(powers) - 1
263 12116706 : p = powers(i)
264 16533876 : IF (ABS(origin(axis)) <= 0.0_dp) THEN
265 7135766 : values(:, :, i, axis) = source(p:p + la_max, 0:lb_max, axis)
266 : ELSE
267 114026149 : values(:, :, i, axis) = 0.0_dp
268 33125099 : DO k = 0, p
269 21753051 : coefficient = binomial_gen(REAL(p, dp), k)*origin(axis)**(p - k)
270 : values(:, :, i, axis) = values(:, :, i, axis) + &
271 233480105 : coefficient*source(k:k + la_max, 0:lb_max, axis)
272 : END DO
273 : END IF
274 : END DO
275 : END DO
276 1472390 : END SUBROUTINE cartesian_polynomial_1d
277 :
278 : END MODULE ai_cartesian_operators
|
