(git:e7e05ae)
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-2024 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  SUBROUTINE eri_mme_2c_integrate(param, la_min, la_max, lb_min, lb_max, zeta, zetb, rab, &
71  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: s_g_1, s_g_2, s_g_no, s_g_no_h
93  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  REAL(kind=dp), DIMENSION(:), POINTER :: aw
97 
98  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  g_rad=g_rad)
103  cutoff = g_rad**2/2
104  CALL get_minimax_from_cutoff(param%minimax_grid, cutoff, n_aw, aw, grid)
105 
106  cpassert(grid .GT. 0)
107 
108  ! cell info
109  h_inv = param%h_inv
110  hmat = param%hmat
111  vol = param%vol
112 
113  IF (PRESENT(normalize)) THEN
114  norm = normalize
115  ELSE
116  norm = .false.
117  END IF
118 
119  l_max = la_max + lb_max
120 
121  IF (PRESENT(potential)) THEN
122  exact = .true.
123  ELSE
124  exact = .false.
125  END IF
126 
127  IF (exact) THEN
128  is_ortho = .false.
129  ELSE
130  is_ortho = param%is_ortho
131  END IF
132 
133  IF (is_ortho) THEN
134  ALLOCATE (s_g(0:l_max, 3, n_aw))
135  s_g = 0.0_dp
136 
137  IF (param%debug) THEN
138  ALLOCATE (s_g_1(0:l_max))
139  ALLOCATE (s_g_2(0:l_max))
140  END IF
141  ELSE
142  ALLOCATE (s_g_no(ncoset(l_max)))
143  s_g_no(:) = 0.0_dp
144  ALLOCATE (s_g_no_h(ncoset(l_max)))
145  END IF
146 
147  IF (exact) THEN
148  alpha_g = 0.25_dp/zeta + 0.25_dp/zetb
149  ! resolution for Gaussian width
150  g_res = 0.5_dp*param%G_min
151  r_res = 0.5_dp*param%R_min
152 
153  g_rad = exp_radius(l_max, alpha_g, 0.01*param%sum_precision, 1.0_dp, epsabs=g_res)
154  g_bounds(:) = ellipsoid_bounds(g_rad, transpose(hmat)/(2.0_dp*pi))
155  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  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  g_bounds, g_rad, r_bounds, r_rad)
163  alpha_g = aw(i_aw) + 0.25_dp/zeta + 0.25_dp/zetb
164  alpha_r = 0.25_dp/alpha_g
165  IF (is_ortho) THEN ! orthorhombic cell
166 
167  ! 1) precompute Ewald-like sum
168 
169  DO i_xyz = 1, 3
170  lgth = abs(hmat(i_xyz, i_xyz))
171  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  do_g_sum = n_sum_1d(1, i_xyz) < n_sum_1d(2, i_xyz) !G_bounds < R_bounds
175 
176  IF (do_g_sum) THEN
177  CALL pgf_sum_2c_gspace_1d(s_g(:, i_xyz, i_aw), -rab(i_xyz), alpha_g, inv_lgth, g_bounds(i_xyz))
178  IF (PRESENT(g_count)) g_count = g_count + 1
179  ELSE
180  CALL pgf_sum_2c_rspace_1d(s_g(:, i_xyz, i_aw), -rab(i_xyz), alpha_r, lgth, r_bounds(i_xyz))
181  IF (PRESENT(r_count)) r_count = r_count + 1
182  END IF
183 
184  IF (param%debug) THEN
185  ! check consistency of summation methods
186  CALL pgf_sum_2c_gspace_1d(s_g_1, -rab(i_xyz), alpha_g, inv_lgth, g_bounds(i_xyz))
187  CALL pgf_sum_2c_rspace_1d(s_g_2, -rab(i_xyz), alpha_r, lgth, r_bounds(i_xyz))
188  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  cpassert(max_error .LE. param%debug_delta)
191  END IF
192  END DO
193 
194  ELSE ! general cell
195 
196  do_g_sum = n_sum_3d(1) < n_sum_3d(2) !PRODUCT(2*R_bounds+1) .GT. PRODUCT(2*G_bounds+1)
197 
198  IF (do_g_sum) THEN
199  CALL pgf_sum_2c_gspace_3d(s_g_no_h, l_max, -rab, alpha_g, h_inv, g_bounds, g_rad, vol)
200  IF (PRESENT(g_count)) g_count = g_count + 1
201  ELSE
202  CALL pgf_sum_2c_rspace_3d(s_g_no_h, l_max, -rab, alpha_r, hmat, h_inv, r_bounds, r_rad)
203  IF (PRESENT(r_count)) r_count = r_count + 1
204  END IF
205  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  prefac = sqrt(1.0_dp/(zeta*zetb))
212 
213  ! 2) Assemble integral values from Ewald sums
214  DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
215  CALL get_l(jco, lb, bx, by, bz)
216  lb_xyz = [bx, by, bz]
217  DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
218  CALL get_l(ico, la, ax, ay, az)
219  la_xyz = [ax, ay, az]
220  IF (is_ortho) THEN
221  imm = 0.0_dp
222  DO i_aw = 1, n_aw
223  ixyz = 1.0_dp
224  DO i_xyz = 1, 3
225  ixyz = ixyz*s_g(la_xyz(i_xyz) + lb_xyz(i_xyz), i_xyz, i_aw)*prefac
226  END DO
227  imm = imm + aw(n_aw + i_aw)*ixyz
228  END DO
229  ELSE
230  imm = s_g_no(coset(ax + bx, ay + by, az + bz))*prefac**3
231  END IF
232  IF (la + lb .EQ. 0 .AND. .NOT. exact) THEN
233  imm = imm - sum(aw(n_aw + 1:2*n_aw))*prefac**3/vol ! subtracting G = 0 term
234  END IF
235  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  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  hermite_gauss_norm(zetb, lb_xyz)
243  END IF
244  END DO ! la
245  END DO ! lb
246 
247  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  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  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  cpassert(param%is_valid)
284  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  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  habc, o1, o2, o3, gg_count, gr_count, rr_count)
291 
292  END IF
293  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  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: s_g_0
333  REAL(kind=dp), ALLOCATABLE, &
334  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  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  param%sum_precision, g_rads_1=g_rads_1)
345 
346  cutoff = (min(g_rads_1(1), g_rads_1(2) + g_rads_1(3)))**2/2
347 
348  CALL get_minimax_from_cutoff(param%minimax_grid, cutoff, n_aw, aw, grid)
349 
350  cpassert(grid .GT. 0)
351 
352  ! minimax coeffs
353  n_aw = param%minimax_grid(grid)%n_minimax
354  aw => param%minimax_grid(grid)%minimax_aw
355 
356  ! cell info
357  h_inv = param%h_inv
358  hmat = param%hmat
359  vol = param%vol
360 
361  ! prefactor for integral values (unnormalized Hermite Gaussians)
362  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  g_res = 0.5_dp*param%G_min
366  r_res = 0.5_dp*param%R_min
367 
368  ALLOCATE (s_g(n_aw, 3, 0:la_max, 0:lb_max, 0:lc_max))
369 
370  ! G= 0 components
371  IF (lc_min == 0) THEN
372  ALLOCATE (s_g_0(0:la_max + lb_max, 3))
373 
374  alpha_r_0 = zeta*zetb/(zeta + zetb)
375  lmax_0 = la_max + lb_max
376  r_rad_0 = exp_radius(lmax_0, alpha_r_0, param%sum_precision, 1.0_dp, epsabs=r_res)
377 
378  lgth = abs(hmat(1, 1))
379  CALL pgf_sum_2c_rspace_1d(s_g_0(:, 1), rb(1) - ra(1), alpha_r_0, lgth, r_rad_0/lgth)
380  lgth = abs(hmat(2, 2))
381  CALL pgf_sum_2c_rspace_1d(s_g_0(:, 2), rb(2) - ra(2), alpha_r_0, lgth, r_rad_0/lgth)
382  lgth = abs(hmat(3, 3))
383  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  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  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  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  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  r_bounds_3(:, 3))
401 
402  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  zeta, zetb, zetc, n_aw, aw, s_g, s_g_0, prefac, habc, o1, o2, o3)
407 
408  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  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: results_no, s_g_0_no
456  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: s_g_no, s_g_no_1, s_g_no_2, s_g_no_3, &
457  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  REAL(kind=dp), DIMENSION(:), POINTER :: aw
464 
465  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  param%sum_precision, g_rads_1=g_rads_1)
471 
472  cutoff = (min(g_rads_1(1), g_rads_1(2) + g_rads_1(3)))**2/2
473 
474  CALL get_minimax_from_cutoff(param%minimax_grid, cutoff, n_aw, aw, grid)
475 
476  cpassert(grid .GT. 0)
477 
478  ! minimax coeffs
479  n_aw = param%minimax_grid(grid)%n_minimax
480  aw => param%minimax_grid(grid)%minimax_aw
481 
482  ! cell info
483  h_inv = param%h_inv
484  hmat = param%hmat
485  vol = param%vol
486 
487  ! prefactor for integral values (unnormalized Hermite Gaussians)
488  prefac = (zeta*zetb*zetc)**(-1.5_dp)*pi**(11.0_dp/2.0_dp)*4.0_dp
489 
490  IF (param%debug) THEN
491  ALLOCATE (s_g_no_1(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
492  ALLOCATE (s_g_no_2(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
493  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  g_res = 0.5_dp*param%G_min
498  r_res = 0.5_dp*param%R_min
499 
500  ALLOCATE (s_g_no(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
501 
502  s_g_no(:, :, :) = 0.0_dp
503  IF (param%debug) THEN
504  s_g_no_1(:, :, :) = -1.0_dp
505  s_g_no_2(:, :, :) = -1.0_dp
506  s_g_no_3(:, :, :) = -1.0_dp
507  END IF
508  ALLOCATE (s_g_no_h(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
509 
510  ! G= 0 components
511  IF (lc_min == 0) THEN
512  ALLOCATE (s_g_0_no(ncoset(la_max + lb_max)))
513  alpha_g_0 = 0.25_dp/zetb + 0.25_dp/zeta
514  alpha_r_0 = 0.25_dp/alpha_g_0
515  lmax_0 = la_max + lb_max
516  r_0 = rb - ra
517  g_rad_0 = exp_radius(lmax_0, alpha_g_0, param%sum_precision, 1.0_dp, epsabs=g_res)
518  r_rad_0 = exp_radius(lmax_0, alpha_r_0, param%sum_precision, 1.0_dp, epsabs=r_res)
519  g_bound_0 = ellipsoid_bounds(g_rad_0, transpose(hmat)/(2.0_dp*pi))
520  r_bound_0 = ellipsoid_bounds(r_rad_0, h_inv)
521  do_g_sum_0 = product(2*r_bound_0 + 1) .GT. product(2*g_bound_0 + 1)
522  IF (do_g_sum_0) THEN
523  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  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  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  g_bounds_1, g_rads_1, r_bounds_2, r_rads_2, r_bounds_3, r_rads_3)
533  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  method=sum_method, method_out=method)
540  s_g_no(:, :, :) = s_g_no(:, :, :) + aw(n_aw + i_aw)*s_g_no_h(:, :, :)
541 
542  SELECT CASE (method)
543  CASE (1)
544  IF (PRESENT(gg_count)) gg_count = gg_count + 1
545  CASE (2)
546  IF (PRESENT(gr_count)) gr_count = gr_count + 1
547  CASE (3)
548  IF (PRESENT(rr_count)) rr_count = rr_count + 1
549  CASE DEFAULT
550  cpabort("")
551  END SELECT
552 
553  IF (param%debug) THEN
554  nresults = 0
555  ! check consistency of summation methods
556 
557  db_sum1 = (n_sum_3d(1)) .LT. int(param%debug_nsum, kind=int_8)
558  db_sum2 = (n_sum_3d(2)) .LT. int(param%debug_nsum, kind=int_8)
559  db_sum3 = (n_sum_3d(3)) .LT. int(param%debug_nsum, kind=int_8)
560 
561  IF (param%unit_nr > 0) THEN
562  WRITE (param%unit_nr, *) "ERI_MME DEBUG | number of summands (GG / GR / RR)", n_sum_3d
563  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  s_g_no_1(:, :, :) = 0.0_dp
567  s_g_no_2(:, :, :) = 0.0_dp
568  s_g_no_3(:, :, :) = 0.0_dp
569 
570  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  method=1)
576  nresults = nresults + 1
577  END IF
578 
579  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  method=2)
585  nresults = nresults + 1
586  END IF
587 
588  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  method=3)
594  nresults = nresults + 1
595  END IF
596 
597  max_error = 0.0_dp
598  ALLOCATE (results_no(nresults))
599 
600  DO kco = ncoset(lc_min - 1) + 1, ncoset(lc_max)
601  CALL get_l(kco, lc, cx, cy, cz)
602  DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
603  CALL get_l(jco, lb, bx, by, bz)
604  DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
605  CALL get_l(ico, la, ax, ay, az)
606 
607  max_error = 0.0_dp
608  ir = 0
609  IF (db_sum1) THEN
610  ir = ir + 1
611  results_no(ir) = s_g_no_1(ico, jco, kco)
612  END IF
613 
614  IF (db_sum2) THEN
615  ir = ir + 1
616  results_no(ir) = s_g_no_2(ico, jco, kco)
617  END IF
618 
619  IF (db_sum3) THEN
620  ir = ir + 1
621  results_no(ir) = s_g_no_3(ico, jco, kco)
622  END IF
623 
624  max_result = maxval(results_no)
625  min_result = minval(results_no)
626  IF (nresults > 0) max_error = max(max_error, &
627  (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  cpassert(max_error .LE. param%debug_delta)
633  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  zeta, zetb, zetc, n_aw, aw, s_g_no, s_g_0_no, prefac, habc, o1, o2, o3)
641 
642  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  SUBROUTINE eri_mme_3c_collect_ortho(vol, la_min, la_max, lb_min, lb_max, lc_min, lc_max, &
667  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  resc_a_init = (2.0_dp*zeta)**la_min
691  resc_b_init = (2.0_dp*zetb)**lb_min
692  resc_c_init = (2.0_dp*zetc)**lc_min
693 
694  resc_c = resc_c_init
695  lc_prev = lc_min
696  DO kco = ncoset(lc_min - 1) + 1, ncoset(lc_max)
697  CALL get_l(kco, lc, cx, cy, cz)
698  IF (lc_prev < lc) resc_c = resc_c*(2.0_dp*zetc)
699 
700  resc_b = resc_b_init
701  lb_prev = lb_min
702  DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
703  CALL get_l(jco, lb, bx, by, bz)
704  mone_lb = (-1.0_dp)**lb
705  IF (lb_prev < lb) resc_b = resc_b*(2.0_dp*zetb)
706 
707  resc_a = resc_a_init
708  la_prev = la_min
709  DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
710  CALL get_l(ico, la, ax, ay, az)
711 
712  IF (la_prev < la) resc_a = resc_a*(2.0_dp*zeta)
713  ixyz_0 = 0.0_dp
714  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  /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  s_g(1:n_aw, 3, az, bz, cz)) - ixyz_0)
724 
725  ! rescaling needed due to Hermite Gaussians
726  habc(o1 + ico, o2 + jco, o3 + kco) = imm*prefac/(resc_a*resc_b*resc_c)
727  la_prev = la
728  END DO ! la
729  lb_prev = lb
730  END DO ! lb
731  lc_prev = lc
732  END DO ! lc
733 
734  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  PURE SUBROUTINE eri_mme_3c_collect_nonortho(vol, la_min, la_max, lb_min, lb_max, lc_min, lc_max, &
759  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  resc_a_init = (2.0_dp*zeta)**la_min
782  resc_b_init = (2.0_dp*zetb)**lb_min
783  resc_c_init = (2.0_dp*zetc)**lc_min
784 
785  resc_c = resc_c_init
786  lc_prev = lc_min
787  DO kco = ncoset(lc_min - 1) + 1, ncoset(lc_max)
788  CALL get_l(kco, lc, cx, cy, cz)
789  IF (lc_prev < lc) resc_c = resc_c*(2.0_dp*zetc)
790 
791  resc_b = resc_b_init
792  lb_prev = lb_min
793  DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
794  CALL get_l(jco, lb, bx, by, bz)
795  mone_lb = (-1.0_dp)**lb
796  IF (lb_prev < lb) resc_b = resc_b*(2.0_dp*zetb)
797 
798  resc_a = resc_a_init
799  la_prev = la_min
800  DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
801  CALL get_l(ico, la, ax, ay, az)
802 
803  IF (la_prev < la) resc_a = resc_a*(2.0_dp*zeta)
804  IF (lc .GT. 0) THEN
805  imm = s_g(ico, jco, kco)
806  ELSE
807  ijco = coset(ax + bx, ay + by, az + bz)
808  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  habc(o1 + ico, o2 + jco, o3 + kco) = imm*prefac/(resc_a*resc_b*resc_c)
813  la_prev = la
814  END DO ! la
815  lb_prev = lb
816  END DO ! lb
817  lc_prev = lc
818  END DO ! lc
819 
820  END SUBROUTINE
821 
822 END 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_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:72
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:274
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_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:736
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_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:120272
eri_mme_lattice_summation::get_l
pure elemental subroutine, public get_l(lco, l, lx, ly, lz)
...
Definition: eri_mme_lattice_summation.F:121093
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::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::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_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::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:121050
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::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_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:421
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
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