(git:4978e50)
Loading...
Searching...
No Matches
eri_mme_integrate.F
Go to the documentation of this file.
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
23MODULE 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
50CONTAINS
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 SUBROUTINE eri_mme_2c_integrate(param, la_min, la_max, lb_min, lb_max, zeta, zetb, rab, &
72 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: s_g_1, s_g_2, s_g_no, s_g_no_h
96 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 REAL(kind=dp), DIMENSION(:), POINTER :: aw
100
101 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 g_rad=g_rad)
106 cutoff = g_rad**2/2
107 CALL get_minimax_from_cutoff(param%minimax_grid, cutoff, n_aw, aw, grid)
108
109 cpassert(grid > 0)
110
111 ! cell info
112 h_inv = param%h_inv
113 hmat = param%hmat
114 vol = param%vol
115
116 IF (PRESENT(normalize)) THEN
117 norm = normalize
118 ELSE
119 norm = .false.
120 END IF
121
122 l_max = la_max + lb_max
123
124 IF (PRESENT(potential)) THEN
125 exact = .true.
126 ELSE
127 exact = .false.
128 END IF
129
130 IF (exact) THEN
131 is_ortho = .false.
132 ELSE
133 is_ortho = param%is_ortho
134 END IF
135
136 IF (is_ortho) THEN
137 ALLOCATE (s_g(0:l_max, 3, n_aw))
138 s_g = 0.0_dp
139
140 IF (param%debug) THEN
141 ALLOCATE (s_g_1(0:l_max))
142 ALLOCATE (s_g_2(0:l_max))
143 END IF
144 ELSE
145 ALLOCATE (s_g_no(ncoset(l_max)))
146 s_g_no(:) = 0.0_dp
147 ALLOCATE (s_g_no_h(ncoset(l_max)))
148 END IF
149
150 IF (exact) THEN
151 alpha_g = 0.25_dp/zeta + 0.25_dp/zetb
152 ! resolution for Gaussian width
153 g_res = 0.5_dp*param%G_min
154 r_res = 0.5_dp*param%R_min
155
156 g_rad = exp_radius(l_max, alpha_g, rel_prec*param%sum_precision, 1.0_dp, epsabs=g_res)
157 g_bounds(:) = ellipsoid_bounds(g_rad, transpose(hmat)/twopi)
158 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 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 g_bounds, g_rad, r_bounds, r_rad)
166 alpha_g = aw(i_aw) + 0.25_dp/zeta + 0.25_dp/zetb
167 alpha_r = 0.25_dp/alpha_g
168 IF (is_ortho) THEN ! orthorhombic cell
169
170 ! 1) precompute Ewald-like sum
171
172 DO i_xyz = 1, 3
173 lgth = abs(hmat(i_xyz, i_xyz))
174 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 do_g_sum = n_sum_1d(1, i_xyz) < n_sum_1d(2, i_xyz) !G_bounds < R_bounds
178
179 IF (do_g_sum) THEN
180 CALL pgf_sum_2c_gspace_1d(s_g(:, i_xyz, i_aw), -rab(i_xyz), alpha_g, inv_lgth, g_bounds(i_xyz))
181 IF (PRESENT(g_count)) g_count = g_count + 1
182 ELSE
183 CALL pgf_sum_2c_rspace_1d(s_g(:, i_xyz, i_aw), -rab(i_xyz), alpha_r, lgth, r_bounds(i_xyz))
184 IF (PRESENT(r_count)) r_count = r_count + 1
185 END IF
186
187 IF (param%debug) THEN
188 ! check consistency of summation methods
189 CALL pgf_sum_2c_gspace_1d(s_g_1, -rab(i_xyz), alpha_g, inv_lgth, g_bounds(i_xyz))
190 CALL pgf_sum_2c_rspace_1d(s_g_2, -rab(i_xyz), alpha_r, lgth, r_bounds(i_xyz))
191 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 cpassert(max_error <= param%debug_delta)
194 END IF
195 END DO
196
197 ELSE ! general cell
198
199 do_g_sum = n_sum_3d(1) < n_sum_3d(2) !PRODUCT(2*R_bounds+1) > PRODUCT(2*G_bounds+1)
200
201 IF (do_g_sum) THEN
202 CALL pgf_sum_2c_gspace_3d(s_g_no_h, l_max, -rab, alpha_g, h_inv, g_bounds, g_rad, vol)
203 IF (PRESENT(g_count)) g_count = g_count + 1
204 ELSE
205 CALL pgf_sum_2c_rspace_3d(s_g_no_h, l_max, -rab, alpha_r, hmat, h_inv, r_bounds, r_rad)
206 IF (PRESENT(r_count)) r_count = r_count + 1
207 END IF
208 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 prefac = sqrt(1.0_dp/(zeta*zetb))
215
216 ! 2) Assemble integral values from Ewald sums
217 DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
218 CALL get_l(jco, lb, bx, by, bz)
219 lb_xyz = [bx, by, bz]
220 DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
221 CALL get_l(ico, la, ax, ay, az)
222 la_xyz = [ax, ay, az]
223 IF (is_ortho) THEN
224 imm = 0.0_dp
225 DO i_aw = 1, n_aw
226 ixyz = 1.0_dp
227 DO i_xyz = 1, 3
228 ixyz = ixyz*s_g(la_xyz(i_xyz) + lb_xyz(i_xyz), i_xyz, i_aw)*prefac
229 END DO
230 imm = imm + aw(n_aw + i_aw)*ixyz
231 END DO
232 ELSE
233 imm = s_g_no(coset(ax + bx, ay + by, az + bz))*prefac**3
234 END IF
235 IF (la + lb == 0 .AND. .NOT. exact) THEN
236 imm = imm - sum(aw(n_aw + 1:2*n_aw))*prefac**3/vol ! subtracting G = 0 term
237 END IF
238 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 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 hermite_gauss_norm(zetb, lb_xyz)
246 END IF
247 END DO ! la
248 END DO ! lb
249
250 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 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 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 cpassert(param%is_valid)
287 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 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 habc, o1, o2, o3, gg_count, gr_count, rr_count)
294
295 END IF
296 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 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: s_g_0
336 REAL(kind=dp), ALLOCATABLE, &
337 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 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 param%sum_precision, g_rads_1=g_rads_1)
348
349 cutoff = (min(g_rads_1(1), g_rads_1(2) + g_rads_1(3)))**2/2
350
351 CALL get_minimax_from_cutoff(param%minimax_grid, cutoff, n_aw, aw, grid)
352
353 cpassert(grid > 0)
354
355 ! minimax coeffs
356 n_aw = param%minimax_grid(grid)%n_minimax
357 aw => param%minimax_grid(grid)%minimax_aw
358
359 ! cell info
360 h_inv = param%h_inv
361 hmat = param%hmat
362 vol = param%vol
363
364 ! prefactor for integral values (unnormalized Hermite Gaussians)
365 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 g_res = 0.5_dp*param%G_min
369 r_res = 0.5_dp*param%R_min
370
371 ALLOCATE (s_g(n_aw, 3, 0:la_max, 0:lb_max, 0:lc_max))
372
373 ! G= 0 components
374 IF (lc_min == 0) THEN
375 ALLOCATE (s_g_0(0:la_max + lb_max, 3))
376
377 alpha_r_0 = zeta*zetb/(zeta + zetb)
378 lmax_0 = la_max + lb_max
379 r_rad_0 = exp_radius(lmax_0, alpha_r_0, param%sum_precision, 1.0_dp, epsabs=r_res)
380
381 lgth = abs(hmat(1, 1))
382 CALL pgf_sum_2c_rspace_1d(s_g_0(:, 1), rb(1) - ra(1), alpha_r_0, lgth, r_rad_0/lgth)
383 lgth = abs(hmat(2, 2))
384 CALL pgf_sum_2c_rspace_1d(s_g_0(:, 2), rb(2) - ra(2), alpha_r_0, lgth, r_rad_0/lgth)
385 lgth = abs(hmat(3, 3))
386 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 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 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 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 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 r_bounds_3(:, 3))
404
405 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 zeta, zetb, zetc, n_aw, aw, s_g, s_g_0, prefac, habc, o1, o2, o3)
410
411 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 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: results_no, s_g_0_no
459 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: s_g_no, s_g_no_1, s_g_no_2, s_g_no_3, &
460 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 REAL(kind=dp), DIMENSION(:), POINTER :: aw
467
468 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 param%sum_precision, g_rads_1=g_rads_1)
474
475 cutoff = (min(g_rads_1(1), g_rads_1(2) + g_rads_1(3)))**2/2
476
477 CALL get_minimax_from_cutoff(param%minimax_grid, cutoff, n_aw, aw, grid)
478
479 cpassert(grid > 0)
480
481 ! minimax coeffs
482 n_aw = param%minimax_grid(grid)%n_minimax
483 aw => param%minimax_grid(grid)%minimax_aw
484
485 ! cell info
486 h_inv = param%h_inv
487 hmat = param%hmat
488 vol = param%vol
489
490 ! prefactor for integral values (unnormalized Hermite Gaussians)
491 prefac = (zeta*zetb*zetc)**(-1.5_dp)*pi**(11.0_dp/2.0_dp)*4.0_dp
492
493 IF (param%debug) THEN
494 ALLOCATE (s_g_no_1(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
495 ALLOCATE (s_g_no_2(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
496 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 g_res = 0.5_dp*param%G_min
501 r_res = 0.5_dp*param%R_min
502
503 ALLOCATE (s_g_no(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
504
505 s_g_no(:, :, :) = 0.0_dp
506 IF (param%debug) THEN
507 s_g_no_1(:, :, :) = -1.0_dp
508 s_g_no_2(:, :, :) = -1.0_dp
509 s_g_no_3(:, :, :) = -1.0_dp
510 END IF
511 ALLOCATE (s_g_no_h(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
512
513 ! G= 0 components
514 IF (lc_min == 0) THEN
515 ALLOCATE (s_g_0_no(ncoset(la_max + lb_max)))
516 alpha_g_0 = 0.25_dp/zetb + 0.25_dp/zeta
517 alpha_r_0 = 0.25_dp/alpha_g_0
518 lmax_0 = la_max + lb_max
519 r_0 = rb - ra
520 g_rad_0 = exp_radius(lmax_0, alpha_g_0, param%sum_precision, 1.0_dp, epsabs=g_res)
521 r_rad_0 = exp_radius(lmax_0, alpha_r_0, param%sum_precision, 1.0_dp, epsabs=r_res)
522 g_bound_0 = ellipsoid_bounds(g_rad_0, transpose(hmat)/(2.0_dp*pi))
523 r_bound_0 = ellipsoid_bounds(r_rad_0, h_inv)
524 do_g_sum_0 = product(2*r_bound_0 + 1) > product(2*g_bound_0 + 1)
525 IF (do_g_sum_0) THEN
526 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 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 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 g_bounds_1, g_rads_1, r_bounds_2, r_rads_2, r_bounds_3, r_rads_3)
536 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 method=sum_method, method_out=method)
543 s_g_no(:, :, :) = s_g_no(:, :, :) + aw(n_aw + i_aw)*s_g_no_h(:, :, :)
544
545 SELECT CASE (method)
546 CASE (1)
547 IF (PRESENT(gg_count)) gg_count = gg_count + 1
548 CASE (2)
549 IF (PRESENT(gr_count)) gr_count = gr_count + 1
550 CASE (3)
551 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 "for eri_mme_3c_integrate_nonortho")
556 END SELECT
557
558 IF (param%debug) THEN
559 nresults = 0
560 ! check consistency of summation methods
561
562 db_sum1 = (n_sum_3d(1)) < int(param%debug_nsum, kind=int_8)
563 db_sum2 = (n_sum_3d(2)) < int(param%debug_nsum, kind=int_8)
564 db_sum3 = (n_sum_3d(3)) < int(param%debug_nsum, kind=int_8)
565
566 IF (param%unit_nr > 0) THEN
567 WRITE (param%unit_nr, *) "ERI_MME DEBUG | number of summands (GG / GR / RR)", n_sum_3d
568 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 s_g_no_1(:, :, :) = 0.0_dp
572 s_g_no_2(:, :, :) = 0.0_dp
573 s_g_no_3(:, :, :) = 0.0_dp
574
575 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 method=1)
581 nresults = nresults + 1
582 END IF
583
584 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 method=2)
590 nresults = nresults + 1
591 END IF
592
593 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 method=3)
599 nresults = nresults + 1
600 END IF
601
602 max_error = 0.0_dp
603 ALLOCATE (results_no(nresults))
604
605 DO kco = ncoset(lc_min - 1) + 1, ncoset(lc_max)
606 CALL get_l(kco, lc, cx, cy, cz)
607 DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
608 CALL get_l(jco, lb, bx, by, bz)
609 DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
610 CALL get_l(ico, la, ax, ay, az)
611
612 max_error = 0.0_dp
613 ir = 0
614 IF (db_sum1) THEN
615 ir = ir + 1
616 results_no(ir) = s_g_no_1(ico, jco, kco)
617 END IF
618
619 IF (db_sum2) THEN
620 ir = ir + 1
621 results_no(ir) = s_g_no_2(ico, jco, kco)
622 END IF
623
624 IF (db_sum3) THEN
625 ir = ir + 1
626 results_no(ir) = s_g_no_3(ico, jco, kco)
627 END IF
628
629 max_result = maxval(results_no)
630 min_result = minval(results_no)
631 IF (nresults > 0) max_error = max(max_error, &
632 (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 cpassert(max_error <= param%debug_delta)
638 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 zeta, zetb, zetc, n_aw, aw, s_g_no, s_g_0_no, prefac, habc, o1, o2, o3)
646
647 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 SUBROUTINE eri_mme_3c_collect_ortho(vol, la_min, la_max, lb_min, lb_max, lc_min, lc_max, &
672 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 resc_a_init = (2.0_dp*zeta)**la_min
696 resc_b_init = (2.0_dp*zetb)**lb_min
697 resc_c_init = (2.0_dp*zetc)**lc_min
698
699 resc_c = resc_c_init
700 lc_prev = lc_min
701 DO kco = ncoset(lc_min - 1) + 1, ncoset(lc_max)
702 CALL get_l(kco, lc, cx, cy, cz)
703 IF (lc_prev < lc) resc_c = resc_c*(2.0_dp*zetc)
704
705 resc_b = resc_b_init
706 lb_prev = lb_min
707 DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
708 CALL get_l(jco, lb, bx, by, bz)
709 mone_lb = (-1.0_dp)**lb
710 IF (lb_prev < lb) resc_b = resc_b*(2.0_dp*zetb)
711
712 resc_a = resc_a_init
713 la_prev = la_min
714 DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
715 CALL get_l(ico, la, ax, ay, az)
716
717 IF (la_prev < la) resc_a = resc_a*(2.0_dp*zeta)
718 ixyz_0 = 0.0_dp
719 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 /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 s_g(1:n_aw, 3, az, bz, cz)) - ixyz_0)
729
730 ! rescaling needed due to Hermite Gaussians
731 habc(o1 + ico, o2 + jco, o3 + kco) = imm*prefac/(resc_a*resc_b*resc_c)
732 la_prev = la
733 END DO ! la
734 lb_prev = lb
735 END DO ! lb
736 lc_prev = lc
737 END DO ! lc
738
739 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 PURE SUBROUTINE eri_mme_3c_collect_nonortho(vol, la_min, la_max, lb_min, lb_max, lc_min, lc_max, &
764 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 resc_a_init = (2.0_dp*zeta)**la_min
787 resc_b_init = (2.0_dp*zetb)**lb_min
788 resc_c_init = (2.0_dp*zetc)**lc_min
789
790 resc_c = resc_c_init
791 lc_prev = lc_min
792 DO kco = ncoset(lc_min - 1) + 1, ncoset(lc_max)
793 CALL get_l(kco, lc, cx, cy, cz)
794 IF (lc_prev < lc) resc_c = resc_c*(2.0_dp*zetc)
795
796 resc_b = resc_b_init
797 lb_prev = lb_min
798 DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
799 CALL get_l(jco, lb, bx, by, bz)
800 mone_lb = (-1.0_dp)**lb
801 IF (lb_prev < lb) resc_b = resc_b*(2.0_dp*zetb)
802
803 resc_a = resc_a_init
804 la_prev = la_min
805 DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
806 CALL get_l(ico, la, ax, ay, az)
807
808 IF (la_prev < la) resc_a = resc_a*(2.0_dp*zeta)
809 IF (lc > 0) THEN
810 imm = s_g(ico, jco, kco)
811 ELSE
812 ijco = coset(ax + bx, ay + by, az + bz)
813 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 habc(o1 + ico, o2 + jco, o3 + kco) = imm*prefac/(resc_a*resc_b*resc_c)
818 la_prev = la
819 END DO ! la
820 lb_prev = lb
821 END DO ! lb
822 lc_prev = lc
823 END DO ! lc
824
825 END SUBROUTINE eri_mme_3c_collect_nonortho
826
827END MODULE eri_mme_integrate
ao_util
All kind of helpful little routines.
Definition ao_util.F:14
ao_util::exp_radius
real(kind=dp) function, public exp_radius(l, alpha, threshold, prefactor, epsabs, epsrel, rlow)
The radius of a primitive Gaussian function for a given threshold is calculated. g(r) = prefactor*r**...
Definition ao_util.F:96
eri_mme_gaussian
Methods related to properties of Hermite and Cartesian Gaussian functions.
Definition eri_mme_gaussian.F:15
eri_mme_gaussian::hermite_gauss_norm
pure real(kind=dp) function, public hermite_gauss_norm(zet, l)
Norm of 1d Hermite-Gauss functions.
Definition eri_mme_gaussian.F:76
eri_mme_integrate
Minimax-Ewald (MME) method for calculating 2-center and 3-center electron repulsion integrals (ERI) o...
Definition eri_mme_integrate.F:23
eri_mme_integrate::eri_mme_3c_integrate
subroutine, public eri_mme_3c_integrate(param, la_min, la_max, lb_min, lb_max, lc_min, lc_max, zeta, zetb, zetc, ra, rb, rc, habc, o1, o2, o3, gg_count, gr_count, rr_count)
Low-level integration routine for 3-center ERIs.
Definition eri_mme_integrate.F:277
eri_mme_integrate::eri_mme_2c_integrate
subroutine, public eri_mme_2c_integrate(param, la_min, la_max, lb_min, lb_max, zeta, zetb, rab, hab, o1, o2, g_count, r_count, normalize, potential, pot_par)
Low-level integration routine for 2-center ERIs.
Definition eri_mme_integrate.F:73
eri_mme_lattice_summation
Ewald sums to represent integrals in direct and reciprocal lattice.
Definition eri_mme_lattice_summation.F:15
eri_mme_lattice_summation::pgf_sum_2c_gspace_1d
pure subroutine, public pgf_sum_2c_gspace_1d(s_g, r, alpha, inv_lgth, g_c)
Compute Ewald-like sum for 2-center ERIs in G space in 1 dimension S_G(l, alpha) = (-i)^l*inv_lgth*su...
Definition eri_mme_lattice_summation.F:359
eri_mme_lattice_summation::eri_mme_2c_get_bounds
subroutine, public eri_mme_2c_get_bounds(hmat, h_inv, vol, is_ortho, g_min, r_min, la_max, lb_max, zeta, zetb, a_mm, sum_precision, n_sum_1d, n_sum_3d, g_bounds, g_rad, r_bounds, r_rad)
Get summation bounds for 2c integrals.
Definition eri_mme_lattice_summation.F:177
eri_mme_lattice_summation::get_l
pure elemental subroutine, public get_l(lco, l, lx, ly, lz)
...
Definition eri_mme_lattice_summation.F:121094
eri_mme_lattice_summation::pgf_sum_3c_3d
subroutine, public pgf_sum_3c_3d(s_g, la_max, lb_max, lc_max, ra, rb, rc, zeta, zetb, zetc, a_mm, hmat, h_inv, vol, g_bounds_1, r_bounds_2, r_bounds_3, g_rads_1, r_rads_2, r_rads_3, method, method_out, order)
As pgf_sum_3c_1d but 3d sum required for non-orthorhombic cells.
Definition eri_mme_lattice_summation.F:120273
eri_mme_lattice_summation::pgf_sum_2c_rspace_3d
pure subroutine, public pgf_sum_2c_rspace_3d(s_r, l_max, r, alpha, hmat, h_inv, r_c, r_rad)
As pgf_sum_2c_rspace_1d but 3d sum required for non-orthorhombic cells.
Definition eri_mme_lattice_summation.F:570
eri_mme_lattice_summation::pgf_sum_2c_gspace_3d
pure subroutine, public pgf_sum_2c_gspace_3d(s_g, l_max, r, alpha, h_inv, g_c, g_rad, vol, potential, pot_par)
As pgf_sum_2c_gspace_1d but 3d sum required for non-orthorhombic cells.
Definition eri_mme_lattice_summation.F:414
eri_mme_lattice_summation::eri_mme_3c_get_rads
subroutine, public eri_mme_3c_get_rads(la_max, lb_max, lc_max, zeta, zetb, zetc, a_mm, g_min, r_min, sum_precision, g_rads_1, r_rads_2, r_rads_3)
Get summation radii for 3c integrals.
Definition eri_mme_lattice_summation.F:111
eri_mme_lattice_summation::pgf_sum_3c_1d
subroutine, public pgf_sum_3c_1d(s_g, ra, rb, rc, zeta, zetb, zetc, a_mm, lgth, r_bounds_3)
Compute Ewald-like sum for 3-center integrals in 1 dimension S_G(l, m, n, alpha, beta,...
Definition eri_mme_lattice_summation.F:737
eri_mme_lattice_summation::eri_mme_2c_get_rads
subroutine, public eri_mme_2c_get_rads(la_max, lb_max, zeta, zetb, a_mm, g_min, r_min, sum_precision, g_rad, r_rad)
Get summation radii for 2c integrals.
Definition eri_mme_lattice_summation.F:73
eri_mme_lattice_summation::ellipsoid_bounds
pure real(kind=dp) function, dimension(3), public ellipsoid_bounds(s_rad, s_to_e)
Compute bounding box for ellipsoid. This is needed in order to find summation bounds for sphere for s...
Definition eri_mme_lattice_summation.F:121051
eri_mme_lattice_summation::eri_mme_3c_get_bounds
subroutine, public eri_mme_3c_get_bounds(hmat, h_inv, vol, is_ortho, g_min, r_min, la_max, lb_max, lc_max, zeta, zetb, zetc, a_mm, sum_precision, n_sum_1d, n_sum_3d, g_bounds_1, g_rads_1, r_bounds_2, r_rads_2, r_bounds_3, r_rads_3)
Get summation bounds for 3c integrals.
Definition eri_mme_lattice_summation.F:244
eri_mme_lattice_summation::pgf_sum_2c_rspace_1d
pure subroutine, public pgf_sum_2c_rspace_1d(s_r, r, alpha, lgth, r_c)
Compute Ewald-like sum for 2-center ERIs in R space in 1 dimension S_R(l, alpha) = SQRT(alpha/pi) sum...
Definition eri_mme_lattice_summation.F:514
eri_mme_types
Types and initialization / release routines for Minimax-Ewald method for electron repulsion integrals...
Definition eri_mme_types.F:16
eri_mme_types::get_minimax_from_cutoff
subroutine, public get_minimax_from_cutoff(grids, cutoff, n_minimax, minimax_aw, grid_no)
...
Definition eri_mme_types.F:426
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::coset
integer, dimension(:, :, :), allocatable, public coset
Definition orbital_pointers.F:45
eri_mme_types::eri_mme_param
Definition eri_mme_types.F:63