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