Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Methods aiming for error estimate and automatic cutoff calibration.
10 : !> integrals.
11 : !> \par History
12 : !> 2015 09 created
13 : !> \author Patrick Seewald
14 : ! **************************************************************************************************
15 :
16 : MODULE eri_mme_error_control
17 : USE ao_util, ONLY: exp_radius
18 : USE eri_mme_gaussian, ONLY: get_minimax_coeff_v_gspace,&
19 : hermite_gauss_norm
20 : USE eri_mme_lattice_summation, ONLY: pgf_sum_2c_gspace_1d_deltal
21 : USE kinds, ONLY: dp
22 : USE mathconstants, ONLY: pi,&
23 : twopi
24 : USE message_passing, ONLY: mp_para_env_type
25 : #include "../base/base_uses.f90"
26 :
27 : IMPLICIT NONE
28 :
29 : PRIVATE
30 :
31 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
32 :
33 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eri_mme_error_control'
34 :
35 : PUBLIC :: calibrate_cutoff, cutoff_minimax_error, minimax_error, cutoff_error
36 : CONTAINS
37 :
38 : ! **************************************************************************************************
39 : !> \brief Find optimal cutoff minimizing errors due to minimax approximation and
40 : !> due to finite cutoff using bisection on the difference of the errors
41 : !> \param hmat ...
42 : !> \param h_inv ...
43 : !> \param G_min ...
44 : !> \param vol ...
45 : !> \param zet_min Minimum exponent
46 : !> \param l_mm Total ang. mom. quantum number
47 : !> \param zet_max Max. exponents to estimate cutoff error
48 : !> \param l_max_zet Max. total ang. mom. quantum numbers to estimate cutoff error
49 : !> \param n_minimax Number of terms in minimax approximation
50 : !> \param cutoff_l Initial guess of lower bound for cutoff
51 : !> \param cutoff_r Initial guess of upper bound for cutoff
52 : !> \param tol Tolerance (cutoff precision)
53 : !> \param delta to modify initial guess interval
54 : !> \param cutoff Best cutoff
55 : !> \param err_mm Minimax error
56 : !> \param err_c Cutoff error
57 : !> \param C_mm Scaling constant to generalize AM-GM upper bound estimate to
58 : !> minimax approx.
59 : !> \param para_env ...
60 : !> \param print_calib ...
61 : !> \param unit_nr ...
62 : ! **************************************************************************************************
63 84 : SUBROUTINE calibrate_cutoff(hmat, h_inv, G_min, vol, zet_min, l_mm, zet_max, l_max_zet, &
64 : n_minimax, cutoff_l, cutoff_r, tol, delta, &
65 : cutoff, err_mm, err_c, C_mm, para_env, print_calib, unit_nr)
66 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv
67 : REAL(KIND=dp), INTENT(IN) :: G_min
68 : REAL(KIND=dp) :: vol
69 : REAL(KIND=dp), INTENT(IN) :: zet_min
70 : INTEGER, INTENT(IN) :: l_mm
71 : REAL(KIND=dp), INTENT(IN) :: zet_max
72 : INTEGER, INTENT(IN) :: l_max_zet, n_minimax
73 : REAL(KIND=dp), INTENT(IN) :: cutoff_l, cutoff_r, tol, delta
74 : REAL(KIND=dp), INTENT(OUT) :: cutoff, err_mm, err_c, C_mm
75 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
76 : LOGICAL, INTENT(IN) :: print_calib
77 : INTEGER, INTENT(IN) :: unit_nr
78 :
79 : INTEGER :: i, iter1, iter2, max_iter
80 : LOGICAL :: do_print, valid_initial
81 : REAL(KIND=dp) :: cutoff_mid, delta_c_mid, delta_mm_mid
82 84 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: minimax_aw
83 : REAL(KIND=dp), DIMENSION(2) :: cutoff_lr, delta_c, delta_mm
84 :
85 84 : do_print = unit_nr > 0 .AND. print_calib
86 : IF (do_print) THEN
87 0 : WRITE (unit_nr, '(/T2, A)') "ERI_MME| Basis set parameters for estimating minimax error"
88 0 : WRITE (unit_nr, '(T2, A, T67, ES12.2, 1X, I1)') "ERI_MME| exp, l:", zet_min, l_mm
89 0 : WRITE (unit_nr, '(T2, A)') "ERI_MME| Basis set parameters for estimating cutoff error"
90 0 : WRITE (unit_nr, '(T2, A, T67, ES12.2, 1X, I1)') "ERI_MME| exp, l:", zet_max, l_max_zet
91 : END IF
92 :
93 84 : max_iter = 100
94 :
95 84 : IF ((cutoff_r - cutoff_l)/(0.5_dp*(cutoff_r + cutoff_l)) <= tol) THEN
96 : CALL cp_abort(__LOCATION__, "difference of boundaries for cutoff "// &
97 0 : "(MAX - MIN) must be greater than cutoff precision.")
98 : END IF
99 :
100 84 : IF ((delta >= 1.0_dp) .OR. (delta <= 0.0_dp)) THEN
101 : CALL cp_abort(__LOCATION__, &
102 0 : "relative delta to modify initial cutoff interval (DELTA) must be in (0, 1)")
103 : END IF
104 :
105 84 : cutoff_lr(1) = cutoff_l
106 84 : cutoff_lr(2) = cutoff_r
107 :
108 252 : ALLOCATE (minimax_aw(2*n_minimax))
109 :
110 84 : IF (do_print) THEN
111 0 : WRITE (unit_nr, '(/T2, A)') "ERI_MME| Calibrating cutoff by bisecting error(minimax) - error(cutoff)"
112 0 : WRITE (unit_nr, '(T2, A, T72, ES9.2)') "ERI_MME| Rel. cutoff precision", tol
113 0 : WRITE (unit_nr, '(T2, A, T77, F4.1)') "ERI_MME| Rel. cutoff delta to modify initial interval", delta
114 : END IF
115 :
116 : ! 1) find valid initial values for bisection
117 84 : DO iter1 = 1, max_iter + 1
118 84 : IF (iter1 > max_iter) THEN
119 : CALL cp_abort(__LOCATION__, &
120 : "Maximum number of iterations in bisection to determine initial "// &
121 0 : "cutoff interval has been exceeded.")
122 : END IF
123 :
124 84 : cutoff_lr(1) = MAX(cutoff_lr(1), 0.5_dp*G_min**2)
125 : ! approx.) is hit
126 :
127 252 : DO i = 1, 2
128 : CALL cutoff_minimax_error(cutoff_lr(i), hmat, h_inv, vol, G_min, zet_min, l_mm, zet_max, l_max_zet, &
129 252 : n_minimax, minimax_aw, delta_mm(i), delta_c(i), C_mm, para_env)
130 : END DO
131 :
132 84 : valid_initial = .TRUE.
133 84 : IF ((delta_mm(1) - delta_c(1)) > 0) THEN
134 0 : cutoff_lr(1) = cutoff_lr(1)*(1.0_dp - ABS(delta))
135 0 : valid_initial = .FALSE.
136 : END IF
137 84 : IF ((delta_mm(2) - delta_c(2)) < 0) THEN
138 0 : cutoff_lr(2) = cutoff_lr(2)*(1.0_dp + ABS(delta))
139 : valid_initial = .FALSE.
140 : END IF
141 :
142 84 : IF (valid_initial) EXIT
143 : END DO
144 :
145 : ! 2) bisection to find cutoff s.t. err_minimax(cutoff) - err_cutoff(cutoff) = 0
146 84 : IF (do_print) WRITE (unit_nr, '(/T2, A)') &
147 0 : "ERI_MME| Step, cutoff (min, max, mid), err(minimax), err(cutoff), err diff"
148 :
149 1190 : DO iter2 = 1, max_iter + 1
150 1190 : IF (iter2 > max_iter) THEN
151 : CALL cp_abort(__LOCATION__, &
152 0 : "Maximum number of iterations in bisection to determine cutoff has been exceeded")
153 : END IF
154 :
155 1190 : cutoff_mid = 0.5_dp*(cutoff_lr(1) + cutoff_lr(2))
156 : CALL cutoff_minimax_error(cutoff_mid, hmat, h_inv, vol, G_min, zet_min, l_mm, zet_max, l_max_zet, &
157 1190 : n_minimax, minimax_aw, delta_mm_mid, delta_c_mid, C_mm, para_env)
158 1190 : IF (do_print) WRITE (unit_nr, '(T11, I2, F11.1, F11.1, F11.1, 3X, ES9.2, 3X, ES9.2, 3X, ES9.2)') &
159 0 : iter2, cutoff_lr(1), cutoff_lr(2), cutoff_mid, &
160 0 : delta_mm_mid, delta_c_mid, delta_mm_mid - delta_c_mid
161 :
162 1190 : IF ((cutoff_lr(2) - cutoff_lr(1))/cutoff_mid < tol) EXIT
163 2380 : IF (delta_mm_mid - delta_c_mid > 0) THEN
164 776 : cutoff_lr(2) = cutoff_mid
165 : delta_mm(2) = delta_mm_mid
166 : delta_c(2) = delta_c_mid
167 : ELSE
168 330 : cutoff_lr(1) = cutoff_mid
169 : delta_mm(1) = delta_mm_mid
170 : delta_c(1) = delta_c_mid
171 : END IF
172 : END DO
173 84 : err_mm = delta_mm_mid
174 84 : err_c = delta_c_mid
175 84 : cutoff = cutoff_mid
176 :
177 84 : IF (do_print) THEN
178 0 : WRITE (unit_nr, '(/T2, A)') "ERI_MME| Cutoff calibration number of steps:"
179 0 : WRITE (unit_nr, '(T2, A, T79, I2)') "ERI_MME| Steps for initial interval", iter1 - 1
180 0 : WRITE (unit_nr, '(T2, A, T79, I2/)') "ERI_MME| Bisection iteration steps", iter2 - 1
181 : END IF
182 :
183 84 : END SUBROUTINE calibrate_cutoff
184 :
185 : ! **************************************************************************************************
186 : !> \brief Compute upper bounds for the errors of 2-center ERI's (P|P) due
187 : !> to minimax approximation and due to finite cutoff, where P is a
188 : !> normalized Hermite Gaussian.
189 : !> \param cutoff ...
190 : !> \param hmat ...
191 : !> \param h_inv ...
192 : !> \param vol ...
193 : !> \param G_min ...
194 : !> \param zet_min Exponent of P to estimate minimax error
195 : !> \param l_mm total ang. mom. quantum number of P to estimate minimax error
196 : !> \param zet_max Max. exponents of P to estimate cutoff error
197 : !> \param l_max_zet Max. total ang. mom. quantum numbers of P to estimate cutoff error
198 : !> \param n_minimax Number of terms in minimax approximation
199 : !> \param minimax_aw Minimax coefficients
200 : !> \param err_mm Minimax error
201 : !> \param err_ctff Cutoff error
202 : !> \param C_mm Scaling constant to generalize AM-GM upper bound estimate to
203 : !> minimax approx.
204 : !> \param para_env ...
205 : ! **************************************************************************************************
206 1396 : SUBROUTINE cutoff_minimax_error(cutoff, hmat, h_inv, vol, G_min, zet_min, l_mm, zet_max, l_max_zet, &
207 1396 : n_minimax, minimax_aw, err_mm, err_ctff, C_mm, para_env)
208 : REAL(KIND=dp), INTENT(IN) :: cutoff
209 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv
210 : REAL(KIND=dp), INTENT(IN) :: vol, G_min, zet_min
211 : INTEGER, INTENT(IN) :: l_mm
212 : REAL(KIND=dp), INTENT(IN) :: zet_max
213 : INTEGER, INTENT(IN) :: l_max_zet, n_minimax
214 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: minimax_aw
215 : REAL(KIND=dp), INTENT(OUT) :: err_mm, err_ctff, C_mm
216 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
217 :
218 : REAL(KIND=dp) :: delta_mm
219 :
220 : CALL minimax_error(cutoff, hmat, vol, G_min, zet_min, l_mm, &
221 1396 : n_minimax, minimax_aw, err_mm, delta_mm)
222 : CALL cutoff_error(cutoff, h_inv, G_min, zet_max, l_max_zet, &
223 1396 : n_minimax, minimax_aw, err_ctff, C_mm, para_env)
224 :
225 1396 : END SUBROUTINE cutoff_minimax_error
226 :
227 : ! **************************************************************************************************
228 : !> \brief Minimax error, simple analytical formula
229 : !> Note minimax error may blow up for small exponents. This is also observed numerically,
230 : !> but in this case, error estimate is no upper bound.
231 : !> \param cutoff ...
232 : !> \param hmat ...
233 : !> \param vol ...
234 : !> \param G_min ...
235 : !> \param zet_min Exponent of P to estimate minimax error
236 : !> \param l_mm total ang. mom. quantum number of P to estimate minimax error
237 : !> \param n_minimax Number of terms in minimax approximation
238 : !> \param minimax_aw Minimax coefficients
239 : !> \param err_mm Minimax error
240 : !> \param delta_mm ...
241 : !> \param potential ...
242 : !> \param pot_par ...
243 : ! **************************************************************************************************
244 86888 : SUBROUTINE minimax_error(cutoff, hmat, vol, G_min, zet_min, l_mm, &
245 86888 : n_minimax, minimax_aw, err_mm, delta_mm, potential, pot_par)
246 : REAL(KIND=dp), INTENT(IN) :: cutoff
247 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat
248 : REAL(KIND=dp), INTENT(IN) :: vol, G_min, zet_min
249 : INTEGER, INTENT(IN) :: l_mm, n_minimax
250 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: minimax_aw
251 : REAL(KIND=dp), INTENT(OUT) :: err_mm, delta_mm
252 : INTEGER, INTENT(IN), OPTIONAL :: potential
253 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par
254 :
255 : INTEGER :: i_xyz
256 : REAL(KIND=dp) :: prod_mm_k
257 :
258 : CALL get_minimax_coeff_v_gspace(n_minimax, cutoff, G_min, minimax_aw(:), &
259 86888 : potential=potential, pot_par=pot_par, err_minimax=delta_mm)
260 :
261 86888 : prod_mm_k = 1.0_dp
262 347552 : DO i_xyz = 1, 3
263 : prod_mm_k = prod_mm_k*(ABS(hmat(i_xyz, i_xyz))/twopi + &
264 348032 : MERGE(SQRT(2.0_dp/(zet_min*pi))*EXP(-1.0_dp), 0.0_dp, l_mm > 0))
265 : END DO
266 86888 : err_mm = 32*pi**4/vol*delta_mm*prod_mm_k
267 :
268 86888 : END SUBROUTINE minimax_error
269 :
270 : ! **************************************************************************************************
271 : !> \brief Cutoff error, estimating G > G_c part of Ewald sum by using C/3 * 1/(Gx^2*Gy^2*Gz^2)^1/3 as an
272 : !> upper bound for 1/G^2 (AM-GM inequality) and its minimax approximation (factor C).
273 : !>
274 : !> Note: usually, minimax approx. falls off faster than 1/G**2, so C should be approximately 1.
275 : !> The error is calculated for all l up to l_max and golden section search algorithm is
276 : !> applied to find the exponent that maximizes cutoff error.
277 : !> \param cutoff ...
278 : !> \param h_inv ...
279 : !> \param G_min ...
280 : !> \param zet_max Max. exponents of P to estimate cutoff error
281 : !> \param l_max_zet Max. total ang. mom. quantum numbers of P to estimate cutoff error
282 : !> \param n_minimax Number of terms in minimax approximation
283 : !> \param minimax_aw Minimax coefficients
284 : !> \param err_ctff Cutoff error
285 : !> \param C_mm Scaling constant to generalize AM-GM upper bound estimate to
286 : !> minimax approx.
287 : !> \param para_env ...
288 : ! **************************************************************************************************
289 1396 : SUBROUTINE cutoff_error(cutoff, h_inv, G_min, zet_max, l_max_zet, &
290 1396 : n_minimax, minimax_aw, err_ctff, C_mm, para_env)
291 : REAL(KIND=dp), INTENT(IN) :: cutoff
292 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: h_inv
293 : REAL(KIND=dp), INTENT(IN) :: G_min, zet_max
294 : INTEGER, INTENT(IN) :: l_max_zet, n_minimax
295 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: minimax_aw
296 : REAL(KIND=dp), INTENT(OUT) :: err_ctff, C_mm
297 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
298 :
299 : INTEGER :: i_aw, iG, iter, max_iter, nG
300 : REAL(KIND=dp) :: C, dG, eps_zet, err0, err1, err_c, err_ctff_curr, err_ctff_prev, err_d, G, &
301 : G_1, G_c, gr, zet_a, zet_b, zet_c, zet_d, zet_div, zet_max_tmp
302 :
303 : ! parameters for finding exponent maximizing cutoff error
304 :
305 1396 : eps_zet = 1.0E-05_dp ! tolerance for exponent
306 1396 : zet_div = 2.0_dp ! sampling constant for finding initial values of exponents
307 1396 : max_iter = 100 ! maximum number of iterations in golden section search
308 1396 : G_c = SQRT(2.0*cutoff)
309 :
310 1396 : zet_max_tmp = zet_max
311 :
312 : ! 2) Cutoff error, estimating G > G_c part of Ewald sum by using C/3 * 1/(Gx^2*Gy^2*Gz^2)^1/3 as an
313 : ! upper bound for 1/G^2 (AM-GM inequality) and its minimax approximation (factor C).
314 : ! Note: usually, minimax approx. falls off faster than 1/G**2, so C should be approximately 1.
315 : ! The error is calculated for all l up to l_max and golden section search algorithm is
316 : ! applied to find the exponent that maximizes cutoff error.
317 24190 : G_1 = SQRT(1.0_dp/(3.0_dp*MINVAL(minimax_aw(1:n_minimax))))
318 :
319 1396 : C_mm = 0.0_dp
320 1396 : IF (G_1 > G_c) THEN
321 762 : nG = 1000
322 762 : dG = (G_1 - G_c)/nG
323 762 : G = G_c
324 762762 : DO iG = 1, nG
325 762000 : G = MIN(G, G_c)
326 762000 : C = 0.0_dp
327 20478000 : DO i_aw = 1, n_minimax
328 20478000 : C = C + 3.0_dp*minimax_aw(n_minimax + i_aw)*EXP(-3.0_dp*minimax_aw(i_aw)*G**2)*G**2
329 : END DO
330 762000 : C_mm = MAX(C, C_mm)
331 762762 : G = G + dG
332 : END DO
333 : ELSE
334 3712 : DO i_aw = 1, n_minimax
335 3712 : C_mm = C_mm + 3.0_dp*minimax_aw(n_minimax + i_aw)*EXP(-3.0_dp*minimax_aw(i_aw)*G_c**2)*G_c**2
336 : END DO
337 : END IF
338 1396 : C = MAX(1.0_dp, C_mm)
339 :
340 1396 : err_ctff_prev = 0.0_dp
341 1396 : gr = 0.5_dp*(SQRT(5.0_dp) - 1.0_dp) ! golden ratio
342 : ! Find valid starting values for golden section search
343 2754 : DO iter = 1, max_iter + 1
344 2754 : IF (iter > max_iter) THEN
345 : CALL cp_abort(__LOCATION__, "Maximum number of iterations for finding "// &
346 0 : "exponent maximizing cutoff error has been exceeded.")
347 : END IF
348 :
349 2754 : CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_max_tmp, C, err_ctff_curr, para_env)
350 2754 : IF (err_ctff_prev >= err_ctff_curr) THEN
351 1396 : zet_a = zet_max_tmp
352 1396 : zet_b = MIN(zet_max_tmp*zet_div**2, zet_max)
353 1396 : EXIT
354 : ELSE
355 1358 : err_ctff_prev = err_ctff_curr
356 : END IF
357 4112 : zet_max_tmp = zet_max_tmp/zet_div
358 : END DO
359 :
360 : ! Golden section search
361 1396 : zet_c = zet_b - gr*(zet_b - zet_a)
362 1396 : zet_d = zet_a + gr*(zet_b - zet_a)
363 23210 : DO iter = 1, max_iter + 1
364 23210 : IF (ABS(zet_c - zet_d) < eps_zet*(zet_a + zet_b)) THEN
365 1396 : CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_a, C, err0, para_env)
366 1396 : CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_b, C, err1, para_env)
367 1396 : err_ctff_curr = MAX(err0, err1)
368 1396 : EXIT
369 : END IF
370 21814 : CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_c, C, err_c, para_env)
371 21814 : CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_d, C, err_d, para_env)
372 21814 : IF (err_c > err_d) THEN
373 2132 : zet_b = zet_d
374 2132 : zet_d = zet_c
375 2132 : zet_c = zet_b - gr*(zet_b - zet_a)
376 : ELSE
377 19682 : zet_a = zet_c
378 19682 : zet_c = zet_d
379 19682 : zet_d = zet_a + gr*(zet_b - zet_a)
380 : END IF
381 : END DO
382 1396 : err_ctff = err_ctff_curr
383 :
384 1396 : END SUBROUTINE cutoff_error
385 :
386 : ! **************************************************************************************************
387 : !> \brief Calculate cutoff error estimate by using C_mm/3 * 1/(Gx^2*Gy^2*Gz^2)^1/3
388 : !> as an upper bound for 1/G^2 (and its minimax approximation) for |G| > G_c.
389 : !> Error is referring to a basis function P with fixed exponent zet_max and
390 : !> max. angular momentum l_max_zet.
391 : !> \param cutoff ...
392 : !> \param h_inv ...
393 : !> \param G_min ...
394 : !> \param l_max_zet ...
395 : !> \param zet_max ...
396 : !> \param C_mm ...
397 : !> \param err_c ...
398 : !> \param para_env ...
399 : ! **************************************************************************************************
400 49174 : SUBROUTINE cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_max, C_mm, err_c, para_env)
401 : REAL(KIND=dp), INTENT(IN) :: cutoff
402 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: h_inv
403 : REAL(KIND=dp), INTENT(IN) :: G_min
404 : INTEGER, INTENT(IN) :: l_max_zet
405 : REAL(KIND=dp), INTENT(IN) :: zet_max, C_mm
406 : REAL(KIND=dp), INTENT(OUT) :: err_c
407 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
408 :
409 : INTEGER :: ax, ay, az, G_l, G_u, Gl_first, Gl_last, &
410 : Gu_first, Gu_last, i_xyz, l, my_p, &
411 : n_Gl, n_Gl_left, n_Gl_p, n_Gu, &
412 : n_Gu_left, n_Gu_p, n_p
413 : REAL(KIND=dp) :: alpha_G, eps_G, err_c_l, G_c, G_rad, &
414 : G_res, inv_lgth, prefactor, sum_G_diff
415 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: S_G_l, S_G_u
416 :
417 49174 : G_c = SQRT(2.0_dp*cutoff)
418 49174 : eps_G = TINY(eps_G) ! sum up to machine precision
419 49174 : G_res = 0.5_dp*G_min ! resolution for screening
420 :
421 49174 : err_c = 0.0_dp
422 49174 : alpha_G = 1.0_dp/(2.0_dp*zet_max)
423 49174 : prefactor = 1.0_dp/zet_max
424 :
425 196696 : ALLOCATE (S_G_l(0:2*l_max_zet, 3))
426 98348 : ALLOCATE (S_G_u(0:2*l_max_zet, 3))
427 :
428 49174 : G_rad = exp_radius(2*l_max_zet, alpha_G, eps_G, prefactor, epsabs=G_res)
429 :
430 : ! Parallelization of sum over G vectors
431 49174 : my_p = para_env%mepos ! mpi rank
432 49174 : n_p = para_env%num_pe ! total number of processes
433 :
434 196696 : DO i_xyz = 1, 3
435 147522 : inv_lgth = ABS(h_inv(i_xyz, i_xyz))
436 :
437 147522 : G_l = FLOOR(G_c/(inv_lgth*twopi))
438 147522 : G_u = FLOOR(G_rad/(inv_lgth*twopi))
439 :
440 147522 : IF (G_u < G_l) G_u = G_l
441 :
442 : ! Serial code:
443 : ! !Sum |G| <= G_c
444 : ! CALL pgf_sum_2c_gspace_1d_deltal(S_G_l(:, i_xyz), alpha_G, inv_lgth, -G_l, G_l, &
445 : ! 2.0_dp/3.0_dp, prefactor)
446 : ! !Sum |G| > G_c
447 : ! CALL pgf_sum_2c_gspace_1d_deltal(S_G_u(:, i_xyz), alpha_G, inv_lgth, G_l + 1, G_u, &
448 : ! 2.0_dp/3.0_dp, prefactor)
449 :
450 : ! Parallel code:
451 147522 : n_Gu = MAX((G_u - G_l), 0)
452 147522 : n_Gl = 2*G_l + 1
453 147522 : n_Gu_p = n_Gu/n_p
454 147522 : n_Gl_p = n_Gl/n_p
455 147522 : n_Gu_left = MOD(n_Gu, n_p)
456 147522 : n_Gl_left = MOD(n_Gl, n_p)
457 :
458 147522 : IF (my_p < n_Gu_left) THEN
459 34831 : Gu_first = G_l + 1 + (n_Gu_p + 1)*my_p
460 34831 : Gu_last = G_l + 1 + (n_Gu_p + 1)*(my_p + 1) - 1
461 : ELSE
462 112691 : Gu_first = G_l + 1 + n_Gu_left + n_Gu_p*my_p
463 112691 : Gu_last = G_l + 1 + n_Gu_left + n_Gu_p*(my_p + 1) - 1
464 : END IF
465 :
466 147522 : IF (my_p < n_Gl_left) THEN
467 73761 : Gl_first = -G_l + (n_Gl_p + 1)*my_p
468 73761 : Gl_last = -G_l + (n_Gl_p + 1)*(my_p + 1) - 1
469 : ELSE
470 73761 : Gl_first = -G_l + n_Gl_left + n_Gl_p*my_p
471 73761 : Gl_last = -G_l + n_Gl_left + n_Gl_p*(my_p + 1) - 1
472 : END IF
473 :
474 : ! Sum |G| <= G_c
475 : CALL pgf_sum_2c_gspace_1d_deltal(S_G_l(:, i_xyz), alpha_G, inv_lgth, Gl_first, Gl_last, &
476 147522 : 2.0_dp/3.0_dp, prefactor)
477 :
478 : ! Sum |G| > G_c
479 : CALL pgf_sum_2c_gspace_1d_deltal(S_G_u(:, i_xyz), alpha_G, inv_lgth, Gu_first, Gu_last, &
480 196696 : 2.0_dp/3.0_dp, prefactor)
481 : END DO
482 :
483 49174 : CALL para_env%sum(S_G_l)
484 49174 : CALL para_env%sum(S_G_u)
485 :
486 879226 : S_G_u = S_G_u*2.0_dp ! to include negative values of G
487 :
488 187516 : DO l = 0, l_max_zet
489 482956 : DO ax = 0, l
490 994330 : DO ay = 0, l - ax
491 560548 : az = l - ax - ay
492 :
493 : ! Compute prod_k (S_G_l(l_k,k) + S_G_u(l_k,k)) - prod_k (S_G_l(l_k,k)) with k in {x, y, z}
494 : ! Note: term by term multiplication to avoid subtraction for numerical stability
495 : sum_G_diff = S_G_u(2*ax, 1)*S_G_u(2*ay, 2)*S_G_u(2*az, 3) + &
496 : S_G_u(2*ax, 1)*S_G_u(2*ay, 2)*S_G_l(2*az, 3) + &
497 : S_G_u(2*ax, 1)*S_G_l(2*ay, 2)*S_G_u(2*az, 3) + &
498 : S_G_l(2*ax, 1)*S_G_u(2*ay, 2)*S_G_u(2*az, 3) + &
499 : S_G_u(2*ax, 1)*S_G_l(2*ay, 2)*S_G_l(2*az, 3) + &
500 : S_G_l(2*ax, 1)*S_G_u(2*ay, 2)*S_G_l(2*az, 3) + &
501 560548 : S_G_l(2*ax, 1)*S_G_l(2*ay, 2)*S_G_u(2*az, 3)
502 :
503 : err_c_l = 4.0_dp*pi**4*hermite_gauss_norm(zet_max, [ax, ay, az])**2* &
504 2242192 : C_mm/3.0_dp*sum_G_diff
505 :
506 855988 : err_c = MAX(err_c, err_c_l)
507 : END DO
508 : END DO
509 : END DO
510 :
511 49174 : DEALLOCATE (S_G_u, S_G_l)
512 :
513 49174 : END SUBROUTINE cutoff_error_fixed_exp
514 :
515 : END MODULE eri_mme_error_control
|
