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