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