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 Minimax-Ewald (MME) method for calculating 2-center and 3-center
10 : !> electron repulsion integrals (ERI) of periodic systems using a
11 : !> Hermite Gaussian basis.
12 : !> The method relies on analytical Fourier transforms of Cartesian and
13 : !> Hermite Gaussian functions and Poisson summation formula to represent
14 : !> ERIs as a discrete sum over direct lattice vectors or reciprocal
15 : !> lattice vectors. The reciprocal space potential 1/G^2 is approximated
16 : !> by a linear combination of Gaussians employing minimax approximation.
17 : !> Not yet implemented: 3c ERIs for nonorthogonal cells.
18 : !> \par History
19 : !> 2015 09 created
20 : !> \author Patrick Seewald
21 : ! **************************************************************************************************
22 :
23 : MODULE eri_mme_integrate
24 : USE ao_util, ONLY: exp_radius
25 : USE eri_mme_gaussian, ONLY: hermite_gauss_norm
26 : USE eri_mme_lattice_summation, ONLY: &
27 : ellipsoid_bounds, eri_mme_2c_get_bounds, eri_mme_2c_get_rads, eri_mme_3c_get_bounds, &
28 : eri_mme_3c_get_rads, get_l, pgf_sum_2c_gspace_1d, pgf_sum_2c_gspace_3d, &
29 : pgf_sum_2c_rspace_1d, pgf_sum_2c_rspace_3d, pgf_sum_3c_1d, pgf_sum_3c_3d
30 : USE eri_mme_types, ONLY: eri_mme_param,&
31 : get_minimax_from_cutoff
32 : USE kinds, ONLY: dp,&
33 : int_8
34 : USE mathconstants, ONLY: pi,&
35 : twopi
36 : USE orbital_pointers, ONLY: coset,&
37 : ncoset
38 : #include "../base/base_uses.f90"
39 :
40 : IMPLICIT NONE
41 :
42 : PRIVATE
43 :
44 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
45 :
46 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eri_mme_integrate'
47 :
48 : PUBLIC :: eri_mme_2c_integrate, eri_mme_3c_integrate
49 :
50 : CONTAINS
51 :
52 : ! **************************************************************************************************
53 : !> \brief Low-level integration routine for 2-center ERIs.
54 : !> \param param ...
55 : !> \param la_min ...
56 : !> \param la_max ...
57 : !> \param lb_min ...
58 : !> \param lb_max ...
59 : !> \param zeta ...
60 : !> \param zetb ...
61 : !> \param rab ...
62 : !> \param hab ...
63 : !> \param o1 ...
64 : !> \param o2 ...
65 : !> \param G_count ...
66 : !> \param R_count ...
67 : !> \param normalize calculate integrals w.r.t. normalized Hermite-Gaussians
68 : !> \param potential use exact potential instead of minimax approx. (for testing only)
69 : !> \param pot_par potential parameter
70 : ! **************************************************************************************************
71 55934 : SUBROUTINE eri_mme_2c_integrate(param, la_min, la_max, lb_min, lb_max, zeta, zetb, rab, &
72 55934 : hab, o1, o2, G_count, R_count, normalize, potential, pot_par)
73 : TYPE(eri_mme_param), INTENT(IN) :: param
74 : INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max
75 : REAL(KIND=dp), INTENT(IN) :: zeta, zetb
76 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
77 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: hab
78 : INTEGER, INTENT(IN) :: o1, o2
79 : INTEGER, INTENT(INOUT), OPTIONAL :: G_count, R_count
80 : LOGICAL, INTENT(IN), OPTIONAL :: normalize
81 : INTEGER, INTENT(IN), OPTIONAL :: potential
82 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par
83 :
84 : REAL(KIND=dp), PARAMETER :: rel_prec = 1.0E-02_dp
85 :
86 : INTEGER :: ax, ay, az, bx, by, bz, grid, i_aw, &
87 : i_xyz, ico, jco, l_max, la, lb, n_aw
88 : INTEGER(KIND=int_8), DIMENSION(2) :: n_sum_3d
89 : INTEGER(KIND=int_8), DIMENSION(2, 3) :: n_sum_1d
90 : INTEGER, DIMENSION(3) :: la_xyz, lb_xyz
91 : LOGICAL :: do_g_sum, exact, is_ortho, norm
92 : REAL(KIND=dp) :: alpha_G, alpha_R, cutoff, G_rad, G_res, &
93 : Imm, inv_lgth, Ixyz, lgth, max_error, &
94 : prefac, R_rad, R_res, vol
95 55934 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: S_G_1, S_G_2, S_G_no, S_G_no_H
96 55934 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: S_G
97 : REAL(KIND=dp), DIMENSION(3) :: G_bounds, R_bounds
98 : REAL(KIND=dp), DIMENSION(3, 3) :: h_inv, hmat
99 55934 : REAL(KIND=dp), DIMENSION(:), POINTER :: aw
100 :
101 0 : CPASSERT(param%is_valid)
102 :
103 : grid = 0
104 : CALL eri_mme_2c_get_rads(la_max, lb_max, zeta, zetb, 0.0_dp, param%G_min, param%R_min, param%sum_precision, &
105 55934 : G_rad=G_rad)
106 55934 : cutoff = G_rad**2/2
107 55934 : CALL get_minimax_from_cutoff(param%minimax_grid, cutoff, n_aw, aw, grid)
108 :
109 55934 : CPASSERT(grid > 0)
110 :
111 : ! cell info
112 727142 : h_inv = param%h_inv
113 727142 : hmat = param%hmat
114 55934 : vol = param%vol
115 :
116 55934 : IF (PRESENT(normalize)) THEN
117 2304 : norm = normalize
118 : ELSE
119 : norm = .FALSE.
120 : END IF
121 :
122 55934 : l_max = la_max + lb_max
123 :
124 55934 : IF (PRESENT(potential)) THEN
125 : exact = .TRUE.
126 : ELSE
127 : exact = .FALSE.
128 : END IF
129 :
130 : IF (exact) THEN
131 192 : is_ortho = .FALSE.
132 : ELSE
133 55742 : is_ortho = param%is_ortho
134 : END IF
135 :
136 55934 : IF (is_ortho) THEN
137 225320 : ALLOCATE (S_G(0:l_max, 3, n_aw))
138 45064 : S_G = 0.0_dp
139 :
140 45064 : IF (param%debug) THEN
141 0 : ALLOCATE (S_G_1(0:l_max))
142 0 : ALLOCATE (S_G_2(0:l_max))
143 : END IF
144 : ELSE
145 32610 : ALLOCATE (S_G_no(ncoset(l_max)))
146 10870 : S_G_no(:) = 0.0_dp
147 21740 : ALLOCATE (S_G_no_H(ncoset(l_max)))
148 : END IF
149 :
150 55934 : IF (exact) THEN
151 192 : alpha_G = 0.25_dp/zeta + 0.25_dp/zetb
152 : ! resolution for Gaussian width
153 192 : G_res = 0.5_dp*param%G_min
154 192 : R_res = 0.5_dp*param%R_min
155 :
156 192 : G_rad = exp_radius(l_max, alpha_G, rel_prec*param%sum_precision, 1.0_dp, epsabs=G_res)
157 2496 : G_bounds(:) = ellipsoid_bounds(G_rad, TRANSPOSE(hmat)/twopi)
158 768 : CALL pgf_sum_2c_gspace_3d(S_G_no, l_max, -rab, alpha_G, h_inv, G_bounds, G_rad, vol, potential, pot_par)
159 : ELSE
160 :
161 947846 : DO i_aw = 1, n_aw
162 :
163 : CALL eri_mme_2c_get_bounds(hmat, h_inv, vol, is_ortho, param%G_min, param%R_min, la_max, lb_max, &
164 : zeta, zetb, aw(i_aw), param%sum_precision, n_sum_1d, n_sum_3d, &
165 892104 : G_bounds, G_rad, R_bounds, R_rad)
166 892104 : alpha_G = aw(i_aw) + 0.25_dp/zeta + 0.25_dp/zetb
167 892104 : alpha_R = 0.25_dp/alpha_G
168 947846 : IF (is_ortho) THEN ! orthorhombic cell
169 :
170 : ! 1) precompute Ewald-like sum
171 :
172 2981984 : DO i_xyz = 1, 3
173 2236488 : lgth = ABS(hmat(i_xyz, i_xyz))
174 2236488 : inv_lgth = ABS(h_inv(i_xyz, i_xyz))
175 :
176 : ! perform sum in R or G space. Choose the space in which less summands are required for convergence
177 2236488 : do_g_sum = n_sum_1d(1, i_xyz) < n_sum_1d(2, i_xyz) !G_bounds < R_bounds
178 :
179 2236488 : IF (do_g_sum) THEN
180 202238 : CALL pgf_sum_2c_gspace_1d(S_G(:, i_xyz, i_aw), -rab(i_xyz), alpha_G, inv_lgth, G_bounds(i_xyz))
181 202238 : IF (PRESENT(G_count)) G_count = G_count + 1
182 : ELSE
183 2034250 : CALL pgf_sum_2c_rspace_1d(S_G(:, i_xyz, i_aw), -rab(i_xyz), alpha_R, lgth, R_bounds(i_xyz))
184 2034250 : IF (PRESENT(R_count)) R_count = R_count + 1
185 : END IF
186 :
187 2981984 : IF (param%debug) THEN
188 : ! check consistency of summation methods
189 0 : CALL pgf_sum_2c_gspace_1d(S_G_1, -rab(i_xyz), alpha_G, inv_lgth, G_bounds(i_xyz))
190 0 : CALL pgf_sum_2c_rspace_1d(S_G_2, -rab(i_xyz), alpha_R, lgth, R_bounds(i_xyz))
191 0 : max_error = MAXVAL(ABS(S_G_1 - S_G_2)/(0.5_dp*(ABS(S_G_1) + ABS(S_G_2)) + 1.0_dp))
192 :
193 0 : CPASSERT(max_error <= param%debug_delta)
194 : END IF
195 : END DO
196 :
197 : ELSE ! general cell
198 :
199 146608 : do_g_sum = n_sum_3d(1) < n_sum_3d(2) !PRODUCT(2*R_bounds+1) > PRODUCT(2*G_bounds+1)
200 :
201 146608 : IF (do_g_sum) THEN
202 109556 : CALL pgf_sum_2c_gspace_3d(S_G_no_H, l_max, -rab, alpha_G, h_inv, G_bounds, G_rad, vol)
203 27389 : IF (PRESENT(G_count)) G_count = G_count + 1
204 : ELSE
205 476876 : CALL pgf_sum_2c_rspace_3d(S_G_no_H, l_max, -rab, alpha_R, hmat, h_inv, R_bounds, R_rad)
206 119219 : IF (PRESENT(R_count)) R_count = R_count + 1
207 : END IF
208 2579036 : S_G_no(:) = S_G_no(:) + aw(n_aw + i_aw)*S_G_no_H
209 : END IF
210 : END DO
211 : END IF
212 :
213 : ! prefactor for integral values (unnormalized Hermite Gaussians)
214 55934 : prefac = SQRT(1.0_dp/(zeta*zetb))
215 :
216 : ! 2) Assemble integral values from Ewald sums
217 294428 : DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
218 238494 : CALL get_l(jco, lb, bx, by, bz)
219 953976 : lb_xyz = [bx, by, bz]
220 3577034 : DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
221 3282606 : CALL get_l(ico, la, ax, ay, az)
222 13130424 : la_xyz = [ax, ay, az]
223 3282606 : IF (is_ortho) THEN
224 2481319 : Imm = 0.0_dp
225 41625051 : DO i_aw = 1, n_aw
226 : Ixyz = 1.0_dp
227 156574928 : DO i_xyz = 1, 3
228 156574928 : Ixyz = Ixyz*S_G(la_xyz(i_xyz) + lb_xyz(i_xyz), i_xyz, i_aw)*prefac
229 : END DO
230 41625051 : Imm = Imm + aw(n_aw + i_aw)*Ixyz
231 : END DO
232 : ELSE
233 801287 : Imm = S_G_no(coset(ax + bx, ay + by, az + bz))*prefac**3
234 : END IF
235 3282606 : IF (la + lb == 0 .AND. .NOT. exact) THEN
236 226050 : Imm = Imm - SUM(aw(n_aw + 1:2*n_aw))*prefac**3/vol ! subtracting G = 0 term
237 : END IF
238 3521100 : IF (.NOT. norm) THEN
239 : ! rescaling needed due to Hermite Gaussians (such that they can be contracted same way as Cartesian Gaussians)
240 : ! and factor of 4 pi**4 (-1)**lb
241 963438 : hab(o1 + ico, o2 + jco) = Imm*4.0_dp*pi**4/((2.0_dp*zeta)**la*(-2.0_dp*zetb)**lb)
242 : ELSE
243 : ! same thing for normalized Hermite Gaussians
244 : hab(o1 + ico, o2 + jco) = Imm*4.0_dp*pi**4*(-1.0_dp)**lb*hermite_gauss_norm(zeta, la_xyz)* &
245 2319168 : hermite_gauss_norm(zetb, lb_xyz)
246 : END IF
247 : END DO ! la
248 : END DO ! lb
249 :
250 111868 : END SUBROUTINE eri_mme_2c_integrate
251 :
252 : ! **************************************************************************************************
253 : !> \brief Low-level integration routine for 3-center ERIs
254 : !> \param param ...
255 : !> \param la_min ...
256 : !> \param la_max ...
257 : !> \param lb_min ...
258 : !> \param lb_max ...
259 : !> \param lc_min ...
260 : !> \param lc_max ...
261 : !> \param zeta ...
262 : !> \param zetb ...
263 : !> \param zetc ...
264 : !> \param RA ...
265 : !> \param RB ...
266 : !> \param RC ...
267 : !> \param habc ...
268 : !> \param o1 ...
269 : !> \param o2 ...
270 : !> \param o3 ...
271 : !> \param GG_count ...
272 : !> \param GR_count ...
273 : !> \param RR_count ...
274 : ! **************************************************************************************************
275 154554 : SUBROUTINE eri_mme_3c_integrate(param, la_min, la_max, lb_min, lb_max, lc_min, lc_max, zeta, zetb, zetc, RA, RB, RC, &
276 154554 : habc, o1, o2, o3, GG_count, GR_count, RR_count)
277 : TYPE(eri_mme_param), INTENT(IN) :: param
278 : INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max, lc_min, &
279 : lc_max
280 : REAL(KIND=dp), INTENT(IN) :: zeta, zetb, zetc
281 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: RA, RB, RC
282 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: habc
283 : INTEGER, INTENT(IN) :: o1, o2, o3
284 : INTEGER, INTENT(INOUT), OPTIONAL :: GG_count, GR_count, RR_count
285 :
286 154554 : CPASSERT(param%is_valid)
287 154554 : IF (param%is_ortho) THEN
288 : CALL eri_mme_3c_integrate_ortho(param, la_min, la_max, lb_min, lb_max, lc_min, lc_max, zeta, zetb, zetc, RA, RB, RC, &
289 134154 : habc, o1, o2, o3, RR_count)
290 :
291 : ELSE
292 : CALL eri_mme_3c_integrate_nonortho(param, la_min, la_max, lb_min, lb_max, lc_min, lc_max, zeta, zetb, zetc, RA, RB, RC, &
293 20400 : habc, o1, o2, o3, GG_count, GR_count, RR_count)
294 :
295 : END IF
296 154554 : END SUBROUTINE eri_mme_3c_integrate
297 :
298 : ! **************************************************************************************************
299 : !> \brief ...
300 : !> \param param ...
301 : !> \param la_min ...
302 : !> \param la_max ...
303 : !> \param lb_min ...
304 : !> \param lb_max ...
305 : !> \param lc_min ...
306 : !> \param lc_max ...
307 : !> \param zeta ...
308 : !> \param zetb ...
309 : !> \param zetc ...
310 : !> \param RA ...
311 : !> \param RB ...
312 : !> \param RC ...
313 : !> \param habc ...
314 : !> \param o1 ...
315 : !> \param o2 ...
316 : !> \param o3 ...
317 : !> \param RR_count ...
318 : ! **************************************************************************************************
319 134154 : SUBROUTINE eri_mme_3c_integrate_ortho(param, la_min, la_max, lb_min, lb_max, lc_min, lc_max, zeta, zetb, zetc, RA, RB, RC, &
320 134154 : habc, o1, o2, o3, RR_count)
321 : TYPE(eri_mme_param), INTENT(IN) :: param
322 : INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max, lc_min, &
323 : lc_max
324 : REAL(KIND=dp), INTENT(IN) :: zeta, zetb, zetc
325 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: RA, RB, RC
326 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: habc
327 : INTEGER, INTENT(IN) :: o1, o2, o3
328 : INTEGER, INTENT(INOUT), OPTIONAL :: RR_count
329 :
330 : INTEGER :: grid, i_aw, lmax_0, n_aw
331 : INTEGER(KIND=int_8), DIMENSION(3) :: n_sum_3d
332 : INTEGER(KIND=int_8), DIMENSION(3, 3) :: n_sum_1d
333 : REAL(KIND=dp) :: alpha_R_0, cutoff, G_res, lgth, prefac, &
334 : R_rad_0, R_res, vol
335 134154 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: S_G_0
336 : REAL(KIND=dp), ALLOCATABLE, &
337 134154 : DIMENSION(:, :, :, :, :) :: S_G
338 : REAL(KIND=dp), DIMENSION(2) :: R_rads_3
339 : REAL(KIND=dp), DIMENSION(2, 3) :: R_bounds_3
340 : REAL(KIND=dp), DIMENSION(3) :: G_rads_1, R_rads_2
341 : REAL(KIND=dp), DIMENSION(3, 3) :: G_bounds_1, h_inv, hmat, R_bounds_2
342 134154 : REAL(KIND=dp), DIMENSION(:), POINTER :: aw
343 :
344 : grid = 0
345 :
346 : CALL eri_mme_3c_get_rads(la_max, lb_max, lc_max, zeta, zetb, zetc, 0.0_dp, param%G_min, param%R_min, &
347 134154 : param%sum_precision, G_rads_1=G_rads_1)
348 :
349 134154 : cutoff = (MIN(G_rads_1(1), G_rads_1(2) + G_rads_1(3)))**2/2
350 :
351 134154 : CALL get_minimax_from_cutoff(param%minimax_grid, cutoff, n_aw, aw, grid)
352 :
353 134154 : CPASSERT(grid > 0)
354 :
355 : ! minimax coeffs
356 134154 : n_aw = param%minimax_grid(grid)%n_minimax
357 134154 : aw => param%minimax_grid(grid)%minimax_aw
358 :
359 : ! cell info
360 1744002 : h_inv = param%h_inv
361 1744002 : hmat = param%hmat
362 134154 : vol = param%vol
363 :
364 : ! prefactor for integral values (unnormalized Hermite Gaussians)
365 134154 : prefac = (zeta*zetb*zetc)**(-1.5_dp)*pi**(11.0_dp/2.0_dp)*4.0_dp
366 :
367 : ! Preparations for G=0 component
368 134154 : G_res = 0.5_dp*param%G_min
369 134154 : R_res = 0.5_dp*param%R_min
370 :
371 804924 : ALLOCATE (S_G(n_aw, 3, 0:la_max, 0:lb_max, 0:lc_max))
372 :
373 : ! G= 0 components
374 134154 : IF (lc_min == 0) THEN
375 237912 : ALLOCATE (S_G_0(0:la_max + lb_max, 3))
376 :
377 59478 : alpha_R_0 = zeta*zetb/(zeta + zetb)
378 59478 : lmax_0 = la_max + lb_max
379 59478 : R_rad_0 = exp_radius(lmax_0, alpha_R_0, param%sum_precision, 1.0_dp, epsabs=R_res)
380 :
381 59478 : lgth = ABS(hmat(1, 1))
382 59478 : CALL pgf_sum_2c_rspace_1d(S_G_0(:, 1), RB(1) - RA(1), alpha_R_0, lgth, R_rad_0/lgth)
383 59478 : lgth = ABS(hmat(2, 2))
384 59478 : CALL pgf_sum_2c_rspace_1d(S_G_0(:, 2), RB(2) - RA(2), alpha_R_0, lgth, R_rad_0/lgth)
385 59478 : lgth = ABS(hmat(3, 3))
386 59478 : CALL pgf_sum_2c_rspace_1d(S_G_0(:, 3), RB(3) - RA(3), alpha_R_0, lgth, R_rad_0/lgth)
387 : END IF
388 :
389 2242194 : DO i_aw = 1, n_aw
390 : CALL eri_mme_3c_get_bounds(hmat, h_inv, vol, .TRUE., param%G_min, param%R_min, la_max, lb_max, lc_max, &
391 : zeta, zetb, zetc, aw(i_aw), param%sum_precision, n_sum_1d, n_sum_3d, &
392 2108040 : G_bounds_1, G_rads_1, R_bounds_2, R_rads_2, R_bounds_3, R_rads_3)
393 : CALL pgf_sum_3c_1d(S_G(i_aw, 1, :, :, :), RA(1), RB(1), RC(1), &
394 : zeta, zetb, zetc, aw(i_aw), ABS(hmat(1, 1)), &
395 2108040 : R_bounds_3(:, 1))
396 :
397 : CALL pgf_sum_3c_1d(S_G(i_aw, 2, :, :, :), RA(2), RB(2), RC(2), &
398 : zeta, zetb, zetc, aw(i_aw), ABS(hmat(2, 2)), &
399 2108040 : R_bounds_3(:, 2))
400 :
401 : CALL pgf_sum_3c_1d(S_G(i_aw, 3, :, :, :), RA(3), RB(3), RC(3), &
402 : zeta, zetb, zetc, aw(i_aw), ABS(hmat(3, 3)), &
403 2108040 : R_bounds_3(:, 3))
404 :
405 2242194 : IF (PRESENT(RR_count)) RR_count = RR_count + 3
406 : END DO
407 :
408 : CALL eri_mme_3c_collect_ortho(vol, la_min, la_max, lb_min, lb_max, lc_min, lc_max, &
409 134154 : zeta, zetb, zetc, n_aw, aw, S_G, S_G_0, prefac, habc, o1, o2, o3)
410 :
411 268308 : END SUBROUTINE eri_mme_3c_integrate_ortho
412 :
413 : ! **************************************************************************************************
414 : !> \brief ...
415 : !> \param param ...
416 : !> \param la_min ...
417 : !> \param la_max ...
418 : !> \param lb_min ...
419 : !> \param lb_max ...
420 : !> \param lc_min ...
421 : !> \param lc_max ...
422 : !> \param zeta ...
423 : !> \param zetb ...
424 : !> \param zetc ...
425 : !> \param RA ...
426 : !> \param RB ...
427 : !> \param RC ...
428 : !> \param habc ...
429 : !> \param o1 ...
430 : !> \param o2 ...
431 : !> \param o3 ...
432 : !> \param GG_count ...
433 : !> \param GR_count ...
434 : !> \param RR_count ...
435 : ! **************************************************************************************************
436 20400 : SUBROUTINE eri_mme_3c_integrate_nonortho(param, la_min, la_max, lb_min, lb_max, lc_min, lc_max, zeta, zetb, zetc, RA, RB, RC, &
437 20400 : habc, o1, o2, o3, GG_count, GR_count, RR_count)
438 :
439 : TYPE(eri_mme_param), INTENT(IN) :: param
440 : INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max, lc_min, &
441 : lc_max
442 : REAL(KIND=dp), INTENT(IN) :: zeta, zetb, zetc
443 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: RA, RB, RC
444 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: habc
445 : INTEGER, INTENT(IN) :: o1, o2, o3
446 : INTEGER, INTENT(INOUT), OPTIONAL :: GG_count, GR_count, RR_count
447 :
448 : INTEGER :: ax, ay, az, bx, by, bz, cx, cy, cz, &
449 : grid, i_aw, ico, ir, jco, kco, la, lb, &
450 : lc, lmax_0, method, n_aw, nresults, &
451 : sum_method
452 : INTEGER(KIND=int_8), DIMENSION(3) :: n_sum_3d
453 : INTEGER(KIND=int_8), DIMENSION(3, 3) :: n_sum_1d
454 : LOGICAL :: db_sum1, db_sum2, db_sum3, do_g_sum_0
455 : REAL(KIND=dp) :: alpha_G_0, alpha_R_0, cutoff, G_rad_0, &
456 : G_res, max_error, max_result, &
457 : min_result, prefac, R_rad_0, R_res, vol
458 20400 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: results_no, S_G_0_no
459 20400 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: S_G_no, S_G_no_1, S_G_no_2, S_G_no_3, &
460 20400 : S_G_no_H
461 : REAL(KIND=dp), DIMENSION(2) :: R_rads_3
462 : REAL(KIND=dp), DIMENSION(2, 3) :: R_bounds_3
463 : REAL(KIND=dp), DIMENSION(3) :: G_bound_0, G_rads_1, R_0, R_bound_0, &
464 : R_rads_2
465 : REAL(KIND=dp), DIMENSION(3, 3) :: G_bounds_1, h_inv, hmat, R_bounds_2
466 20400 : REAL(KIND=dp), DIMENSION(:), POINTER :: aw
467 :
468 0 : CPASSERT(param%is_valid)
469 :
470 : grid = 0
471 :
472 : CALL eri_mme_3c_get_rads(la_max, lb_max, lc_max, zeta, zetb, zetc, 0.0_dp, param%G_min, param%R_min, &
473 20400 : param%sum_precision, G_rads_1=G_rads_1)
474 :
475 20400 : cutoff = (MIN(G_rads_1(1), G_rads_1(2) + G_rads_1(3)))**2/2
476 :
477 20400 : CALL get_minimax_from_cutoff(param%minimax_grid, cutoff, n_aw, aw, grid)
478 :
479 20400 : CPASSERT(grid > 0)
480 :
481 : ! minimax coeffs
482 20400 : n_aw = param%minimax_grid(grid)%n_minimax
483 20400 : aw => param%minimax_grid(grid)%minimax_aw
484 :
485 : ! cell info
486 265200 : h_inv = param%h_inv
487 265200 : hmat = param%hmat
488 20400 : vol = param%vol
489 :
490 : ! prefactor for integral values (unnormalized Hermite Gaussians)
491 20400 : prefac = (zeta*zetb*zetc)**(-1.5_dp)*pi**(11.0_dp/2.0_dp)*4.0_dp
492 :
493 20400 : IF (param%debug) THEN
494 0 : ALLOCATE (S_G_no_1(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
495 0 : ALLOCATE (S_G_no_2(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
496 0 : ALLOCATE (S_G_no_3(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
497 : END IF
498 :
499 : ! Preparations for G=0 component
500 20400 : G_res = 0.5_dp*param%G_min
501 20400 : R_res = 0.5_dp*param%R_min
502 :
503 102000 : ALLOCATE (S_G_no(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
504 :
505 20400 : S_G_no(:, :, :) = 0.0_dp
506 20400 : IF (param%debug) THEN
507 0 : S_G_no_1(:, :, :) = -1.0_dp
508 0 : S_G_no_2(:, :, :) = -1.0_dp
509 0 : S_G_no_3(:, :, :) = -1.0_dp
510 : END IF
511 102000 : ALLOCATE (S_G_no_H(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
512 :
513 : ! G= 0 components
514 20400 : IF (lc_min == 0) THEN
515 27000 : ALLOCATE (S_G_0_no(ncoset(la_max + lb_max)))
516 9000 : alpha_G_0 = 0.25_dp/zetb + 0.25_dp/zeta
517 9000 : alpha_R_0 = 0.25_dp/alpha_G_0
518 9000 : lmax_0 = la_max + lb_max
519 36000 : R_0 = RB - RA
520 9000 : G_rad_0 = exp_radius(lmax_0, alpha_G_0, param%sum_precision, 1.0_dp, epsabs=G_res)
521 9000 : R_rad_0 = exp_radius(lmax_0, alpha_R_0, param%sum_precision, 1.0_dp, epsabs=R_res)
522 117000 : G_bound_0 = ellipsoid_bounds(G_rad_0, TRANSPOSE(hmat)/(2.0_dp*pi))
523 9000 : R_bound_0 = ellipsoid_bounds(R_rad_0, h_inv)
524 63000 : do_g_sum_0 = PRODUCT(2*R_bound_0 + 1) > PRODUCT(2*G_bound_0 + 1)
525 9000 : IF (do_g_sum_0) THEN
526 0 : CALL pgf_sum_2c_gspace_3d(S_G_0_no, lmax_0, R_0, alpha_G_0, h_inv, G_bound_0, G_rad_0, vol)
527 : ELSE
528 9000 : CALL pgf_sum_2c_rspace_3d(S_G_0_no, lmax_0, R_0, alpha_R_0, hmat, h_inv, R_bound_0, R_rad_0)
529 : END IF
530 : END IF
531 :
532 168900 : DO i_aw = 1, n_aw
533 : CALL eri_mme_3c_get_bounds(hmat, h_inv, vol, .FALSE., param%G_min, param%R_min, la_max, lb_max, lc_max, &
534 : zeta, zetb, zetc, aw(i_aw), param%sum_precision, n_sum_1d, n_sum_3d, &
535 148500 : G_bounds_1, G_rads_1, R_bounds_2, R_rads_2, R_bounds_3, R_rads_3)
536 594000 : sum_method = MINLOC(n_sum_3d, DIM=1)
537 :
538 : CALL pgf_sum_3c_3d(S_G_no_H, la_max, lb_max, lc_max, RA, RB, RC, &
539 : zeta, zetb, zetc, aw(i_aw), hmat, h_inv, vol, &
540 : G_bounds_1, R_bounds_2, R_bounds_3, &
541 : G_rads_1, R_rads_2, R_rads_3, &
542 148500 : method=sum_method, method_out=method)
543 24569280 : S_G_no(:, :, :) = S_G_no(:, :, :) + aw(n_aw + i_aw)*S_G_no_H(:, :, :)
544 :
545 0 : SELECT CASE (method)
546 : CASE (1)
547 0 : IF (PRESENT(GG_count)) GG_count = GG_count + 1
548 : CASE (2)
549 5846 : IF (PRESENT(GR_count)) GR_count = GR_count + 1
550 : CASE (3)
551 142654 : IF (PRESENT(RR_count)) RR_count = RR_count + 1
552 : CASE DEFAULT
553 : CALL cp_abort(__LOCATION__, &
554 : "Only 1, 2, 3 are supported as method "// &
555 148500 : "for eri_mme_3c_integrate_nonortho")
556 : END SELECT
557 :
558 317400 : IF (param%debug) THEN
559 0 : nresults = 0
560 : ! check consistency of summation methods
561 :
562 0 : db_sum1 = (n_sum_3d(1)) < INT(param%debug_nsum, KIND=int_8)
563 0 : db_sum2 = (n_sum_3d(2)) < INT(param%debug_nsum, KIND=int_8)
564 0 : db_sum3 = (n_sum_3d(3)) < INT(param%debug_nsum, KIND=int_8)
565 :
566 0 : IF (param%unit_nr > 0) THEN
567 0 : WRITE (param%unit_nr, *) "ERI_MME DEBUG | number of summands (GG / GR / RR)", n_sum_3d
568 0 : WRITE (param%unit_nr, *) "ERI_MME DEBUG | sum methods to be compared (GG / GR / RR)", db_sum1, db_sum2, db_sum3
569 : END IF
570 :
571 0 : S_G_no_1(:, :, :) = 0.0_dp
572 0 : S_G_no_2(:, :, :) = 0.0_dp
573 0 : S_G_no_3(:, :, :) = 0.0_dp
574 :
575 0 : IF (db_sum1) THEN
576 : CALL pgf_sum_3c_3d(S_G_no_1, la_max, lb_max, lc_max, RA, RB, RC, &
577 : zeta, zetb, zetc, aw(i_aw), hmat, h_inv, vol, &
578 : G_bounds_1, R_bounds_2, R_bounds_3, &
579 : G_rads_1, R_rads_2, R_rads_3, &
580 0 : method=1)
581 0 : nresults = nresults + 1
582 : END IF
583 :
584 0 : IF (db_sum2) THEN
585 : CALL pgf_sum_3c_3d(S_G_no_2, la_max, lb_max, lc_max, RA, RB, RC, &
586 : zeta, zetb, zetc, aw(i_aw), hmat, h_inv, vol, &
587 : G_bounds_1, R_bounds_2, R_bounds_3, &
588 : G_rads_1, R_rads_2, R_rads_3, &
589 0 : method=2)
590 0 : nresults = nresults + 1
591 : END IF
592 :
593 0 : IF (db_sum3) THEN
594 : CALL pgf_sum_3c_3d(S_G_no_3, la_max, lb_max, lc_max, RA, RB, RC, &
595 : zeta, zetb, zetc, aw(i_aw), hmat, h_inv, vol, &
596 : G_bounds_1, R_bounds_2, R_bounds_3, &
597 : G_rads_1, R_rads_2, R_rads_3, &
598 0 : method=3)
599 0 : nresults = nresults + 1
600 : END IF
601 :
602 0 : max_error = 0.0_dp
603 0 : ALLOCATE (results_no(nresults))
604 :
605 0 : DO kco = ncoset(lc_min - 1) + 1, ncoset(lc_max)
606 0 : CALL get_l(kco, lc, cx, cy, cz)
607 0 : DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
608 0 : CALL get_l(jco, lb, bx, by, bz)
609 0 : DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
610 0 : CALL get_l(ico, la, ax, ay, az)
611 :
612 0 : max_error = 0.0_dp
613 0 : ir = 0
614 0 : IF (db_sum1) THEN
615 0 : ir = ir + 1
616 0 : results_no(ir) = S_G_no_1(ico, jco, kco)
617 : END IF
618 :
619 0 : IF (db_sum2) THEN
620 0 : ir = ir + 1
621 0 : results_no(ir) = S_G_no_2(ico, jco, kco)
622 : END IF
623 :
624 0 : IF (db_sum3) THEN
625 0 : ir = ir + 1
626 0 : results_no(ir) = S_G_no_3(ico, jco, kco)
627 : END IF
628 :
629 0 : max_result = MAXVAL(results_no)
630 0 : min_result = MINVAL(results_no)
631 0 : IF (nresults > 0) max_error = MAX(max_error, &
632 0 : (max_result - min_result)/(0.5_dp*(ABS(max_result) + ABS(min_result)) + 1.0_dp))
633 : END DO
634 : END DO
635 : END DO
636 :
637 0 : CPASSERT(max_error <= param%debug_delta)
638 0 : DEALLOCATE (results_no)
639 : END IF
640 : END DO
641 :
642 : ! assemble integral values
643 :
644 : CALL eri_mme_3c_collect_nonortho(vol, la_min, la_max, lb_min, lb_max, lc_min, lc_max, &
645 20400 : zeta, zetb, zetc, n_aw, aw, S_G_no, S_G_0_no, prefac, habc, o1, o2, o3)
646 :
647 61200 : END SUBROUTINE eri_mme_3c_integrate_nonortho
648 :
649 : ! **************************************************************************************************
650 : !> \brief ...
651 : !> \param vol ...
652 : !> \param la_min ...
653 : !> \param la_max ...
654 : !> \param lb_min ...
655 : !> \param lb_max ...
656 : !> \param lc_min ...
657 : !> \param lc_max ...
658 : !> \param zeta ...
659 : !> \param zetb ...
660 : !> \param zetc ...
661 : !> \param n_aw ...
662 : !> \param aw ...
663 : !> \param S_G ...
664 : !> \param S_G_0 ...
665 : !> \param prefac ...
666 : !> \param habc ...
667 : !> \param o1 ...
668 : !> \param o2 ...
669 : !> \param o3 ...
670 : ! **************************************************************************************************
671 134154 : SUBROUTINE eri_mme_3c_collect_ortho(vol, la_min, la_max, lb_min, lb_max, lc_min, lc_max, &
672 134154 : zeta, zetb, zetc, n_aw, aw, S_G, S_G_0, prefac, habc, o1, o2, o3)
673 : REAL(KIND=dp), INTENT(IN) :: vol
674 : INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max, lc_min, &
675 : lc_max
676 : REAL(KIND=dp), INTENT(IN) :: zeta, zetb, zetc
677 : INTEGER, INTENT(IN) :: n_aw
678 : REAL(KIND=dp), DIMENSION(2*n_aw), INTENT(IN) :: aw
679 : REAL(KIND=dp), DIMENSION(:, :, 0:, 0:, 0:), &
680 : INTENT(IN) :: S_G
681 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
682 : INTENT(IN) :: S_G_0
683 : REAL(KIND=dp), INTENT(IN) :: prefac
684 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: habc
685 : INTEGER, INTENT(IN) :: o1, o2, o3
686 :
687 : INTEGER :: ax, ay, az, bx, by, bz, cx, cy, cz, ico, &
688 : jco, kco, la, la_prev, lb, lb_prev, &
689 : lc, lc_prev
690 : REAL(KIND=dp) :: Imm, Ixyz_0, mone_lb, resc_a, &
691 : resc_a_init, resc_b, resc_b_init, &
692 : resc_c, resc_c_init
693 :
694 : ! Initialization of rescaling factors due to Hermite Gaussians
695 134154 : resc_a_init = (2.0_dp*zeta)**la_min
696 134154 : resc_b_init = (2.0_dp*zetb)**lb_min
697 134154 : resc_c_init = (2.0_dp*zetc)**lc_min
698 :
699 134154 : resc_c = resc_c_init
700 134154 : lc_prev = lc_min
701 903872 : DO kco = ncoset(lc_min - 1) + 1, ncoset(lc_max)
702 769718 : CALL get_l(kco, lc, cx, cy, cz)
703 769718 : IF (lc_prev < lc) resc_c = resc_c*(2.0_dp*zetc)
704 :
705 769718 : resc_b = resc_b_init
706 769718 : lb_prev = lb_min
707 6336708 : DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
708 5566990 : CALL get_l(jco, lb, bx, by, bz)
709 5566990 : mone_lb = (-1.0_dp)**lb
710 5566990 : IF (lb_prev < lb) resc_b = resc_b*(2.0_dp*zetb)
711 :
712 5566990 : resc_a = resc_a_init
713 5566990 : la_prev = la_min
714 47679670 : DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
715 42112680 : CALL get_l(ico, la, ax, ay, az)
716 :
717 42112680 : IF (la_prev < la) resc_a = resc_a*(2.0_dp*zeta)
718 42112680 : Ixyz_0 = 0.0_dp
719 42112680 : IF (lc == 0) THEN
720 : Ixyz_0 = S_G_0(ax + bx, 1)* &
721 : S_G_0(ay + by, 2)* &
722 : S_G_0(az + bz, 3) &
723 2781458 : /vol*mone_lb
724 : END IF
725 : Imm = SUM(aw(n_aw + 1:2*n_aw)*( &
726 : S_G(1:n_aw, 1, ax, bx, cx)* &
727 : S_G(1:n_aw, 2, ay, by, cy)* &
728 741910920 : S_G(1:n_aw, 3, az, bz, cz)) - Ixyz_0)
729 :
730 : ! rescaling needed due to Hermite Gaussians
731 42112680 : habc(o1 + ico, o2 + jco, o3 + kco) = Imm*prefac/(resc_a*resc_b*resc_c)
732 47679670 : la_prev = la
733 : END DO ! la
734 6336708 : lb_prev = lb
735 : END DO ! lb
736 903872 : lc_prev = lc
737 : END DO ! lc
738 :
739 134154 : END SUBROUTINE eri_mme_3c_collect_ortho
740 :
741 : ! **************************************************************************************************
742 : !> \brief ...
743 : !> \param vol ...
744 : !> \param la_min ...
745 : !> \param la_max ...
746 : !> \param lb_min ...
747 : !> \param lb_max ...
748 : !> \param lc_min ...
749 : !> \param lc_max ...
750 : !> \param zeta ...
751 : !> \param zetb ...
752 : !> \param zetc ...
753 : !> \param n_aw ...
754 : !> \param aw ...
755 : !> \param S_G ...
756 : !> \param S_G_0 ...
757 : !> \param prefac ...
758 : !> \param habc ...
759 : !> \param o1 ...
760 : !> \param o2 ...
761 : !> \param o3 ...
762 : ! **************************************************************************************************
763 20400 : PURE SUBROUTINE eri_mme_3c_collect_nonortho(vol, la_min, la_max, lb_min, lb_max, lc_min, lc_max, &
764 20400 : zeta, zetb, zetc, n_aw, aw, S_G, S_G_0, prefac, habc, o1, o2, o3)
765 : REAL(KIND=dp), INTENT(IN) :: vol
766 : INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max, lc_min, &
767 : lc_max
768 : REAL(KIND=dp), INTENT(IN) :: zeta, zetb, zetc
769 : INTEGER, INTENT(IN) :: n_aw
770 : REAL(KIND=dp), DIMENSION(2*n_aw), INTENT(IN) :: aw
771 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: S_G
772 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
773 : INTENT(IN) :: S_G_0
774 : REAL(KIND=dp), INTENT(IN) :: prefac
775 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: habc
776 : INTEGER, INTENT(IN) :: o1, o2, o3
777 :
778 : INTEGER :: ax, ay, az, bx, by, bz, cx, cy, cz, ico, &
779 : ijco, jco, kco, la, la_prev, lb, &
780 : lb_prev, lc, lc_prev
781 : REAL(KIND=dp) :: Imm, mone_lb, resc_a, resc_a_init, &
782 : resc_b, resc_b_init, resc_c, &
783 : resc_c_init
784 :
785 : ! Initialization of rescaling factors due to Hermite Gaussians
786 20400 : resc_a_init = (2.0_dp*zeta)**la_min
787 20400 : resc_b_init = (2.0_dp*zetb)**lb_min
788 20400 : resc_c_init = (2.0_dp*zetc)**lc_min
789 :
790 20400 : resc_c = resc_c_init
791 20400 : lc_prev = lc_min
792 100350 : DO kco = ncoset(lc_min - 1) + 1, ncoset(lc_max)
793 79950 : CALL get_l(kco, lc, cx, cy, cz)
794 79950 : IF (lc_prev < lc) resc_c = resc_c*(2.0_dp*zetc)
795 :
796 79950 : resc_b = resc_b_init
797 79950 : lb_prev = lb_min
798 505890 : DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
799 425940 : CALL get_l(jco, lb, bx, by, bz)
800 425940 : mone_lb = (-1.0_dp)**lb
801 425940 : IF (lb_prev < lb) resc_b = resc_b*(2.0_dp*zetb)
802 :
803 425940 : resc_a = resc_a_init
804 425940 : la_prev = la_min
805 4209543 : DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
806 3783603 : CALL get_l(ico, la, ax, ay, az)
807 :
808 3783603 : IF (la_prev < la) resc_a = resc_a*(2.0_dp*zeta)
809 3783603 : IF (lc > 0) THEN
810 3631431 : Imm = S_G(ico, jco, kco)
811 : ELSE
812 152172 : ijco = coset(ax + bx, ay + by, az + bz)
813 810828 : Imm = S_G(ico, jco, kco) - SUM(aw(n_aw + 1:2*n_aw))*S_G_0(ijco)/vol*mone_lb
814 : END IF
815 :
816 : ! rescaling needed due to Hermite Gaussians
817 3783603 : habc(o1 + ico, o2 + jco, o3 + kco) = Imm*prefac/(resc_a*resc_b*resc_c)
818 4209543 : la_prev = la
819 : END DO ! la
820 505890 : lb_prev = lb
821 : END DO ! lb
822 100350 : lc_prev = lc
823 : END DO ! lc
824 :
825 20400 : END SUBROUTINE eri_mme_3c_collect_nonortho
826 :
827 : END MODULE eri_mme_integrate
|
