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 Ewald sums to represent integrals in direct and reciprocal lattice.
10 : !> \par History
11 : !> 2015 09 created
12 : !> \author Patrick Seewald
13 : ! **************************************************************************************************
14 :
15 : MODULE eri_mme_lattice_summation
16 : #:include "eri_mme_unroll.fypp"
17 :
18 : USE ao_util, ONLY: exp_radius
19 : USE eri_mme_gaussian, ONLY: create_gaussian_overlap_dist_to_hermite, &
20 : create_hermite_to_cartesian, &
21 : eri_mme_coulomb, eri_mme_yukawa, &
22 : eri_mme_longrange
23 : USE kinds, ONLY: dp, &
24 : int_8
25 : USE mathconstants, ONLY: gaussi, &
26 : pi, &
27 : twopi
28 : USE orbital_pointers, ONLY: coset, &
29 : indco, &
30 : ncoset
31 : #include "../base/base_uses.f90"
32 :
33 : IMPLICIT NONE
34 :
35 : PRIVATE
36 :
37 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
38 :
39 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eri_mme_lattice_summation'
40 :
41 : PUBLIC :: &
42 : ellipsoid_bounds, &
43 : eri_mme_2c_get_bounds, &
44 : eri_mme_2c_get_rads, &
45 : eri_mme_3c_get_bounds, &
46 : eri_mme_3c_get_rads, &
47 : get_l, &
48 : pgf_sum_2c_gspace_1d, &
49 : pgf_sum_2c_gspace_1d_deltal, &
50 : pgf_sum_2c_gspace_3d, &
51 : pgf_sum_2c_rspace_1d, &
52 : pgf_sum_2c_rspace_3d, &
53 : pgf_sum_3c_1d, &
54 : pgf_sum_3c_3d
55 :
56 : INTEGER, PARAMETER, PRIVATE :: exp_w = 50, div_w = 10
57 :
58 : CONTAINS
59 :
60 : ! **************************************************************************************************
61 : !> \brief Get summation radii for 2c integrals
62 : !> \param la_max ...
63 : !> \param lb_max ...
64 : !> \param zeta ...
65 : !> \param zetb ...
66 : !> \param a_mm ...
67 : !> \param G_min ...
68 : !> \param R_min ...
69 : !> \param sum_precision ...
70 : !> \param G_rad ...
71 : !> \param R_rad ...
72 : ! **************************************************************************************************
73 948038 : SUBROUTINE eri_mme_2c_get_rads(la_max, lb_max, zeta, zetb, a_mm, G_min, R_min, sum_precision, G_rad, R_rad)
74 : INTEGER, INTENT(IN) :: la_max, lb_max
75 : REAL(KIND=dp), INTENT(IN) :: zeta, zetb, a_mm, G_min, R_min, &
76 : sum_precision
77 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: G_rad, R_rad
78 :
79 : INTEGER :: l_max
80 : REAL(KIND=dp) :: alpha_G, alpha_R, G_res, R_res
81 :
82 948038 : l_max = la_max + lb_max
83 948038 : alpha_G = a_mm + 0.25_dp/zeta + 0.25_dp/zetb
84 948038 : alpha_R = 0.25_dp/alpha_G
85 :
86 948038 : G_res = 0.5_dp*G_min
87 948038 : R_res = 0.5_dp*R_min
88 :
89 948038 : IF (PRESENT(G_rad)) G_rad = exp_radius(l_max, alpha_G, sum_precision, 1.0_dp, epsabs=G_res)
90 948038 : IF (PRESENT(R_rad)) R_rad = exp_radius(l_max, alpha_R, sum_precision, 1.0_dp, epsabs=R_res)
91 :
92 948038 : END SUBROUTINE eri_mme_2c_get_rads
93 :
94 : ! **************************************************************************************************
95 : !> \brief Get summation radii for 3c integrals
96 : !> \param la_max ...
97 : !> \param lb_max ...
98 : !> \param lc_max ...
99 : !> \param zeta ...
100 : !> \param zetb ...
101 : !> \param zetc ...
102 : !> \param a_mm ...
103 : !> \param G_min ...
104 : !> \param R_min ...
105 : !> \param sum_precision ...
106 : !> \param G_rads_1 ...
107 : !> \param R_rads_2 ...
108 : !> \param R_rads_3 ...
109 : ! **************************************************************************************************
110 2411094 : SUBROUTINE eri_mme_3c_get_rads(la_max, lb_max, lc_max, zeta, zetb, zetc, a_mm, G_min, R_min, &
111 : sum_precision, G_rads_1, R_rads_2, R_rads_3)
112 : INTEGER, INTENT(IN) :: la_max, lb_max, lc_max
113 : REAL(KIND=dp), INTENT(IN) :: zeta, zetb, zetc, a_mm
114 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: G_min, R_min
115 : REAL(KIND=dp), INTENT(IN) :: sum_precision
116 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: G_rads_1, R_rads_2
117 : REAL(KIND=dp), DIMENSION(2), INTENT(OUT), OPTIONAL :: R_rads_3
118 :
119 : REAL(KIND=dp) :: alpha, alpha_R, beta, G_res, gamma, R_res
120 :
121 : ! exponents in G space
122 2411094 : alpha = 0.25_dp/zeta
123 2411094 : beta = 0.25_dp/zetb
124 2411094 : gamma = 0.25_dp/zetc + a_mm
125 :
126 : ! Summation radii and number of summands for all lattice summation methods
127 : ! sum method 1
128 2411094 : IF (PRESENT(G_rads_1)) THEN
129 2411094 : G_res = 0.5_dp*G_min
130 2411094 : G_rads_1(1) = exp_radius(la_max, alpha, sum_precision, 1.0_dp, G_res)
131 2411094 : G_rads_1(2) = exp_radius(lb_max, beta, sum_precision, 1.0_dp, G_res)
132 2411094 : G_rads_1(3) = exp_radius(lc_max, gamma, sum_precision, 1.0_dp, G_res)
133 : END IF
134 :
135 : ! sum method 2
136 2411094 : IF (PRESENT(R_rads_2)) THEN
137 2256540 : R_res = 0.5_dp*R_min
138 2256540 : R_rads_2(1) = exp_radius(lb_max + lc_max, 0.25_dp/(beta + gamma), sum_precision, 1.0_dp, epsabs=R_res)
139 2256540 : R_rads_2(2) = exp_radius(lc_max + la_max, 0.25_dp/(alpha + gamma), sum_precision, 1.0_dp, epsabs=R_res)
140 2256540 : R_rads_2(3) = exp_radius(lb_max + la_max, 0.25_dp/(alpha + beta), sum_precision, 1.0_dp, epsabs=R_res)
141 : END IF
142 :
143 : ! sum method 3
144 :
145 2411094 : IF (PRESENT(R_rads_3)) THEN
146 2256540 : R_res = 0.5_dp*R_min
147 2256540 : alpha_R = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm)
148 2256540 : R_rads_3(1) = exp_radius(la_max + lb_max, zeta*zetb/(zeta + zetb), sum_precision, 1.0_dp, R_res)
149 2256540 : R_rads_3(2) = exp_radius(la_max + lb_max + lc_max, alpha_R, sum_precision, 1.0_dp, R_res)
150 : END IF
151 :
152 2411094 : END SUBROUTINE eri_mme_3c_get_rads
153 :
154 : ! **************************************************************************************************
155 : !> \brief Get summation bounds for 2c integrals
156 : !> \param hmat ...
157 : !> \param h_inv ...
158 : !> \param vol ...
159 : !> \param is_ortho ...
160 : !> \param G_min ...
161 : !> \param R_min ...
162 : !> \param la_max ...
163 : !> \param lb_max ...
164 : !> \param zeta ...
165 : !> \param zetb ...
166 : !> \param a_mm ...
167 : !> \param sum_precision ...
168 : !> \param n_sum_1d ...
169 : !> \param n_sum_3d ...
170 : !> \param G_bounds ...
171 : !> \param G_rad ...
172 : !> \param R_bounds ...
173 : !> \param R_rad ...
174 : ! **************************************************************************************************
175 892104 : SUBROUTINE eri_mme_2c_get_bounds(hmat, h_inv, vol, is_ortho, G_min, R_min, la_max, lb_max, &
176 : zeta, zetb, a_mm, sum_precision, n_sum_1d, n_sum_3d, &
177 : G_bounds, G_rad, R_bounds, R_rad)
178 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv
179 : REAL(KIND=dp), INTENT(IN) :: vol
180 : LOGICAL, INTENT(IN) :: is_ortho
181 : REAL(KIND=dp), INTENT(IN) :: G_min, R_min
182 : INTEGER, INTENT(IN) :: la_max, lb_max
183 : REAL(KIND=dp), INTENT(IN) :: zeta, zetb, a_mm, sum_precision
184 : INTEGER(KIND=int_8), DIMENSION(2, 3), INTENT(OUT) :: n_sum_1d
185 : INTEGER(KIND=int_8), DIMENSION(2), INTENT(OUT) :: n_sum_3d
186 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: G_bounds
187 : REAL(KIND=dp), INTENT(OUT) :: G_rad
188 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: R_bounds
189 : REAL(KIND=dp), INTENT(OUT) :: R_rad
190 :
191 : INTEGER :: i_xyz
192 : REAL(KIND=dp) :: ns_G, ns_R, vol_G
193 :
194 892104 : CALL eri_mme_2c_get_rads(la_max, lb_max, zeta, zetb, a_mm, G_min, R_min, sum_precision, G_rad, R_rad)
195 :
196 11597352 : G_bounds = ellipsoid_bounds(G_rad, TRANSPOSE(hmat)/(2.0_dp*pi))
197 892104 : R_bounds = ellipsoid_bounds(R_rad, h_inv)
198 :
199 892104 : vol_G = twopi**3/vol
200 :
201 892104 : IF (is_ortho) THEN
202 2981984 : DO i_xyz = 1, 3
203 2236488 : ns_G = 2.0_dp*G_bounds(i_xyz)
204 2236488 : ns_R = 2.0_dp*R_bounds(i_xyz)
205 2236488 : n_sum_1d(1, i_xyz) = nsum_2c_gspace_1d(ns_G, la_max, lb_max)
206 2981984 : n_sum_1d(2, i_xyz) = nsum_2c_rspace_1d(ns_R, la_max, lb_max)
207 : END DO
208 : ELSE
209 146608 : ns_G = 4._dp/3._dp*pi*G_rad**3/vol_G
210 146608 : ns_R = 4._dp/3._dp*pi*R_rad**3/vol
211 146608 : n_sum_3d(1) = nsum_2c_gspace_3d(ns_G, la_max, lb_max)
212 146608 : n_sum_3d(2) = nsum_2c_rspace_3d(ns_R, la_max, lb_max)
213 : END IF
214 :
215 892104 : END SUBROUTINE eri_mme_2c_get_bounds
216 :
217 : ! **************************************************************************************************
218 : !> \brief Get summation bounds for 3c integrals
219 : !> \param hmat ...
220 : !> \param h_inv ...
221 : !> \param vol ...
222 : !> \param is_ortho ...
223 : !> \param G_min ...
224 : !> \param R_min ...
225 : !> \param la_max ...
226 : !> \param lb_max ...
227 : !> \param lc_max ...
228 : !> \param zeta ...
229 : !> \param zetb ...
230 : !> \param zetc ...
231 : !> \param a_mm ...
232 : !> \param sum_precision ...
233 : !> \param n_sum_1d ...
234 : !> \param n_sum_3d ...
235 : !> \param G_bounds_1 ...
236 : !> \param G_rads_1 ...
237 : !> \param R_bounds_2 ...
238 : !> \param R_rads_2 ...
239 : !> \param R_bounds_3 ...
240 : !> \param R_rads_3 ...
241 : ! **************************************************************************************************
242 2256540 : SUBROUTINE eri_mme_3c_get_bounds(hmat, h_inv, vol, is_ortho, G_min, R_min, la_max, lb_max, lc_max, &
243 : zeta, zetb, zetc, a_mm, sum_precision, n_sum_1d, n_sum_3d, &
244 : G_bounds_1, G_rads_1, R_bounds_2, R_rads_2, R_bounds_3, R_rads_3)
245 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv
246 : REAL(KIND=dp), INTENT(IN) :: vol
247 : LOGICAL, INTENT(IN) :: is_ortho
248 : REAL(KIND=dp), INTENT(IN) :: G_min, R_min
249 : INTEGER, INTENT(IN) :: la_max, lb_max, lc_max
250 : REAL(KIND=dp), INTENT(IN) :: zeta, zetb, zetc, a_mm, sum_precision
251 : INTEGER(KIND=int_8), DIMENSION(3, 3), INTENT(OUT) :: n_sum_1d
252 : INTEGER(KIND=int_8), DIMENSION(3), INTENT(OUT) :: n_sum_3d
253 : REAL(KIND=dp), DIMENSION(3, 3) :: G_bounds_1
254 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: G_rads_1
255 : REAL(KIND=dp), DIMENSION(3, 3) :: R_bounds_2
256 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: R_rads_2
257 : REAL(KIND=dp), DIMENSION(2, 3) :: R_bounds_3
258 : REAL(KIND=dp), DIMENSION(2), INTENT(OUT) :: R_rads_3
259 :
260 : INTEGER :: i, i_xyz, order_1, order_2
261 : REAL(KIND=dp) :: ns1_G1, ns1_G2, ns2_G, ns2_R, ns3_R1, &
262 : ns3_R2, vol_G
263 : CALL eri_mme_3c_get_rads(la_max, lb_max, lc_max, zeta, zetb, zetc, a_mm, G_min, R_min, sum_precision, &
264 2256540 : G_rads_1=G_rads_1, R_rads_2=R_rads_2, R_rads_3=R_rads_3)
265 :
266 2256540 : vol_G = twopi**3/vol
267 :
268 9026160 : order_1 = MAXLOC(G_rads_1, DIM=1)
269 9026160 : order_2 = MINLOC(G_rads_1, DIM=1)
270 :
271 9026160 : DO i = 1, 3
272 110570460 : G_bounds_1(i, :) = ellipsoid_bounds(G_rads_1(i), TRANSPOSE(hmat)/(2.0_dp*pi))
273 : END DO
274 :
275 9026160 : DO i = 1, 3
276 29335020 : R_bounds_2(i, :) = ellipsoid_bounds(R_rads_2(i), h_inv)
277 : END DO
278 :
279 6769620 : DO i = 1, 2
280 20308860 : R_bounds_3(i, :) = ellipsoid_bounds(R_rads_3(i), h_inv)
281 : END DO
282 :
283 2256540 : IF (is_ortho) THEN
284 8432160 : DO i_xyz = 1, 3
285 :
286 6324120 : ns3_R1 = 2.0_dp*R_bounds_3(1, i_xyz)
287 6324120 : ns3_R2 = 2.0_dp*R_bounds_3(2, i_xyz)
288 :
289 8432160 : n_sum_1d(3, i_xyz) = nsum_3c_rspace_1d(ns3_R1, ns3_R2)
290 : END DO
291 :
292 : ELSE
293 :
294 594000 : order_1 = MAXLOC(G_rads_1, DIM=1)
295 594000 : order_2 = MINLOC(G_rads_1, DIM=1)
296 :
297 71062 : SELECT CASE (order_1)
298 : CASE (1)
299 71062 : ns1_G1 = 4._dp/3._dp*pi*G_rads_1(2)**3/vol_G
300 71062 : ns1_G2 = 4._dp/3._dp*pi*G_rads_1(3)**3/vol_G
301 : CASE (2)
302 53482 : ns1_G1 = 4._dp/3._dp*pi*G_rads_1(1)**3/vol_G
303 53482 : ns1_G2 = 4._dp/3._dp*pi*G_rads_1(3)**3/vol_G
304 : CASE (3)
305 23956 : ns1_G1 = 4._dp/3._dp*pi*G_rads_1(1)**3/vol_G
306 148500 : ns1_G2 = 4._dp/3._dp*pi*G_rads_1(2)**3/vol_G
307 : END SELECT
308 :
309 148500 : n_sum_3d(1) = nsum_3c_gspace_3d(ns1_G1, ns1_G2, la_max, lb_max, lc_max)
310 :
311 148500 : ns2_G = 4._dp/3._dp*pi*G_rads_1(order_2)**3/vol_G
312 148500 : ns2_R = 4._dp/3._dp*pi*R_rads_2(order_2)**3/vol
313 :
314 148500 : ns3_R1 = 4._dp/3._dp*pi*R_rads_3(1)**3/vol
315 148500 : ns3_R2 = 4._dp/3._dp*pi*R_rads_3(2)**3/vol
316 :
317 30492 : SELECT CASE (order_2)
318 : CASE (1)
319 30492 : n_sum_3d(2) = nsum_product_3c_gspace_3d(ns2_G, ns2_R, la_max, lb_max, lc_max)
320 : CASE (2)
321 24372 : n_sum_3d(2) = nsum_product_3c_gspace_3d(ns2_G, ns2_R, lb_max, la_max, lc_max)
322 : CASE (3)
323 148500 : n_sum_3d(2) = nsum_product_3c_gspace_3d(ns2_G, ns2_R, lc_max, lb_max, la_max)
324 : END SELECT
325 :
326 148500 : n_sum_3d(3) = nsum_3c_rspace_3d(ns3_R1, ns3_R2, la_max, lb_max, lc_max)
327 : END IF
328 :
329 2256540 : END SUBROUTINE eri_mme_3c_get_bounds
330 :
331 : ! **************************************************************************************************
332 : !> \brief Roughly estimated number of floating point operations
333 : !> \param ns_G ...
334 : !> \param l ...
335 : !> \param m ...
336 : !> \return ...
337 : ! **************************************************************************************************
338 2236488 : PURE FUNCTION nsum_2c_gspace_1d(ns_G, l, m)
339 : REAL(KIND=dp), INTENT(IN) :: ns_G
340 : INTEGER, INTENT(IN) :: l, m
341 : INTEGER(KIND=int_8) :: nsum_2c_gspace_1d
342 :
343 2236488 : nsum_2c_gspace_1d = NINT(ns_G*(2*exp_w + (l + m + 1)*5), KIND=int_8)
344 2236488 : END FUNCTION nsum_2c_gspace_1d
345 :
346 : ! **************************************************************************************************
347 : !> \brief Compute Ewald-like sum for 2-center ERIs in G space in 1 dimension
348 : !> S_G(l, alpha) = (-i)^l*inv_lgth*sum_G( C(l, alpha, G) exp(iGR) ), with
349 : !> C(l, alpha, r) = r^l exp(-alpha*r^2),
350 : !> dG = inv_lgth*twopi and G = -G_bound*dG, (-G_bound + 1)*dG, ..., G_bound*dG
351 : !> for all l < = l_max.
352 : !> \param S_G ...
353 : !> \param R ...
354 : !> \param alpha ...
355 : !> \param inv_lgth ...
356 : !> \param G_c ...
357 : !> \note S_G is real.
358 : ! **************************************************************************************************
359 202238 : PURE SUBROUTINE pgf_sum_2c_gspace_1d(S_G, R, alpha, inv_lgth, G_c)
360 : REAL(KIND=dp), DIMENSION(0:), INTENT(OUT) :: S_G
361 : REAL(KIND=dp), INTENT(IN) :: R, alpha, inv_lgth, G_c
362 :
363 : COMPLEX(KIND=dp) :: exp_tot
364 202238 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: S_G_c
365 : INTEGER :: gg, l, l_max
366 : REAL(KIND=dp) :: dG, G, G_pow_l
367 :
368 202238 : dG = inv_lgth*twopi
369 202238 : l_max = UBOUND(S_G, 1)
370 :
371 606714 : ALLOCATE (S_G_c(0:l_max))
372 202238 : S_G_c(:) = 0.0_dp
373 1079898 : DO gg = -FLOOR(G_c), FLOOR(G_c)
374 877660 : G = gg*dG
375 877660 : exp_tot = EXP(-alpha*G**2)*EXP(gaussi*G*R) ! cost: 2 exp_w flops
376 877660 : G_pow_l = 1.0_dp
377 4490891 : DO l = 0, l_max
378 3410993 : S_G_c(l) = S_G_c(l) + G_pow_l*(-1.0_dp)**l*exp_tot ! cost: 4 flops
379 4288653 : G_pow_l = G_pow_l*G ! cost: 1 flop
380 : END DO
381 : END DO
382 :
383 1848744 : S_G(:) = REAL(S_G_c(0:l_max)*i_pow([(l, l=0, l_max)]))*inv_lgth
384 202238 : END SUBROUTINE pgf_sum_2c_gspace_1d
385 :
386 : ! **************************************************************************************************
387 : !> \brief Roughly estimated number of floating point operations
388 : !> \param ns_G ...
389 : !> \param l ...
390 : !> \param m ...
391 : !> \return ...
392 : ! **************************************************************************************************
393 146608 : PURE FUNCTION nsum_2c_gspace_3d(ns_G, l, m)
394 : REAL(KIND=dp), INTENT(IN) :: ns_G
395 : INTEGER, INTENT(IN) :: l, m
396 : INTEGER(KIND=int_8) :: nsum_2c_gspace_3d
397 :
398 146608 : nsum_2c_gspace_3d = NINT(ns_G*(2*exp_w + ncoset(l + m)*7), KIND=int_8)
399 146608 : END FUNCTION nsum_2c_gspace_3d
400 :
401 : ! **************************************************************************************************
402 : !> \brief As pgf_sum_2c_gspace_1d but 3d sum required for non-orthorhombic cells
403 : !> \param S_G ...
404 : !> \param l_max ...
405 : !> \param R ...
406 : !> \param alpha ...
407 : !> \param h_inv ...
408 : !> \param G_c ...
409 : !> \param G_rad ...
410 : !> \param vol ...
411 : !> \param coulomb ...
412 : !> \note MMME Method is not very efficient for non-orthorhombic cells
413 : ! **************************************************************************************************
414 27581 : PURE SUBROUTINE pgf_sum_2c_gspace_3d(S_G, l_max, R, alpha, h_inv, G_c, G_rad, vol, potential, pot_par)
415 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: S_G
416 : INTEGER, INTENT(IN) :: l_max
417 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R
418 : REAL(KIND=dp), INTENT(IN) :: alpha
419 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: h_inv
420 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G_c
421 : REAL(KIND=dp), INTENT(IN) :: G_rad, vol
422 : INTEGER, INTENT(IN), OPTIONAL :: potential
423 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par
424 :
425 : COMPLEX(KIND=dp) :: exp_tot
426 : INTEGER, DIMENSION(3) :: l_xyz
427 : INTEGER :: gx, gy, gz, k, l, lco, lx, ly, lz
428 55162 : COMPLEX(KIND=dp), DIMENSION(ncoset(l_max)) :: Ig
429 : REAL(KIND=dp) :: G_rads_sq, G_sq
430 : REAL(KIND=dp), DIMENSION(3) :: G, G_x, G_y, G_z
431 27581 : REAL(KIND=dp), DIMENSION(3, 0:l_max) :: G_pow_l
432 : REAL(KIND=dp), DIMENSION(3, 3) :: ht
433 :
434 358553 : ht = twopi*TRANSPOSE(h_inv)
435 660268 : Ig(:) = 0.0_dp
436 :
437 27581 : G_rads_sq = G_rad**2
438 :
439 231084 : DO gx = -FLOOR(G_c(1)), FLOOR(G_c(1))
440 814012 : G_x = ht(:, 1)*gx
441 1913121 : DO gy = -FLOOR(G_c(2)), FLOOR(G_c(2))
442 6728148 : G_y = ht(:, 2)*gy
443 18506843 : DO gz = -FLOOR(G_c(3)), FLOOR(G_c(3))
444 66485212 : G_z = ht(:, 3)*gz
445 :
446 66485212 : G = G_x + G_y + G_z
447 16621303 : G_sq = G(1)**2 + G(2)**2 + G(3)**2
448 16621303 : IF (G_sq > G_rads_sq) CYCLE
449 :
450 6396767 : IF (PRESENT(potential)) THEN
451 1033920 : IF (gx == 0 .AND. gy == 0 .AND. gz == 0) CYCLE
452 : END IF
453 :
454 25586300 : exp_tot = EXP(-alpha*G_sq)*EXP(gaussi*DOT_PRODUCT(G, R)) ! cost: 2 exp_w flops
455 6396575 : IF (PRESENT(potential)) THEN
456 1378304 : SELECT CASE (potential)
457 : CASE (eri_mme_coulomb)
458 344576 : exp_tot = exp_tot/G_sq
459 : CASE (eri_mme_yukawa)
460 344576 : exp_tot = exp_tot/(G_sq + pot_par**2)
461 : !exp_tot = exp_tot/G_sq
462 : CASE (eri_mme_longrange)
463 1033728 : exp_tot = exp_tot/G_sq*EXP(-G_sq/pot_par**2)
464 : END SELECT
465 : END IF
466 25586300 : DO k = 1, 3
467 19189725 : G_pow_l(k, 0) = 1.0_dp
468 110227076 : DO l = 1, l_max
469 103830501 : G_pow_l(k, l) = G_pow_l(k, l - 1)*G(k)
470 : END DO
471 : END DO
472 466764837 : DO lco = 1, ncoset(l_max)
473 458686225 : CALL get_l(lco, l, lx, ly, lz)
474 1834744900 : l_xyz = [lx, ly, lz]
475 : Ig(coset(lx, ly, lz)) = Ig(coset(lx, ly, lz)) + & ! cost: 7 flops
476 : G_pow_l(1, lx)*G_pow_l(2, ly)*G_pow_l(3, lz)* &
477 475307528 : exp_tot*(-1.0_dp)**l*i_pow(l)
478 : END DO
479 : END DO
480 : END DO
481 : END DO
482 660268 : S_G(:) = REAL(Ig(:), KIND=dp)/vol
483 27581 : END SUBROUTINE pgf_sum_2c_gspace_3d
484 :
485 : ! **************************************************************************************************
486 : !> \brief Roughly estimated number of floating point operations
487 : !> \param ns_R ...
488 : !> \param l ...
489 : !> \param m ...
490 : !> \return ...
491 : ! **************************************************************************************************
492 2236488 : PURE FUNCTION nsum_2c_rspace_1d(ns_R, l, m)
493 : REAL(KIND=dp), INTENT(IN) :: ns_R
494 : INTEGER, INTENT(IN) :: l, m
495 : INTEGER(KIND=int_8) :: nsum_2c_rspace_1d
496 :
497 2236488 : nsum_2c_rspace_1d = NINT(ns_R*(exp_w + (l + m + 1)*3), KIND=int_8)
498 2236488 : END FUNCTION nsum_2c_rspace_1d
499 :
500 : ! **************************************************************************************************
501 : !> \brief Compute Ewald-like sum for 2-center ERIs in R space in 1 dimension
502 : !> S_R(l, alpha) = SQRT(alpha/pi) sum_R'( H(l, alpha, R-R') ),
503 : !> with H(l, alpha, R) = (-d/dR)^l exp(-alpha*R^2),
504 : !> dR = lgth and R' = -R_min*dR, (-R_min + 1)*dR, ..., R_max*dR,
505 : !> for all l < = l_max.
506 : !> \param S_R ...
507 : !> \param R ...
508 : !> \param alpha ...
509 : !> \param lgth ...
510 : !> \param R_c ...
511 : !> \note result is equivalent to pgf_sum_2c_gspace_1d with
512 : !> S_R(l, alpha) = S_G(l, 1/(4*alpha))
513 : ! **************************************************************************************************
514 2212684 : PURE SUBROUTINE pgf_sum_2c_rspace_1d(S_R, R, alpha, lgth, R_c)
515 : REAL(KIND=dp), DIMENSION(0:), INTENT(OUT) :: S_R
516 : REAL(KIND=dp), INTENT(IN) :: R, alpha, lgth, R_c
517 :
518 : INTEGER :: l, l_max, rr
519 : REAL(KIND=dp) :: dR, exp_tot, R_pow_l, Rp
520 2212684 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: h_to_c
521 :
522 2212684 : dR = lgth
523 2212684 : l_max = UBOUND(S_R, 1)
524 :
525 : ! 1) compute sum over C(l, alpha, R - R') instead of H(l, alpha, R - R')
526 9453605 : S_R(:) = 0.0_dp
527 2212684 : Rp = R - R_c*dR
528 9268511 : DO rr = CEILING(-R_c - R/dR), FLOOR(R_c - R/dR)
529 7055827 : Rp = R + rr*dR
530 7055827 : exp_tot = EXP(-alpha*Rp**2) ! cost: exp_w flops
531 7055827 : R_pow_l = 1.0_dp
532 32695273 : DO l = 0, l_max
533 23426762 : S_R(l) = S_R(l) + R_pow_l*exp_tot ! cost: 2 flops
534 30482589 : R_pow_l = R_pow_l*Rp ! cost: 1 flop
535 : END DO
536 : END DO
537 :
538 : ! 2) C --> H
539 2212684 : CALL create_hermite_to_cartesian(alpha, l_max, h_to_c)
540 9453605 : S_R = MATMUL(TRANSPOSE(h_to_c(0:l_max, 0:l_max)), S_R)*SQRT(alpha/pi)
541 2212684 : END SUBROUTINE pgf_sum_2c_rspace_1d
542 :
543 : ! **************************************************************************************************
544 : !> \brief Roughly estimated number of floating point operations
545 : !> \param ns_R ...
546 : !> \param l ...
547 : !> \param m ...
548 : !> \return ...
549 : ! **************************************************************************************************
550 146608 : PURE FUNCTION nsum_2c_rspace_3d(ns_R, l, m)
551 : REAL(KIND=dp), INTENT(IN) :: ns_R
552 : INTEGER, INTENT(IN) :: l, m
553 : INTEGER(KIND=int_8) :: nsum_2c_rspace_3d
554 :
555 146608 : nsum_2c_rspace_3d = NINT(ns_R*(exp_w + ncoset(l + m)*(4 + ncoset(l + m)*4)), KIND=int_8)
556 146608 : END FUNCTION nsum_2c_rspace_3d
557 :
558 : ! **************************************************************************************************
559 : !> \brief As pgf_sum_2c_rspace_1d but 3d sum required for non-orthorhombic cells
560 : !> \param S_R ...
561 : !> \param l_max ...
562 : !> \param R ...
563 : !> \param alpha ...
564 : !> \param hmat ...
565 : !> \param h_inv ...
566 : !> \param R_c ...
567 : !> \param R_rad ...
568 : !> \note MMME Method is not very efficient for non-orthorhombic cells
569 : ! **************************************************************************************************
570 128219 : PURE SUBROUTINE pgf_sum_2c_rspace_3d(S_R, l_max, R, alpha, hmat, h_inv, R_c, R_rad)
571 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: S_R
572 : INTEGER, INTENT(IN) :: l_max
573 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R
574 : REAL(KIND=dp), INTENT(IN) :: alpha
575 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv
576 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R_c
577 : REAL(KIND=dp), INTENT(IN) :: R_rad
578 :
579 : INTEGER :: k, l, lco, llx, lly, llz, lx, ly, lz, &
580 : sx, sy, sz
581 : REAL(KIND=dp) :: exp_tot, R_rad_sq, R_sq
582 128219 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: h_to_c
583 : REAL(KIND=dp), DIMENSION(3) :: R_l, R_r, Rp, Rx, Ry, Rz, s_shift
584 256438 : REAL(KIND=dp), DIMENSION(3, 0:l_max) :: R_pow_l
585 128219 : REAL(KIND=dp), DIMENSION(ncoset(l_max)) :: S_R_C
586 :
587 2057776 : S_R(:) = 0.0_dp
588 2057776 : S_R_C(:) = 0.0_dp
589 :
590 2051504 : s_shift = MATMUL(h_inv, -R)
591 512876 : R_l = -R_c + s_shift
592 512876 : R_r = R_c + s_shift
593 :
594 128219 : R_rad_sq = R_rad**2
595 :
596 432234 : DO sx = CEILING(R_l(1)), FLOOR(R_r(1))
597 1216060 : Rx = hmat(:, 1)*sx
598 1470931 : DO sy = CEILING(R_l(2)), FLOOR(R_r(2))
599 4154788 : Ry = hmat(:, 2)*sy
600 5870059 : DO sz = CEILING(R_l(3)), FLOOR(R_r(3))
601 18109388 : Rz = hmat(:, 3)*sz
602 18109388 : Rp = Rx + Ry + Rz
603 4527347 : R_sq = (Rp(1) + R(1))**2 + (Rp(2) + R(2))**2 + (Rp(3) + R(3))**2
604 4527347 : IF (R_sq > R_rad_sq) CYCLE
605 1814156 : exp_tot = EXP(-alpha*R_sq) ! cost: exp_w flops
606 7256624 : DO k = 1, 3
607 5442468 : R_pow_l(k, 0) = 1.0_dp
608 16096739 : DO l = 1, l_max
609 14282583 : R_pow_l(k, l) = R_pow_l(k, l - 1)*(Rp(k) + R(k))
610 : END DO
611 : END DO
612 22199618 : DO lco = 1, ncoset(l_max)
613 19346765 : CALL get_l(lco, l, lx, ly, lz)
614 : ! cost: 4 flops
615 23874112 : S_R_C(coset(lx, ly, lz)) = S_R_C(coset(lx, ly, lz)) + R_pow_l(1, lx)*R_pow_l(2, ly)*R_pow_l(3, lz)*exp_tot
616 : END DO
617 : END DO
618 : END DO
619 : END DO
620 :
621 128219 : CALL create_hermite_to_cartesian(alpha, l_max, h_to_c)
622 :
623 2057776 : DO lco = 1, ncoset(l_max)
624 1929557 : CALL get_l(lco, l, lx, ly, lz)
625 :
626 5692372 : DO llx = 0, lx
627 11963703 : DO lly = 0, ly
628 20729884 : DO llz = 0, lz
629 : S_R(coset(lx, ly, lz)) = S_R(coset(lx, ly, lz)) + & ! cost: 4 flops
630 : h_to_c(llx, lx)*h_to_c(lly, ly)*h_to_c(llz, lz)* &
631 17095288 : S_R_C(coset(llx, lly, llz))
632 : END DO
633 : END DO
634 : END DO
635 : END DO
636 2057776 : S_R(:) = S_R(:)*(alpha/pi)**1.5_dp
637 :
638 128219 : END SUBROUTINE pgf_sum_2c_rspace_3d
639 :
640 : ! **************************************************************************************************
641 : !> \brief Compute 1d sum
642 : !> S_G(l, alpha) = inv_lgth*sum_G( C(l, alpha, delta_l, G) ) with
643 : !> C(l, alpha, delta_l, G) = prefactor*|G|^(l-delta_l) exp(-alpha*G^2)
644 : !> if G not equal 0
645 : !> C(l = 0, alpha, delta_l, 0) = 1, C(l>0, alpha, delta_l, 0) = 0
646 : !> dG = inv_lgth*twopi and G = -G_bound*dG, (-G_bound + 1)*dG, ..., G_bound*dG
647 : !> for all l < = l_max.
648 : !> \param S_G ...
649 : !> \param alpha ...
650 : !> \param inv_lgth ...
651 : !> \param G_min ...
652 : !> \param G_c ...
653 : !> \param delta_l ...
654 : !> \param prefactor ...
655 : !> \note needed for cutoff error estimate
656 : ! **************************************************************************************************
657 295044 : PURE SUBROUTINE pgf_sum_2c_gspace_1d_deltal(S_G, alpha, inv_lgth, G_min, G_c, delta_l, prefactor)
658 : REAL(KIND=dp), DIMENSION(0:), INTENT(OUT) :: S_G
659 : REAL(KIND=dp), INTENT(IN) :: alpha, inv_lgth
660 : INTEGER, INTENT(IN) :: G_min, G_c
661 : REAL(KIND=dp), INTENT(IN) :: delta_l, prefactor
662 :
663 : INTEGER :: k, k0, k1, l, l_max
664 : REAL(KIND=dp) :: dG, exp_tot, G, prefac
665 :
666 295044 : prefac = prefactor*inv_lgth
667 295044 : dG = inv_lgth*twopi ! positive
668 295044 : l_max = UBOUND(S_G, 1)
669 :
670 1660104 : S_G(:) = 0.0_dp
671 295044 : IF (0 < G_min) THEN
672 : k0 = G_min; k1 = 0
673 73761 : ELSE IF (G_c < 0) THEN
674 : k0 = 0; k1 = G_c
675 : ELSE ! Gmin <= 0 /\ 0 <= Gc
676 73761 : S_G(0) = prefac
677 : k0 = 1; k1 = -1
678 : END IF
679 : ! positive range
680 : IF (0 < k0) THEN
681 14185303 : DO k = k0, G_c
682 13890259 : G = k*dG; exp_tot = EXP(-alpha*G**2)*prefac
683 83043248 : DO l = 0, l_max
684 82748204 : S_G(l) = S_G(l) + G**(l - delta_l)*exp_tot
685 : END DO
686 : END DO
687 : END IF
688 : ! negative range
689 295044 : IF (k1 < 0) THEN
690 4816941 : DO k = G_min, k1
691 4743180 : G = k*dG; exp_tot = EXP(-alpha*G**2)*prefac
692 27823869 : DO l = 0, l_max
693 27750108 : S_G(l) = S_G(l) + ABS(G)**(l - delta_l)*exp_tot
694 : END DO
695 : END DO
696 : END IF
697 295044 : END SUBROUTINE pgf_sum_2c_gspace_1d_deltal
698 :
699 : ! **************************************************************************************************
700 : !> \brief Compute Ewald-like sum for 3-center integrals in 1 dimension
701 : !> S_G(l, m, n, alpha, beta, gamma) = i^(l+m+n)*(-1)^(l+m)*inv_lgth^2*
702 : !> sum_G sum_G'( exp(i G R1)
703 : !> C(l,alpha,G) C(m,beta,G'-G) C(n,gamma,G') exp(i G' R2) )
704 : !> for all l < = l_max, m <= m_max, n <= n_max.
705 : !> a_mm is the minimax exponent.
706 : !> alpha = 1/(4 zeta), beta = 1/(4 zetb), gamma = 1/(4 zetc) + a_mm
707 : !> R1 = RB-RA; R2 = RC-RB
708 : !> Note on method / order arguments:
709 : !> Three equivalent methods (Poisson summation) to compute this sum over
710 : !> Cartesian Gaussians C or Hermite Gaussians H and
711 : !> reciprocal lattice vectors G or direct lattice vectors R:
712 : !> - method 1: sum_G sum_G' C(G) C(G,G') C(G')
713 : !> - method 2: sum_G sum_R C(G) C(R)
714 : !> - method 3: sum_R sum_R' H(R, R')
715 : !> The order parameter selects the Gaussian functions over which the sum is performed
716 : !> method 1: order = 1, 2, 3
717 : !> method 2: order = 1, 2, 3
718 : !> method 3: order = 1
719 : !> If method and order are not present, the method / order that converges fastest is
720 : !> automatically chosen.
721 : !> \param S_G ...
722 : !> \param RA ...
723 : !> \param RB ...
724 : !> \param RC ...
725 : !> \param zeta ...
726 : !> \param zetb ...
727 : !> \param zetc ...
728 : !> \param a_mm ...
729 : !> \param lgth ...
730 : !> \param G_bounds_1 ...
731 : !> \param R_bounds_2 ...
732 : !> \param R_bounds_3 ...
733 : !> \param method ...
734 : !> \param method_out ...
735 : !> \param order ...
736 : ! **************************************************************************************************
737 6324120 : SUBROUTINE pgf_sum_3c_1d(S_G, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_bounds_3)
738 : REAL(KIND=dp), DIMENSION(0:, 0:, 0:), INTENT(OUT) :: S_G
739 : REAL(KIND=dp), INTENT(IN) :: RA, RB, RC, zeta, zetb, zetc, a_mm, lgth
740 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: R_bounds_3
741 :
742 : INTEGER :: l_max, m_max, n_max
743 : INTEGER :: nR1, nR2
744 : INTEGER :: prop_exp
745 :
746 6324120 : l_max = UBOUND(S_G, 1)
747 6324120 : m_max = UBOUND(S_G, 2)
748 6324120 : n_max = UBOUND(S_G, 3)
749 :
750 6324120 : nR1 = 2*FLOOR(R_bounds_3(1)) + 1
751 6324120 : nR2 = 2*FLOOR(R_bounds_3(2)) + 1
752 :
753 6324120 : IF (nR1*nR2 > 1 + nR1*2) THEN
754 : prop_exp = 1
755 : ELSE
756 2997690 : prop_exp = 0
757 : END IF
758 :
759 25296480 : IF (MAXVAL([l_max, m_max, n_max]) > ${lmax_unroll}$) THEN
760 92160 : CALL pgf_sum_3c_rspace_1d_generic(S_G, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_bounds_3)
761 : ELSE
762 : #:for l_max in range(0,lmax_unroll+1)
763 6231960 : IF (l_max == ${l_max}$) THEN
764 : #:for m_max in range(0,lmax_unroll+1)
765 6231960 : IF (m_max == ${m_max}$) THEN
766 : #:for n_max in range(0, lmax_unroll+1)
767 6231960 : IF (n_max == ${n_max}$) THEN
768 : #:for prop_exp in range(0,2)
769 6231960 : IF (prop_exp == ${prop_exp}$) THEN
770 : CALL pgf_sum_3c_rspace_1d_${l_max}$_${m_max}$_${n_max}$_exp_${prop_exp}$ (S_G, RA, RB, RC, &
771 6231960 : zeta, zetb, zetc, a_mm, lgth, R_bounds_3)
772 6231960 : RETURN
773 : END IF
774 : #:endfor
775 : END IF
776 : #:endfor
777 : END IF
778 : #:endfor
779 : END IF
780 : #:endfor
781 : END IF
782 :
783 : END SUBROUTINE pgf_sum_3c_1d
784 :
785 : ! **************************************************************************************************
786 : !> \brief Roughly estimated number of floating point operations
787 : !> \param ns_G1 ...
788 : !> \param ns_G2 ...
789 : !> \param l ...
790 : !> \param m ...
791 : !> \param n ...
792 : !> \return ...
793 : ! **************************************************************************************************
794 0 : PURE FUNCTION nsum_3c_gspace_1d()
795 : INTEGER(KIND=int_8) :: nsum_3c_gspace_1d
796 :
797 0 : nsum_3c_gspace_1d = 15
798 0 : END FUNCTION nsum_3c_gspace_1d
799 :
800 : ! **************************************************************************************************
801 : !> \brief Roughly estimated number of floating point operations
802 : !> \param ns_G ...
803 : !> \param ns_R ...
804 : !> \param l ...
805 : !> \param m ...
806 : !> \param n ...
807 : !> \return ...
808 : ! **************************************************************************************************
809 0 : PURE FUNCTION nsum_product_3c_gspace_1d(ns_G, ns_R)
810 : REAL(KIND=dp), INTENT(IN) :: ns_G, ns_R
811 : INTEGER(KIND=int_8) :: nsum_product_3c_gspace_1d
812 :
813 0 : nsum_product_3c_gspace_1d = MIN(19, NINT(ns_G*(3 + ns_R*2)))
814 0 : END FUNCTION nsum_product_3c_gspace_1d
815 :
816 : ! **************************************************************************************************
817 : !> \brief Roughly estimated number of floating point operations
818 : !> \param ns_R1 ...
819 : !> \param ns_R2 ...
820 : !> \param l ...
821 : !> \param m ...
822 : !> \param n ...
823 : !> \return ...
824 : ! **************************************************************************************************
825 6324120 : PURE FUNCTION nsum_3c_rspace_1d(ns_R1, ns_R2)
826 : REAL(KIND=dp), INTENT(IN) :: ns_R1, ns_R2
827 : INTEGER(KIND=int_8) :: nsum_3c_rspace_1d
828 :
829 6324120 : nsum_3c_rspace_1d = NINT(MIN((4 + ns_R1*2), ns_R1*(ns_R2 + 1)), KIND=int_8)
830 6324120 : END FUNCTION nsum_3c_rspace_1d
831 :
832 : ! **************************************************************************************************
833 : !> \brief Helper routine: compute SQRT(alpha/pi) (-1)^n sum_(R, R') sum_{t=0}^{l+m} E(t,l,m) H(RC - P(R) - R', t + n, alpha)
834 : !> with alpha = 1.0_dp/((a + b + c)/((a + b)*c) + 4.0_dp*a_mm),
835 : !> P(R) = (a*(RA + R) + b*RB)/(a + b)
836 : !> \param S_R ...
837 : !> \param RA ...
838 : !> \param RB ...
839 : !> \param RC ...
840 : !> \param zeta ...
841 : !> \param zetb ...
842 : !> \param zetc ...
843 : !> \param a_mm ...
844 : !> \param lgth ...
845 : !> \param R_c ...
846 : ! **************************************************************************************************
847 92160 : PURE SUBROUTINE pgf_sum_3c_rspace_1d_generic(S_R, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_c)
848 : REAL(KIND=dp), DIMENSION(0:, 0:, 0:), INTENT(OUT) :: S_R
849 : REAL(KIND=dp), INTENT(IN) :: RA, RB, RC, zeta, zetb, zetc, a_mm, lgth
850 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: R_c
851 :
852 : INTEGER :: ll, mm, l, k, l_max, m, m_max, n, n_max, rr1, rr2, t, l_tot_max
853 : REAL(KIND=dp) :: alpha, dR, exp_tot, R, R1, R2, R_offset, &
854 : R_pow_t, R_tmp, rr1_delta, rr2_delta, c1, c2, c3
855 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: S_R_t
856 92160 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: h_to_c
857 92160 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: E
858 :
859 92160 : dR = lgth
860 92160 : alpha = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm)
861 92160 : l_max = UBOUND(S_R, 1)
862 92160 : m_max = UBOUND(S_R, 2)
863 92160 : n_max = UBOUND(S_R, 3)
864 92160 : l_tot_max = l_max + m_max + n_max
865 :
866 276480 : ALLOCATE (S_R_t(0:l_max + m_max + n_max))
867 460800 : ALLOCATE (E(-1:l_max + m_max + 1, -1:l_max, -1:m_max))
868 :
869 2903040 : S_R(:, :, :) = 0.0_dp
870 :
871 92160 : R_offset = RC - (zeta*RA + zetb*RB)/(zeta + zetb)
872 :
873 : ! inline CALL create_hermite_to_cartesian(alpha, l_tot_max, h_to_c)
874 368640 : ALLOCATE (h_to_c(-1:l_tot_max + 1, 0:l_tot_max))
875 92160 : h_to_c(:, :) = 0.0_dp
876 92160 : h_to_c(0, 0) = 1.0_dp
877 737280 : DO l = 0, l_tot_max - 1
878 3962880 : DO k = 0, l + 1
879 3870720 : h_to_c(k, l + 1) = -(k + 1)*h_to_c(k + 1, l) + 2.0_dp*alpha*h_to_c(k - 1, l)
880 : END DO
881 : END DO
882 :
883 92160 : rr1_delta = (RA - RB)/dR
884 1428480 : DO rr1 = CEILING(-R_c(1) + rr1_delta), FLOOR(R_c(1) + rr1_delta)
885 1336320 : S_R_t(:) = 0.0_dp
886 1336320 : R1 = rr1*dR
887 1336320 : R_tmp = R_offset + R1*zeta/(zeta + zetb)
888 1336320 : rr2_delta = -R_tmp/dR
889 15398400 : DO rr2 = CEILING(-R_c(2) + rr2_delta), FLOOR(R_c(2) + rr2_delta)
890 14062080 : R2 = rr2*dR
891 14062080 : R = R_tmp + R2
892 14062080 : exp_tot = EXP(-alpha*R**2) ! cost: exp_w flops
893 14062080 : R_pow_t = 1.0_dp
894 127895040 : DO t = 0, l_max + m_max + n_max
895 112496640 : S_R_t(t) = S_R_t(t) + R_pow_t*exp_tot ! cost: 2 flops
896 126558720 : R_pow_t = R_pow_t*R ! cost: 1 flop
897 : END DO
898 : END DO
899 :
900 : ! C --> H
901 12026880 : S_R_t(:) = MATMUL(TRANSPOSE(h_to_c(0:l_max + m_max + n_max, 0:l_max + m_max + n_max)), S_R_t)*SQRT(alpha/pi)
902 :
903 : ! H --> HH
904 : !inline CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, RA-R1, RB, 2, E)
905 :
906 1336320 : E(:, :, :) = 0.0_dp
907 1336320 : E(0, 0, 0) = EXP(-zeta*zetb/(zeta + zetb)*(RA - R1 - RB)**2)
908 :
909 1336320 : c1 = 1.0_dp/(zeta + zetb)
910 1336320 : c2 = 2.0_dp*(zetb/(zeta + zetb))*(RB - (RA - R1))
911 1336320 : c3 = 2.0_dp*(zeta/(zeta + zetb))*(RA - R1 - RB)
912 :
913 7948800 : DO mm = 0, m_max - 1
914 8179200 : DO ll = 0, l_max - 1
915 7879680 : DO t = 0, ll + mm + 1
916 : E(t, ll + 1, mm) = zeta*(c1*E(t - 1, ll, mm) + &
917 : c2*E(t, ll, mm) + &
918 : 2*(t + 1)*E(t + 1, ll, mm) - &
919 1036800 : 2*ll*E(t, ll - 1, mm))
920 : E(t, ll, mm + 1) = zetb*(c1*E(t - 1, ll, mm) + &
921 : c3*E(t, ll, mm) + &
922 : 2*(t + 1)*E(t + 1, ll, mm) - &
923 1267200 : 2*mm*E(t, ll, mm - 1))
924 : END DO
925 : END DO
926 : END DO
927 :
928 1451520 : DO ll = 0, l_max - 1
929 2142720 : DO t = 0, ll + m_max + 1
930 : E(t, ll + 1, m_max) = zeta*(c1*E(t - 1, ll, m_max) + &
931 : c2*E(t, ll, m_max) + &
932 : 2*(t + 1)*E(t + 1, ll, m_max) - &
933 806400 : 2*ll*E(t, ll - 1, m_max))
934 : END DO
935 : END DO
936 :
937 7948800 : DO mm = 0, m_max - 1
938 34560000 : DO t = 0, l_max + mm + 1
939 : E(t, l_max, mm + 1) = zetb*(c1*E(t - 1, l_max, mm) + &
940 : c3*E(t, l_max, mm) + &
941 : 2*(t + 1)*E(t + 1, l_max, mm) - &
942 33223680 : 2*mm*E(t, l_max, mm - 1))
943 : END DO
944 : END DO
945 :
946 5391360 : DO n = 0, n_max
947 29007360 : DO m = 0, m_max
948 51724800 : DO l = 0, l_max
949 132364800 : DO t = 0, l + m
950 108656640 : S_R(l, m, n) = S_R(l, m, n) + E(t, l, m)*(-1)**n*S_R_t(t + n) ! cost: 5 flops
951 : END DO
952 : END DO
953 : END DO
954 : END DO
955 : END DO
956 :
957 2903040 : S_R = S_R*pi**(-0.5_dp)*((zeta + zetb)/(zeta*zetb))**(-0.5_dp)
958 92160 : END SUBROUTINE pgf_sum_3c_rspace_1d_generic
959 :
960 : ! **************************************************************************************************
961 : !> \brief Helper routine: compute SQRT(alpha/pi) (-1)^n sum_(R, R') sum_{t=0}^{l+m} E(t,l,m) H(RC - P(R) - R', t + n, alpha)
962 : !> with alpha = 1.0_dp/((a + b + c)/((a + b)*c) + 4.0_dp*a_mm),
963 : !> P(R) = (a*(RA + R) + b*RB)/(a + b)
964 : !> \param S_R ...
965 : !> \param RA ...
966 : !> \param RB ...
967 : !> \param RC ...
968 : !> \param zeta ...
969 : !> \param zetb ...
970 : !> \param zetc ...
971 : !> \param a_mm ...
972 : !> \param lgth ...
973 : !> \param R_c ...
974 : ! **************************************************************************************************
975 : #:for prop_exp in range(0,2)
976 : #:for l_max in range(0, lmax_unroll+1)
977 : #:for m_max in range(0, lmax_unroll+1)
978 : #:for n_max in range(0, lmax_unroll+1)
979 : #:set l_tot_max = l_max + m_max + n_max
980 6231960 : PURE SUBROUTINE pgf_sum_3c_rspace_1d_${l_max}$_${m_max}$_${n_max}$_exp_${prop_exp}$ ( &
981 6231960 : S_R, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_c)
982 : REAL(KIND=dp), DIMENSION(0:, 0:, 0:), INTENT(OUT) :: S_R
983 : REAL(KIND=dp), INTENT(IN) :: RA, RB, RC, zeta, zetb, zetc, a_mm, lgth
984 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: R_c
985 : INTEGER :: rr1, rr2, rr2_l, rr2_r, rr1_l, rr1_r
986 : REAL(KIND=dp) :: alpha, alpha_E, dR, exp2_Rsq, R, R1, R_offset, &
987 : R_pow_t, R_tmp, rr1_delta, rr2_delta
988 :
989 : #:if l_tot_max > 0
990 : REAL(KIND=dp) :: c1, c2, c3
991 : #:endif
992 : REAL(KIND=dp) :: ${", ".join(["S_R_t_{}".format(t) for t in range(0,l_tot_max+1)])}$
993 : REAL(KIND=dp) :: ${", ".join(["S_R_t2_{}".format(t) for t in range(0,l_tot_max+1)])}$
994 : REAL(KIND=dp) :: ${", ".join([", ".join(["h_to_c_{}_{}".format(l1,l2) for l1 in range(0,l2+1)]) for l2 in range(0,l_tot_max+1)])}$
995 : REAL(KIND=dp) :: ${", ".join([", ".join([", ".join(["E_{}_{}_{}".format(t,l,m) for t in range(0,l+m+1)]) for l in range(0,l_max+1)]) for m in range(0,m_max+1)])}$
996 :
997 : #:if prop_exp
998 : REAL(KIND=dp) :: exp2_2RdR, exp_dRsq, exp_2dRsq !exp_E_dRsq, exp_E_2dRsq, exp_E_2RdR, exp_E_Rsq
999 : #:endif
1000 :
1001 6231960 : dR = lgth
1002 6231960 : alpha = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm)
1003 :
1004 113325552 : S_R(:, :, :) = 0.0_dp
1005 :
1006 6231960 : R_offset = RC - (zeta*RA + zetb*RB)/(zeta + zetb)
1007 :
1008 6231960 : h_to_c_0_0 = SQRT(alpha/pi)
1009 :
1010 : #:for l in range(0, l_tot_max)
1011 : #:for k in range(0, l+2)
1012 : #:if k<l or k>0
1013 5509656 : h_to_c_${k}$_${l+1}$ = #{if k<l}#-${k+1}$*h_to_c_${k+1}$_${l}$#{endif}##{if k > 0}#+2*alpha*h_to_c_${k-1}$_${l}$#{endif}#
1014 : #:else
1015 5509656 : h_to_c_${k}$_${l+1}$ = 0.0_dp
1016 : #:endif
1017 : #:endfor
1018 : #:endfor
1019 :
1020 : #:if prop_exp
1021 3268062 : exp_dRsq = exp(-alpha*dR*dR)
1022 3268062 : exp_2dRsq = exp_dRsq*exp_dRsq
1023 : #:endif
1024 :
1025 6231960 : rr1_delta = (RA - RB)/dR
1026 :
1027 6231960 : rr1_l = CEILING(-R_c(1) + rr1_delta)
1028 6231960 : rr1_r = FLOOR(R_c(1) + rr1_delta)
1029 :
1030 6231960 : R1 = rr1_l*dR
1031 :
1032 6231960 : alpha_E = zeta*zetb/(zeta + zetb)
1033 :
1034 20231352 : DO rr1 = rr1_l, rr1_r
1035 : #:for t in range(0, l_tot_max + 1)
1036 13999392 : S_R_t_${t}$ = 0.0_dp
1037 13999392 : S_R_t2_${t}$ = 0.0_dp
1038 : #:endfor
1039 13999392 : R_tmp = R_offset + R1*zeta/(zeta + zetb)
1040 13999392 : rr2_delta = -R_tmp/dR
1041 :
1042 13999392 : rr2_l = CEILING(-R_c(2) + rr2_delta)
1043 13999392 : rr2_r = FLOOR(R_c(2) + rr2_delta)
1044 :
1045 13999392 : R = R_tmp + (rr2_l)*dR
1046 :
1047 : #:if prop_exp
1048 9119092 : exp2_2RdR = exp(-2*alpha*R*dR)
1049 9119092 : exp2_Rsq = exp(-alpha*R*R)
1050 : #:endif
1051 :
1052 71598485 : DO rr2 = rr2_l, rr2_r
1053 57599093 : R_pow_t = 1.0_dp
1054 : #:if not prop_exp
1055 7179779 : exp2_Rsq = exp(-alpha*R*R)
1056 : #:endif
1057 : #:for t in range(0, l_tot_max + 1)
1058 57599093 : S_R_t_${t}$ = S_R_t_${t}$+R_pow_t*exp2_Rsq
1059 : #:if t<l_tot_max
1060 51468255 : R_pow_t = R_pow_t*R
1061 : #:endif
1062 : #:endfor
1063 :
1064 : #:if prop_exp
1065 50419314 : exp2_Rsq = exp2_Rsq*exp_dRsq*exp2_2RdR
1066 50419314 : exp2_2RdR = exp2_2RdR*exp_2dRsq
1067 : #:endif
1068 71598485 : R = R + dR
1069 : END DO
1070 :
1071 : ! C --> H
1072 : #:for l in range(0, l_tot_max+1)
1073 : #:for k in range(0, l+1)
1074 13999392 : S_R_t2_${l}$ = S_R_t2_${l}$+h_to_c_${k}$_${l}$*S_R_t_${k}$
1075 : #:endfor
1076 : #:endfor
1077 :
1078 : ! H --> HH
1079 13999392 : E_0_0_0 = exp(-alpha_E*(RA - RB - R1)*(RA - RB - R1))
1080 :
1081 : #:if l_tot_max > 0
1082 12382032 : c1 = 1.0_dp/(zeta + zetb)
1083 12382032 : c2 = 2.0_dp*(zetb/(zeta + zetb))*(RB - (RA - R1))
1084 12382032 : c3 = 2.0_dp*(zeta/(zeta + zetb))*(RA - R1 - RB)
1085 : #:endif
1086 :
1087 : #:for m in range(0,m_max+1)
1088 : #:for l in range(0,l_max+1)
1089 : #:for t in range(0,l+m+2)
1090 : #:if l < l_max
1091 : E_${t}$_${l+1}$_${m}$ = zeta*(#{if t>0}#c1*E_${t-1}$_${l}$_${m}$#{endif}# &
1092 : #{if t<=l+m}# +c2*E_${t}$_${l}$_${m}$&#{endif}#
1093 : #{if t<l+m}# +${2*(t+1)}$*E_${t+1}$_${l}$_${m}$ &#{endif}#
1094 8760456 : #{if l>0 and t<=l-1+m}#-${2*l}$*E_${t}$_${l-1}$_${m}$#{endif}#)
1095 : #:endif
1096 : #:if m < m_max
1097 : E_${t}$_${l}$_${m+1}$ = zetb*(#{if t>0}#c1*E_${t-1}$_${l}$_${m}$#{endif}# &
1098 : #{if t<=l+m}#+c3*E_${t}$_${l}$_${m}$&#{endif}#
1099 : #{if t<l+m}#+${2*(t+1)}$*E_${t+1}$_${l}$_${m}$&#{endif}#
1100 8761656 : #{if m>0 and t<=m-1+l}#-${2*m}$*E_${t}$_${l}$_${m-1}$#{endif}#)
1101 : #:endif
1102 : #:endfor
1103 : #:endfor
1104 : #:endfor
1105 :
1106 : #:for n in range(0, n_max+1)
1107 : #:for m in range(0, m_max+1)
1108 : #:for l in range(0, l_max+1)
1109 : #:for t in range(0, l+m+1)
1110 13999392 : S_R(${l}$, ${m}$, ${n}$) = S_R(${l}$, ${m}$, ${n}$) + E_${t}$_${l}$_${m}$*(${(-1)**n}$)*S_R_t2_${t+n}$ ! cost: 5 flops
1111 : #:endfor
1112 : #:endfor
1113 : #:endfor
1114 : #:endfor
1115 20231352 : R1 = R1 + dR
1116 : END DO
1117 :
1118 113325552 : S_R = S_R*pi**(-0.5_dp)*((zeta + zetb)/(zeta*zetb))**(-0.5_dp)
1119 6231960 : END SUBROUTINE pgf_sum_3c_rspace_1d_${l_max}$_${m_max}$_${n_max}$_exp_${prop_exp}$
1120 : #:endfor
1121 : #:endfor
1122 : #:endfor
1123 : #:endfor
1124 :
1125 : ! **************************************************************************************************
1126 : !> \brief As pgf_sum_3c_1d but 3d sum required for non-orthorhombic cells
1127 : !> \param S_G ...
1128 : !> \param la_max ...
1129 : !> \param lb_max ...
1130 : !> \param lc_max ...
1131 : !> \param RA ...
1132 : !> \param RB ...
1133 : !> \param RC ...
1134 : !> \param zeta ...
1135 : !> \param zetb ...
1136 : !> \param zetc ...
1137 : !> \param a_mm ...
1138 : !> \param hmat ...
1139 : !> \param h_inv ...
1140 : !> \param vol ...
1141 : !> \param G_bounds_1 ...
1142 : !> \param R_bounds_2 ...
1143 : !> \param R_bounds_3 ...
1144 : !> \param G_rads_1 ...
1145 : !> \param R_rads_2 ...
1146 : !> \param R_rads_3 ...
1147 : !> \param method ...
1148 : !> \param method_out ...
1149 : !> \param order ...
1150 : ! **************************************************************************************************
1151 148500 : SUBROUTINE pgf_sum_3c_3d(S_G, la_max, lb_max, lc_max, RA, RB, RC, &
1152 : zeta, zetb, zetc, a_mm, hmat, h_inv, vol, &
1153 : G_bounds_1, R_bounds_2, R_bounds_3, &
1154 : G_rads_1, R_rads_2, R_rads_3, &
1155 : method, method_out, order)
1156 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: S_G
1157 : INTEGER, INTENT(IN) :: la_max, lb_max, lc_max
1158 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: RA, RB, RC
1159 : REAL(KIND=dp), INTENT(IN) :: zeta, zetb, zetc, a_mm
1160 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv
1161 : REAL(KIND=dp), INTENT(IN) :: vol
1162 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: G_bounds_1, R_bounds_2
1163 : REAL(KIND=dp), DIMENSION(2, 3), INTENT(IN) :: R_bounds_3
1164 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G_rads_1, R_rads_2
1165 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: R_rads_3
1166 : INTEGER, INTENT(IN) :: method
1167 : INTEGER, INTENT(OUT), OPTIONAL :: method_out
1168 : INTEGER, INTENT(IN), OPTIONAL :: order
1169 :
1170 : INTEGER :: l_max, m_max, n_max, sum_method, &
1171 : sum_order
1172 : LOGICAL :: assert_stm
1173 : REAL(KIND=dp) :: alpha, beta, G_rad, gamma, R_rad
1174 148500 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: S_G_tmp
1175 : REAL(KIND=dp), DIMENSION(3) :: G_bound, R1, R2, R_bound
1176 : REAL(KIND=dp), DIMENSION(3, 3) :: ht
1177 :
1178 148500 : IF (PRESENT(order)) THEN
1179 0 : sum_order = order
1180 : ELSE
1181 148500 : sum_order = 0
1182 : END IF
1183 :
1184 148500 : sum_method = method
1185 :
1186 148500 : alpha = 0.25_dp/zeta
1187 148500 : beta = 0.25_dp/zetb
1188 148500 : gamma = 0.25_dp/zetc + a_mm
1189 :
1190 148500 : l_max = la_max
1191 148500 : m_max = lb_max
1192 148500 : n_max = lc_max
1193 :
1194 594000 : R1 = RB - RA
1195 594000 : R2 = RC - RB
1196 :
1197 1930500 : ht = twopi*TRANSPOSE(h_inv)
1198 :
1199 0 : SELECT CASE (sum_method)
1200 : CASE (1) ! sum_G sum_G' C(G) C(G,G') C(G')
1201 :
1202 0 : IF (sum_order == 0) THEN
1203 0 : sum_order = MAXLOC(G_bounds_1(:, 1), DIM=1)
1204 : assert_stm = MINLOC(G_bounds_1(:, 1), DIM=1) == &
1205 : MINLOC(G_bounds_1(:, 2), DIM=1) .AND. &
1206 : MINLOC(G_bounds_1(:, 1), DIM=1) == &
1207 0 : MINLOC(G_bounds_1(:, 3), DIM=1)
1208 0 : CPASSERT(assert_stm)
1209 : END IF
1210 :
1211 0 : CALL pgf_sum_3c_gspace_3d(S_G, l_max, m_max, n_max, R1, R2, alpha, beta, gamma, ht, vol, G_bounds_1, G_rads_1, sum_order)
1212 :
1213 : CASE (2) ! sum_G sum_R C(G) C(R)
1214 5846 : IF (sum_order == 0) THEN
1215 23384 : sum_order = MINLOC(G_bounds_1(:, 1), DIM=1)
1216 : assert_stm = MINLOC(G_bounds_1(:, 1), DIM=1) == &
1217 : MINLOC(G_bounds_1(:, 2), DIM=1) .AND. &
1218 : MINLOC(G_bounds_1(:, 1), DIM=1) == &
1219 93536 : MINLOC(G_bounds_1(:, 3), DIM=1)
1220 0 : CPASSERT(assert_stm)
1221 : END IF
1222 5846 : R_rad = R_rads_2(sum_order)
1223 5846 : G_rad = G_rads_1(sum_order)
1224 23384 : R_bound = R_bounds_2(sum_order, :)
1225 23384 : G_bound = G_bounds_1(sum_order, :)
1226 142654 : SELECT CASE (sum_order)
1227 : CASE (1)
1228 0 : ALLOCATE (S_G_tmp(ncoset(l_max), ncoset(m_max), ncoset(n_max)))
1229 : CALL pgf_sum_product_3c_gspace_3d(S_G_tmp, l_max, m_max, n_max, R1, R2, alpha, beta, gamma, hmat, h_inv, vol, &
1230 0 : R_bound, G_bound, R_rad, G_rad)
1231 0 : S_G = RESHAPE(S_G_tmp, SHAPE(S_G), order=[1, 2, 3])
1232 : CASE (2)
1233 0 : ALLOCATE (S_G_tmp(ncoset(m_max), ncoset(l_max), ncoset(n_max)))
1234 : CALL pgf_sum_product_3c_gspace_3d(S_G_tmp, m_max, l_max, n_max, -R1, R1 + R2, beta, alpha, gamma, hmat, h_inv, vol, &
1235 0 : R_bound, G_bound, R_rad, G_rad)
1236 0 : S_G = RESHAPE(S_G_tmp, SHAPE(S_G), order=[2, 1, 3])
1237 : CASE (3)
1238 29230 : ALLOCATE (S_G_tmp(ncoset(n_max), ncoset(m_max), ncoset(l_max)))
1239 : CALL pgf_sum_product_3c_gspace_3d(S_G_tmp, n_max, m_max, l_max, -R2, -R1, gamma, beta, alpha, hmat, h_inv, vol, &
1240 40922 : R_bound, G_bound, R_rad, G_rad)
1241 29230 : S_G = RESHAPE(S_G_tmp, SHAPE(S_G), order=[3, 2, 1])
1242 : END SELECT
1243 : CASE (3) ! sum_R sum_R' H(R, R')
1244 : CALL pgf_sum_3c_rspace_3d(S_G, l_max, m_max, n_max, RA, RB, RC, zeta, zetb, zetc, a_mm, hmat, h_inv, &
1245 142654 : R_bounds_3, R_rads_3)
1246 24412738 : S_G = S_G*pi**(-1.5_dp)*((zeta + zetb)/(zeta*zetb))**(-1.5_dp)
1247 : END SELECT
1248 :
1249 148500 : IF (PRESENT(method_out)) method_out = sum_method
1250 :
1251 148500 : END SUBROUTINE pgf_sum_3c_3d
1252 :
1253 : ! **************************************************************************************************
1254 : !> \brief Roughly estimated number of floating point operations
1255 : !> \param ns_G1 ...
1256 : !> \param ns_G2 ...
1257 : !> \param l ...
1258 : !> \param m ...
1259 : !> \param n ...
1260 : !> \return ...
1261 : ! **************************************************************************************************
1262 148500 : PURE FUNCTION nsum_3c_gspace_3d(ns_G1, ns_G2, l, m, n)
1263 : REAL(KIND=dp), INTENT(IN) :: ns_G1, ns_G2
1264 : INTEGER, INTENT(IN) :: l, m, n
1265 : INTEGER(KIND=int_8) :: nsum_3c_gspace_3d
1266 :
1267 148500 : nsum_3c_gspace_3d = NINT(ns_G1*ns_G2*(5*exp_w + ncoset(l)*ncoset(m)*ncoset(n)*4), KIND=int_8)
1268 :
1269 148500 : END FUNCTION nsum_3c_gspace_3d
1270 :
1271 : ! **************************************************************************************************
1272 : !> \brief ...
1273 : !> \param S_G ...
1274 : !> \param l_max ...
1275 : !> \param m_max ...
1276 : !> \param n_max ...
1277 : !> \param R1 ...
1278 : !> \param R2 ...
1279 : !> \param alpha ...
1280 : !> \param beta ...
1281 : !> \param gamma ...
1282 : !> \param ht ...
1283 : !> \param vol ...
1284 : !> \param G_c ...
1285 : !> \param G_rad ...
1286 : !> \param sum_order ...
1287 : !> \param coulomb ...
1288 : ! **************************************************************************************************
1289 0 : PURE SUBROUTINE pgf_sum_3c_gspace_3d(S_G, l_max, m_max, n_max, R1, R2, alpha, beta, gamma, &
1290 : ht, vol, G_c, G_rad, sum_order, coulomb)
1291 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: S_G
1292 : INTEGER, INTENT(IN) :: l_max, m_max, n_max
1293 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R1, R2
1294 : REAL(KIND=dp), INTENT(IN) :: alpha, beta, gamma
1295 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: ht
1296 : REAL(KIND=dp), INTENT(IN) :: vol
1297 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: G_c
1298 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G_rad
1299 : INTEGER, INTENT(IN) :: sum_order
1300 : LOGICAL, INTENT(IN), OPTIONAL :: coulomb
1301 :
1302 : INTEGER, DIMENSION(3) :: G1c, G2c, G3c
1303 : INTEGER :: g1x, g1y, g1z, g2x, g2y, g2z, g3x, g3y, &
1304 : g3z
1305 : COMPLEX(KIND=dp), DIMENSION(ncoset(l_max), ncoset( &
1306 0 : m_max), ncoset(n_max)) :: S_G_c
1307 : LOGICAL :: use_coulomb
1308 : REAL(KIND=dp) :: G1_sq, G2_sq, G3_sq
1309 : REAL(KIND=dp), DIMENSION(3) :: G1, G1_x, G1_y, G1_z, G2, G2_x, G2_y, &
1310 : G2_z, G3, G3_x, G3_y, G3_z, G_rads_sq
1311 :
1312 0 : S_G_c(:, :, :) = 0.0_dp
1313 :
1314 0 : G1c = FLOOR(G_c(1, :))
1315 0 : G2c = FLOOR(G_c(2, :))
1316 0 : G3c = FLOOR(G_c(3, :))
1317 :
1318 : ! we can not precompute exponentials for general cell
1319 : ! Perform double G sum
1320 0 : G_rads_sq = G_rad**2
1321 :
1322 0 : IF (PRESENT(coulomb)) THEN
1323 0 : use_coulomb = coulomb
1324 : ELSE
1325 0 : use_coulomb = .FALSE.
1326 : END IF
1327 :
1328 0 : SELECT CASE (sum_order)
1329 : CASE (1)
1330 0 : DO g2x = -G2c(1), G2c(1)
1331 0 : G2_x = ht(:, 1)*g2x
1332 0 : DO g2y = -G2c(2), G2c(2)
1333 0 : G2_y = ht(:, 2)*g2y
1334 0 : DO g2z = -G2c(3), G2c(3)
1335 0 : G2_z = ht(:, 3)*g2z
1336 0 : G2 = G2_x + G2_y + G2_z
1337 0 : G2_sq = G2(1)**2 + G2(2)**2 + G2(3)**2
1338 0 : IF (G2_sq > G_rads_sq(2)) CYCLE
1339 0 : DO g3x = -G3c(1), G3c(1)
1340 0 : G3_x = ht(:, 1)*g3x
1341 0 : DO g3y = -G3c(2), G3c(2)
1342 0 : G3_y = ht(:, 2)*g3y
1343 0 : DO g3z = -G3c(3), G3c(3)
1344 0 : G3_z = ht(:, 3)*g3z
1345 0 : G3 = G3_x + G3_y + G3_z
1346 0 : G3_sq = G3(1)**2 + G3(2)**2 + G3(3)**2
1347 0 : IF (G3_sq > G_rads_sq(3)) CYCLE
1348 0 : G1 = G3 - G2
1349 0 : G1_sq = G1(1)**2 + G1(2)**2 + G1(3)**2
1350 0 : IF (G1_sq > G_rads_sq(1)) CYCLE
1351 0 : IF (use_coulomb) THEN
1352 0 : IF (g3x == 0 .AND. g3y == 0 .AND. g3z == 0) CYCLE
1353 : END IF
1354 : CALL pgf_product_3c_gspace_3d(S_G_c, G1, G1_sq, G2, G2_sq, G3, G3_sq, l_max, m_max, n_max, &
1355 0 : alpha, beta, gamma, R1, R2, use_coulomb)
1356 : END DO
1357 : END DO
1358 : END DO
1359 : END DO
1360 : END DO
1361 : END DO
1362 : CASE (2)
1363 0 : DO g1x = -G1c(1), G1c(1)
1364 0 : G1_x = ht(:, 1)*g1x
1365 0 : DO g1y = -G1c(2), G1c(2)
1366 0 : G1_y = ht(:, 2)*g1y
1367 0 : DO g1z = -G1c(3), G1c(3)
1368 0 : G1_z = ht(:, 3)*g1z
1369 0 : G1 = G1_x + G1_y + G1_z
1370 0 : G1_sq = G1(1)**2 + G1(2)**2 + G1(3)**2
1371 0 : IF (G1_sq > G_rads_sq(1)) CYCLE
1372 0 : DO g3x = -G3c(1), G3c(1)
1373 0 : G3_x = ht(:, 1)*g3x
1374 0 : DO g3y = -G3c(2), G3c(2)
1375 0 : G3_y = ht(:, 2)*g3y
1376 0 : DO g3z = -G3c(3), G3c(3)
1377 0 : G3_z = ht(:, 3)*g3z
1378 0 : G3 = G3_x + G3_y + G3_z
1379 0 : G3_sq = G3(1)**2 + G3(2)**2 + G3(3)**2
1380 0 : IF (G3_sq > G_rads_sq(3)) CYCLE
1381 0 : G2 = G3 - G1
1382 0 : G2_sq = G2(1)**2 + G2(2)**2 + G2(3)**2
1383 0 : IF (G2_sq > G_rads_sq(2)) CYCLE
1384 0 : IF (use_coulomb) THEN
1385 0 : IF (g3x == 0 .AND. g3y == 0 .AND. g3z == 0) CYCLE
1386 : END IF
1387 : CALL pgf_product_3c_gspace_3d(S_G_c, G1, G1_sq, G2, G2_sq, G3, G3_sq, l_max, m_max, n_max, &
1388 0 : alpha, beta, gamma, R1, R2, use_coulomb)
1389 : END DO
1390 : END DO
1391 : END DO
1392 : END DO
1393 : END DO
1394 : END DO
1395 : CASE (3)
1396 0 : DO g1x = -G1c(1), G1c(1)
1397 0 : G1_x = ht(:, 1)*g1x
1398 0 : DO g1y = -G1c(2), G1c(2)
1399 0 : G1_y = ht(:, 2)*g1y
1400 0 : DO g1z = -G1c(3), G1c(3)
1401 0 : G1_z = ht(:, 3)*g1z
1402 0 : G1 = G1_x + G1_y + G1_z
1403 0 : G1_sq = G1(1)**2 + G1(2)**2 + G1(3)**2
1404 0 : IF (G1_sq > G_rads_sq(1)) CYCLE
1405 0 : DO g2x = -G2c(1), G2c(1)
1406 0 : G2_x = ht(:, 1)*g2x
1407 0 : DO g2y = -G2c(2), G2c(2)
1408 0 : G2_y = ht(:, 2)*g2y
1409 0 : DO g2z = -G2c(3), G2c(3)
1410 0 : G2_z = ht(:, 3)*g2z
1411 0 : G2 = G2_x + G2_y + G2_z
1412 0 : G2_sq = G2(1)**2 + G2(2)**2 + G2(3)**2
1413 0 : IF (G2_sq > G_rads_sq(2)) CYCLE
1414 0 : G3 = G1 + G2
1415 0 : G3_sq = G3(1)**2 + G3(2)**2 + G3(3)**2
1416 0 : IF (G3_sq > G_rads_sq(3)) CYCLE
1417 0 : IF (use_coulomb) THEN
1418 0 : IF (g1x + g2x == 0 .AND. g1y + g2y == 0 .AND. g1z + g2z == 0) CYCLE
1419 : END IF
1420 : CALL pgf_product_3c_gspace_3d(S_G_c, G1, G1_sq, G2, G2_sq, G3, G3_sq, l_max, m_max, n_max, &
1421 0 : alpha, beta, gamma, R1, R2, use_coulomb)
1422 : END DO
1423 : END DO
1424 : END DO
1425 : END DO
1426 : END DO
1427 : END DO
1428 : END SELECT
1429 0 : S_G = REAL(S_G_c, KIND=dp)/vol**2
1430 :
1431 0 : END SUBROUTINE pgf_sum_3c_gspace_3d
1432 :
1433 : ! **************************************************************************************************
1434 : !> \brief ...
1435 : !> \param S_G ...
1436 : !> \param G1 ...
1437 : !> \param G1_sq ...
1438 : !> \param G2 ...
1439 : !> \param G2_sq ...
1440 : !> \param G3 ...
1441 : !> \param G3_sq ...
1442 : !> \param l_max ...
1443 : !> \param m_max ...
1444 : !> \param n_max ...
1445 : !> \param alpha ...
1446 : !> \param beta ...
1447 : !> \param gamma ...
1448 : !> \param R1 ...
1449 : !> \param R2 ...
1450 : !> \param coulomb ...
1451 : ! **************************************************************************************************
1452 0 : PURE SUBROUTINE pgf_product_3c_gspace_3d(S_G, G1, G1_sq, G2, G2_sq, G3, G3_sq, l_max, m_max, n_max, &
1453 : alpha, beta, gamma, R1, R2, coulomb)
1454 :
1455 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: S_G
1456 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G1
1457 : REAL(KIND=dp), INTENT(IN) :: G1_sq
1458 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G2
1459 : REAL(KIND=dp), INTENT(IN) :: G2_sq
1460 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G3
1461 : REAL(KIND=dp), INTENT(IN) :: G3_sq
1462 : INTEGER, INTENT(IN) :: l_max, m_max, n_max
1463 : REAL(KIND=dp), INTENT(IN) :: alpha, beta, gamma
1464 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R1, R2
1465 : LOGICAL, INTENT(IN) :: coulomb
1466 :
1467 : COMPLEX(KIND=dp) :: exp_tot
1468 : INTEGER :: k, l, lco, lx, ly, lz, m, mco, mx, my, &
1469 : mz, n, nco, nx, ny, nz
1470 0 : COMPLEX(KIND=dp), DIMENSION(ncoset(n_max)) :: S_G_n
1471 0 : COMPLEX(KIND=dp), DIMENSION(ncoset(m_max)) :: S_G_m
1472 0 : COMPLEX(KIND=dp), DIMENSION(ncoset(l_max)) :: S_G_l
1473 0 : REAL(KIND=dp), DIMENSION(3, 0:l_max) :: G1_pow_l
1474 0 : REAL(KIND=dp), DIMENSION(3, 0:m_max) :: G2_pow_m
1475 0 : REAL(KIND=dp), DIMENSION(3, 0:n_max) :: G3_pow_n
1476 :
1477 : exp_tot = EXP(gaussi*DOT_PRODUCT(G1, R1))*EXP(-alpha*G1_sq)* & ! cost: 5 exp_w flops
1478 : EXP(-beta*G2_sq)* &
1479 0 : EXP(-gamma*G3_sq)*EXP(gaussi*DOT_PRODUCT(G3, R2))
1480 :
1481 0 : IF (coulomb) exp_tot = exp_tot/G3_sq
1482 :
1483 0 : DO k = 1, 3
1484 0 : G1_pow_l(k, 0) = 1.0_dp
1485 0 : DO l = 1, l_max
1486 0 : G1_pow_l(k, l) = G1_pow_l(k, l - 1)*G1(k)
1487 : END DO
1488 0 : G2_pow_m(k, 0) = 1.0_dp
1489 0 : DO m = 1, m_max
1490 0 : G2_pow_m(k, m) = G2_pow_m(k, m - 1)*G2(k)
1491 : END DO
1492 0 : G3_pow_n(k, 0) = 1.0_dp
1493 0 : DO n = 1, n_max
1494 0 : G3_pow_n(k, n) = G3_pow_n(k, n - 1)*G3(k)
1495 : END DO
1496 : END DO
1497 :
1498 0 : DO lco = 1, ncoset(l_max)
1499 0 : CALL get_l(lco, l, lx, ly, lz)
1500 0 : S_G_l(lco) = G1_pow_l(1, lx)*G1_pow_l(2, ly)*G1_pow_l(3, lz)*(-1.0_dp)**l*i_pow(l)
1501 : END DO
1502 :
1503 0 : DO mco = 1, ncoset(m_max)
1504 0 : CALL get_l(mco, m, mx, my, mz)
1505 0 : S_G_m(mco) = G2_pow_m(1, mx)*G2_pow_m(2, my)*G2_pow_m(3, mz)*(-1.0_dp)**m*i_pow(m)
1506 : END DO
1507 :
1508 0 : DO nco = 1, ncoset(n_max)
1509 0 : CALL get_l(nco, n, nx, ny, nz)
1510 0 : S_G_n(nco) = G3_pow_n(1, nx)*G3_pow_n(2, ny)*G3_pow_n(3, nz)*i_pow(n)
1511 : END DO
1512 :
1513 0 : DO nco = 1, ncoset(n_max)
1514 0 : DO mco = 1, ncoset(m_max)
1515 0 : DO lco = 1, ncoset(l_max)
1516 : S_G(lco, mco, nco) = S_G(lco, mco, nco) + & ! cost: 4 flops
1517 : S_G_l(lco)*S_G_m(mco)*S_G_n(nco)* &
1518 0 : exp_tot
1519 : END DO
1520 : END DO
1521 : END DO
1522 :
1523 0 : END SUBROUTINE pgf_product_3c_gspace_3d
1524 :
1525 : ! **************************************************************************************************
1526 : !> \brief Roughly estimated number of floating point operations
1527 : !> \param ns_G ...
1528 : !> \param ns_R ...
1529 : !> \param l ...
1530 : !> \param m ...
1531 : !> \param n ...
1532 : !> \return ...
1533 : ! **************************************************************************************************
1534 148500 : PURE FUNCTION nsum_product_3c_gspace_3d(ns_G, ns_R, l, m, n)
1535 : REAL(KIND=dp), INTENT(IN) :: ns_G, ns_R
1536 : INTEGER, INTENT(IN) :: l, m, n
1537 : INTEGER(KIND=int_8) :: nsum_product_3c_gspace_3d
1538 :
1539 : nsum_product_3c_gspace_3d = &
1540 : NINT( &
1541 : ns_G*( &
1542 : (exp_w*2) + &
1543 : ns_R*(exp_w*2 + ncoset(l + m)*7) + &
1544 : 3*nsum_gaussian_overlap(l, m, 1) + &
1545 : ncoset(l)*ncoset(m)*(ncoset(l + m)*4 + ncoset(n)*8)), &
1546 148500 : KIND=int_8)
1547 148500 : END FUNCTION nsum_product_3c_gspace_3d
1548 :
1549 : ! **************************************************************************************************
1550 : !> \brief ...
1551 : !> \param S_G ...
1552 : !> \param l_max ...
1553 : !> \param m_max ...
1554 : !> \param n_max ...
1555 : !> \param R1 ...
1556 : !> \param R2 ...
1557 : !> \param alpha ...
1558 : !> \param beta ...
1559 : !> \param gamma ...
1560 : !> \param hmat ...
1561 : !> \param h_inv ...
1562 : !> \param vol ...
1563 : !> \param R_c ...
1564 : !> \param G_c ...
1565 : !> \param R_rad ...
1566 : !> \param G_rad ...
1567 : ! **************************************************************************************************
1568 5846 : PURE SUBROUTINE pgf_sum_product_3c_gspace_3d(S_G, l_max, m_max, n_max, R1, R2, alpha, beta, gamma, &
1569 : hmat, h_inv, vol, R_c, G_c, R_rad, G_rad)
1570 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: S_G
1571 : INTEGER, INTENT(IN) :: l_max, m_max, n_max
1572 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R1, R2
1573 : REAL(KIND=dp), INTENT(IN) :: alpha, beta, gamma
1574 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv
1575 : REAL(KIND=dp), INTENT(IN) :: vol
1576 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R_c, G_c
1577 : REAL(KIND=dp), INTENT(IN) :: R_rad, G_rad
1578 :
1579 : COMPLEX(KIND=dp) :: exp_tot
1580 : INTEGER :: gx, gy, gz, k, l, lco, lnco, lx, ly, lz, &
1581 : m, mco, n, nco
1582 : COMPLEX(KIND=dp), &
1583 11692 : DIMENSION(ncoset(m_max), ncoset(n_max)) :: S_R
1584 : COMPLEX(KIND=dp), DIMENSION(ncoset(l_max), ncoset( &
1585 11692 : m_max), ncoset(n_max)) :: S_G_c
1586 : REAL(KIND=dp) :: G_rad_sq, G_sq, R_rad_sq
1587 : REAL(KIND=dp), DIMENSION(3) :: G, G_x, G_y, G_z
1588 11692 : REAL(KIND=dp), DIMENSION(3, 0:l_max) :: G_pow_l
1589 : REAL(KIND=dp), DIMENSION(3, 3) :: ht
1590 5846 : REAL(KIND=dp), DIMENSION(ncoset(l_max)) :: S_G_c_l
1591 :
1592 276845 : S_G_c(:, :, :) = 0.0_dp
1593 :
1594 5846 : G_rad_sq = G_rad**2
1595 5846 : R_rad_sq = R_rad**2
1596 :
1597 5846 : lnco = ncoset(l_max)
1598 :
1599 75998 : ht = twopi*TRANSPOSE(h_inv)
1600 36078 : DO gx = -FLOOR(G_c(1)), FLOOR(G_c(1))
1601 120928 : G_x = ht(:, 1)*gx
1602 194252 : DO gy = -FLOOR(G_c(2)), FLOOR(G_c(2))
1603 632696 : G_y = ht(:, 2)*gy
1604 1028374 : DO gz = -FLOOR(G_c(3)), FLOOR(G_c(3))
1605 3359872 : G_z = ht(:, 3)*gz
1606 3359872 : G = G_x + G_y + G_z
1607 839968 : G_sq = G(1)**2 + G(2)**2 + G(3)**2
1608 839968 : IF (G_sq > G_rad_sq) CYCLE
1609 :
1610 1350880 : exp_tot = EXP(-alpha*G_sq)*EXP(gaussi*DOT_PRODUCT(G, R1)) ! cost: exp_w*2 flops
1611 :
1612 1350880 : DO k = 1, 3
1613 1013160 : G_pow_l(k, 0) = 1.0_dp
1614 2236258 : DO l = 1, l_max
1615 1898538 : G_pow_l(k, l) = G_pow_l(k, l - 1)*G(k)
1616 : END DO
1617 : END DO
1618 :
1619 337720 : CALL pgf_sum_product_3c_rspace_3d(S_R, m_max, n_max, G, R2, beta, gamma, hmat, h_inv, vol, R_c, R_rad_sq)
1620 :
1621 1947011 : DO lco = 1, ncoset(l_max)
1622 1609291 : CALL get_l(lco, l, lx, ly, lz)
1623 1947011 : S_G_c_l(lco) = G_pow_l(1, lx)*G_pow_l(2, ly)*G_pow_l(3, lz)*(-1.0_dp)**l
1624 : END DO
1625 :
1626 1519228 : DO nco = 1, ncoset(n_max)
1627 1023334 : CALL get_l(nco, n)
1628 4612179 : DO mco = 1, ncoset(m_max)
1629 2748877 : CALL get_l(mco, m)
1630 15176745 : DO lco = 1, ncoset(l_max)
1631 11404534 : CALL get_l(lco, l)
1632 : S_G_c(lco, mco, nco) = & ! cost: 8 flops
1633 : S_G_c(lco, mco, nco) + &
1634 : S_G_c_l(lco)* &
1635 14153411 : exp_tot*i_pow(l + m + n)*(-1.0_dp)**m*S_R(mco, nco)
1636 : END DO
1637 : END DO
1638 : END DO
1639 : END DO
1640 : END DO
1641 : END DO
1642 276845 : S_G = REAL(S_G_c, KIND=dp)/vol**2
1643 :
1644 5846 : END SUBROUTINE pgf_sum_product_3c_gspace_3d
1645 :
1646 : ! **************************************************************************************************
1647 : !> \brief ...
1648 : !> \param S_R ...
1649 : !> \param l_max ...
1650 : !> \param m_max ...
1651 : !> \param G ...
1652 : !> \param R ...
1653 : !> \param alpha ...
1654 : !> \param beta ...
1655 : !> \param hmat ...
1656 : !> \param h_inv ...
1657 : !> \param vol ...
1658 : !> \param R_c ...
1659 : !> \param R_rad_sq ...
1660 : ! **************************************************************************************************
1661 337720 : PURE SUBROUTINE pgf_sum_product_3c_rspace_3d(S_R, l_max, m_max, G, R, alpha, beta, hmat, h_inv, vol, R_c, R_rad_sq)
1662 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: S_R
1663 : INTEGER, INTENT(IN) :: l_max, m_max
1664 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G, R
1665 : REAL(KIND=dp), INTENT(IN) :: alpha, beta
1666 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv
1667 : REAL(KIND=dp), INTENT(IN) :: vol
1668 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R_c
1669 : REAL(KIND=dp), INTENT(IN) :: R_rad_sq
1670 :
1671 : COMPLEX(KIND=dp) :: exp_tot
1672 : INTEGER :: k, l, lco, lx, ly, lz, m, mco, mx, my, &
1673 : mz, sx, sy, sz, t, tco, tx, ty, tz
1674 675440 : COMPLEX(KIND=dp), DIMENSION(ncoset(l_max + m_max)) :: S_R_t
1675 : REAL(KIND=dp) :: c1, c2, Rp_sq
1676 : REAL(KIND=dp), &
1677 675440 : DIMENSION(-1:l_max + m_max + 1, -1:l_max, -1:m_max) :: E1, E2, E3
1678 : REAL(KIND=dp), DIMENSION(3) :: R_l, R_r, Rp, Rx, Ry, Rz, s_shift
1679 675440 : REAL(KIND=dp), DIMENSION(3, 0:l_max + m_max) :: R_pow_t
1680 :
1681 337720 : c1 = 0.25_dp/(alpha + beta)
1682 337720 : c2 = alpha/(alpha + beta)
1683 :
1684 2198096 : S_R_t(:) = 0.0_dp
1685 4109931 : S_R(:, :) = 0.0_dp
1686 :
1687 4390360 : s_shift = MATMUL(h_inv, R)
1688 1350880 : R_l = -R_c + s_shift
1689 1350880 : R_r = R_c + s_shift
1690 :
1691 1084696 : DO sx = CEILING(R_l(1)), FLOOR(R_r(1))
1692 2987904 : Rx = hmat(:, 1)*sx
1693 3046393 : DO sy = CEILING(R_l(2)), FLOOR(R_r(2))
1694 7846788 : Ry = hmat(:, 2)*sy
1695 8205714 : DO sz = CEILING(R_l(3)), FLOOR(R_r(3))
1696 21988164 : Rz = hmat(:, 3)*sz
1697 21988164 : Rp = Rx + Ry + Rz - R
1698 5497041 : Rp_sq = Rp(1)**2 + Rp(2)**2 + Rp(3)**2
1699 5497041 : IF (Rp_sq > R_rad_sq) CYCLE
1700 :
1701 8361484 : exp_tot = EXP(-c1*Rp_sq)*EXP(-gaussi*c2*DOT_PRODUCT(Rp, G)) ! cost: exp_w*2 flops
1702 8361484 : DO k = 1, 3
1703 6271113 : R_pow_t(k, 0) = 1.0_dp
1704 13292668 : DO t = 1, l_max + m_max
1705 11202297 : R_pow_t(k, t) = R_pow_t(k, t - 1)*Rp(k)
1706 : END DO
1707 : END DO
1708 :
1709 13088657 : DO tco = 1, ncoset(l_max + m_max)
1710 9036589 : CALL get_l(tco, t, tx, ty, tz)
1711 : S_R_t(tco) = S_R_t(tco) + & ! cost: 7 flops
1712 : R_pow_t(1, tx)*R_pow_t(2, ty)*R_pow_t(3, tz)* &
1713 14533630 : (-1.0_dp)**t*i_pow(t)*exp_tot
1714 : END DO
1715 :
1716 : END DO
1717 : END DO
1718 : END DO
1719 :
1720 337720 : CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, alpha, beta, G(1), 0.0_dp, 1, E1)
1721 337720 : CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, alpha, beta, G(2), 0.0_dp, 1, E2)
1722 337720 : CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, alpha, beta, G(3), 0.0_dp, 1, E3)
1723 :
1724 1361054 : DO mco = 1, ncoset(m_max)
1725 1023334 : CALL get_l(mco, m, mx, my, mz)
1726 4109931 : DO lco = 1, ncoset(l_max)
1727 2748877 : CALL get_l(lco, l, lx, ly, lz)
1728 8010112 : DO tx = 0, lx + mx
1729 13267665 : DO ty = 0, ly + my
1730 19515063 : DO tz = 0, lz + mz
1731 8996275 : tco = coset(tx, ty, tz)
1732 : S_R(lco, mco) = S_R(lco, mco) + & ! cost: 4 flops
1733 15277162 : E1(tx, lx, mx)*E2(ty, ly, my)*E3(tz, lz, mz)*S_R_t(tco)
1734 :
1735 : END DO
1736 : END DO
1737 : END DO
1738 : END DO
1739 : END DO
1740 :
1741 4109931 : S_R(:, :) = S_R(:, :)*vol/(twopi)**3*(pi/(alpha + beta))**1.5_dp
1742 :
1743 337720 : END SUBROUTINE pgf_sum_product_3c_rspace_3d
1744 :
1745 : ! **************************************************************************************************
1746 : !> \brief Roughly estimated number of floating point operations
1747 : !> \param ns_R1 ...
1748 : !> \param ns_R2 ...
1749 : !> \param l ...
1750 : !> \param m ...
1751 : !> \param n ...
1752 : !> \return ...
1753 : ! **************************************************************************************************
1754 148500 : PURE FUNCTION nsum_3c_rspace_3d(ns_R1, ns_R2, l, m, n)
1755 : REAL(KIND=dp), INTENT(IN) :: ns_R1, ns_R2
1756 : INTEGER, INTENT(IN) :: l, m, n
1757 : INTEGER(KIND=int_8) :: nsum_3c_rspace_3d
1758 :
1759 : nsum_3c_rspace_3d = &
1760 : NINT( &
1761 : ns_R1*( &
1762 : ns_R2*(exp_w + ncoset(l + m + n)*4) + &
1763 : 3*nsum_gaussian_overlap(l, m, 2) + &
1764 : ncoset(m)*ncoset(l)*( &
1765 : ncoset(l + m)*2 + ncoset(n)*ncoset(l + m)*4)), &
1766 148500 : KIND=int_8)
1767 :
1768 148500 : END FUNCTION nsum_3c_rspace_3d
1769 :
1770 : ! **************************************************************************************************
1771 : !> \brief ...
1772 : !> \param S_R ...
1773 : !> \param l_max ...
1774 : !> \param m_max ...
1775 : !> \param n_max ...
1776 : !> \param RA ...
1777 : !> \param RB ...
1778 : !> \param RC ...
1779 : !> \param zeta ...
1780 : !> \param zetb ...
1781 : !> \param zetc ...
1782 : !> \param a_mm ...
1783 : !> \param hmat ...
1784 : !> \param h_inv ...
1785 : !> \param R_c ...
1786 : !> \param R_rad ...
1787 : ! **************************************************************************************************
1788 142654 : SUBROUTINE pgf_sum_3c_rspace_3d(S_R, l_max, m_max, n_max, RA, RB, RC, zeta, zetb, zetc, a_mm, &
1789 : hmat, h_inv, R_c, R_rad)
1790 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: S_R
1791 : INTEGER, INTENT(IN) :: l_max, m_max, n_max
1792 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: RA, RB, RC
1793 : REAL(KIND=dp), INTENT(IN) :: zeta, zetb, zetc, a_mm
1794 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv
1795 : REAL(KIND=dp), DIMENSION(2, 3), INTENT(IN) :: R_c
1796 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: R_rad
1797 :
1798 : INTEGER :: k, l, lco, lx, ly, lz, m, mco, mx, my, &
1799 : mz, n, nco, nx, ny, nz, s1x, s1y, s1z, &
1800 : s2x, s2y, s2z, t, tco, tnco, ttco, &
1801 : ttx, tty, ttz, tx, ty, tz
1802 : REAL(KIND=dp) :: alpha, exp_tot, R1_sq, R_sq
1803 142654 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: h_to_c
1804 : REAL(KIND=dp), &
1805 285308 : DIMENSION(-1:l_max + m_max + 1, -1:l_max, -1:m_max) :: E1, E2, E3
1806 : REAL(KIND=dp), DIMENSION(2) :: R_rad_sq
1807 : REAL(KIND=dp), DIMENSION(3) :: R, R1, R1_l, R1_r, R1_tmp, R1x, R1y, &
1808 : R1z, R2_l, R2_r, R2x, R2y, R2z, &
1809 : R_offset, R_tmp, s1_shift, s2_shift
1810 285308 : REAL(KIND=dp), DIMENSION(3, 0:l_max + m_max + n_max) :: R_pow_t
1811 : REAL(KIND=dp), DIMENSION(ncoset(l_max + m_max), &
1812 285308 : ncoset(l_max), ncoset(m_max)) :: Eco
1813 : REAL(KIND=dp), &
1814 285308 : DIMENSION(ncoset(l_max + m_max + n_max)) :: S_R_t
1815 : REAL(KIND=dp), DIMENSION(ncoset(l_max + m_max + n_max) &
1816 285308 : , ncoset(l_max + m_max + n_max)) :: h_to_c_tco
1817 :
1818 142654 : alpha = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm)
1819 :
1820 71800162 : Eco(:, :, :) = 0.0_dp
1821 24264238 : S_R(:, :, :) = 0.0_dp
1822 157194456 : h_to_c_tco(:, :) = 0.0_dp
1823 :
1824 570616 : R_offset = RC - (zeta*RA + zetb*RB)/(zeta + zetb)
1825 :
1826 142654 : CALL create_hermite_to_cartesian(alpha, l_max + m_max + n_max, h_to_c)
1827 :
1828 2682070 : DO tco = 1, ncoset(l_max + m_max + n_max)
1829 2539416 : CALL get_l(tco, t, tx, ty, tz)
1830 8078511 : DO ttx = 0, tx
1831 18802080 : DO tty = 0, ty
1832 37176724 : DO ttz = 0, tz
1833 20914060 : ttco = coset(ttx, tty, ttz)
1834 31780283 : h_to_c_tco(ttco, tco) = h_to_c(ttx, tx)*h_to_c(tty, ty)*h_to_c(ttz, tz)
1835 : END DO
1836 : END DO
1837 : END DO
1838 : END DO
1839 :
1840 2282464 : s1_shift = MATMUL(h_inv, RA - RB)
1841 570616 : R1_l = -R_c(1, :) + s1_shift
1842 570616 : R1_r = R_c(1, :) + s1_shift
1843 :
1844 427962 : R_rad_sq = R_rad**2
1845 :
1846 348524 : DO s1x = CEILING(R1_l(1)), FLOOR(R1_r(1))
1847 823480 : R1x = hmat(:, 1)*s1x
1848 741859 : DO s1y = CEILING(R1_l(2)), FLOOR(R1_r(2))
1849 1573340 : R1y = hmat(:, 2)*s1y
1850 1532744 : DO s1z = CEILING(R1_l(3)), FLOOR(R1_r(3))
1851 3734156 : R1z = hmat(:, 3)*s1z
1852 3734156 : R1 = R1x + R1y + R1z
1853 13798961 : S_R_t(:) = 0.0_dp
1854 3734156 : R1_tmp = R1 - (RA - RB)
1855 933539 : R1_sq = R1_tmp(1)**2 + R1_tmp(2)**2 + R1_tmp(3)**2
1856 :
1857 933539 : IF (R1_sq > R_rad_sq(1)) CYCLE
1858 :
1859 1625580 : R_tmp = R_offset + R1*zeta/(zeta + zetb)
1860 6502320 : s2_shift = -MATMUL(h_inv, R_tmp)
1861 1625580 : R2_l = -R_c(2, :) + s2_shift
1862 1625580 : R2_r = R_c(2, :) + s2_shift
1863 1468530 : DO s2x = CEILING(R2_l(1)), FLOOR(R2_r(1))
1864 4248540 : R2x = hmat(:, 1)*s2x
1865 6145785 : DO s2y = CEILING(R2_l(2)), FLOOR(R2_r(2))
1866 18709020 : R2y = hmat(:, 2)*s2y
1867 34911377 : DO s2z = CEILING(R2_l(3)), FLOOR(R2_r(3))
1868 116687948 : R2z = hmat(:, 3)*s2z
1869 116687948 : R = R_tmp + R2x + R2y + R2z
1870 29171987 : R_sq = R(1)**2 + R(2)**2 + R(3)**2
1871 :
1872 29171987 : IF (R_sq > R_rad_sq(2)) CYCLE
1873 :
1874 14475027 : exp_tot = EXP(-alpha*R_sq) ! cost: exp_w flops
1875 57900108 : DO k = 1, 3
1876 43425081 : R_pow_t(k, 0) = 1.0_dp
1877 156067857 : DO t = 1, l_max + m_max + n_max
1878 141592830 : R_pow_t(k, t) = R_pow_t(k, t - 1)*R(k)
1879 : END DO
1880 : END DO
1881 272529341 : DO tco = 1, ncoset(l_max + m_max + n_max)
1882 253377059 : CALL get_l(tco, t, tx, ty, tz)
1883 282549046 : S_R_t(tco) = S_R_t(tco) + R_pow_t(1, tx)*R_pow_t(2, ty)*R_pow_t(3, tz)*exp_tot ! cost: 4 flops
1884 : END DO
1885 : END DO
1886 : END DO
1887 : END DO
1888 :
1889 5981328 : S_R_t(:) = MATMUL(TRANSPOSE(h_to_c_tco), S_R_t)*(alpha/pi)**1.5_dp
1890 :
1891 406395 : CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, RA(1) - R1(1), RB(1), 2, E1)
1892 406395 : CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, RA(2) - R1(2), RB(2), 2, E2)
1893 406395 : CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, RA(3) - R1(3), RB(3), 2, E3)
1894 :
1895 1416942 : DO mco = 1, ncoset(m_max)
1896 1010547 : CALL get_l(mco, m, mx, my, mz)
1897 5054628 : DO lco = 1, ncoset(l_max)
1898 3637686 : CALL get_l(lco, l, lx, ly, lz)
1899 10890282 : DO tx = 0, lx + mx
1900 20345458 : DO ty = 0, ly + my
1901 33814276 : DO tz = 0, lz + mz
1902 17106504 : tco = coset(tx, ty, tz)
1903 27572227 : Eco(tco, lco, mco) = E1(tx, lx, mx)*E2(ty, ly, my)*E3(tz, lz, mz) ! cost: 2 flops
1904 : END DO
1905 : END DO
1906 : END DO
1907 : END DO
1908 : END DO
1909 :
1910 3040511 : DO nco = 1, ncoset(n_max)
1911 2240781 : CALL get_l(nco, n, nx, ny, nz)
1912 17172306 : DO tco = 1, ncoset(l_max + m_max)
1913 14525130 : CALL get_l(tco, t, tx, ty, tz)
1914 14525130 : tnco = coset(tx + nx, ty + ny, tz + nz)
1915 : S_R(:, :, nco) = S_R(:, :, nco) + & ! cost: 4 flops
1916 : Eco(tco, :, :)* &
1917 1070960721 : (-1)**n*S_R_t(tnco)
1918 :
1919 : END DO
1920 : END DO
1921 : END DO
1922 : END DO
1923 : END DO
1924 142654 : END SUBROUTINE pgf_sum_3c_rspace_3d
1925 :
1926 : ! **************************************************************************************************
1927 : !> \brief Compute bounding box for ellipsoid. This is needed in order to find summation bounds for
1928 : !> sphere for sums over non-orthogonal lattice vectors.
1929 : !> \param s_rad sphere radius
1930 : !> \param s_to_e sphere to ellipsoid trafo
1931 : !> \return ...
1932 : ! **************************************************************************************************
1933 19854720 : PURE FUNCTION ellipsoid_bounds(s_rad, s_to_e)
1934 : REAL(KIND=dp), INTENT(IN) :: s_rad
1935 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: s_to_e
1936 : REAL(KIND=dp), DIMENSION(3) :: ellipsoid_bounds
1937 :
1938 : INTEGER :: i_xyz
1939 :
1940 79418880 : DO i_xyz = 1, 3
1941 79418880 : ellipsoid_bounds(i_xyz) = SQRT(s_to_e(i_xyz, 1)**2 + s_to_e(i_xyz, 2)**2 + s_to_e(i_xyz, 3)**2)*s_rad
1942 : END DO
1943 :
1944 : END FUNCTION ellipsoid_bounds
1945 :
1946 : ! **************************************************************************************************
1947 : !> \brief Roughly estimated number of floating point operations
1948 : !> \param l ...
1949 : !> \param m ...
1950 : !> \param H_or_C_product ...
1951 : !> \return ...
1952 : ! **************************************************************************************************
1953 297000 : PURE FUNCTION nsum_gaussian_overlap(l, m, H_or_C_product)
1954 : INTEGER, INTENT(IN) :: l, m, H_or_C_product
1955 : INTEGER :: nsum_gaussian_overlap
1956 :
1957 : INTEGER :: loop
1958 :
1959 297000 : nsum_gaussian_overlap = exp_w
1960 297000 : loop = (m + 1)*(l + 1)*(m + l + 2)
1961 297000 : IF (H_or_C_product == 1) THEN
1962 148500 : nsum_gaussian_overlap = nsum_gaussian_overlap + loop*16
1963 : ELSE
1964 148500 : nsum_gaussian_overlap = nsum_gaussian_overlap + loop*32
1965 : END IF
1966 297000 : END FUNCTION nsum_gaussian_overlap
1967 :
1968 : ! **************************************************************************************************
1969 : !> \brief ...
1970 : !> \param lco ...
1971 : !> \param l ...
1972 : !> \param lx ...
1973 : !> \param ly ...
1974 : !> \param lz ...
1975 : ! **************************************************************************************************
1976 843147983 : PURE ELEMENTAL SUBROUTINE get_l(lco, l, lx, ly, lz)
1977 : INTEGER, INTENT(IN) :: lco
1978 : INTEGER, INTENT(OUT) :: l
1979 : INTEGER, INTENT(OUT), OPTIONAL :: lx, ly, lz
1980 :
1981 3372591932 : l = SUM(indco(:, lco))
1982 843147983 : IF (PRESENT(lx)) lx = indco(1, lco)
1983 843147983 : IF (PRESENT(ly)) ly = indco(2, lco)
1984 843147983 : IF (PRESENT(lz)) lz = indco(3, lco)
1985 843147983 : END SUBROUTINE get_l
1986 :
1987 : ! **************************************************************************************************
1988 : !> \brief ...
1989 : !> \param i ...
1990 : !> \return ...
1991 : ! **************************************************************************************************
1992 479950601 : PURE ELEMENTAL FUNCTION i_pow(i)
1993 : INTEGER, INTENT(IN) :: i
1994 : COMPLEX(KIND=dp) :: i_pow
1995 :
1996 : COMPLEX(KIND=dp), DIMENSION(0:3), PARAMETER :: &
1997 : ip = [(1.0_dp, 0.0_dp), (0.0_dp, 1.0_dp), (-1.0_dp, 0.0_dp), (0.0_dp, -1.0_dp)]
1998 :
1999 479950601 : i_pow = ip(MOD(i, 4))
2000 :
2001 479950601 : END FUNCTION i_pow
2002 :
2003 3955399 : END MODULE eri_mme_lattice_summation
|
