Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Interface to Minimax-Ewald method for periodic ERI's to be used in CP2K.
10 : ! **************************************************************************************************
11 :
12 : MODULE cp_eri_mme_interface
13 : USE basis_set_types, ONLY: gto_basis_set_type
14 : USE cell_methods, ONLY: cell_create,&
15 : init_cell
16 : USE cell_types, ONLY: cell_release,&
17 : cell_type,&
18 : pbc
19 : USE cp_log_handling, ONLY: cp_get_default_logger,&
20 : cp_logger_type
21 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
22 : cp_print_key_section_create,&
23 : cp_print_key_unit_nr,&
24 : medium_print_level
25 : USE eri_mme_test, ONLY: eri_mme_2c_perf_acc_test,&
26 : eri_mme_3c_perf_acc_test
27 : USE eri_mme_types, ONLY: eri_mme_coulomb,&
28 : eri_mme_init,&
29 : eri_mme_longrange,&
30 : eri_mme_param,&
31 : eri_mme_print_grid_info,&
32 : eri_mme_release,&
33 : eri_mme_set_params,&
34 : eri_mme_yukawa
35 : USE input_keyword_types, ONLY: keyword_create,&
36 : keyword_release,&
37 : keyword_type
38 : USE input_section_types, ONLY: section_add_keyword,&
39 : section_add_subsection,&
40 : section_create,&
41 : section_release,&
42 : section_type,&
43 : section_vals_get_subs_vals,&
44 : section_vals_type,&
45 : section_vals_val_get
46 : USE input_val_types, ONLY: real_t
47 : USE kinds, ONLY: default_string_length,&
48 : dp
49 : USE message_passing, ONLY: mp_para_env_type
50 : USE orbital_pointers, ONLY: init_orbital_pointers
51 : USE qs_kind_types, ONLY: get_qs_kind,&
52 : qs_kind_type
53 : USE string_utilities, ONLY: s2a
54 : #include "./base/base_uses.f90"
55 :
56 : IMPLICIT NONE
57 :
58 : PRIVATE
59 :
60 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
61 :
62 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_eri_mme_interface'
63 :
64 : PUBLIC :: &
65 : cp_eri_mme_finalize, &
66 : cp_eri_mme_init_read_input, &
67 : cp_eri_mme_param, &
68 : cp_eri_mme_perf_acc_test, &
69 : cp_eri_mme_set_params, &
70 : cp_eri_mme_update_local_counts, &
71 : create_eri_mme_section, &
72 : create_eri_mme_test_section
73 :
74 : INTERFACE cp_eri_mme_set_params
75 : MODULE PROCEDURE eri_mme_set_params_from_basis
76 : MODULE PROCEDURE eri_mme_set_params_custom
77 : END INTERFACE
78 :
79 : TYPE cp_eri_mme_param
80 : TYPE(cp_logger_type), POINTER :: logger => NULL()
81 : ! There is a bug with some older compilers preventing a default initialization of derived types with allocatable components
82 : #if defined(FTN_NO_DEFAULT_INIT)
83 : TYPE(eri_mme_param) :: par = eri_mme_param(minimax_grid=NULL())
84 : #else
85 : TYPE(eri_mme_param) :: par = eri_mme_param()
86 : #endif
87 : TYPE(section_vals_type), &
88 : POINTER :: mme_section => NULL()
89 : INTEGER :: G_count_2c = 0, R_count_2c = 0
90 : INTEGER :: GG_count_3c = 0, GR_count_3c = 0, RR_count_3c = 0
91 : LOGICAL :: do_calib = .FALSE.
92 : END TYPE cp_eri_mme_param
93 :
94 : CONTAINS
95 : ! **************************************************************************************************
96 : !> \brief Create main input section
97 : !> \param section ...
98 : !> \param default_n_minimax ...
99 : ! **************************************************************************************************
100 83190 : SUBROUTINE create_eri_mme_section(section, default_n_minimax)
101 : TYPE(section_type), POINTER :: section
102 : INTEGER, INTENT(IN), OPTIONAL :: default_n_minimax
103 :
104 : INTEGER :: my_default_n_minimax
105 : TYPE(keyword_type), POINTER :: keyword
106 : TYPE(section_type), POINTER :: print_key, subsection
107 :
108 83190 : NULLIFY (keyword, print_key, subsection)
109 83190 : CPASSERT(.NOT. ASSOCIATED(section))
110 :
111 83190 : IF (PRESENT(default_n_minimax)) THEN
112 9238 : my_default_n_minimax = default_n_minimax
113 : ELSE
114 73952 : my_default_n_minimax = 20
115 : END IF
116 :
117 : CALL section_create(section, __LOCATION__, name="ERI_MME", &
118 : description="Parameters for the calculation of periodic electron repulsion "// &
119 : "integrals (ERI) using the Minimax-Ewald (MME) method. "// &
120 : "Note: N_MINIMAX is the only parameter to be tuned for accuracy, "// &
121 : "all other parameters can be left to default. MME method is faster "// &
122 : "than numerical GPW.", &
123 83190 : n_keywords=5, n_subsections=1)
124 :
125 : CALL keyword_create(keyword, __LOCATION__, &
126 : name="N_MINIMAX", &
127 : description="Number of terms in minimax approximation of "// &
128 : "reciprocal space potential. ", &
129 83190 : default_i_val=my_default_n_minimax)
130 83190 : CALL section_add_keyword(section, keyword)
131 83190 : CALL keyword_release(keyword)
132 :
133 : CALL keyword_create(keyword, __LOCATION__, &
134 : name="CUTOFF", &
135 : description="User-defined energy cutoff to be used only if "// &
136 : "DO_CALIBRATE_CUTOFF is set to .FALSE. ", &
137 83190 : default_r_val=300.0_dp)
138 83190 : CALL section_add_keyword(section, keyword)
139 83190 : CALL keyword_release(keyword)
140 :
141 : CALL keyword_create(keyword, __LOCATION__, &
142 : name="SUM_PRECISION", &
143 : description="Terms in lattice sums are ignored if absolute value smaller than this value.", &
144 83190 : default_r_val=1.0E-16_dp)
145 83190 : CALL section_add_keyword(section, keyword)
146 83190 : CALL keyword_release(keyword)
147 :
148 : CALL keyword_create(keyword, __LOCATION__, &
149 : name="DO_CALIBRATE_CUTOFF", &
150 : description="Whether the energy cutoff shall be calibrated to "// &
151 : "minimize upper bound error estimate. ", &
152 : default_l_val=.TRUE., &
153 83190 : lone_keyword_l_val=.TRUE.)
154 83190 : CALL section_add_keyword(section, keyword)
155 83190 : CALL keyword_release(keyword)
156 :
157 : CALL keyword_create(keyword, __LOCATION__, &
158 : name="DO_ERROR_ESTIMATE", &
159 : description="Whether the error due to minimax approx. and cutoff shall be estimated", &
160 : default_l_val=.TRUE., &
161 83190 : lone_keyword_l_val=.TRUE.)
162 83190 : CALL section_add_keyword(section, keyword)
163 83190 : CALL keyword_release(keyword)
164 :
165 : CALL cp_print_key_section_create(print_key, __LOCATION__, "ERI_MME_INFO", &
166 : description="Controls the printing info.", &
167 83190 : print_level=medium_print_level, filename="__STD_OUT__")
168 83190 : CALL section_add_subsection(section, print_key)
169 83190 : CALL section_release(print_key)
170 :
171 : CALL keyword_create(keyword, __LOCATION__, &
172 : name="PRINT_CALIB", &
173 : description="Print detailed info on calibration. ", &
174 : default_l_val=.FALSE., &
175 83190 : lone_keyword_l_val=.TRUE.)
176 83190 : CALL section_add_keyword(section, keyword)
177 83190 : CALL keyword_release(keyword)
178 :
179 : CALL keyword_create(keyword, __LOCATION__, &
180 : name="DEBUG", &
181 : description="debug mode (consistency of summation methods is checked).", &
182 : default_l_val=.FALSE., &
183 83190 : lone_keyword_l_val=.TRUE.)
184 83190 : CALL section_add_keyword(section, keyword)
185 83190 : CALL keyword_release(keyword)
186 :
187 : CALL keyword_create(keyword, __LOCATION__, &
188 : name="DEBUG_TOLERANCE", &
189 : description="tolerance for rel. numerical error in debug mode.", &
190 83190 : default_r_val=1.0E-06_dp)
191 83190 : CALL section_add_keyword(section, keyword)
192 83190 : CALL keyword_release(keyword)
193 :
194 : CALL keyword_create(keyword, __LOCATION__, &
195 : name="DEBUG_NSUM_MAX", &
196 : description="restrict debug mode for non-ortho cells to this number of summands. "// &
197 : "Sums with more terms are not checked.", &
198 83190 : default_i_val=1000000)
199 83190 : CALL section_add_keyword(section, keyword)
200 83190 : CALL keyword_release(keyword)
201 :
202 : CALL section_create(subsection, __LOCATION__, name="CUTOFF_CALIB", &
203 : description="Parameters for the calibration of the energy cutoff by "// &
204 : "minimizing the errors due to finite cutoff and minimax approximation. "// &
205 : "Implemented as bisection of error(minimax) - error(cutoff). Not "// &
206 : "implemented for non-orthorhombic cells. ", &
207 83190 : n_keywords=5, n_subsections=0)
208 :
209 : CALL keyword_create(keyword, __LOCATION__, &
210 : name="MIN", &
211 : description="Initial guess of lower bound for cutoff. ", &
212 83190 : default_r_val=10.0_dp)
213 83190 : CALL section_add_keyword(subsection, keyword)
214 83190 : CALL keyword_release(keyword)
215 :
216 : CALL keyword_create(keyword, __LOCATION__, &
217 : name="MAX", &
218 : description="Initial guess of upper bound for cutoff. ", &
219 83190 : default_r_val=10000.0_dp)
220 83190 : CALL section_add_keyword(subsection, keyword)
221 83190 : CALL keyword_release(keyword)
222 :
223 : CALL keyword_create(keyword, __LOCATION__, &
224 : name="DELTA", &
225 : description="Relative widening of cutoff interval in case starting "// &
226 : "values are not valid. ", &
227 83190 : default_r_val=0.9_dp)
228 83190 : CALL section_add_keyword(subsection, keyword)
229 83190 : CALL keyword_release(keyword)
230 :
231 : CALL keyword_create(keyword, __LOCATION__, &
232 : name="EPS", &
233 : description="Relative cutoff precision required to stop calibration. ", &
234 83190 : default_r_val=0.01_dp)
235 83190 : CALL section_add_keyword(subsection, keyword)
236 83190 : CALL keyword_release(keyword)
237 :
238 83190 : CALL section_add_subsection(section, subsection)
239 83190 : CALL section_release(subsection)
240 :
241 83190 : END SUBROUTINE create_eri_mme_section
242 :
243 : ! **************************************************************************************************
244 : !> \brief Read input and initialize parameter type
245 : !> \param mme_section ...
246 : !> \param param ...
247 : ! **************************************************************************************************
248 700 : SUBROUTINE cp_eri_mme_init_read_input(mme_section, param)
249 : TYPE(section_vals_type), POINTER :: mme_section
250 : TYPE(cp_eri_mme_param), INTENT(INOUT) :: param
251 :
252 : INTEGER :: debug_nsum, n_minimax, unit_nr
253 : LOGICAL :: debug, do_calib_cutoff, do_error_est, &
254 : print_calib
255 : REAL(KIND=dp) :: cutoff, cutoff_delta, cutoff_eps, &
256 : cutoff_max, cutoff_min, debug_delta, &
257 : sum_precision
258 : TYPE(cp_logger_type), POINTER :: logger
259 : TYPE(section_vals_type), POINTER :: subsection
260 :
261 50 : logger => cp_get_default_logger()
262 : unit_nr = cp_print_key_unit_nr(logger, mme_section, "ERI_MME_INFO", &
263 50 : extension=".eri_mme")
264 :
265 50 : NULLIFY (subsection)
266 50 : CALL section_vals_val_get(mme_section, "N_MINIMAX", i_val=n_minimax)
267 :
268 50 : CALL section_vals_val_get(mme_section, "CUTOFF", r_val=cutoff)
269 50 : CALL section_vals_val_get(mme_section, "SUM_PRECISION", r_val=sum_precision)
270 50 : CALL section_vals_val_get(mme_section, "DO_CALIBRATE_CUTOFF", l_val=do_calib_cutoff)
271 50 : CALL section_vals_val_get(mme_section, "DO_ERROR_ESTIMATE", l_val=do_error_est)
272 50 : CALL section_vals_val_get(mme_section, "PRINT_CALIB", l_val=print_calib)
273 50 : subsection => section_vals_get_subs_vals(mme_section, "CUTOFF_CALIB")
274 50 : CALL section_vals_val_get(subsection, "MIN", r_val=cutoff_min)
275 50 : CALL section_vals_val_get(subsection, "MAX", r_val=cutoff_max)
276 50 : CALL section_vals_val_get(subsection, "EPS", r_val=cutoff_eps)
277 50 : CALL section_vals_val_get(subsection, "DELTA", r_val=cutoff_delta)
278 50 : CALL section_vals_val_get(mme_section, "DEBUG", l_val=debug)
279 50 : CALL section_vals_val_get(mme_section, "DEBUG_TOLERANCE", r_val=debug_delta)
280 50 : CALL section_vals_val_get(mme_section, "DEBUG_NSUM_MAX", i_val=debug_nsum)
281 :
282 50 : param%mme_section => mme_section
283 :
284 : CALL eri_mme_init(param%par, n_minimax, &
285 : cutoff, do_calib_cutoff, do_error_est, cutoff_min, cutoff_max, cutoff_eps, cutoff_delta, &
286 50 : sum_precision, debug, debug_delta, debug_nsum, unit_nr, print_calib)
287 :
288 50 : param%do_calib = do_calib_cutoff
289 :
290 50 : param%G_count_2c = 0
291 50 : param%R_count_2c = 0
292 50 : param%GG_count_3c = 0
293 50 : param%GR_count_3c = 0
294 50 : param%RR_count_3c = 0
295 :
296 50 : param%logger => logger
297 50 : END SUBROUTINE cp_eri_mme_init_read_input
298 :
299 : ! **************************************************************************************************
300 : !> \brief Release eri mme data. Prints some statistics on summation methods chosen.
301 : !> \param param ...
302 : ! **************************************************************************************************
303 50 : SUBROUTINE cp_eri_mme_finalize(param)
304 : TYPE(cp_eri_mme_param), INTENT(INOUT) :: param
305 :
306 : INTEGER :: count_2c, count_3c, unit_nr
307 :
308 50 : count_2c = param%G_count_2c + param%R_count_2c
309 50 : count_3c = param%GG_count_3c + param%GR_count_3c + param%RR_count_3c
310 :
311 50 : unit_nr = param%par%unit_nr
312 :
313 50 : IF (unit_nr > 0) THEN
314 18 : IF (count_2c .GT. 0) THEN
315 18 : WRITE (unit_nr, '(/T2, A)') "ERI_MME| Percentage of 2-center integrals evaluated in"
316 18 : WRITE (unit_nr, '(T2, A, T76, F5.1)') "ERI_MME| G space:", &
317 36 : 100.0_dp*param%G_count_2c/count_2c
318 18 : WRITE (unit_nr, '(T2, A, T76, F5.1/)') "ERI_MME| R space:", &
319 36 : 100.0_dp*param%R_count_2c/count_2c
320 : END IF
321 18 : IF (count_3c .GT. 0) THEN
322 6 : WRITE (unit_nr, '(/T2, A)') "ERI_MME| Percentage of 3-center integrals evaluated in"
323 6 : WRITE (unit_nr, '(T2, A, T76, F5.1)') "ERI_MME| G/G space:", &
324 12 : 100.0_dp*param%GG_count_3c/count_3c
325 6 : WRITE (unit_nr, '(T2, A, T76, F5.1)') "ERI_MME| G/R space:", &
326 12 : 100.0_dp*param%GR_count_3c/count_3c
327 6 : WRITE (unit_nr, '(T2, A, T76, F5.1/)') "ERI_MME| R/R space:", &
328 12 : 100.0_dp*param%RR_count_3c/count_3c
329 : END IF
330 : END IF
331 :
332 50 : CALL eri_mme_release(param%par)
333 50 : CALL cp_print_key_finished_output(unit_nr, param%logger, param%mme_section, "ERI_MME_INFO")
334 50 : END SUBROUTINE cp_eri_mme_finalize
335 :
336 : ! **************************************************************************************************
337 : !> \brief Set parameters for MME method by deriving basis info from basis set.
338 : !> Cutoff can be auto-calibrated to minimize total error.
339 : !> \param param ...
340 : !> \param cell ...
341 : !> \param qs_kind_set ...
342 : !> \param basis_type_1 ...
343 : !> \param basis_type_2 ...
344 : !> \param para_env ...
345 : !> \param potential ...
346 : !> \param pot_par ...
347 : ! **************************************************************************************************
348 142 : SUBROUTINE eri_mme_set_params_from_basis(param, cell, qs_kind_set, basis_type_1, basis_type_2, para_env, &
349 : potential, pot_par)
350 : TYPE(cp_eri_mme_param), INTENT(INOUT) :: param
351 : TYPE(cell_type), POINTER :: cell
352 : TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
353 : POINTER :: qs_kind_set
354 : CHARACTER(len=*), INTENT(IN) :: basis_type_1
355 : CHARACTER(len=*), INTENT(IN), OPTIONAL :: basis_type_2
356 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
357 : INTEGER, INTENT(IN), OPTIONAL :: potential
358 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par
359 :
360 : CHARACTER(LEN=*), PARAMETER :: routineN = 'eri_mme_set_params_from_basis'
361 :
362 : INTEGER :: handle, l_max, l_max_zet
363 : REAL(KIND=dp) :: zet_max, zet_min
364 :
365 142 : CALL timeset(routineN, handle)
366 :
367 : CALL error_est_pgf_params_from_basis(qs_kind_set, basis_type_1, basis_type_2, &
368 158 : zet_min, zet_max, l_max_zet, l_max)
369 :
370 : CALL eri_mme_set_params_custom(param, cell%hmat, cell%orthorhombic, &
371 : zet_min, zet_max, l_max_zet, &
372 : l_max, para_env, &
373 142 : potential, pot_par)
374 :
375 142 : CALL timestop(handle)
376 142 : END SUBROUTINE eri_mme_set_params_from_basis
377 :
378 : ! **************************************************************************************************
379 : !> \brief Wrapper for eri_mme_set_params
380 : !> \param param ...
381 : !> \param hmat ...
382 : !> \param is_ortho ...
383 : !> \param zet_min ...
384 : !> \param zet_max ...
385 : !> \param l_max_zet ...
386 : !> \param l_max ...
387 : !> \param para_env ...
388 : !> \param potential ...
389 : !> \param pot_par ...
390 : ! **************************************************************************************************
391 158 : SUBROUTINE eri_mme_set_params_custom(param, hmat, is_ortho, zet_min, zet_max, l_max_zet, l_max, para_env, &
392 : potential, pot_par)
393 : TYPE(cp_eri_mme_param), INTENT(INOUT) :: param
394 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat
395 : LOGICAL, INTENT(IN) :: is_ortho
396 : REAL(KIND=dp), INTENT(IN) :: zet_min, zet_max
397 : INTEGER, INTENT(IN) :: l_max_zet, l_max
398 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
399 : INTEGER, INTENT(IN), OPTIONAL :: potential
400 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par
401 :
402 : REAL(KIND=dp), PARAMETER :: eps_changed = 1.0E-14_dp
403 :
404 158 : IF (param%do_calib) THEN
405 156 : IF (.NOT. param%par%is_valid) THEN
406 64 : param%par%do_calib_cutoff = .TRUE.
407 : ELSE
408 : ! only calibrate cutoff if parameters (cell, basis coefficients) have changed
409 : IF (ALL(ABS(param%par%hmat - hmat) < eps_changed) .AND. &
410 : ABS(param%par%zet_min - zet_min) < eps_changed .AND. &
411 784 : ABS(param%par%zet_max - zet_max) < eps_changed .AND. &
412 : param%par%l_max_zet == l_max_zet) THEN
413 40 : param%par%do_calib_cutoff = .FALSE.
414 : ELSE
415 52 : param%par%do_calib_cutoff = .TRUE.
416 : END IF
417 : END IF
418 : ELSE
419 2 : param%par%do_calib_cutoff = .FALSE.
420 : END IF
421 :
422 : CALL eri_mme_set_params(param%par, hmat, is_ortho, zet_min, zet_max, l_max_zet, l_max, para_env, &
423 158 : potential, pot_par)
424 :
425 158 : CALL eri_mme_print_info(param)
426 158 : END SUBROUTINE eri_mme_set_params_custom
427 :
428 : ! **************************************************************************************************
429 : !> \brief Get basis parameters for estimating cutoff and minimax error from cp2k basis
430 : !> \param qs_kind_set ...
431 : !> \param basis_type_1 ...
432 : !> \param basis_type_2 ...
433 : !> \param zet_min Smallest exponent, used to estimate error due to minimax approx.
434 : !> \param zet_max contains max. exponent,
435 : !> used to estimate cutoff error
436 : !> \param l_max_zet contains the largest l for max. exponent,
437 : !> used to estimate cutoff error
438 : !> \param l_max ...
439 : ! **************************************************************************************************
440 142 : SUBROUTINE error_est_pgf_params_from_basis(qs_kind_set, basis_type_1, basis_type_2, zet_min, zet_max, l_max_zet, l_max)
441 : TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
442 : POINTER :: qs_kind_set
443 : CHARACTER(len=*), INTENT(IN) :: basis_type_1
444 : CHARACTER(len=*), INTENT(IN), OPTIONAL :: basis_type_2
445 : REAL(KIND=dp), INTENT(OUT) :: zet_min, zet_max
446 : INTEGER, INTENT(OUT) :: l_max_zet, l_max
447 :
448 : CHARACTER(LEN=*), PARAMETER :: routineN = 'error_est_pgf_params_from_basis'
449 :
450 : CHARACTER(len=default_string_length) :: basis_type
451 : INTEGER :: handle, ibasis, ikind, ipgf, iset, l_m, &
452 : l_zet, nbasis, nkind, nset
453 142 : INTEGER, DIMENSION(:), POINTER :: npgf
454 : REAL(KIND=dp) :: zet_m
455 : TYPE(gto_basis_set_type), POINTER :: basis_set
456 :
457 142 : CALL timeset(routineN, handle)
458 :
459 142 : l_m = 0
460 142 : zet_m = 0.0_dp
461 142 : l_zet = -1
462 142 : zet_min = -1.0_dp
463 :
464 142 : nkind = SIZE(qs_kind_set)
465 142 : nbasis = MERGE(2, 1, PRESENT(basis_type_2))
466 :
467 : ! 1) get global max l and max zet
468 : ! (and min zet for minimax error)
469 342 : DO ikind = 1, nkind
470 724 : DO ibasis = 1, nbasis
471 382 : IF (ibasis .EQ. 1) THEN
472 200 : basis_type = basis_type_1
473 : ELSE
474 182 : basis_type = basis_type_2
475 : END IF
476 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set, &
477 382 : basis_type=basis_type)
478 382 : CPASSERT(ASSOCIATED(basis_set))
479 382 : npgf => basis_set%npgf
480 382 : nset = basis_set%nset
481 2572 : l_m = MAX(l_m, MAXVAL(basis_set%lmax(:)))
482 2772 : DO iset = 1, nset
483 7258 : zet_m = MAX(zet_m, MAXVAL(basis_set%zet(1:npgf(iset), iset)))
484 2572 : IF (zet_min .LT. 0.0_dp) THEN
485 872 : zet_min = MINVAL(basis_set%zet(1:npgf(iset), iset))
486 : ELSE
487 6386 : zet_min = MIN(zet_min, MINVAL(basis_set%zet(1:npgf(iset), iset)))
488 : END IF
489 : END DO
490 : END DO
491 : END DO
492 :
493 : ! 2) get largest zet for max l and largest l for max zet
494 342 : DO ikind = 1, nkind
495 724 : DO ibasis = 1, nbasis
496 382 : IF (ibasis .EQ. 1) THEN
497 200 : basis_type = basis_type_1
498 : ELSE
499 182 : basis_type = basis_type_2
500 : END IF
501 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set, &
502 382 : basis_type=basis_type)
503 2772 : DO iset = 1, basis_set%nset
504 5450 : DO ipgf = 1, basis_set%npgf(iset)
505 : IF (ABS(zet_m - basis_set%zet(ipgf, iset)) .LE. (zet_m*1.0E-12_dp) &
506 5068 : .AND. (basis_set%lmax(iset) .GT. l_zet)) THEN
507 : l_zet = basis_set%lmax(iset)
508 : END IF
509 : END DO
510 : END DO
511 : END DO
512 : END DO
513 :
514 142 : CPASSERT(l_zet .GE. 0)
515 :
516 142 : zet_max = zet_m
517 :
518 : ! l + 1 because we may calculate forces
519 : ! this is probably a safe choice also for the case that forces are not needed
520 142 : l_max_zet = l_zet + 1
521 142 : l_max = l_m + 1
522 :
523 142 : CALL timestop(handle)
524 142 : END SUBROUTINE error_est_pgf_params_from_basis
525 :
526 : ! **************************************************************************************************
527 : !> \brief ...
528 : !> \param param ...
529 : ! **************************************************************************************************
530 158 : SUBROUTINE eri_mme_print_info(param)
531 : TYPE(cp_eri_mme_param) :: param
532 :
533 : INTEGER :: igrid, unit_nr
534 : LOGICAL :: print_multigrids
535 : TYPE(cp_logger_type), POINTER :: logger
536 :
537 158 : print_multigrids = .FALSE.
538 :
539 158 : logger => param%logger
540 :
541 158 : unit_nr = param%par%unit_nr
542 :
543 158 : IF (unit_nr > 0) THEN
544 44 : SELECT CASE (param%par%potential)
545 : CASE (eri_mme_coulomb)
546 21 : WRITE (unit_nr, '(/T2, A)') "ERI_MME| Potential: Coulomb"
547 : CASE (eri_mme_yukawa)
548 1 : WRITE (unit_nr, '(/T2, A, ES9.2)') "ERI_MME| Potential: Yukawa with a=", param%par%pot_par
549 : CASE (eri_mme_longrange)
550 23 : WRITE (unit_nr, '(/T2, A, ES9.2)') "ERI_MME| Potential: long-range Coulomb with a=", param%par%pot_par
551 : END SELECT
552 : END IF
553 :
554 : IF (unit_nr > 0) THEN
555 23 : WRITE (unit_nr, '(/T2, A, T71, F10.1)') "ERI_MME| Cutoff for ERIs [a.u.]:", param%par%cutoff
556 23 : WRITE (unit_nr, '(/T2, A, T78, I3/)') "ERI_MME| Number of terms in minimax approximation:", param%par%n_minimax
557 : END IF
558 158 : IF (param%par%is_ortho) THEN
559 122 : IF (unit_nr > 0) THEN
560 18 : IF (param%par%do_error_est) THEN
561 18 : WRITE (unit_nr, '(T2, A)') "ERI_MME| Estimated absolute error for normalized Hermite-Gaussian basis"
562 18 : WRITE (unit_nr, '(T2, A, T72, ES9.2)') "ERI_MME| Minimax error:", param%par%err_mm
563 18 : WRITE (unit_nr, '(T2, A, T72, ES9.2)') "ERI_MME| Cutoff error:", param%par%err_c
564 18 : WRITE (unit_nr, '(T2, A, T72, ES9.2)') "ERI_MME| Total error (minimax + cutoff):", param%par%err_mm + param%par%err_c
565 : END IF
566 18 : IF (param%par%print_calib) &
567 0 : WRITE (unit_nr, '(T2, A, T68, F13.10)') "ERI_MME| Minimax scaling constant in AM-GM estimate:", param%par%C_mm
568 : END IF
569 : END IF
570 :
571 : IF (print_multigrids) THEN
572 : DO igrid = 1, param%par%n_grids
573 : CALL eri_mme_print_grid_info(param%par%minimax_grid(igrid), igrid, unit_nr)
574 : END DO
575 : END IF
576 :
577 54 : IF (unit_nr > 0) WRITE (unit_nr, *)
578 :
579 158 : END SUBROUTINE eri_mme_print_info
580 :
581 : ! **************************************************************************************************
582 : !> \brief Create input section for unit testing
583 : !> \param section ...
584 : ! **************************************************************************************************
585 9238 : SUBROUTINE create_eri_mme_test_section(section)
586 : TYPE(section_type), INTENT(INOUT), POINTER :: section
587 :
588 : TYPE(keyword_type), POINTER :: keyword
589 : TYPE(section_type), POINTER :: subsection
590 :
591 9238 : NULLIFY (keyword, subsection)
592 :
593 9238 : CPASSERT(.NOT. ASSOCIATED(section))
594 : CALL section_create(section, __LOCATION__, name="ERI_MME_TEST", &
595 : description="Parameters for testing the ERI_MME method for electron repulsion integrals. "// &
596 : "Testing w.r.t. performance and accuracy. ", &
597 9238 : n_keywords=5, n_subsections=1)
598 :
599 9238 : CALL create_eri_mme_section(subsection)
600 9238 : CALL section_add_subsection(section, subsection)
601 9238 : CALL section_release(subsection)
602 :
603 : CALL keyword_create(keyword, __LOCATION__, &
604 : name="_SECTION_PARAMETERS_", &
605 : description="Controls the activation the ERI_MME test. ", &
606 : default_l_val=.FALSE., &
607 9238 : lone_keyword_l_val=.TRUE.)
608 9238 : CALL section_add_keyword(section, keyword)
609 9238 : CALL keyword_release(keyword)
610 :
611 : CALL keyword_create(keyword, __LOCATION__, name="TEST_3C", &
612 : description="Whether to test 3-center integrals.", &
613 : default_l_val=.TRUE., &
614 9238 : lone_keyword_l_val=.TRUE.)
615 9238 : CALL section_add_keyword(section, keyword)
616 9238 : CALL keyword_release(keyword)
617 :
618 : CALL keyword_create(keyword, __LOCATION__, name="TEST_2C", &
619 : description="Whether to test 2-center integrals.", &
620 : default_l_val=.TRUE., &
621 9238 : lone_keyword_l_val=.TRUE.)
622 9238 : CALL section_add_keyword(section, keyword)
623 9238 : CALL keyword_release(keyword)
624 :
625 : CALL keyword_create(keyword, __LOCATION__, name="ABC", &
626 : description="Specify the lengths of the cell vectors A, B, and C. ", &
627 : usage="ABC 10.000 10.000 10.000", unit_str="angstrom", &
628 9238 : n_var=3, type_of_var=real_t)
629 9238 : CALL section_add_keyword(section, keyword)
630 9238 : CALL keyword_release(keyword)
631 :
632 : CALL keyword_create(keyword, __LOCATION__, name="MIN_NPOS", &
633 : description="Minimum number of atomic distances to consider. ", &
634 9238 : default_i_val=8)
635 9238 : CALL section_add_keyword(section, keyword)
636 9238 : CALL keyword_release(keyword)
637 :
638 : CALL keyword_create(keyword, __LOCATION__, name="NREP", &
639 : description="Number of repeated calculation of each integral. ", &
640 9238 : default_i_val=1)
641 9238 : CALL section_add_keyword(section, keyword)
642 9238 : CALL keyword_release(keyword)
643 :
644 : CALL keyword_create(keyword, __LOCATION__, name="CHECK_2C_ACCURACY", &
645 : description="Whether integrals should be compared to reference values, "// &
646 : "created on the fly by exact method (G-space sum on grid without "// &
647 : "minimax approximation). Note: only feasible for not too many "// &
648 : "integrals and maximum exponent around 10.0. ", &
649 : default_l_val=.FALSE., &
650 9238 : lone_keyword_l_val=.TRUE.)
651 :
652 9238 : CALL section_add_keyword(section, keyword)
653 9238 : CALL keyword_release(keyword)
654 :
655 : CALL keyword_create(keyword, __LOCATION__, name="LMAX", &
656 : description="Maximum total angular momentum quantum number. ", &
657 9238 : default_i_val=6)
658 9238 : CALL section_add_keyword(section, keyword)
659 9238 : CALL keyword_release(keyword)
660 :
661 : CALL keyword_create(keyword, __LOCATION__, name="ZET_MIN", &
662 : description="Minimum exponent. ", &
663 9238 : default_r_val=0.001_dp)
664 9238 : CALL section_add_keyword(section, keyword)
665 9238 : CALL keyword_release(keyword)
666 :
667 : CALL keyword_create(keyword, __LOCATION__, name="ZET_MAX", &
668 : description="Maximum exponent. ", &
669 9238 : default_r_val=1.0_dp)
670 9238 : CALL section_add_keyword(section, keyword)
671 9238 : CALL keyword_release(keyword)
672 :
673 : CALL keyword_create(keyword, __LOCATION__, name="NZET", &
674 : description="Number of exponents (logarithmic partition between ZET_MIN and ZET_MAX). ", &
675 9238 : default_i_val=4)
676 9238 : CALL section_add_keyword(section, keyword)
677 9238 : CALL keyword_release(keyword)
678 :
679 : CALL keyword_create(keyword, __LOCATION__, name="NSAMPLE_3C", &
680 : description="If NSAMPLE_3C = N, only calculate every Nth 3-center integral.", &
681 9238 : default_i_val=1)
682 9238 : CALL section_add_keyword(section, keyword)
683 9238 : CALL keyword_release(keyword)
684 :
685 : CALL keyword_create(keyword, __LOCATION__, name="POTENTIAL", &
686 : description="Operator to test", &
687 : default_i_val=eri_mme_coulomb, &
688 : enum_i_vals=(/eri_mme_coulomb, eri_mme_yukawa, eri_mme_longrange/), &
689 : enum_c_vals=s2a("COULOMB", "YUKAWA", "LONGRANGE"), &
690 9238 : enum_desc=s2a("1/r", "exp(-a*r)/r", "erf(a*r)/r"))
691 9238 : CALL section_add_keyword(section, keyword)
692 9238 : CALL keyword_release(keyword)
693 :
694 : CALL keyword_create(keyword, __LOCATION__, name="POTENTIAL_PARAM", &
695 : description="Parameter 'a' for chosen potential, ignored for Coulomb", &
696 9238 : default_r_val=0.0_dp)
697 9238 : CALL section_add_keyword(section, keyword)
698 9238 : CALL keyword_release(keyword)
699 :
700 9238 : END SUBROUTINE create_eri_mme_test_section
701 :
702 : ! **************************************************************************************************
703 : !> \brief Update local counters to gather statistics on different paths taken in MME algorithm
704 : !> (each Ewald sum can be performed over direct or reciprocal lattice vectors)
705 : !> \param param ...
706 : !> \param para_env ...
707 : !> \param G_count_2c ...
708 : !> \param R_count_2c ...
709 : !> \param GG_count_3c ...
710 : !> \param GR_count_3c ...
711 : !> \param RR_count_3c ...
712 : ! **************************************************************************************************
713 388 : SUBROUTINE cp_eri_mme_update_local_counts(param, para_env, G_count_2c, R_count_2c, GG_count_3c, GR_count_3c, RR_count_3c)
714 : TYPE(cp_eri_mme_param), INTENT(INOUT) :: param
715 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
716 : INTEGER, INTENT(INOUT), OPTIONAL :: G_count_2c, R_count_2c, GG_count_3c, &
717 : GR_count_3c, RR_count_3c
718 :
719 388 : IF (PRESENT(G_count_2c)) THEN
720 192 : CALL para_env%sum(G_count_2c)
721 192 : param%G_count_2c = param%G_count_2c + G_count_2c
722 : END IF
723 :
724 388 : IF (PRESENT(R_count_2c)) THEN
725 192 : CALL para_env%sum(R_count_2c)
726 192 : param%R_count_2c = param%R_count_2c + R_count_2c
727 : END IF
728 :
729 388 : IF (PRESENT(GG_count_3c)) THEN
730 204 : CALL para_env%sum(GG_count_3c)
731 204 : param%GG_count_3c = param%GG_count_3c + GG_count_3c
732 : END IF
733 :
734 388 : IF (PRESENT(GR_count_3c)) THEN
735 204 : CALL para_env%sum(GR_count_3c)
736 204 : param%GR_count_3c = param%GR_count_3c + GR_count_3c
737 : END IF
738 :
739 388 : IF (PRESENT(RR_count_3c)) THEN
740 204 : CALL para_env%sum(RR_count_3c)
741 204 : param%RR_count_3c = param%RR_count_3c + RR_count_3c
742 : END IF
743 :
744 388 : END SUBROUTINE cp_eri_mme_update_local_counts
745 :
746 : ! **************************************************************************************************
747 : !> \brief ...
748 : !> \param para_env ...
749 : !> \param iw ...
750 : !> \param eri_mme_test_section ...
751 : ! **************************************************************************************************
752 8 : SUBROUTINE cp_eri_mme_perf_acc_test(para_env, iw, eri_mme_test_section)
753 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
754 : INTEGER, INTENT(IN) :: iw
755 : TYPE(section_vals_type), POINTER :: eri_mme_test_section
756 :
757 : INTEGER :: count_r, G_count, GG_count, GR_count, i, &
758 : ix, iy, iz, l_max, min_nR, nR, nR_xyz, &
759 : nrep, nsample, nzet, potential, &
760 : R_count, RR_count
761 : LOGICAL :: test_2c, test_3c, test_accuracy
762 : REAL(KIND=dp) :: pot_par, zet_fac, zetmax, zetmin
763 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: zet
764 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rabc
765 8 : REAL(KIND=dp), DIMENSION(:), POINTER :: cell_par
766 : TYPE(cell_type), POINTER :: box
767 200 : TYPE(cp_eri_mme_param) :: param
768 : TYPE(section_vals_type), POINTER :: eri_mme_section
769 :
770 8 : NULLIFY (box, eri_mme_section, cell_par)
771 :
772 8 : eri_mme_section => section_vals_get_subs_vals(eri_mme_test_section, "ERI_MME")
773 8 : CALL cp_eri_mme_init_read_input(eri_mme_section, param)
774 8 : CALL section_vals_val_get(eri_mme_test_section, "TEST_3C", l_val=test_3c)
775 8 : CALL section_vals_val_get(eri_mme_test_section, "TEST_2C", l_val=test_2c)
776 :
777 8 : CALL section_vals_val_get(eri_mme_test_section, "ABC", r_vals=cell_par)
778 8 : CALL section_vals_val_get(eri_mme_test_section, "MIN_NPOS", i_val=min_nR)
779 8 : CALL section_vals_val_get(eri_mme_test_section, "NREP", i_val=nrep)
780 8 : CALL section_vals_val_get(eri_mme_test_section, "CHECK_2C_ACCURACY", l_val=test_accuracy)
781 8 : CALL section_vals_val_get(eri_mme_test_section, "LMAX", i_val=l_max)
782 8 : CALL section_vals_val_get(eri_mme_test_section, "ZET_MIN", r_val=zetmin)
783 8 : CALL section_vals_val_get(eri_mme_test_section, "ZET_MAX", r_val=zetmax)
784 8 : CALL section_vals_val_get(eri_mme_test_section, "NZET", i_val=nzet)
785 8 : CALL section_vals_val_get(eri_mme_test_section, "NSAMPLE_3C", i_val=nsample)
786 8 : CALL section_vals_val_get(eri_mme_test_section, "POTENTIAL", i_val=potential)
787 8 : CALL section_vals_val_get(eri_mme_test_section, "POTENTIAL_PARAM", r_val=pot_par)
788 :
789 8 : IF (nzet .LE. 0) &
790 0 : CPABORT("Number of exponents NZET must be greater than 0.")
791 :
792 8 : CALL init_orbital_pointers(l_max)
793 :
794 : ! Create ranges of zet to be tested
795 24 : ALLOCATE (zet(nzet))
796 :
797 8 : zet(1) = zetmin
798 8 : IF (nzet .GT. 1) THEN
799 8 : zet_fac = (zetmax/zetmin)**(1.0_dp/(nzet - 1))
800 44 : DO i = 1, nzet - 1
801 44 : zet(i + 1) = zet(i)*zet_fac
802 : END DO
803 : END IF
804 :
805 : ! initialize cell
806 8 : CALL cell_create(box)
807 104 : box%hmat = 0.0_dp
808 8 : box%hmat(1, 1) = cell_par(1)
809 8 : box%hmat(2, 2) = cell_par(2)
810 8 : box%hmat(3, 3) = cell_par(3)
811 8 : CALL init_cell(box)
812 :
813 : ! Create range of rab (atomic distances) to be tested
814 8 : nR_xyz = CEILING(REAL(min_nR, KIND=dp)**(1.0_dp/3.0_dp) - 1.0E-06)
815 8 : nR = nR_xyz**3
816 :
817 24 : ALLOCATE (rabc(3, nR))
818 8 : count_r = 0
819 24 : DO ix = 1, nR_xyz
820 56 : DO iy = 1, nR_xyz
821 112 : DO iz = 1, nR_xyz
822 64 : count_r = count_r + 1
823 : ! adding 10% of cell size to positions to avoid atoms exactly at boundary or center of a cell
824 : rabc(:, count_r) = pbc([ix*ABS(cell_par(1)), &
825 : iy*ABS(cell_par(2)), &
826 : iz*ABS(cell_par(3))]/nR_xyz + &
827 288 : 0.1_dp*ABS(cell_par(:)), box)
828 : END DO
829 : END DO
830 : END DO
831 :
832 : ! initialize MME method
833 : CALL cp_eri_mme_set_params(param, box%hmat, box%orthorhombic, MINVAL(zet), MAXVAL(zet), l_max, l_max, para_env, &
834 112 : potential, pot_par)
835 :
836 8 : IF (iw > 0) WRITE (iw, '(T2, A, T61, I20)') "ERI_MME| Number of atomic distances:", nR
837 :
838 8 : G_count = 0; R_count = 0
839 8 : GG_count = 0; GR_count = 0; RR_count = 0
840 :
841 8 : IF (test_2c) CALL eri_mme_2c_perf_acc_test(param%par, l_max, zet, rabc, nrep, test_accuracy, para_env, iw, &
842 8 : potential=potential, pot_par=pot_par, G_count=G_count, R_count=R_count)
843 8 : IF (test_3c) CALL eri_mme_3c_perf_acc_test(param%par, l_max, zet, rabc, nrep, nsample, &
844 : para_env, iw, potential=potential, pot_par=pot_par, &
845 2 : GG_count=GG_count, GR_count=GR_count, RR_count=RR_count)
846 8 : CALL cp_eri_mme_update_local_counts(param, para_env, G_count, R_count, GG_count, GR_count, RR_count)
847 8 : CALL cp_eri_mme_finalize(param)
848 8 : CALL cell_release(box)
849 32 : END SUBROUTINE cp_eri_mme_perf_acc_test
850 :
851 0 : END MODULE cp_eri_mme_interface
|
