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 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 106316 : 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 106316 : NULLIFY (keyword, print_key, subsection)
109 106316 : CPASSERT(.NOT. ASSOCIATED(section))
110 :
111 106316 : IF (PRESENT(default_n_minimax)) THEN
112 10624 : my_default_n_minimax = default_n_minimax
113 : ELSE
114 95692 : 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 106316 : 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 106316 : default_i_val=my_default_n_minimax)
130 106316 : CALL section_add_keyword(section, keyword)
131 106316 : 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 106316 : default_r_val=300.0_dp)
138 106316 : CALL section_add_keyword(section, keyword)
139 106316 : 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 106316 : default_r_val=1.0E-16_dp)
145 106316 : CALL section_add_keyword(section, keyword)
146 106316 : 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 106316 : lone_keyword_l_val=.TRUE.)
154 106316 : CALL section_add_keyword(section, keyword)
155 106316 : 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 106316 : lone_keyword_l_val=.TRUE.)
162 106316 : CALL section_add_keyword(section, keyword)
163 106316 : 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 106316 : print_level=medium_print_level, filename="__STD_OUT__")
168 106316 : CALL section_add_subsection(section, print_key)
169 106316 : 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 106316 : lone_keyword_l_val=.TRUE.)
176 106316 : CALL section_add_keyword(section, keyword)
177 106316 : 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 106316 : lone_keyword_l_val=.TRUE.)
184 106316 : CALL section_add_keyword(section, keyword)
185 106316 : 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 106316 : default_r_val=1.0E-06_dp)
191 106316 : CALL section_add_keyword(section, keyword)
192 106316 : 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 106316 : default_i_val=1000000)
199 106316 : CALL section_add_keyword(section, keyword)
200 106316 : 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 106316 : 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 106316 : default_r_val=10.0_dp)
213 106316 : CALL section_add_keyword(subsection, keyword)
214 106316 : CALL keyword_release(keyword)
215 :
216 : CALL keyword_create(keyword, __LOCATION__, &
217 : name="MAX", &
218 : description="Initial guess of upper bound for cutoff. ", &
219 106316 : default_r_val=10000.0_dp)
220 106316 : CALL section_add_keyword(subsection, keyword)
221 106316 : 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 106316 : default_r_val=0.9_dp)
228 106316 : CALL section_add_keyword(subsection, keyword)
229 106316 : CALL keyword_release(keyword)
230 :
231 : CALL keyword_create(keyword, __LOCATION__, &
232 : name="EPS", &
233 : description="Relative cutoff precision required to stop calibration. ", &
234 106316 : default_r_val=0.01_dp)
235 106316 : CALL section_add_keyword(subsection, keyword)
236 106316 : CALL keyword_release(keyword)
237 :
238 106316 : CALL section_add_subsection(section, subsection)
239 106316 : CALL section_release(subsection)
240 :
241 106316 : 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 > 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 > 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 == 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 5068 : zet_m = MAX(zet_m, MAXVAL(basis_set%zet(1:npgf(iset), iset)))
484 2572 : IF (zet_min < 0.0_dp) THEN
485 730 : zet_min = MINVAL(basis_set%zet(1:npgf(iset), iset))
486 : ELSE
487 4338 : 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 == 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)) <= (zet_m*1.0E-12_dp) &
506 5068 : .AND. (basis_set%lmax(iset) > 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 >= 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) THEN
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 : END IF
571 :
572 : IF (print_multigrids) THEN
573 : DO igrid = 1, param%par%n_grids
574 : CALL eri_mme_print_grid_info(param%par%minimax_grid(igrid), igrid, unit_nr)
575 : END DO
576 : END IF
577 :
578 54 : IF (unit_nr > 0) WRITE (unit_nr, *)
579 :
580 158 : END SUBROUTINE eri_mme_print_info
581 :
582 : ! **************************************************************************************************
583 : !> \brief Create input section for unit testing
584 : !> \param section ...
585 : ! **************************************************************************************************
586 10624 : SUBROUTINE create_eri_mme_test_section(section)
587 : TYPE(section_type), INTENT(INOUT), POINTER :: section
588 :
589 : TYPE(keyword_type), POINTER :: keyword
590 : TYPE(section_type), POINTER :: subsection
591 :
592 10624 : NULLIFY (keyword, subsection)
593 :
594 10624 : CPASSERT(.NOT. ASSOCIATED(section))
595 : CALL section_create(section, __LOCATION__, name="ERI_MME_TEST", &
596 : description="Parameters for testing the ERI_MME method for electron repulsion integrals. "// &
597 : "Testing w.r.t. performance and accuracy. ", &
598 10624 : n_keywords=5, n_subsections=1)
599 :
600 10624 : CALL create_eri_mme_section(subsection)
601 10624 : CALL section_add_subsection(section, subsection)
602 10624 : CALL section_release(subsection)
603 :
604 : CALL keyword_create(keyword, __LOCATION__, &
605 : name="_SECTION_PARAMETERS_", &
606 : description="Controls the activation the ERI_MME test. ", &
607 : default_l_val=.FALSE., &
608 10624 : lone_keyword_l_val=.TRUE.)
609 10624 : CALL section_add_keyword(section, keyword)
610 10624 : CALL keyword_release(keyword)
611 :
612 : CALL keyword_create(keyword, __LOCATION__, name="TEST_3C", &
613 : description="Whether to test 3-center integrals.", &
614 : default_l_val=.TRUE., &
615 10624 : lone_keyword_l_val=.TRUE.)
616 10624 : CALL section_add_keyword(section, keyword)
617 10624 : CALL keyword_release(keyword)
618 :
619 : CALL keyword_create(keyword, __LOCATION__, name="TEST_2C", &
620 : description="Whether to test 2-center integrals.", &
621 : default_l_val=.TRUE., &
622 10624 : lone_keyword_l_val=.TRUE.)
623 10624 : CALL section_add_keyword(section, keyword)
624 10624 : CALL keyword_release(keyword)
625 :
626 : CALL keyword_create(keyword, __LOCATION__, name="ABC", &
627 : description="Specify the lengths of the cell vectors A, B, and C. ", &
628 : usage="ABC 10.000 10.000 10.000", unit_str="angstrom", &
629 10624 : n_var=3, type_of_var=real_t)
630 10624 : CALL section_add_keyword(section, keyword)
631 10624 : CALL keyword_release(keyword)
632 :
633 : CALL keyword_create(keyword, __LOCATION__, name="MIN_NPOS", &
634 : description="Minimum number of atomic distances to consider. ", &
635 10624 : default_i_val=8)
636 10624 : CALL section_add_keyword(section, keyword)
637 10624 : CALL keyword_release(keyword)
638 :
639 : CALL keyword_create(keyword, __LOCATION__, name="NREP", &
640 : description="Number of repeated calculation of each integral. ", &
641 10624 : default_i_val=1)
642 10624 : CALL section_add_keyword(section, keyword)
643 10624 : CALL keyword_release(keyword)
644 :
645 : CALL keyword_create(keyword, __LOCATION__, name="CHECK_2C_ACCURACY", &
646 : description="Whether integrals should be compared to reference values, "// &
647 : "created on the fly by exact method (G-space sum on grid without "// &
648 : "minimax approximation). Note: only feasible for not too many "// &
649 : "integrals and maximum exponent around 10.0. ", &
650 : default_l_val=.FALSE., &
651 10624 : lone_keyword_l_val=.TRUE.)
652 :
653 10624 : CALL section_add_keyword(section, keyword)
654 10624 : CALL keyword_release(keyword)
655 :
656 : CALL keyword_create(keyword, __LOCATION__, name="LMAX", &
657 : description="Maximum total angular momentum quantum number. ", &
658 10624 : default_i_val=6)
659 10624 : CALL section_add_keyword(section, keyword)
660 10624 : CALL keyword_release(keyword)
661 :
662 : CALL keyword_create(keyword, __LOCATION__, name="ZET_MIN", &
663 : description="Minimum exponent. ", &
664 10624 : default_r_val=0.001_dp)
665 10624 : CALL section_add_keyword(section, keyword)
666 10624 : CALL keyword_release(keyword)
667 :
668 : CALL keyword_create(keyword, __LOCATION__, name="ZET_MAX", &
669 : description="Maximum exponent. ", &
670 10624 : default_r_val=1.0_dp)
671 10624 : CALL section_add_keyword(section, keyword)
672 10624 : CALL keyword_release(keyword)
673 :
674 : CALL keyword_create(keyword, __LOCATION__, name="NZET", &
675 : description="Number of exponents (logarithmic partition between ZET_MIN and ZET_MAX). ", &
676 10624 : default_i_val=4)
677 10624 : CALL section_add_keyword(section, keyword)
678 10624 : CALL keyword_release(keyword)
679 :
680 : CALL keyword_create(keyword, __LOCATION__, name="NSAMPLE_3C", &
681 : description="If NSAMPLE_3C = N, only calculate every Nth 3-center integral.", &
682 10624 : default_i_val=1)
683 10624 : CALL section_add_keyword(section, keyword)
684 10624 : CALL keyword_release(keyword)
685 :
686 : CALL keyword_create(keyword, __LOCATION__, name="POTENTIAL", &
687 : description="Operator to test", &
688 : default_i_val=eri_mme_coulomb, &
689 : enum_i_vals=[eri_mme_coulomb, eri_mme_yukawa, eri_mme_longrange], &
690 : enum_c_vals=s2a("COULOMB", "YUKAWA", "LONGRANGE"), &
691 10624 : enum_desc=s2a("1/r", "exp(-a*r)/r", "erf(a*r)/r"))
692 10624 : CALL section_add_keyword(section, keyword)
693 10624 : CALL keyword_release(keyword)
694 :
695 : CALL keyword_create(keyword, __LOCATION__, name="POTENTIAL_PARAM", &
696 : description="Parameter 'a' for chosen potential, ignored for Coulomb", &
697 10624 : default_r_val=0.0_dp)
698 10624 : CALL section_add_keyword(section, keyword)
699 10624 : CALL keyword_release(keyword)
700 :
701 10624 : END SUBROUTINE create_eri_mme_test_section
702 :
703 : ! **************************************************************************************************
704 : !> \brief Update local counters to gather statistics on different paths taken in MME algorithm
705 : !> (each Ewald sum can be performed over direct or reciprocal lattice vectors)
706 : !> \param param ...
707 : !> \param para_env ...
708 : !> \param G_count_2c ...
709 : !> \param R_count_2c ...
710 : !> \param GG_count_3c ...
711 : !> \param GR_count_3c ...
712 : !> \param RR_count_3c ...
713 : ! **************************************************************************************************
714 400 : 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)
715 : TYPE(cp_eri_mme_param), INTENT(INOUT) :: param
716 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
717 : INTEGER, INTENT(INOUT), OPTIONAL :: G_count_2c, R_count_2c, GG_count_3c, &
718 : GR_count_3c, RR_count_3c
719 :
720 400 : IF (PRESENT(G_count_2c)) THEN
721 192 : CALL para_env%sum(G_count_2c)
722 192 : param%G_count_2c = param%G_count_2c + G_count_2c
723 : END IF
724 :
725 400 : IF (PRESENT(R_count_2c)) THEN
726 192 : CALL para_env%sum(R_count_2c)
727 192 : param%R_count_2c = param%R_count_2c + R_count_2c
728 : END IF
729 :
730 400 : IF (PRESENT(GG_count_3c)) THEN
731 216 : CALL para_env%sum(GG_count_3c)
732 216 : param%GG_count_3c = param%GG_count_3c + GG_count_3c
733 : END IF
734 :
735 400 : IF (PRESENT(GR_count_3c)) THEN
736 216 : CALL para_env%sum(GR_count_3c)
737 216 : param%GR_count_3c = param%GR_count_3c + GR_count_3c
738 : END IF
739 :
740 400 : IF (PRESENT(RR_count_3c)) THEN
741 216 : CALL para_env%sum(RR_count_3c)
742 216 : param%RR_count_3c = param%RR_count_3c + RR_count_3c
743 : END IF
744 :
745 400 : END SUBROUTINE cp_eri_mme_update_local_counts
746 :
747 : ! **************************************************************************************************
748 : !> \brief ...
749 : !> \param para_env ...
750 : !> \param iw ...
751 : !> \param eri_mme_test_section ...
752 : ! **************************************************************************************************
753 8 : SUBROUTINE cp_eri_mme_perf_acc_test(para_env, iw, eri_mme_test_section)
754 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
755 : INTEGER, INTENT(IN) :: iw
756 : TYPE(section_vals_type), POINTER :: eri_mme_test_section
757 :
758 : INTEGER :: count_r, G_count, GG_count, GR_count, i, &
759 : ix, iy, iz, l_max, min_nR, nR, nR_xyz, &
760 : nrep, nsample, nzet, potential, &
761 : R_count, RR_count
762 : LOGICAL :: test_2c, test_3c, test_accuracy
763 : REAL(KIND=dp) :: pot_par, zet_fac, zetmax, zetmin
764 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: zet
765 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rabc
766 8 : REAL(KIND=dp), DIMENSION(:), POINTER :: cell_par
767 : TYPE(cell_type), POINTER :: box
768 200 : TYPE(cp_eri_mme_param) :: param
769 : TYPE(section_vals_type), POINTER :: eri_mme_section
770 :
771 8 : NULLIFY (box, eri_mme_section, cell_par)
772 :
773 8 : eri_mme_section => section_vals_get_subs_vals(eri_mme_test_section, "ERI_MME")
774 8 : CALL cp_eri_mme_init_read_input(eri_mme_section, param)
775 8 : CALL section_vals_val_get(eri_mme_test_section, "TEST_3C", l_val=test_3c)
776 8 : CALL section_vals_val_get(eri_mme_test_section, "TEST_2C", l_val=test_2c)
777 :
778 8 : CALL section_vals_val_get(eri_mme_test_section, "ABC", r_vals=cell_par)
779 8 : CALL section_vals_val_get(eri_mme_test_section, "MIN_NPOS", i_val=min_nR)
780 8 : CALL section_vals_val_get(eri_mme_test_section, "NREP", i_val=nrep)
781 8 : CALL section_vals_val_get(eri_mme_test_section, "CHECK_2C_ACCURACY", l_val=test_accuracy)
782 8 : CALL section_vals_val_get(eri_mme_test_section, "LMAX", i_val=l_max)
783 8 : CALL section_vals_val_get(eri_mme_test_section, "ZET_MIN", r_val=zetmin)
784 8 : CALL section_vals_val_get(eri_mme_test_section, "ZET_MAX", r_val=zetmax)
785 8 : CALL section_vals_val_get(eri_mme_test_section, "NZET", i_val=nzet)
786 8 : CALL section_vals_val_get(eri_mme_test_section, "NSAMPLE_3C", i_val=nsample)
787 8 : CALL section_vals_val_get(eri_mme_test_section, "POTENTIAL", i_val=potential)
788 8 : CALL section_vals_val_get(eri_mme_test_section, "POTENTIAL_PARAM", r_val=pot_par)
789 :
790 8 : IF (nzet <= 0) THEN
791 0 : CPABORT("Number of exponents NZET must be greater than 0.")
792 : END IF
793 :
794 8 : CALL init_orbital_pointers(l_max)
795 :
796 : ! Create ranges of zet to be tested
797 24 : ALLOCATE (zet(nzet))
798 :
799 8 : zet(1) = zetmin
800 8 : IF (nzet > 1) THEN
801 8 : zet_fac = (zetmax/zetmin)**(1.0_dp/(nzet - 1))
802 44 : DO i = 1, nzet - 1
803 44 : zet(i + 1) = zet(i)*zet_fac
804 : END DO
805 : END IF
806 :
807 : ! initialize cell
808 8 : CALL cell_create(box)
809 104 : box%hmat = 0.0_dp
810 8 : box%hmat(1, 1) = cell_par(1)
811 8 : box%hmat(2, 2) = cell_par(2)
812 8 : box%hmat(3, 3) = cell_par(3)
813 8 : CALL init_cell(box)
814 :
815 : ! Create range of rab (atomic distances) to be tested
816 8 : nR_xyz = CEILING(REAL(min_nR, KIND=dp)**(1.0_dp/3.0_dp) - 1.0E-06)
817 8 : nR = nR_xyz**3
818 :
819 24 : ALLOCATE (rabc(3, nR))
820 8 : count_r = 0
821 24 : DO ix = 1, nR_xyz
822 56 : DO iy = 1, nR_xyz
823 112 : DO iz = 1, nR_xyz
824 64 : count_r = count_r + 1
825 : ! adding 10% of cell size to positions to avoid atoms exactly at boundary or center of a cell
826 : rabc(:, count_r) = pbc([ix*ABS(cell_par(1)), &
827 : iy*ABS(cell_par(2)), &
828 : iz*ABS(cell_par(3))]/nR_xyz + &
829 288 : 0.1_dp*ABS(cell_par(:)), box)
830 : END DO
831 : END DO
832 : END DO
833 :
834 : ! initialize MME method
835 : CALL cp_eri_mme_set_params(param, box%hmat, box%orthorhombic, MINVAL(zet), MAXVAL(zet), l_max, l_max, para_env, &
836 96 : potential, pot_par)
837 :
838 8 : IF (iw > 0) WRITE (iw, '(T2, A, T61, I20)') "ERI_MME| Number of atomic distances:", nR
839 :
840 8 : G_count = 0; R_count = 0
841 8 : GG_count = 0; GR_count = 0; RR_count = 0
842 :
843 8 : IF (test_2c) CALL eri_mme_2c_perf_acc_test(param%par, l_max, zet, rabc, nrep, test_accuracy, para_env, iw, &
844 8 : potential=potential, pot_par=pot_par, G_count=G_count, R_count=R_count)
845 8 : IF (test_3c) CALL eri_mme_3c_perf_acc_test(param%par, l_max, zet, rabc, nrep, nsample, &
846 : para_env, iw, potential=potential, pot_par=pot_par, &
847 2 : GG_count=GG_count, GR_count=GR_count, RR_count=RR_count)
848 8 : CALL cp_eri_mme_update_local_counts(param, para_env, G_count, R_count, GG_count, GR_count, RR_count)
849 8 : CALL cp_eri_mme_finalize(param)
850 8 : CALL cell_release(box)
851 32 : END SUBROUTINE cp_eri_mme_perf_acc_test
852 :
853 0 : END MODULE cp_eri_mme_interface
|
