Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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
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 2404794 : 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 2404794 : alpha = 0.25_dp/zeta
123 2404794 : beta = 0.25_dp/zetb
124 2404794 : 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 2404794 : IF (PRESENT(G_rads_1)) THEN
129 2404794 : G_res = 0.5_dp*G_min
130 2404794 : G_rads_1(1) = exp_radius(la_max, alpha, sum_precision, 1.0_dp, G_res)
131 2404794 : G_rads_1(2) = exp_radius(lb_max, beta, sum_precision, 1.0_dp, G_res)
132 2404794 : 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 2404794 : IF (PRESENT(R_rads_2)) THEN
137 2250540 : R_res = 0.5_dp*R_min
138 2250540 : R_rads_2(1) = exp_radius(lb_max + lc_max, 0.25_dp/(beta + gamma), sum_precision, 1.0_dp, epsabs=R_res)
139 2250540 : R_rads_2(2) = exp_radius(lc_max + la_max, 0.25_dp/(alpha + gamma), sum_precision, 1.0_dp, epsabs=R_res)
140 2250540 : 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 2404794 : IF (PRESENT(R_rads_3)) THEN
146 2250540 : R_res = 0.5_dp*R_min
147 2250540 : alpha_R = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm)
148 2250540 : R_rads_3(1) = exp_radius(la_max + lb_max, zeta*zetb/(zeta + zetb), sum_precision, 1.0_dp, R_res)
149 2250540 : 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 2404794 : END SUBROUTINE
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
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 2250540 : 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 2250540 : G_rads_1=G_rads_1, R_rads_2=R_rads_2, R_rads_3=R_rads_3)
265 :
266 2250540 : vol_G = twopi**3/vol
267 :
268 9002160 : order_1 = MAXLOC(G_rads_1, DIM=1)
269 9002160 : order_2 = MINLOC(G_rads_1, DIM=1)
270 :
271 9002160 : DO i = 1, 3
272 110276460 : G_bounds_1(i, :) = ellipsoid_bounds(G_rads_1(i), TRANSPOSE(hmat)/(2.0_dp*pi))
273 : END DO
274 :
275 9002160 : DO i = 1, 3
276 29257020 : R_bounds_2(i, :) = ellipsoid_bounds(R_rads_2(i), h_inv)
277 : END DO
278 :
279 6751620 : DO i = 1, 2
280 20254860 : R_bounds_3(i, :) = ellipsoid_bounds(R_rads_3(i), h_inv)
281 : END DO
282 :
283 2250540 : IF (is_ortho) THEN
284 8408160 : DO i_xyz = 1, 3
285 :
286 6306120 : ns3_R1 = 2.0_dp*R_bounds_3(1, i_xyz)
287 6306120 : ns3_R2 = 2.0_dp*R_bounds_3(2, i_xyz)
288 :
289 8408160 : 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 2250540 : END SUBROUTINE
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
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 1025491 : 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
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 .GT. G_rads_sq) CYCLE
449 :
450 6396767 : IF (PRESENT(potential)) THEN
451 1033920 : IF (gx .EQ. 0 .AND. gy .EQ. 0 .AND. gz .EQ. 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
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 : 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 4425368 : 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
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 .GT. 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 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 ! cost: 4 flops
615 : END DO
616 : END DO
617 : END DO
618 : END DO
619 :
620 : CALL create_hermite_to_cartesian(alpha, l_max, h_to_c)
621 :
622 2057776 : DO lco = 1, ncoset(l_max)
623 1929557 : CALL get_l(lco, l, lx, ly, lz)
624 :
625 5692372 : DO llx = 0, lx
626 11963703 : DO lly = 0, ly
627 20729884 : DO llz = 0, lz
628 : S_R(coset(lx, ly, lz)) = S_R(coset(lx, ly, lz)) + & ! cost: 4 flops
629 : h_to_c(llx, lx)*h_to_c(lly, ly)*h_to_c(llz, lz)* &
630 17095288 : S_R_C(coset(llx, lly, llz))
631 : END DO
632 : END DO
633 : END DO
634 : END DO
635 2057776 : S_R(:) = S_R(:)*(alpha/pi)**1.5_dp
636 :
637 256438 : END SUBROUTINE pgf_sum_2c_rspace_3d
638 :
639 : ! **************************************************************************************************
640 : !> \brief Compute 1d sum
641 : !> S_G(l, alpha) = inv_lgth*sum_G( C(l, alpha, delta_l, G) ) with
642 : !> C(l, alpha, delta_l, G) = prefactor*|G|^(l-delta_l) exp(-alpha*G^2)
643 : !> if G not equal 0
644 : !> C(l = 0, alpha, delta_l, 0) = 1, C(l>0, alpha, delta_l, 0) = 0
645 : !> dG = inv_lgth*twopi and G = -G_bound*dG, (-G_bound + 1)*dG, ..., G_bound*dG
646 : !> for all l < = l_max.
647 : !> \param S_G ...
648 : !> \param alpha ...
649 : !> \param inv_lgth ...
650 : !> \param G_min ...
651 : !> \param G_c ...
652 : !> \param delta_l ...
653 : !> \param prefactor ...
654 : !> \note needed for cutoff error estimate
655 : ! **************************************************************************************************
656 295044 : PURE SUBROUTINE pgf_sum_2c_gspace_1d_deltal(S_G, alpha, inv_lgth, G_min, G_c, delta_l, prefactor)
657 : REAL(KIND=dp), DIMENSION(0:), INTENT(OUT) :: S_G
658 : REAL(KIND=dp), INTENT(IN) :: alpha, inv_lgth
659 : INTEGER, INTENT(IN) :: G_min, G_c
660 : REAL(KIND=dp), INTENT(IN) :: delta_l, prefactor
661 :
662 : INTEGER :: k, k0, k1, l, l_max
663 : REAL(KIND=dp) :: dG, exp_tot, G, prefac
664 :
665 295044 : prefac = prefactor*inv_lgth
666 295044 : dG = inv_lgth*twopi ! positive
667 295044 : l_max = UBOUND(S_G, 1)
668 :
669 1660104 : S_G(:) = 0.0_dp
670 295044 : IF (0 .LT. G_min) THEN
671 : k0 = G_min; k1 = 0
672 73761 : ELSE IF (G_c .LT. 0) THEN
673 : k0 = 0; k1 = G_c
674 : ELSE ! Gmin <= 0 /\ 0 <= Gc
675 73761 : S_G(0) = prefac
676 : k0 = 1; k1 = -1
677 : END IF
678 : ! positive range
679 : IF (0 .LT. k0) THEN
680 14185303 : DO k = k0, G_c
681 13890259 : G = k*dG; exp_tot = EXP(-alpha*G**2)*prefac
682 83043248 : DO l = 0, l_max
683 82748204 : S_G(l) = S_G(l) + G**(l - delta_l)*exp_tot
684 : END DO
685 : END DO
686 : END IF
687 : ! negative range
688 295044 : IF (k1 .LT. 0) THEN
689 4816941 : DO k = G_min, k1
690 4743180 : G = k*dG; exp_tot = EXP(-alpha*G**2)*prefac
691 27823869 : DO l = 0, l_max
692 27750108 : S_G(l) = S_G(l) + ABS(G)**(l - delta_l)*exp_tot
693 : END DO
694 : END DO
695 : END IF
696 295044 : END SUBROUTINE pgf_sum_2c_gspace_1d_deltal
697 :
698 : ! **************************************************************************************************
699 : !> \brief Compute Ewald-like sum for 3-center integrals in 1 dimension
700 : !> S_G(l, m, n, alpha, beta, gamma) = i^(l+m+n)*(-1)^(l+m)*inv_lgth^2*
701 : !> sum_G sum_G'( exp(i G R1)
702 : !> C(l,alpha,G) C(m,beta,G'-G) C(n,gamma,G') exp(i G' R2) )
703 : !> for all l < = l_max, m <= m_max, n <= n_max.
704 : !> a_mm is the minimax exponent.
705 : !> alpha = 1/(4 zeta), beta = 1/(4 zetb), gamma = 1/(4 zetc) + a_mm
706 : !> R1 = RB-RA; R2 = RC-RB
707 : !> Note on method / order arguments:
708 : !> Three equivalent methods (Poisson summation) to compute this sum over
709 : !> Cartesian Gaussians C or Hermite Gaussians H and
710 : !> reciprocal lattice vectors G or direct lattice vectors R:
711 : !> - method 1: sum_G sum_G' C(G) C(G,G') C(G')
712 : !> - method 2: sum_G sum_R C(G) C(R)
713 : !> - method 3: sum_R sum_R' H(R, R')
714 : !> The order parameter selects the Gaussian functions over which the sum is performed
715 : !> method 1: order = 1, 2, 3
716 : !> method 2: order = 1, 2, 3
717 : !> method 3: order = 1
718 : !> If method and order are not present, the method / order that converges fastest is
719 : !> automatically chosen.
720 : !> \param S_G ...
721 : !> \param RA ...
722 : !> \param RB ...
723 : !> \param RC ...
724 : !> \param zeta ...
725 : !> \param zetb ...
726 : !> \param zetc ...
727 : !> \param a_mm ...
728 : !> \param lgth ...
729 : !> \param G_bounds_1 ...
730 : !> \param R_bounds_2 ...
731 : !> \param R_bounds_3 ...
732 : !> \param method ...
733 : !> \param method_out ...
734 : !> \param order ...
735 : ! **************************************************************************************************
736 6306120 : SUBROUTINE pgf_sum_3c_1d(S_G, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_bounds_3)
737 : REAL(KIND=dp), DIMENSION(0:, 0:, 0:), INTENT(OUT) :: S_G
738 : REAL(KIND=dp), INTENT(IN) :: RA, RB, RC, zeta, zetb, zetc, a_mm, lgth
739 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: R_bounds_3
740 :
741 : INTEGER :: l_max, m_max, n_max
742 : INTEGER :: nR1, nR2
743 : INTEGER :: prop_exp
744 :
745 6306120 : l_max = UBOUND(S_G, 1)
746 6306120 : m_max = UBOUND(S_G, 2)
747 6306120 : n_max = UBOUND(S_G, 3)
748 :
749 6306120 : nR1 = 2*FLOOR(R_bounds_3(1)) + 1
750 6306120 : nR2 = 2*FLOOR(R_bounds_3(2)) + 1
751 :
752 6306120 : IF (nR1*nR2 > 1 + nR1*2) THEN
753 : prop_exp = 1
754 : ELSE
755 2986110 : prop_exp = 0
756 : END IF
757 :
758 25224480 : IF (MAXVAL([l_max, m_max, n_max]) .GT. ${lmax_unroll}$) THEN
759 92160 : CALL pgf_sum_3c_rspace_1d_generic(S_G, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_bounds_3)
760 : ELSE
761 : #:for l_max in range(0,lmax_unroll+1)
762 6213960 : IF (l_max == ${l_max}$) THEN
763 : #:for m_max in range(0,lmax_unroll+1)
764 6213960 : IF (m_max == ${m_max}$) THEN
765 : #:for n_max in range(0, lmax_unroll+1)
766 6213960 : IF (n_max == ${n_max}$) THEN
767 : #:for prop_exp in range(0,2)
768 6213960 : IF (prop_exp == ${prop_exp}$) THEN
769 : CALL pgf_sum_3c_rspace_1d_${l_max}$_${m_max}$_${n_max}$_exp_${prop_exp}$ (S_G, RA, RB, RC, &
770 6213960 : zeta, zetb, zetc, a_mm, lgth, R_bounds_3)
771 6213960 : RETURN
772 : END IF
773 : #:endfor
774 : END IF
775 : #:endfor
776 : END IF
777 : #:endfor
778 : END IF
779 : #:endfor
780 : END IF
781 :
782 : END SUBROUTINE pgf_sum_3c_1d
783 :
784 : ! **************************************************************************************************
785 : !> \brief Roughly estimated number of floating point operations
786 : !> \param ns_G1 ...
787 : !> \param ns_G2 ...
788 : !> \param l ...
789 : !> \param m ...
790 : !> \param n ...
791 : !> \return ...
792 : ! **************************************************************************************************
793 0 : PURE FUNCTION nsum_3c_gspace_1d()
794 : INTEGER(KIND=int_8) :: nsum_3c_gspace_1d
795 :
796 0 : nsum_3c_gspace_1d = 15
797 0 : END FUNCTION
798 :
799 : ! **************************************************************************************************
800 : !> \brief Roughly estimated number of floating point operations
801 : !> \param ns_G ...
802 : !> \param ns_R ...
803 : !> \param l ...
804 : !> \param m ...
805 : !> \param n ...
806 : !> \return ...
807 : ! **************************************************************************************************
808 0 : PURE FUNCTION nsum_product_3c_gspace_1d(ns_G, ns_R)
809 : REAL(KIND=dp), INTENT(IN) :: ns_G, ns_R
810 : INTEGER(KIND=int_8) :: nsum_product_3c_gspace_1d
811 :
812 0 : nsum_product_3c_gspace_1d = MIN(19, NINT(ns_G*(3 + ns_R*2)))
813 0 : END FUNCTION
814 :
815 : ! **************************************************************************************************
816 : !> \brief Roughly estimated number of floating point operations
817 : !> \param ns_R1 ...
818 : !> \param ns_R2 ...
819 : !> \param l ...
820 : !> \param m ...
821 : !> \param n ...
822 : !> \return ...
823 : ! **************************************************************************************************
824 6306120 : PURE FUNCTION nsum_3c_rspace_1d(ns_R1, ns_R2)
825 : REAL(KIND=dp), INTENT(IN) :: ns_R1, ns_R2
826 : INTEGER(KIND=int_8) :: nsum_3c_rspace_1d
827 :
828 6306120 : nsum_3c_rspace_1d = NINT(MIN((4 + ns_R1*2), ns_R1*(ns_R2 + 1)), KIND=int_8)
829 6306120 : END FUNCTION
830 :
831 : ! **************************************************************************************************
832 : !> \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)
833 : !> with alpha = 1.0_dp/((a + b + c)/((a + b)*c) + 4.0_dp*a_mm),
834 : !> P(R) = (a*(RA + R) + b*RB)/(a + b)
835 : !> \param S_R ...
836 : !> \param RA ...
837 : !> \param RB ...
838 : !> \param RC ...
839 : !> \param zeta ...
840 : !> \param zetb ...
841 : !> \param zetc ...
842 : !> \param a_mm ...
843 : !> \param lgth ...
844 : !> \param R_c ...
845 : ! **************************************************************************************************
846 92160 : PURE SUBROUTINE pgf_sum_3c_rspace_1d_generic(S_R, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_c)
847 : REAL(KIND=dp), DIMENSION(0:, 0:, 0:), INTENT(OUT) :: S_R
848 : REAL(KIND=dp), INTENT(IN) :: RA, RB, RC, zeta, zetb, zetc, a_mm, lgth
849 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: R_c
850 :
851 : INTEGER :: ll, mm, l, k, l_max, m, m_max, n, n_max, rr1, rr2, t, l_tot_max
852 : REAL(KIND=dp) :: alpha, dR, exp_tot, R, R1, R2, R_offset, &
853 : R_pow_t, R_tmp, rr1_delta, rr2_delta, c1, c2, c3
854 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: S_R_t
855 92160 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: h_to_c
856 92160 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: E
857 :
858 92160 : dR = lgth
859 92160 : alpha = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm)
860 92160 : l_max = UBOUND(S_R, 1)
861 92160 : m_max = UBOUND(S_R, 2)
862 92160 : n_max = UBOUND(S_R, 3)
863 92160 : l_tot_max = l_max + m_max + n_max
864 :
865 276480 : ALLOCATE (S_R_t(0:l_max + m_max + n_max))
866 460800 : ALLOCATE (E(-1:l_max + m_max + 1, -1:l_max, -1:m_max))
867 :
868 2903040 : S_R(:, :, :) = 0.0_dp
869 :
870 92160 : R_offset = RC - (zeta*RA + zetb*RB)/(zeta + zetb)
871 :
872 : ! inline CALL create_hermite_to_cartesian(alpha, l_tot_max, h_to_c)
873 368640 : ALLOCATE (h_to_c(-1:l_tot_max + 1, 0:l_tot_max))
874 8202240 : h_to_c(:, :) = 0.0_dp
875 92160 : h_to_c(0, 0) = 1.0_dp
876 737280 : DO l = 0, l_tot_max - 1
877 3962880 : DO k = 0, l + 1
878 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)
879 : END DO
880 : END DO
881 :
882 92160 : rr1_delta = (RA - RB)/dR
883 1428480 : DO rr1 = CEILING(-R_c(1) + rr1_delta), FLOOR(R_c(1) + rr1_delta)
884 12026880 : S_R_t(:) = 0.0_dp
885 1336320 : R1 = rr1*dR
886 1336320 : R_tmp = R_offset + R1*zeta/(zeta + zetb)
887 1336320 : rr2_delta = -R_tmp/dR
888 15398400 : DO rr2 = CEILING(-R_c(2) + rr2_delta), FLOOR(R_c(2) + rr2_delta)
889 14062080 : R2 = rr2*dR
890 14062080 : R = R_tmp + R2
891 14062080 : exp_tot = EXP(-alpha*R**2) ! cost: exp_w flops
892 14062080 : R_pow_t = 1.0_dp
893 127895040 : DO t = 0, l_max + m_max + n_max
894 112496640 : S_R_t(t) = S_R_t(t) + R_pow_t*exp_tot ! cost: 2 flops
895 126558720 : R_pow_t = R_pow_t*R ! cost: 1 flop
896 : END DO
897 : END DO
898 :
899 : ! C --> H
900 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)
901 :
902 : ! H --> HH
903 : !inline CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, RA-R1, RB, 2, E)
904 :
905 183191040 : E(:, :, :) = 0.0_dp
906 1336320 : E(0, 0, 0) = EXP(-zeta*zetb/(zeta + zetb)*(RA - R1 - RB)**2)
907 :
908 1336320 : c1 = 1.0_dp/(zeta + zetb)
909 1336320 : c2 = 2.0_dp*(zetb/(zeta + zetb))*(RB - (RA - R1))
910 1336320 : c3 = 2.0_dp*(zeta/(zeta + zetb))*(RA - R1 - RB)
911 :
912 7948800 : DO mm = 0, m_max - 1
913 8179200 : DO ll = 0, l_max - 1
914 7879680 : DO t = 0, ll + mm + 1
915 : E(t, ll + 1, mm) = zeta*(c1*E(t - 1, ll, mm) + &
916 : c2*E(t, ll, mm) + &
917 : 2*(t + 1)*E(t + 1, ll, mm) - &
918 1036800 : 2*ll*E(t, ll - 1, mm))
919 : E(t, ll, mm + 1) = zetb*(c1*E(t - 1, ll, mm) + &
920 : c3*E(t, ll, mm) + &
921 : 2*(t + 1)*E(t + 1, ll, mm) - &
922 1267200 : 2*mm*E(t, ll, mm - 1))
923 : END DO
924 : END DO
925 : END DO
926 :
927 1451520 : DO ll = 0, l_max - 1
928 2142720 : DO t = 0, ll + m_max + 1
929 : E(t, ll + 1, m_max) = zeta*(c1*E(t - 1, ll, m_max) + &
930 : c2*E(t, ll, m_max) + &
931 : 2*(t + 1)*E(t + 1, ll, m_max) - &
932 806400 : 2*ll*E(t, ll - 1, m_max))
933 : END DO
934 : END DO
935 :
936 7948800 : DO mm = 0, m_max - 1
937 34560000 : DO t = 0, l_max + mm + 1
938 : E(t, l_max, mm + 1) = zetb*(c1*E(t - 1, l_max, mm) + &
939 : c3*E(t, l_max, mm) + &
940 : 2*(t + 1)*E(t + 1, l_max, mm) - &
941 33223680 : 2*mm*E(t, l_max, mm - 1))
942 : END DO
943 : END DO
944 :
945 5391360 : DO n = 0, n_max
946 29007360 : DO m = 0, m_max
947 51724800 : DO l = 0, l_max
948 132364800 : DO t = 0, l + m
949 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
950 : END DO
951 : END DO
952 : END DO
953 : END DO
954 : END DO
955 :
956 2903040 : S_R = S_R*pi**(-0.5_dp)*((zeta + zetb)/(zeta*zetb))**(-0.5_dp)
957 92160 : END SUBROUTINE
958 :
959 : ! **************************************************************************************************
960 : !> \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)
961 : !> with alpha = 1.0_dp/((a + b + c)/((a + b)*c) + 4.0_dp*a_mm),
962 : !> P(R) = (a*(RA + R) + b*RB)/(a + b)
963 : !> \param S_R ...
964 : !> \param RA ...
965 : !> \param RB ...
966 : !> \param RC ...
967 : !> \param zeta ...
968 : !> \param zetb ...
969 : !> \param zetc ...
970 : !> \param a_mm ...
971 : !> \param lgth ...
972 : !> \param R_c ...
973 : ! **************************************************************************************************
974 : #:for prop_exp in range(0,2)
975 : #:for l_max in range(0, lmax_unroll+1)
976 : #:for m_max in range(0, lmax_unroll+1)
977 : #:for n_max in range(0, lmax_unroll+1)
978 : #:set l_tot_max = l_max + m_max + n_max
979 6213960 : PURE SUBROUTINE pgf_sum_3c_rspace_1d_${l_max}$_${m_max}$_${n_max}$_exp_${prop_exp}$ ( &
980 6213960 : S_R, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_c)
981 : REAL(KIND=dp), DIMENSION(0:, 0:, 0:), INTENT(OUT) :: S_R
982 : REAL(KIND=dp), INTENT(IN) :: RA, RB, RC, zeta, zetb, zetc, a_mm, lgth
983 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: R_c
984 : INTEGER :: rr1, rr2, rr2_l, rr2_r, rr1_l, rr1_r
985 : REAL(KIND=dp) :: alpha, alpha_E, dR, exp2_Rsq, R, R1, R_offset, &
986 : R_pow_t, R_tmp, rr1_delta, rr2_delta
987 :
988 : #:if l_tot_max > 0
989 : REAL(KIND=dp) :: c1, c2, c3
990 : #:endif
991 : REAL(KIND=dp) :: ${", ".join(["S_R_t_{}".format(t) for t in range(0,l_tot_max+1)])}$
992 : REAL(KIND=dp) :: ${", ".join(["S_R_t2_{}".format(t) for t in range(0,l_tot_max+1)])}$
993 : 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)])}$
994 : 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)])}$
995 :
996 : #:if prop_exp
997 : REAL(KIND=dp) :: exp2_2RdR, exp_dRsq, exp_2dRsq !exp_E_dRsq, exp_E_2dRsq, exp_E_2RdR, exp_E_Rsq
998 : #:endif
999 :
1000 6213960 : dR = lgth
1001 6213960 : alpha = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm)
1002 :
1003 112951872 : S_R(:, :, :) = 0.0_dp
1004 :
1005 6213960 : R_offset = RC - (zeta*RA + zetb*RB)/(zeta + zetb)
1006 :
1007 6213960 : h_to_c_0_0 = SQRT(alpha/pi)
1008 :
1009 : #:for l in range(0, l_tot_max)
1010 : #:for k in range(0, l+2)
1011 : #:if k<l or k>0
1012 5491656 : 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}#
1013 : #:else
1014 5491656 : h_to_c_${k}$_${l+1}$ = 0.0_dp
1015 : #:endif
1016 : #:endfor
1017 : #:endfor
1018 :
1019 : #:if prop_exp
1020 3261642 : exp_dRsq = exp(-alpha*dR*dR)
1021 3261642 : exp_2dRsq = exp_dRsq*exp_dRsq
1022 : #:endif
1023 :
1024 6213960 : rr1_delta = (RA - RB)/dR
1025 :
1026 6213960 : rr1_l = CEILING(-R_c(1) + rr1_delta)
1027 6213960 : rr1_r = FLOOR(R_c(1) + rr1_delta)
1028 :
1029 6213960 : R1 = rr1_l*dR
1030 :
1031 6213960 : alpha_E = zeta*zetb/(zeta + zetb)
1032 :
1033 20182392 : DO rr1 = rr1_l, rr1_r
1034 : #:for t in range(0, l_tot_max + 1)
1035 13968432 : S_R_t_${t}$ = 0.0_dp
1036 13968432 : S_R_t2_${t}$ = 0.0_dp
1037 : #:endfor
1038 13968432 : R_tmp = R_offset + R1*zeta/(zeta + zetb)
1039 13968432 : rr2_delta = -R_tmp/dR
1040 :
1041 13968432 : rr2_l = CEILING(-R_c(2) + rr2_delta)
1042 13968432 : rr2_r = FLOOR(R_c(2) + rr2_delta)
1043 :
1044 13968432 : R = R_tmp + (rr2_l)*dR
1045 :
1046 : #:if prop_exp
1047 9106454 : exp2_2RdR = exp(-2*alpha*R*dR)
1048 9106454 : exp2_Rsq = exp(-alpha*R*R)
1049 : #:endif
1050 :
1051 71456359 : DO rr2 = rr2_l, rr2_r
1052 57487927 : R_pow_t = 1.0_dp
1053 : #:if not prop_exp
1054 7156859 : exp2_Rsq = exp(-alpha*R*R)
1055 : #:endif
1056 : #:for t in range(0, l_tot_max + 1)
1057 57487927 : S_R_t_${t}$ = S_R_t_${t}$+R_pow_t*exp2_Rsq
1058 : #:if t<l_tot_max
1059 51357089 : R_pow_t = R_pow_t*R
1060 : #:endif
1061 : #:endfor
1062 :
1063 : #:if prop_exp
1064 50331068 : exp2_Rsq = exp2_Rsq*exp_dRsq*exp2_2RdR
1065 50331068 : exp2_2RdR = exp2_2RdR*exp_2dRsq
1066 : #:endif
1067 71456359 : R = R + dR
1068 : END DO
1069 :
1070 : ! C --> H
1071 : #:for l in range(0, l_tot_max+1)
1072 : #:for k in range(0, l+1)
1073 13968432 : S_R_t2_${l}$ = S_R_t2_${l}$+h_to_c_${k}$_${l}$*S_R_t_${k}$
1074 : #:endfor
1075 : #:endfor
1076 :
1077 : ! H --> HH
1078 13968432 : E_0_0_0 = exp(-alpha_E*(RA - RB - R1)*(RA - RB - R1))
1079 :
1080 : #:if l_tot_max > 0
1081 12351072 : c1 = 1.0_dp/(zeta + zetb)
1082 12351072 : c2 = 2.0_dp*(zetb/(zeta + zetb))*(RB - (RA - R1))
1083 12351072 : c3 = 2.0_dp*(zeta/(zeta + zetb))*(RA - R1 - RB)
1084 : #:endif
1085 :
1086 : #:for m in range(0,m_max+1)
1087 : #:for l in range(0,l_max+1)
1088 : #:for t in range(0,l+m+2)
1089 : #:if l < l_max
1090 : E_${t}$_${l+1}$_${m}$ = zeta*(#{if t>0}# c1*E_${t-1}$_${l}$_${m}$#{endif}# &
1091 : #{if t<=l+m}# +c2*E_${t}$_${l}$_${m}$&#{endif}#
1092 : #{if t<l+m}# +${2*(t+1)}$*E_${t+1}$_${l}$_${m}$ &#{endif}#
1093 8746776 : #{if l>0 and t<=l-1+m}#-${2*l}$*E_${t}$_${l-1}$_${m}$#{endif}#)
1094 : #:endif
1095 : #:if m < m_max
1096 : E_${t}$_${l}$_${m+1}$ = zetb*(#{if t>0}# c1*E_${t-1}$_${l}$_${m}$#{endif}# &
1097 : #{if t<=l+m}#+c3*E_${t}$_${l}$_${m}$&#{endif}#
1098 : #{if t<l+m}# +${2*(t+1)}$*E_${t+1}$_${l}$_${m}$ &#{endif}#
1099 8747976 : #{if m>0 and t<=m-1+l}#-${2*m}$*E_${t}$_${l}$_${m-1}$#{endif}#)
1100 : #:endif
1101 : #:endfor
1102 : #:endfor
1103 : #:endfor
1104 :
1105 : #:for n in range(0, n_max+1)
1106 : #:for m in range(0, m_max+1)
1107 : #:for l in range(0, l_max+1)
1108 : #:for t in range(0, l+m+1)
1109 13968432 : S_R(${l}$, ${m}$, ${n}$) = S_R(${l}$, ${m}$, ${n}$) + E_${t}$_${l}$_${m}$*(${(-1)**n}$)*S_R_t2_${t+n}$ ! cost: 5 flops
1110 : #:endfor
1111 : #:endfor
1112 : #:endfor
1113 : #:endfor
1114 20182392 : R1 = R1 + dR
1115 : END DO
1116 :
1117 112951872 : S_R = S_R*pi**(-0.5_dp)*((zeta + zetb)/(zeta*zetb))**(-0.5_dp)
1118 6213960 : END SUBROUTINE
1119 : #:endfor
1120 : #:endfor
1121 : #:endfor
1122 : #:endfor
1123 :
1124 : ! **************************************************************************************************
1125 : !> \brief As pgf_sum_3c_1d but 3d sum required for non-orthorhombic cells
1126 : !> \param S_G ...
1127 : !> \param la_max ...
1128 : !> \param lb_max ...
1129 : !> \param lc_max ...
1130 : !> \param RA ...
1131 : !> \param RB ...
1132 : !> \param RC ...
1133 : !> \param zeta ...
1134 : !> \param zetb ...
1135 : !> \param zetc ...
1136 : !> \param a_mm ...
1137 : !> \param hmat ...
1138 : !> \param h_inv ...
1139 : !> \param vol ...
1140 : !> \param G_bounds_1 ...
1141 : !> \param R_bounds_2 ...
1142 : !> \param R_bounds_3 ...
1143 : !> \param G_rads_1 ...
1144 : !> \param R_rads_2 ...
1145 : !> \param R_rads_3 ...
1146 : !> \param method ...
1147 : !> \param method_out ...
1148 : !> \param order ...
1149 : ! **************************************************************************************************
1150 148500 : SUBROUTINE pgf_sum_3c_3d(S_G, la_max, lb_max, lc_max, RA, RB, RC, &
1151 : zeta, zetb, zetc, a_mm, hmat, h_inv, vol, &
1152 : G_bounds_1, R_bounds_2, R_bounds_3, &
1153 : G_rads_1, R_rads_2, R_rads_3, &
1154 : method, method_out, order)
1155 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: S_G
1156 : INTEGER, INTENT(IN) :: la_max, lb_max, lc_max
1157 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: RA, RB, RC
1158 : REAL(KIND=dp), INTENT(IN) :: zeta, zetb, zetc, a_mm
1159 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv
1160 : REAL(KIND=dp), INTENT(IN) :: vol
1161 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: G_bounds_1, R_bounds_2
1162 : REAL(KIND=dp), DIMENSION(2, 3), INTENT(IN) :: R_bounds_3
1163 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G_rads_1, R_rads_2
1164 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: R_rads_3
1165 : INTEGER, INTENT(IN) :: method
1166 : INTEGER, INTENT(OUT), OPTIONAL :: method_out
1167 : INTEGER, INTENT(IN), OPTIONAL :: order
1168 :
1169 : INTEGER :: l_max, m_max, n_max, sum_method, &
1170 : sum_order
1171 : LOGICAL :: assert_stm
1172 : REAL(KIND=dp) :: alpha, beta, G_rad, gamma, R_rad
1173 148500 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: S_G_tmp
1174 : REAL(KIND=dp), DIMENSION(3) :: G_bound, R1, R2, R_bound
1175 : REAL(KIND=dp), DIMENSION(3, 3) :: ht
1176 :
1177 148500 : IF (PRESENT(order)) THEN
1178 0 : sum_order = order
1179 : ELSE
1180 148500 : sum_order = 0
1181 : END IF
1182 :
1183 148500 : sum_method = method
1184 :
1185 148500 : alpha = 0.25_dp/zeta
1186 148500 : beta = 0.25_dp/zetb
1187 148500 : gamma = 0.25_dp/zetc + a_mm
1188 :
1189 148500 : l_max = la_max
1190 148500 : m_max = lb_max
1191 148500 : n_max = lc_max
1192 :
1193 594000 : R1 = RB - RA
1194 594000 : R2 = RC - RB
1195 :
1196 1930500 : ht = twopi*TRANSPOSE(h_inv)
1197 :
1198 0 : SELECT CASE (sum_method)
1199 : CASE (1) ! sum_G sum_G' C(G) C(G,G') C(G')
1200 :
1201 0 : IF (sum_order .EQ. 0) THEN
1202 0 : sum_order = MAXLOC(G_bounds_1(:, 1), DIM=1)
1203 : assert_stm = MINLOC(G_bounds_1(:, 1), DIM=1) == &
1204 : MINLOC(G_bounds_1(:, 2), DIM=1) .AND. &
1205 : MINLOC(G_bounds_1(:, 1), DIM=1) == &
1206 0 : MINLOC(G_bounds_1(:, 3), DIM=1)
1207 0 : CPASSERT(assert_stm)
1208 : END IF
1209 :
1210 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)
1211 :
1212 : CASE (2) ! sum_G sum_R C(G) C(R)
1213 5846 : IF (sum_order .EQ. 0) THEN
1214 23384 : sum_order = MINLOC(G_bounds_1(:, 1), DIM=1)
1215 : assert_stm = MINLOC(G_bounds_1(:, 1), DIM=1) == &
1216 : MINLOC(G_bounds_1(:, 2), DIM=1) .AND. &
1217 : MINLOC(G_bounds_1(:, 1), DIM=1) == &
1218 93536 : MINLOC(G_bounds_1(:, 3), DIM=1)
1219 0 : CPASSERT(assert_stm)
1220 : END IF
1221 5846 : R_rad = R_rads_2(sum_order)
1222 5846 : G_rad = G_rads_1(sum_order)
1223 23384 : R_bound = R_bounds_2(sum_order, :)
1224 23384 : G_bound = G_bounds_1(sum_order, :)
1225 142654 : SELECT CASE (sum_order)
1226 : CASE (1)
1227 0 : ALLOCATE (S_G_tmp(ncoset(l_max), ncoset(m_max), ncoset(n_max)))
1228 : 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, &
1229 0 : R_bound, G_bound, R_rad, G_rad)
1230 0 : S_G = RESHAPE(S_G_tmp, SHAPE(S_G), order=[1, 2, 3])
1231 : CASE (2)
1232 0 : ALLOCATE (S_G_tmp(ncoset(m_max), ncoset(l_max), ncoset(n_max)))
1233 : 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, &
1234 0 : R_bound, G_bound, R_rad, G_rad)
1235 0 : S_G = RESHAPE(S_G_tmp, SHAPE(S_G), order=[2, 1, 3])
1236 : CASE (3)
1237 29230 : ALLOCATE (S_G_tmp(ncoset(n_max), ncoset(m_max), ncoset(l_max)))
1238 : 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, &
1239 40922 : R_bound, G_bound, R_rad, G_rad)
1240 23384 : S_G = RESHAPE(S_G_tmp, SHAPE(S_G), order=[3, 2, 1])
1241 : END SELECT
1242 : CASE (3) ! sum_R sum_R' H(R, R')
1243 : 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, &
1244 142654 : R_bounds_3, R_rads_3)
1245 24412738 : S_G = S_G*pi**(-1.5_dp)*((zeta + zetb)/(zeta*zetb))**(-1.5_dp)
1246 : END SELECT
1247 :
1248 148500 : IF (PRESENT(method_out)) method_out = sum_method
1249 :
1250 148500 : END SUBROUTINE pgf_sum_3c_3d
1251 :
1252 : ! **************************************************************************************************
1253 : !> \brief Roughly estimated number of floating point operations
1254 : !> \param ns_G1 ...
1255 : !> \param ns_G2 ...
1256 : !> \param l ...
1257 : !> \param m ...
1258 : !> \param n ...
1259 : !> \return ...
1260 : ! **************************************************************************************************
1261 148500 : PURE FUNCTION nsum_3c_gspace_3d(ns_G1, ns_G2, l, m, n)
1262 : REAL(KIND=dp), INTENT(IN) :: ns_G1, ns_G2
1263 : INTEGER, INTENT(IN) :: l, m, n
1264 : INTEGER(KIND=int_8) :: nsum_3c_gspace_3d
1265 :
1266 148500 : nsum_3c_gspace_3d = NINT(ns_G1*ns_G2*(5*exp_w + ncoset(l)*ncoset(m)*ncoset(n)*4), KIND=int_8)
1267 :
1268 148500 : END FUNCTION
1269 :
1270 : ! **************************************************************************************************
1271 : !> \brief ...
1272 : !> \param S_G ...
1273 : !> \param l_max ...
1274 : !> \param m_max ...
1275 : !> \param n_max ...
1276 : !> \param R1 ...
1277 : !> \param R2 ...
1278 : !> \param alpha ...
1279 : !> \param beta ...
1280 : !> \param gamma ...
1281 : !> \param ht ...
1282 : !> \param vol ...
1283 : !> \param G_c ...
1284 : !> \param G_rad ...
1285 : !> \param sum_order ...
1286 : !> \param coulomb ...
1287 : ! **************************************************************************************************
1288 0 : PURE SUBROUTINE pgf_sum_3c_gspace_3d(S_G, l_max, m_max, n_max, R1, R2, alpha, beta, gamma, &
1289 : ht, vol, G_c, G_rad, sum_order, coulomb)
1290 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: S_G
1291 : INTEGER, INTENT(IN) :: l_max, m_max, n_max
1292 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R1, R2
1293 : REAL(KIND=dp), INTENT(IN) :: alpha, beta, gamma
1294 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: ht
1295 : REAL(KIND=dp), INTENT(IN) :: vol
1296 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: G_c
1297 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G_rad
1298 : INTEGER, INTENT(IN) :: sum_order
1299 : LOGICAL, INTENT(IN), OPTIONAL :: coulomb
1300 :
1301 : INTEGER, DIMENSION(3) :: G1c, G2c, G3c
1302 : INTEGER :: g1x, g1y, g1z, g2x, g2y, g2z, g3x, g3y, &
1303 : g3z
1304 : COMPLEX(KIND=dp), DIMENSION(ncoset(l_max), ncoset( &
1305 0 : m_max), ncoset(n_max)) :: S_G_c
1306 : LOGICAL :: use_coulomb
1307 : REAL(KIND=dp) :: G1_sq, G2_sq, G3_sq
1308 : REAL(KIND=dp), DIMENSION(3) :: G1, G1_x, G1_y, G1_z, G2, G2_x, G2_y, &
1309 : G2_z, G3, G3_x, G3_y, G3_z, G_rads_sq
1310 :
1311 0 : S_G_c(:, :, :) = 0.0_dp
1312 :
1313 0 : G1c = FLOOR(G_c(1, :))
1314 0 : G2c = FLOOR(G_c(2, :))
1315 0 : G3c = FLOOR(G_c(3, :))
1316 :
1317 : ! we can not precompute exponentials for general cell
1318 : ! Perform double G sum
1319 0 : G_rads_sq = G_rad**2
1320 :
1321 0 : IF (PRESENT(coulomb)) THEN
1322 0 : use_coulomb = coulomb
1323 : ELSE
1324 0 : use_coulomb = .FALSE.
1325 : END IF
1326 :
1327 0 : SELECT CASE (sum_order)
1328 : CASE (1)
1329 0 : DO g2x = -G2c(1), G2c(1)
1330 0 : G2_x = ht(:, 1)*g2x
1331 0 : DO g2y = -G2c(2), G2c(2)
1332 0 : G2_y = ht(:, 2)*g2y
1333 0 : DO g2z = -G2c(3), G2c(3)
1334 0 : G2_z = ht(:, 3)*g2z
1335 0 : G2 = G2_x + G2_y + G2_z
1336 0 : G2_sq = G2(1)**2 + G2(2)**2 + G2(3)**2
1337 0 : IF (G2_sq .GT. G_rads_sq(2)) CYCLE
1338 0 : DO g3x = -G3c(1), G3c(1)
1339 0 : G3_x = ht(:, 1)*g3x
1340 0 : DO g3y = -G3c(2), G3c(2)
1341 0 : G3_y = ht(:, 2)*g3y
1342 0 : DO g3z = -G3c(3), G3c(3)
1343 0 : G3_z = ht(:, 3)*g3z
1344 0 : G3 = G3_x + G3_y + G3_z
1345 0 : G3_sq = G3(1)**2 + G3(2)**2 + G3(3)**2
1346 0 : IF (G3_sq .GT. G_rads_sq(3)) CYCLE
1347 0 : G1 = G3 - G2
1348 0 : G1_sq = G1(1)**2 + G1(2)**2 + G1(3)**2
1349 0 : IF (G1_sq .GT. G_rads_sq(1)) CYCLE
1350 0 : IF (use_coulomb) THEN
1351 0 : IF (g3x == 0 .AND. g3y == 0 .AND. g3z == 0) CYCLE
1352 : END IF
1353 : CALL pgf_product_3c_gspace_3d(S_G_c, G1, G1_sq, G2, G2_sq, G3, G3_sq, l_max, m_max, n_max, &
1354 0 : alpha, beta, gamma, R1, R2, use_coulomb)
1355 : END DO
1356 : END DO
1357 : END DO
1358 : END DO
1359 : END DO
1360 : END DO
1361 : CASE (2)
1362 0 : DO g1x = -G1c(1), G1c(1)
1363 0 : G1_x = ht(:, 1)*g1x
1364 0 : DO g1y = -G1c(2), G1c(2)
1365 0 : G1_y = ht(:, 2)*g1y
1366 0 : DO g1z = -G1c(3), G1c(3)
1367 0 : G1_z = ht(:, 3)*g1z
1368 0 : G1 = G1_x + G1_y + G1_z
1369 0 : G1_sq = G1(1)**2 + G1(2)**2 + G1(3)**2
1370 0 : IF (G1_sq .GT. G_rads_sq(1)) CYCLE
1371 0 : DO g3x = -G3c(1), G3c(1)
1372 0 : G3_x = ht(:, 1)*g3x
1373 0 : DO g3y = -G3c(2), G3c(2)
1374 0 : G3_y = ht(:, 2)*g3y
1375 0 : DO g3z = -G3c(3), G3c(3)
1376 0 : G3_z = ht(:, 3)*g3z
1377 0 : G3 = G3_x + G3_y + G3_z
1378 0 : G3_sq = G3(1)**2 + G3(2)**2 + G3(3)**2
1379 0 : IF (G3_sq .GT. G_rads_sq(3)) CYCLE
1380 0 : G2 = G3 - G1
1381 0 : G2_sq = G2(1)**2 + G2(2)**2 + G2(3)**2
1382 0 : IF (G2_sq .GT. G_rads_sq(2)) CYCLE
1383 0 : IF (use_coulomb) THEN
1384 0 : IF (g3x == 0 .AND. g3y == 0 .AND. g3z == 0) CYCLE
1385 : END IF
1386 : CALL pgf_product_3c_gspace_3d(S_G_c, G1, G1_sq, G2, G2_sq, G3, G3_sq, l_max, m_max, n_max, &
1387 0 : alpha, beta, gamma, R1, R2, use_coulomb)
1388 : END DO
1389 : END DO
1390 : END DO
1391 : END DO
1392 : END DO
1393 : END DO
1394 : CASE (3)
1395 0 : DO g1x = -G1c(1), G1c(1)
1396 0 : G1_x = ht(:, 1)*g1x
1397 0 : DO g1y = -G1c(2), G1c(2)
1398 0 : G1_y = ht(:, 2)*g1y
1399 0 : DO g1z = -G1c(3), G1c(3)
1400 0 : G1_z = ht(:, 3)*g1z
1401 0 : G1 = G1_x + G1_y + G1_z
1402 0 : G1_sq = G1(1)**2 + G1(2)**2 + G1(3)**2
1403 0 : IF (G1_sq .GT. G_rads_sq(1)) CYCLE
1404 0 : DO g2x = -G2c(1), G2c(1)
1405 0 : G2_x = ht(:, 1)*g2x
1406 0 : DO g2y = -G2c(2), G2c(2)
1407 0 : G2_y = ht(:, 2)*g2y
1408 0 : DO g2z = -G2c(3), G2c(3)
1409 0 : G2_z = ht(:, 3)*g2z
1410 0 : G2 = G2_x + G2_y + G2_z
1411 0 : G2_sq = G2(1)**2 + G2(2)**2 + G2(3)**2
1412 0 : IF (G2_sq .GT. G_rads_sq(2)) CYCLE
1413 0 : G3 = G1 + G2
1414 0 : G3_sq = G3(1)**2 + G3(2)**2 + G3(3)**2
1415 0 : IF (G3_sq .GT. G_rads_sq(3)) CYCLE
1416 0 : IF (use_coulomb) THEN
1417 0 : IF (g1x + g2x == 0 .AND. g1y + g2y == 0 .AND. g1z + g2z == 0) CYCLE
1418 : END IF
1419 : CALL pgf_product_3c_gspace_3d(S_G_c, G1, G1_sq, G2, G2_sq, G3, G3_sq, l_max, m_max, n_max, &
1420 0 : alpha, beta, gamma, R1, R2, use_coulomb)
1421 : END DO
1422 : END DO
1423 : END DO
1424 : END DO
1425 : END DO
1426 : END DO
1427 : END SELECT
1428 0 : S_G = REAL(S_G_c, KIND=dp)/vol**2
1429 :
1430 0 : END SUBROUTINE
1431 :
1432 : ! **************************************************************************************************
1433 : !> \brief ...
1434 : !> \param S_G ...
1435 : !> \param G1 ...
1436 : !> \param G1_sq ...
1437 : !> \param G2 ...
1438 : !> \param G2_sq ...
1439 : !> \param G3 ...
1440 : !> \param G3_sq ...
1441 : !> \param l_max ...
1442 : !> \param m_max ...
1443 : !> \param n_max ...
1444 : !> \param alpha ...
1445 : !> \param beta ...
1446 : !> \param gamma ...
1447 : !> \param R1 ...
1448 : !> \param R2 ...
1449 : !> \param coulomb ...
1450 : ! **************************************************************************************************
1451 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, &
1452 : alpha, beta, gamma, R1, R2, coulomb)
1453 :
1454 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: S_G
1455 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G1
1456 : REAL(KIND=dp), INTENT(IN) :: G1_sq
1457 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G2
1458 : REAL(KIND=dp), INTENT(IN) :: G2_sq
1459 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G3
1460 : REAL(KIND=dp), INTENT(IN) :: G3_sq
1461 : INTEGER, INTENT(IN) :: l_max, m_max, n_max
1462 : REAL(KIND=dp), INTENT(IN) :: alpha, beta, gamma
1463 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R1, R2
1464 : LOGICAL, INTENT(IN) :: coulomb
1465 :
1466 : COMPLEX(KIND=dp) :: exp_tot
1467 : INTEGER :: k, l, lco, lx, ly, lz, m, mco, mx, my, &
1468 : mz, n, nco, nx, ny, nz
1469 0 : COMPLEX(KIND=dp), DIMENSION(ncoset(n_max)) :: S_G_n
1470 0 : COMPLEX(KIND=dp), DIMENSION(ncoset(m_max)) :: S_G_m
1471 0 : COMPLEX(KIND=dp), DIMENSION(ncoset(l_max)) :: S_G_l
1472 0 : REAL(KIND=dp), DIMENSION(3, 0:l_max) :: G1_pow_l
1473 0 : REAL(KIND=dp), DIMENSION(3, 0:m_max) :: G2_pow_m
1474 0 : REAL(KIND=dp), DIMENSION(3, 0:n_max) :: G3_pow_n
1475 :
1476 : exp_tot = EXP(gaussi*DOT_PRODUCT(G1, R1))*EXP(-alpha*G1_sq)* & ! cost: 5 exp_w flops
1477 : EXP(-beta*G2_sq)* &
1478 0 : EXP(-gamma*G3_sq)*EXP(gaussi*DOT_PRODUCT(G3, R2))
1479 :
1480 0 : IF (coulomb) exp_tot = exp_tot/G3_sq
1481 :
1482 0 : DO k = 1, 3
1483 0 : G1_pow_l(k, 0) = 1.0_dp
1484 0 : DO l = 1, l_max
1485 0 : G1_pow_l(k, l) = G1_pow_l(k, l - 1)*G1(k)
1486 : END DO
1487 0 : G2_pow_m(k, 0) = 1.0_dp
1488 0 : DO m = 1, m_max
1489 0 : G2_pow_m(k, m) = G2_pow_m(k, m - 1)*G2(k)
1490 : END DO
1491 0 : G3_pow_n(k, 0) = 1.0_dp
1492 0 : DO n = 1, n_max
1493 0 : G3_pow_n(k, n) = G3_pow_n(k, n - 1)*G3(k)
1494 : END DO
1495 : END DO
1496 :
1497 0 : DO lco = 1, ncoset(l_max)
1498 0 : CALL get_l(lco, l, lx, ly, lz)
1499 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)
1500 : END DO
1501 :
1502 0 : DO mco = 1, ncoset(m_max)
1503 0 : CALL get_l(mco, m, mx, my, mz)
1504 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)
1505 : END DO
1506 :
1507 0 : DO nco = 1, ncoset(n_max)
1508 0 : CALL get_l(nco, n, nx, ny, nz)
1509 0 : S_G_n(nco) = G3_pow_n(1, nx)*G3_pow_n(2, ny)*G3_pow_n(3, nz)*i_pow(n)
1510 : END DO
1511 :
1512 0 : DO nco = 1, ncoset(n_max)
1513 0 : DO mco = 1, ncoset(m_max)
1514 0 : DO lco = 1, ncoset(l_max)
1515 : S_G(lco, mco, nco) = S_G(lco, mco, nco) + & ! cost: 4 flops
1516 : S_G_l(lco)*S_G_m(mco)*S_G_n(nco)* &
1517 0 : exp_tot
1518 : END DO
1519 : END DO
1520 : END DO
1521 :
1522 0 : END SUBROUTINE pgf_product_3c_gspace_3d
1523 :
1524 : ! **************************************************************************************************
1525 : !> \brief Roughly estimated number of floating point operations
1526 : !> \param ns_G ...
1527 : !> \param ns_R ...
1528 : !> \param l ...
1529 : !> \param m ...
1530 : !> \param n ...
1531 : !> \return ...
1532 : ! **************************************************************************************************
1533 148500 : PURE FUNCTION nsum_product_3c_gspace_3d(ns_G, ns_R, l, m, n)
1534 : REAL(KIND=dp), INTENT(IN) :: ns_G, ns_R
1535 : INTEGER, INTENT(IN) :: l, m, n
1536 : INTEGER(KIND=int_8) :: nsum_product_3c_gspace_3d
1537 :
1538 : nsum_product_3c_gspace_3d = &
1539 : NINT( &
1540 : ns_G*( &
1541 : (exp_w*2) + &
1542 : ns_R*(exp_w*2 + ncoset(l + m)*7) + &
1543 : 3*nsum_gaussian_overlap(l, m, 1) + &
1544 : ncoset(l)*ncoset(m)*(ncoset(l + m)*4 + ncoset(n)*8)), &
1545 148500 : KIND=int_8)
1546 148500 : END FUNCTION
1547 :
1548 : ! **************************************************************************************************
1549 : !> \brief ...
1550 : !> \param S_G ...
1551 : !> \param l_max ...
1552 : !> \param m_max ...
1553 : !> \param n_max ...
1554 : !> \param R1 ...
1555 : !> \param R2 ...
1556 : !> \param alpha ...
1557 : !> \param beta ...
1558 : !> \param gamma ...
1559 : !> \param hmat ...
1560 : !> \param h_inv ...
1561 : !> \param vol ...
1562 : !> \param R_c ...
1563 : !> \param G_c ...
1564 : !> \param R_rad ...
1565 : !> \param G_rad ...
1566 : ! **************************************************************************************************
1567 5846 : PURE SUBROUTINE pgf_sum_product_3c_gspace_3d(S_G, l_max, m_max, n_max, R1, R2, alpha, beta, gamma, &
1568 : hmat, h_inv, vol, R_c, G_c, R_rad, G_rad)
1569 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: S_G
1570 : INTEGER, INTENT(IN) :: l_max, m_max, n_max
1571 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R1, R2
1572 : REAL(KIND=dp), INTENT(IN) :: alpha, beta, gamma
1573 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv
1574 : REAL(KIND=dp), INTENT(IN) :: vol
1575 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R_c, G_c
1576 : REAL(KIND=dp), INTENT(IN) :: R_rad, G_rad
1577 :
1578 : COMPLEX(KIND=dp) :: exp_tot
1579 : INTEGER :: gx, gy, gz, k, l, lco, lnco, lx, ly, lz, &
1580 : m, mco, n, nco
1581 : COMPLEX(KIND=dp), &
1582 11692 : DIMENSION(ncoset(m_max), ncoset(n_max)) :: S_R
1583 : COMPLEX(KIND=dp), DIMENSION(ncoset(l_max), ncoset( &
1584 11692 : m_max), ncoset(n_max)) :: S_G_c
1585 : REAL(KIND=dp) :: G_rad_sq, G_sq, R_rad_sq
1586 : REAL(KIND=dp), DIMENSION(3) :: G, G_x, G_y, G_z
1587 11692 : REAL(KIND=dp), DIMENSION(3, 0:l_max) :: G_pow_l
1588 : REAL(KIND=dp), DIMENSION(3, 3) :: ht
1589 5846 : REAL(KIND=dp), DIMENSION(ncoset(l_max)) :: S_G_c_l
1590 :
1591 276845 : S_G_c(:, :, :) = 0.0_dp
1592 :
1593 5846 : G_rad_sq = G_rad**2
1594 5846 : R_rad_sq = R_rad**2
1595 :
1596 5846 : lnco = ncoset(l_max)
1597 :
1598 75998 : ht = twopi*TRANSPOSE(h_inv)
1599 36078 : DO gx = -FLOOR(G_c(1)), FLOOR(G_c(1))
1600 120928 : G_x = ht(:, 1)*gx
1601 194252 : DO gy = -FLOOR(G_c(2)), FLOOR(G_c(2))
1602 632696 : G_y = ht(:, 2)*gy
1603 1028374 : DO gz = -FLOOR(G_c(3)), FLOOR(G_c(3))
1604 3359872 : G_z = ht(:, 3)*gz
1605 3359872 : G = G_x + G_y + G_z
1606 839968 : G_sq = G(1)**2 + G(2)**2 + G(3)**2
1607 839968 : IF (G_sq .GT. G_rad_sq) CYCLE
1608 :
1609 1350880 : exp_tot = EXP(-alpha*G_sq)*EXP(gaussi*DOT_PRODUCT(G, R1)) ! cost: exp_w*2 flops
1610 :
1611 1350880 : DO k = 1, 3
1612 1013160 : G_pow_l(k, 0) = 1.0_dp
1613 2236258 : DO l = 1, l_max
1614 1898538 : G_pow_l(k, l) = G_pow_l(k, l - 1)*G(k)
1615 : END DO
1616 : END DO
1617 :
1618 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)
1619 :
1620 1947011 : DO lco = 1, ncoset(l_max)
1621 1609291 : CALL get_l(lco, l, lx, ly, lz)
1622 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
1623 : END DO
1624 :
1625 1519228 : DO nco = 1, ncoset(n_max)
1626 1023334 : CALL get_l(nco, n)
1627 4612179 : DO mco = 1, ncoset(m_max)
1628 2748877 : CALL get_l(mco, m)
1629 15176745 : DO lco = 1, ncoset(l_max)
1630 11404534 : CALL get_l(lco, l)
1631 : S_G_c(lco, mco, nco) = & ! cost: 8 flops
1632 : S_G_c(lco, mco, nco) + &
1633 : S_G_c_l(lco)* &
1634 14153411 : exp_tot*i_pow(l + m + n)*(-1.0_dp)**m*S_R(mco, nco)
1635 : END DO
1636 : END DO
1637 : END DO
1638 : END DO
1639 : END DO
1640 : END DO
1641 276845 : S_G = REAL(S_G_c, KIND=dp)/vol**2
1642 :
1643 5846 : END SUBROUTINE pgf_sum_product_3c_gspace_3d
1644 :
1645 : ! **************************************************************************************************
1646 : !> \brief ...
1647 : !> \param S_R ...
1648 : !> \param l_max ...
1649 : !> \param m_max ...
1650 : !> \param G ...
1651 : !> \param R ...
1652 : !> \param alpha ...
1653 : !> \param beta ...
1654 : !> \param hmat ...
1655 : !> \param h_inv ...
1656 : !> \param vol ...
1657 : !> \param R_c ...
1658 : !> \param R_rad_sq ...
1659 : ! **************************************************************************************************
1660 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)
1661 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: S_R
1662 : INTEGER, INTENT(IN) :: l_max, m_max
1663 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G, R
1664 : REAL(KIND=dp), INTENT(IN) :: alpha, beta
1665 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv
1666 : REAL(KIND=dp), INTENT(IN) :: vol
1667 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R_c
1668 : REAL(KIND=dp), INTENT(IN) :: R_rad_sq
1669 :
1670 : COMPLEX(KIND=dp) :: exp_tot
1671 : INTEGER :: k, l, lco, lx, ly, lz, m, mco, mx, my, &
1672 : mz, sx, sy, sz, t, tco, tx, ty, tz
1673 675440 : COMPLEX(KIND=dp), DIMENSION(ncoset(l_max + m_max)) :: S_R_t
1674 : REAL(KIND=dp) :: c1, c2, Rp_sq
1675 : REAL(KIND=dp), &
1676 675440 : DIMENSION(-1:l_max + m_max + 1, -1:l_max, -1:m_max) :: E1, E2, E3
1677 : REAL(KIND=dp), DIMENSION(3) :: R_l, R_r, Rp, Rx, Ry, Rz, s_shift
1678 675440 : REAL(KIND=dp), DIMENSION(3, 0:l_max + m_max) :: R_pow_t
1679 :
1680 337720 : c1 = 0.25_dp/(alpha + beta)
1681 337720 : c2 = alpha/(alpha + beta)
1682 :
1683 2198096 : S_R_t(:) = 0.0_dp
1684 4109931 : S_R(:, :) = 0.0_dp
1685 :
1686 4390360 : s_shift = MATMUL(h_inv, R)
1687 1350880 : R_l = -R_c + s_shift
1688 1350880 : R_r = R_c + s_shift
1689 :
1690 1084696 : DO sx = CEILING(R_l(1)), FLOOR(R_r(1))
1691 2987904 : Rx = hmat(:, 1)*sx
1692 3046393 : DO sy = CEILING(R_l(2)), FLOOR(R_r(2))
1693 7846788 : Ry = hmat(:, 2)*sy
1694 8205714 : DO sz = CEILING(R_l(3)), FLOOR(R_r(3))
1695 21988164 : Rz = hmat(:, 3)*sz
1696 21988164 : Rp = Rx + Ry + Rz - R
1697 5497041 : Rp_sq = Rp(1)**2 + Rp(2)**2 + Rp(3)**2
1698 5497041 : IF (Rp_sq .GT. R_rad_sq) CYCLE
1699 :
1700 8361484 : exp_tot = EXP(-c1*Rp_sq)*EXP(-gaussi*c2*DOT_PRODUCT(Rp, G)) ! cost: exp_w*2 flops
1701 8361484 : DO k = 1, 3
1702 6271113 : R_pow_t(k, 0) = 1.0_dp
1703 13292668 : DO t = 1, l_max + m_max
1704 11202297 : R_pow_t(k, t) = R_pow_t(k, t - 1)*Rp(k)
1705 : END DO
1706 : END DO
1707 :
1708 13088657 : DO tco = 1, ncoset(l_max + m_max)
1709 9036589 : CALL get_l(tco, t, tx, ty, tz)
1710 : S_R_t(tco) = S_R_t(tco) + & ! cost: 7 flops
1711 : R_pow_t(1, tx)*R_pow_t(2, ty)*R_pow_t(3, tz)* &
1712 14533630 : (-1.0_dp)**t*i_pow(t)*exp_tot
1713 : END DO
1714 :
1715 : END DO
1716 : END DO
1717 : END DO
1718 :
1719 337720 : CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, alpha, beta, G(1), 0.0_dp, 1, E1)
1720 337720 : CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, alpha, beta, G(2), 0.0_dp, 1, E2)
1721 337720 : CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, alpha, beta, G(3), 0.0_dp, 1, E3)
1722 :
1723 1361054 : DO mco = 1, ncoset(m_max)
1724 1023334 : CALL get_l(mco, m, mx, my, mz)
1725 4109931 : DO lco = 1, ncoset(l_max)
1726 2748877 : CALL get_l(lco, l, lx, ly, lz)
1727 8010112 : DO tx = 0, lx + mx
1728 13267665 : DO ty = 0, ly + my
1729 19515063 : DO tz = 0, lz + mz
1730 8996275 : tco = coset(tx, ty, tz)
1731 : S_R(lco, mco) = S_R(lco, mco) + & ! cost: 4 flops
1732 15277162 : E1(tx, lx, mx)*E2(ty, ly, my)*E3(tz, lz, mz)*S_R_t(tco)
1733 :
1734 : END DO
1735 : END DO
1736 : END DO
1737 : END DO
1738 : END DO
1739 :
1740 4109931 : S_R(:, :) = S_R(:, :)*vol/(twopi)**3*(pi/(alpha + beta))**1.5_dp
1741 :
1742 337720 : END SUBROUTINE pgf_sum_product_3c_rspace_3d
1743 :
1744 : ! **************************************************************************************************
1745 : !> \brief Roughly estimated number of floating point operations
1746 : !> \param ns_R1 ...
1747 : !> \param ns_R2 ...
1748 : !> \param l ...
1749 : !> \param m ...
1750 : !> \param n ...
1751 : !> \return ...
1752 : ! **************************************************************************************************
1753 148500 : PURE FUNCTION nsum_3c_rspace_3d(ns_R1, ns_R2, l, m, n)
1754 : REAL(KIND=dp), INTENT(IN) :: ns_R1, ns_R2
1755 : INTEGER, INTENT(IN) :: l, m, n
1756 : INTEGER(KIND=int_8) :: nsum_3c_rspace_3d
1757 :
1758 : nsum_3c_rspace_3d = &
1759 : NINT( &
1760 : ns_R1*( &
1761 : ns_R2*(exp_w + ncoset(l + m + n)*4) + &
1762 : 3*nsum_gaussian_overlap(l, m, 2) + &
1763 : ncoset(m)*ncoset(l)*( &
1764 : ncoset(l + m)*2 + ncoset(n)*ncoset(l + m)*4)), &
1765 148500 : KIND=int_8)
1766 :
1767 148500 : END FUNCTION
1768 :
1769 : ! **************************************************************************************************
1770 : !> \brief ...
1771 : !> \param S_R ...
1772 : !> \param l_max ...
1773 : !> \param m_max ...
1774 : !> \param n_max ...
1775 : !> \param RA ...
1776 : !> \param RB ...
1777 : !> \param RC ...
1778 : !> \param zeta ...
1779 : !> \param zetb ...
1780 : !> \param zetc ...
1781 : !> \param a_mm ...
1782 : !> \param hmat ...
1783 : !> \param h_inv ...
1784 : !> \param R_c ...
1785 : !> \param R_rad ...
1786 : ! **************************************************************************************************
1787 142654 : SUBROUTINE pgf_sum_3c_rspace_3d(S_R, l_max, m_max, n_max, RA, RB, RC, zeta, zetb, zetc, a_mm, &
1788 : hmat, h_inv, R_c, R_rad)
1789 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: S_R
1790 : INTEGER, INTENT(IN) :: l_max, m_max, n_max
1791 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: RA, RB, RC
1792 : REAL(KIND=dp), INTENT(IN) :: zeta, zetb, zetc, a_mm
1793 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv
1794 : REAL(KIND=dp), DIMENSION(2, 3), INTENT(IN) :: R_c
1795 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: R_rad
1796 :
1797 : INTEGER :: k, l, lco, lx, ly, lz, m, mco, mx, my, &
1798 : mz, n, nco, nx, ny, nz, s1x, s1y, s1z, &
1799 : s2x, s2y, s2z, t, tco, tnco, ttco, &
1800 : ttx, tty, ttz, tx, ty, tz
1801 : REAL(KIND=dp) :: alpha, exp_tot, R1_sq, R_sq
1802 142654 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: h_to_c
1803 : REAL(KIND=dp), &
1804 285308 : DIMENSION(-1:l_max + m_max + 1, -1:l_max, -1:m_max) :: E1, E2, E3
1805 : REAL(KIND=dp), DIMENSION(2) :: R_rad_sq
1806 : REAL(KIND=dp), DIMENSION(3) :: R, R1, R1_l, R1_r, R1_tmp, R1x, R1y, &
1807 : R1z, R2_l, R2_r, R2x, R2y, R2z, &
1808 : R_offset, R_tmp, s1_shift, s2_shift
1809 285308 : REAL(KIND=dp), DIMENSION(3, 0:l_max + m_max + n_max) :: R_pow_t
1810 : REAL(KIND=dp), DIMENSION(ncoset(l_max + m_max), &
1811 285308 : ncoset(l_max), ncoset(m_max)) :: Eco
1812 : REAL(KIND=dp), &
1813 285308 : DIMENSION(ncoset(l_max + m_max + n_max)) :: S_R_t
1814 : REAL(KIND=dp), DIMENSION(ncoset(l_max + m_max + n_max) &
1815 285308 : , ncoset(l_max + m_max + n_max)) :: h_to_c_tco
1816 :
1817 142654 : alpha = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm)
1818 :
1819 71800162 : Eco(:, :, :) = 0.0_dp
1820 24264238 : S_R(:, :, :) = 0.0_dp
1821 157194456 : h_to_c_tco(:, :) = 0.0_dp
1822 :
1823 570616 : R_offset = RC - (zeta*RA + zetb*RB)/(zeta + zetb)
1824 :
1825 142654 : CALL create_hermite_to_cartesian(alpha, l_max + m_max + n_max, h_to_c)
1826 :
1827 2682070 : DO tco = 1, ncoset(l_max + m_max + n_max)
1828 2539416 : CALL get_l(tco, t, tx, ty, tz)
1829 8078511 : DO ttx = 0, tx
1830 18802080 : DO tty = 0, ty
1831 37176724 : DO ttz = 0, tz
1832 20914060 : ttco = coset(ttx, tty, ttz)
1833 31780283 : h_to_c_tco(ttco, tco) = h_to_c(ttx, tx)*h_to_c(tty, ty)*h_to_c(ttz, tz)
1834 : END DO
1835 : END DO
1836 : END DO
1837 : END DO
1838 :
1839 2282464 : s1_shift = MATMUL(h_inv, RA - RB)
1840 570616 : R1_l = -R_c(1, :) + s1_shift
1841 570616 : R1_r = R_c(1, :) + s1_shift
1842 :
1843 427962 : R_rad_sq = R_rad**2
1844 :
1845 348524 : DO s1x = CEILING(R1_l(1)), FLOOR(R1_r(1))
1846 823480 : R1x = hmat(:, 1)*s1x
1847 741859 : DO s1y = CEILING(R1_l(2)), FLOOR(R1_r(2))
1848 1573340 : R1y = hmat(:, 2)*s1y
1849 1532744 : DO s1z = CEILING(R1_l(3)), FLOOR(R1_r(3))
1850 3734156 : R1z = hmat(:, 3)*s1z
1851 3734156 : R1 = R1x + R1y + R1z
1852 13798961 : S_R_t(:) = 0.0_dp
1853 3734156 : R1_tmp = R1 - (RA - RB)
1854 933539 : R1_sq = R1_tmp(1)**2 + R1_tmp(2)**2 + R1_tmp(3)**2
1855 :
1856 933539 : IF (R1_sq .GT. R_rad_sq(1)) CYCLE
1857 :
1858 1625580 : R_tmp = R_offset + R1*zeta/(zeta + zetb)
1859 6502320 : s2_shift = -MATMUL(h_inv, R_tmp)
1860 1625580 : R2_l = -R_c(2, :) + s2_shift
1861 1625580 : R2_r = R_c(2, :) + s2_shift
1862 1468530 : DO s2x = CEILING(R2_l(1)), FLOOR(R2_r(1))
1863 4248540 : R2x = hmat(:, 1)*s2x
1864 6145785 : DO s2y = CEILING(R2_l(2)), FLOOR(R2_r(2))
1865 18709020 : R2y = hmat(:, 2)*s2y
1866 34911377 : DO s2z = CEILING(R2_l(3)), FLOOR(R2_r(3))
1867 116687948 : R2z = hmat(:, 3)*s2z
1868 116687948 : R = R_tmp + R2x + R2y + R2z
1869 29171987 : R_sq = R(1)**2 + R(2)**2 + R(3)**2
1870 :
1871 29171987 : IF (R_sq .GT. R_rad_sq(2)) CYCLE
1872 :
1873 14475027 : exp_tot = EXP(-alpha*R_sq) ! cost: exp_w flops
1874 57900108 : DO k = 1, 3
1875 43425081 : R_pow_t(k, 0) = 1.0_dp
1876 156067857 : DO t = 1, l_max + m_max + n_max
1877 141592830 : R_pow_t(k, t) = R_pow_t(k, t - 1)*R(k)
1878 : END DO
1879 : END DO
1880 272529341 : DO tco = 1, ncoset(l_max + m_max + n_max)
1881 253377059 : CALL get_l(tco, t, tx, ty, tz)
1882 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
1883 : END DO
1884 : END DO
1885 : END DO
1886 : END DO
1887 :
1888 5981328 : S_R_t(:) = MATMUL(TRANSPOSE(h_to_c_tco), S_R_t)*(alpha/pi)**1.5_dp
1889 :
1890 406395 : CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, RA(1) - R1(1), RB(1), 2, E1)
1891 406395 : CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, RA(2) - R1(2), RB(2), 2, E2)
1892 406395 : CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, RA(3) - R1(3), RB(3), 2, E3)
1893 :
1894 1416942 : DO mco = 1, ncoset(m_max)
1895 1010547 : CALL get_l(mco, m, mx, my, mz)
1896 5054628 : DO lco = 1, ncoset(l_max)
1897 3637686 : CALL get_l(lco, l, lx, ly, lz)
1898 10890282 : DO tx = 0, lx + mx
1899 20345458 : DO ty = 0, ly + my
1900 33814276 : DO tz = 0, lz + mz
1901 17106504 : tco = coset(tx, ty, tz)
1902 27572227 : Eco(tco, lco, mco) = E1(tx, lx, mx)*E2(ty, ly, my)*E3(tz, lz, mz) ! cost: 2 flops
1903 : END DO
1904 : END DO
1905 : END DO
1906 : END DO
1907 : END DO
1908 :
1909 3040511 : DO nco = 1, ncoset(n_max)
1910 2240781 : CALL get_l(nco, n, nx, ny, nz)
1911 17172306 : DO tco = 1, ncoset(l_max + m_max)
1912 14525130 : CALL get_l(tco, t, tx, ty, tz)
1913 14525130 : tnco = coset(tx + nx, ty + ny, tz + nz)
1914 : S_R(:, :, nco) = S_R(:, :, nco) + & ! cost: 4 flops
1915 : Eco(tco, :, :)* &
1916 1070960721 : (-1)**n*S_R_t(tnco)
1917 :
1918 : END DO
1919 : END DO
1920 : END DO
1921 : END DO
1922 : END DO
1923 142654 : END SUBROUTINE pgf_sum_3c_rspace_3d
1924 :
1925 : ! **************************************************************************************************
1926 : !> \brief Compute bounding box for ellipsoid. This is needed in order to find summation bounds for
1927 : !> sphere for sums over non-orthogonal lattice vectors.
1928 : !> \param s_rad sphere radius
1929 : !> \param s_to_e sphere to ellipsoid trafo
1930 : !> \return ...
1931 : ! **************************************************************************************************
1932 19806720 : PURE FUNCTION ellipsoid_bounds(s_rad, s_to_e)
1933 : REAL(KIND=dp), INTENT(IN) :: s_rad
1934 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: s_to_e
1935 : REAL(KIND=dp), DIMENSION(3) :: ellipsoid_bounds
1936 :
1937 : INTEGER :: i_xyz
1938 :
1939 79226880 : DO i_xyz = 1, 3
1940 79226880 : 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
1941 : END DO
1942 :
1943 : END FUNCTION ellipsoid_bounds
1944 :
1945 : ! **************************************************************************************************
1946 : !> \brief Roughly estimated number of floating point operations
1947 : !> \param l ...
1948 : !> \param m ...
1949 : !> \param H_or_C_product ...
1950 : !> \return ...
1951 : ! **************************************************************************************************
1952 297000 : PURE FUNCTION nsum_gaussian_overlap(l, m, H_or_C_product)
1953 : INTEGER, INTENT(IN) :: l, m, H_or_C_product
1954 : INTEGER :: nsum_gaussian_overlap
1955 :
1956 : INTEGER :: loop
1957 :
1958 297000 : nsum_gaussian_overlap = exp_w
1959 297000 : loop = (m + 1)*(l + 1)*(m + l + 2)
1960 297000 : IF (H_or_C_product .EQ. 1) THEN
1961 148500 : nsum_gaussian_overlap = nsum_gaussian_overlap + loop*16
1962 : ELSE
1963 148500 : nsum_gaussian_overlap = nsum_gaussian_overlap + loop*32
1964 : END IF
1965 297000 : END FUNCTION
1966 :
1967 : ! **************************************************************************************************
1968 : !> \brief ...
1969 : !> \param lco ...
1970 : !> \param l ...
1971 : !> \param lx ...
1972 : !> \param ly ...
1973 : !> \param lz ...
1974 : ! **************************************************************************************************
1975 843119003 : PURE ELEMENTAL SUBROUTINE get_l(lco, l, lx, ly, lz)
1976 : INTEGER, INTENT(IN) :: lco
1977 : INTEGER, INTENT(OUT) :: l
1978 : INTEGER, INTENT(OUT), OPTIONAL :: lx, ly, lz
1979 :
1980 3372476012 : l = SUM(indco(:, lco))
1981 843119003 : IF (PRESENT(lx)) lx = indco(1, lco)
1982 843119003 : IF (PRESENT(ly)) ly = indco(2, lco)
1983 843119003 : IF (PRESENT(lz)) lz = indco(3, lco)
1984 843119003 : END SUBROUTINE
1985 :
1986 : ! **************************************************************************************************
1987 : !> \brief ...
1988 : !> \param i ...
1989 : !> \return ...
1990 : ! **************************************************************************************************
1991 479950601 : PURE ELEMENTAL FUNCTION i_pow(i)
1992 : INTEGER, INTENT(IN) :: i
1993 : COMPLEX(KIND=dp) :: i_pow
1994 :
1995 : COMPLEX(KIND=dp), DIMENSION(0:3), PARAMETER :: &
1996 : 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)/)
1997 :
1998 479950601 : i_pow = ip(MOD(i, 4))
1999 :
2000 479950601 : END FUNCTION
2001 :
2002 3955399 : END MODULE eri_mme_lattice_summation
|
