(git:374b731)
Loading...
Searching...
No Matches
cp_eri_mme_interface.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Interface to Minimax-Ewald method for periodic ERI's to be used in CP2K.
10! **************************************************************************************************
11
14 USE cell_methods, ONLY: cell_create,&
16 USE cell_types, ONLY: cell_release,&
17 cell_type,&
18 pbc
27 USE eri_mme_types, ONLY: eri_mme_coulomb,&
29 eri_mme_longrange,&
34 eri_mme_yukawa
46 USE input_val_types, ONLY: real_t
47 USE kinds, ONLY: default_string_length,&
48 dp
51 USE qs_kind_types, ONLY: get_qs_kind,&
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 :: &
73
75 MODULE PROCEDURE eri_mme_set_params_from_basis
76 MODULE PROCEDURE eri_mme_set_params_custom
77 END INTERFACE
78
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 __GNUC__ < 9 || (__GNUC__ == 9 && __GNUC_MINOR__ < 5)
83 TYPE(eri_mme_param) :: par
84#else
85 TYPE(eri_mme_param) :: par = eri_mme_param()
86#endif
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
94CONTAINS
95! **************************************************************************************************
96!> \brief Create main input section
97!> \param section ...
98!> \param default_n_minimax ...
99! **************************************************************************************************
100 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 NULLIFY (keyword, print_key, subsection)
109 cpassert(.NOT. ASSOCIATED(section))
110
111 IF (PRESENT(default_n_minimax)) THEN
112 my_default_n_minimax = default_n_minimax
113 ELSE
114 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 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 default_i_val=my_default_n_minimax)
130 CALL section_add_keyword(section, keyword)
131 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 default_r_val=300.0_dp)
138 CALL section_add_keyword(section, keyword)
139 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 default_r_val=1.0e-16_dp)
145 CALL section_add_keyword(section, keyword)
146 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 lone_keyword_l_val=.true.)
154 CALL section_add_keyword(section, keyword)
155 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 lone_keyword_l_val=.true.)
162 CALL section_add_keyword(section, keyword)
163 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 print_level=medium_print_level, filename="__STD_OUT__")
168 CALL section_add_subsection(section, print_key)
169 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 lone_keyword_l_val=.true.)
176 CALL section_add_keyword(section, keyword)
177 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 lone_keyword_l_val=.true.)
184 CALL section_add_keyword(section, keyword)
185 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 default_r_val=1.0e-06_dp)
191 CALL section_add_keyword(section, keyword)
192 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 default_i_val=1000000)
199 CALL section_add_keyword(section, keyword)
200 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 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 default_r_val=10.0_dp)
213 CALL section_add_keyword(subsection, keyword)
214 CALL keyword_release(keyword)
215
216 CALL keyword_create(keyword, __location__, &
217 name="MAX", &
218 description="Initial guess of upper bound for cutoff. ", &
219 default_r_val=10000.0_dp)
220 CALL section_add_keyword(subsection, keyword)
221 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 default_r_val=0.9_dp)
228 CALL section_add_keyword(subsection, keyword)
229 CALL keyword_release(keyword)
230
231 CALL keyword_create(keyword, __location__, &
232 name="EPS", &
233 description="Relative cutoff precision required to stop calibration. ", &
234 default_r_val=0.01_dp)
235 CALL section_add_keyword(subsection, keyword)
236 CALL keyword_release(keyword)
237
238 CALL section_add_subsection(section, subsection)
239 CALL section_release(subsection)
240
241 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 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 logger => cp_get_default_logger()
262 unit_nr = cp_print_key_unit_nr(logger, mme_section, "ERI_MME_INFO", &
263 extension=".eri_mme")
264
265 NULLIFY (subsection)
266 CALL section_vals_val_get(mme_section, "N_MINIMAX", i_val=n_minimax)
267
268 CALL section_vals_val_get(mme_section, "CUTOFF", r_val=cutoff)
269 CALL section_vals_val_get(mme_section, "SUM_PRECISION", r_val=sum_precision)
270 CALL section_vals_val_get(mme_section, "DO_CALIBRATE_CUTOFF", l_val=do_calib_cutoff)
271 CALL section_vals_val_get(mme_section, "DO_ERROR_ESTIMATE", l_val=do_error_est)
272 CALL section_vals_val_get(mme_section, "PRINT_CALIB", l_val=print_calib)
273 subsection => section_vals_get_subs_vals(mme_section, "CUTOFF_CALIB")
274 CALL section_vals_val_get(subsection, "MIN", r_val=cutoff_min)
275 CALL section_vals_val_get(subsection, "MAX", r_val=cutoff_max)
276 CALL section_vals_val_get(subsection, "EPS", r_val=cutoff_eps)
277 CALL section_vals_val_get(subsection, "DELTA", r_val=cutoff_delta)
278 CALL section_vals_val_get(mme_section, "DEBUG", l_val=debug)
279 CALL section_vals_val_get(mme_section, "DEBUG_TOLERANCE", r_val=debug_delta)
280 CALL section_vals_val_get(mme_section, "DEBUG_NSUM_MAX", i_val=debug_nsum)
281
282 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 sum_precision, debug, debug_delta, debug_nsum, unit_nr, print_calib)
287
288 param%do_calib = do_calib_cutoff
289
290 param%G_count_2c = 0
291 param%R_count_2c = 0
292 param%GG_count_3c = 0
293 param%GR_count_3c = 0
294 param%RR_count_3c = 0
295
296 param%logger => logger
297 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 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 count_2c = param%G_count_2c + param%R_count_2c
309 count_3c = param%GG_count_3c + param%GR_count_3c + param%RR_count_3c
310
311 unit_nr = param%par%unit_nr
312
313 IF (unit_nr > 0) THEN
314 IF (count_2c .GT. 0) THEN
315 WRITE (unit_nr, '(/T2, A)') "ERI_MME| Percentage of 2-center integrals evaluated in"
316 WRITE (unit_nr, '(T2, A, T76, F5.1)') "ERI_MME| G space:", &
317 100.0_dp*param%G_count_2c/count_2c
318 WRITE (unit_nr, '(T2, A, T76, F5.1/)') "ERI_MME| R space:", &
319 100.0_dp*param%R_count_2c/count_2c
320 END IF
321 IF (count_3c .GT. 0) THEN
322 WRITE (unit_nr, '(/T2, A)') "ERI_MME| Percentage of 3-center integrals evaluated in"
323 WRITE (unit_nr, '(T2, A, T76, F5.1)') "ERI_MME| G/G space:", &
324 100.0_dp*param%GG_count_3c/count_3c
325 WRITE (unit_nr, '(T2, A, T76, F5.1)') "ERI_MME| G/R space:", &
326 100.0_dp*param%GR_count_3c/count_3c
327 WRITE (unit_nr, '(T2, A, T76, F5.1/)') "ERI_MME| R/R space:", &
328 100.0_dp*param%RR_count_3c/count_3c
329 END IF
330 END IF
331
332 CALL eri_mme_release(param%par)
333 CALL cp_print_key_finished_output(unit_nr, param%logger, param%mme_section, "ERI_MME_INFO")
334 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 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 CALL timeset(routinen, handle)
366
367 CALL error_est_pgf_params_from_basis(qs_kind_set, basis_type_1, basis_type_2, &
368 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 potential, pot_par)
374
375 CALL timestop(handle)
376 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 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 IF (param%do_calib) THEN
405 IF (.NOT. param%par%is_valid) THEN
406 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 abs(param%par%zet_max - zet_max) < eps_changed .AND. &
412 param%par%l_max_zet == l_max_zet) THEN
413 param%par%do_calib_cutoff = .false.
414 ELSE
415 param%par%do_calib_cutoff = .true.
416 END IF
417 END IF
418 ELSE
419 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 potential, pot_par)
424
425 CALL eri_mme_print_info(param)
426 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 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 INTEGER, DIMENSION(:), POINTER :: npgf
454 REAL(kind=dp) :: zet_m
455 TYPE(gto_basis_set_type), POINTER :: basis_set
456
457 CALL timeset(routinen, handle)
458
459 l_m = 0
460 zet_m = 0.0_dp
461 l_zet = -1
462 zet_min = -1.0_dp
463
464 nkind = SIZE(qs_kind_set)
465 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 DO ikind = 1, nkind
470 DO ibasis = 1, nbasis
471 IF (ibasis .EQ. 1) THEN
472 basis_type = basis_type_1
473 ELSE
474 basis_type = basis_type_2
475 END IF
476 CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set, &
477 basis_type=basis_type)
478 cpassert(ASSOCIATED(basis_set))
479 npgf => basis_set%npgf
480 nset = basis_set%nset
481 l_m = max(l_m, maxval(basis_set%lmax(:)))
482 DO iset = 1, nset
483 zet_m = max(zet_m, maxval(basis_set%zet(1:npgf(iset), iset)))
484 IF (zet_min .LT. 0.0_dp) THEN
485 zet_min = minval(basis_set%zet(1:npgf(iset), iset))
486 ELSE
487 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 DO ikind = 1, nkind
495 DO ibasis = 1, nbasis
496 IF (ibasis .EQ. 1) THEN
497 basis_type = basis_type_1
498 ELSE
499 basis_type = basis_type_2
500 END IF
501 CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set, &
502 basis_type=basis_type)
503 DO iset = 1, basis_set%nset
504 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 .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 cpassert(l_zet .GE. 0)
515
516 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 l_max_zet = l_zet + 1
521 l_max = l_m + 1
522
523 CALL timestop(handle)
524 END SUBROUTINE error_est_pgf_params_from_basis
525
526! **************************************************************************************************
527!> \brief ...
528!> \param param ...
529! **************************************************************************************************
530 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 print_multigrids = .false.
538
539 logger => param%logger
540
541 unit_nr = param%par%unit_nr
542
543 IF (unit_nr > 0) THEN
544 SELECT CASE (param%par%potential)
545 CASE (eri_mme_coulomb)
546 WRITE (unit_nr, '(/T2, A)') "ERI_MME| Potential: Coulomb"
547 CASE (eri_mme_yukawa)
548 WRITE (unit_nr, '(/T2, A, ES9.2)') "ERI_MME| Potential: Yukawa with a=", param%par%pot_par
549 CASE (eri_mme_longrange)
550 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 WRITE (unit_nr, '(/T2, A, T71, F10.1)') "ERI_MME| Cutoff for ERIs [a.u.]:", param%par%cutoff
556 WRITE (unit_nr, '(/T2, A, T78, I3/)') "ERI_MME| Number of terms in minimax approximation:", param%par%n_minimax
557 END IF
558 IF (param%par%is_ortho) THEN
559 IF (unit_nr > 0) THEN
560 IF (param%par%do_error_est) THEN
561 WRITE (unit_nr, '(T2, A)') "ERI_MME| Estimated absolute error for normalized Hermite-Gaussian basis"
562 WRITE (unit_nr, '(T2, A, T72, ES9.2)') "ERI_MME| Minimax error:", param%par%err_mm
563 WRITE (unit_nr, '(T2, A, T72, ES9.2)') "ERI_MME| Cutoff error:", param%par%err_c
564 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 IF (param%par%print_calib) &
567 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 IF (unit_nr > 0) WRITE (unit_nr, *)
578
579 END SUBROUTINE eri_mme_print_info
580
581! **************************************************************************************************
582!> \brief Create input section for unit testing
583!> \param section ...
584! **************************************************************************************************
585 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 NULLIFY (keyword, subsection)
592
593 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 n_keywords=5, n_subsections=1)
598
599 CALL create_eri_mme_section(subsection)
600 CALL section_add_subsection(section, subsection)
601 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 lone_keyword_l_val=.true.)
608 CALL section_add_keyword(section, keyword)
609 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 lone_keyword_l_val=.true.)
615 CALL section_add_keyword(section, keyword)
616 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 lone_keyword_l_val=.true.)
622 CALL section_add_keyword(section, keyword)
623 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 n_var=3, type_of_var=real_t)
629 CALL section_add_keyword(section, keyword)
630 CALL keyword_release(keyword)
631
632 CALL keyword_create(keyword, __location__, name="MIN_NPOS", &
633 description="Minimum number of atomic distances to consider. ", &
634 default_i_val=8)
635 CALL section_add_keyword(section, keyword)
636 CALL keyword_release(keyword)
637
638 CALL keyword_create(keyword, __location__, name="NREP", &
639 description="Number of repeated calculation of each integral. ", &
640 default_i_val=1)
641 CALL section_add_keyword(section, keyword)
642 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 lone_keyword_l_val=.true.)
651
652 CALL section_add_keyword(section, keyword)
653 CALL keyword_release(keyword)
654
655 CALL keyword_create(keyword, __location__, name="LMAX", &
656 description="Maximum total angular momentum quantum number. ", &
657 default_i_val=6)
658 CALL section_add_keyword(section, keyword)
659 CALL keyword_release(keyword)
660
661 CALL keyword_create(keyword, __location__, name="ZET_MIN", &
662 description="Minimum exponent. ", &
663 default_r_val=0.001_dp)
664 CALL section_add_keyword(section, keyword)
665 CALL keyword_release(keyword)
666
667 CALL keyword_create(keyword, __location__, name="ZET_MAX", &
668 description="Maximum exponent. ", &
669 default_r_val=1.0_dp)
670 CALL section_add_keyword(section, keyword)
671 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 default_i_val=4)
676 CALL section_add_keyword(section, keyword)
677 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 default_i_val=1)
682 CALL section_add_keyword(section, keyword)
683 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 enum_desc=s2a("1/r", "exp(-a*r)/r", "erf(a*r)/r"))
691 CALL section_add_keyword(section, keyword)
692 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 default_r_val=0.0_dp)
697 CALL section_add_keyword(section, keyword)
698 CALL keyword_release(keyword)
699
700 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 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 IF (PRESENT(g_count_2c)) THEN
720 CALL para_env%sum(g_count_2c)
721 param%G_count_2c = param%G_count_2c + g_count_2c
722 END IF
723
724 IF (PRESENT(r_count_2c)) THEN
725 CALL para_env%sum(r_count_2c)
726 param%R_count_2c = param%R_count_2c + r_count_2c
727 END IF
728
729 IF (PRESENT(gg_count_3c)) THEN
730 CALL para_env%sum(gg_count_3c)
731 param%GG_count_3c = param%GG_count_3c + gg_count_3c
732 END IF
733
734 IF (PRESENT(gr_count_3c)) THEN
735 CALL para_env%sum(gr_count_3c)
736 param%GR_count_3c = param%GR_count_3c + gr_count_3c
737 END IF
738
739 IF (PRESENT(rr_count_3c)) THEN
740 CALL para_env%sum(rr_count_3c)
741 param%RR_count_3c = param%RR_count_3c + rr_count_3c
742 END IF
743
744 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rabc
765 REAL(kind=dp), DIMENSION(:), POINTER :: cell_par
766 TYPE(cell_type), POINTER :: box
767 TYPE(cp_eri_mme_param) :: param
768 TYPE(section_vals_type), POINTER :: eri_mme_section
769
770 NULLIFY (box, eri_mme_section, cell_par)
771
772 eri_mme_section => section_vals_get_subs_vals(eri_mme_test_section, "ERI_MME")
773 CALL cp_eri_mme_init_read_input(eri_mme_section, param)
774 CALL section_vals_val_get(eri_mme_test_section, "TEST_3C", l_val=test_3c)
775 CALL section_vals_val_get(eri_mme_test_section, "TEST_2C", l_val=test_2c)
776
777 CALL section_vals_val_get(eri_mme_test_section, "ABC", r_vals=cell_par)
778 CALL section_vals_val_get(eri_mme_test_section, "MIN_NPOS", i_val=min_nr)
779 CALL section_vals_val_get(eri_mme_test_section, "NREP", i_val=nrep)
780 CALL section_vals_val_get(eri_mme_test_section, "CHECK_2C_ACCURACY", l_val=test_accuracy)
781 CALL section_vals_val_get(eri_mme_test_section, "LMAX", i_val=l_max)
782 CALL section_vals_val_get(eri_mme_test_section, "ZET_MIN", r_val=zetmin)
783 CALL section_vals_val_get(eri_mme_test_section, "ZET_MAX", r_val=zetmax)
784 CALL section_vals_val_get(eri_mme_test_section, "NZET", i_val=nzet)
785 CALL section_vals_val_get(eri_mme_test_section, "NSAMPLE_3C", i_val=nsample)
786 CALL section_vals_val_get(eri_mme_test_section, "POTENTIAL", i_val=potential)
787 CALL section_vals_val_get(eri_mme_test_section, "POTENTIAL_PARAM", r_val=pot_par)
788
789 IF (nzet .LE. 0) &
790 cpabort("Number of exponents NZET must be greater than 0.")
791
792 CALL init_orbital_pointers(l_max)
793
794 ! Create ranges of zet to be tested
795 ALLOCATE (zet(nzet))
796
797 zet(1) = zetmin
798 IF (nzet .GT. 1) THEN
799 zet_fac = (zetmax/zetmin)**(1.0_dp/(nzet - 1))
800 DO i = 1, nzet - 1
801 zet(i + 1) = zet(i)*zet_fac
802 END DO
803 END IF
804
805 ! initialize cell
806 CALL cell_create(box)
807 box%hmat = 0.0_dp
808 box%hmat(1, 1) = cell_par(1)
809 box%hmat(2, 2) = cell_par(2)
810 box%hmat(3, 3) = cell_par(3)
811 CALL init_cell(box)
812
813 ! Create range of rab (atomic distances) to be tested
814 nr_xyz = ceiling(real(min_nr, kind=dp)**(1.0_dp/3.0_dp) - 1.0e-06)
815 nr = nr_xyz**3
816
817 ALLOCATE (rabc(3, nr))
818 count_r = 0
819 DO ix = 1, nr_xyz
820 DO iy = 1, nr_xyz
821 DO iz = 1, nr_xyz
822 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 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 potential, pot_par)
835
836 IF (iw > 0) WRITE (iw, '(T2, A, T61, I20)') "ERI_MME| Number of atomic distances:", nr
837
838 g_count = 0; r_count = 0
839 gg_count = 0; gr_count = 0; rr_count = 0
840
841 IF (test_2c) CALL eri_mme_2c_perf_acc_test(param%par, l_max, zet, rabc, nrep, test_accuracy, para_env, iw, &
842 potential=potential, pot_par=pot_par, g_count=g_count, r_count=r_count)
843 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 gg_count=gg_count, gr_count=gr_count, rr_count=rr_count)
846 CALL cp_eri_mme_update_local_counts(param, para_env, g_count, r_count, gg_count, gr_count, rr_count)
847 CALL cp_eri_mme_finalize(param)
848 CALL cell_release(box)
849 END SUBROUTINE cp_eri_mme_perf_acc_test
850
851END MODULE cp_eri_mme_interface
Handles all functions related to the CELL.
subroutine, public init_cell(cell, hmat, periodic)
Initialise/readjust a simulation cell after hmat has been changed.
subroutine, public cell_create(cell, hmat, periodic, tag)
allocates and initializes a cell
Handles all functions related to the CELL.
Definition cell_types.F:15
subroutine, public cell_release(cell)
releases the given cell (see doc/ReferenceCounting.html)
Definition cell_types.F:559
Interface to Minimax-Ewald method for periodic ERI's to be used in CP2K.
subroutine, public create_eri_mme_test_section(section)
Create input section for unit testing.
subroutine, public create_eri_mme_section(section, default_n_minimax)
Create main input section.
subroutine, public cp_eri_mme_perf_acc_test(para_env, iw, eri_mme_test_section)
...
subroutine, public cp_eri_mme_init_read_input(mme_section, param)
Read input and initialize parameter type.
subroutine, public cp_eri_mme_finalize(param)
Release eri mme data. Prints some statistics on summation methods chosen.
subroutine, public cp_eri_mme_update_local_counts(param, para_env, g_count_2c, r_count_2c, gg_count_3c, gr_count_3c, rr_count_3c)
Update local counters to gather statistics on different paths taken in MME algorithm (each Ewald sum ...
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
integer, parameter, public medium_print_level
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
subroutine, public cp_print_key_section_create(print_key_section, location, name, description, print_level, each_iter_names, each_iter_values, add_last, filename, common_iter_levels, citations, unit_str)
creates a print_key section
Methods for testing / debugging.
subroutine, public eri_mme_3c_perf_acc_test(param, l_max, zet, rabc, nrep, nsample, para_env, iw, potential, pot_par, gg_count, gr_count, rr_count)
...
subroutine, public eri_mme_2c_perf_acc_test(param, l_max, zet, rabc, nrep, test_accuracy, para_env, iw, potential, pot_par, g_count, r_count)
Unit test for performance and accuracy.
Types and initialization / release routines for Minimax-Ewald method for electron repulsion integrals...
subroutine, public eri_mme_release(param)
...
subroutine, public eri_mme_init(param, n_minimax, cutoff, do_calib_cutoff, do_error_est, cutoff_min, cutoff_max, cutoff_eps, cutoff_delta, sum_precision, debug, debug_delta, debug_nsum, unit_nr, print_calib)
...
subroutine, public eri_mme_set_params(param, hmat, is_ortho, zet_min, zet_max, l_max_zet, l_max, para_env, potential, pot_par)
Set parameters for MME method with manual specification of basis parameters. Takes care of cutoff cal...
subroutine, public eri_mme_print_grid_info(grid, grid_no, unit_nr)
...
represents keywords in an input
subroutine, public keyword_release(keyword)
releases the given keyword (see doc/ReferenceCounting.html)
subroutine, public keyword_create(keyword, location, name, description, usage, type_of_var, n_var, repeats, variants, default_val, default_l_val, default_r_val, default_lc_val, default_c_val, default_i_val, default_l_vals, default_r_vals, default_c_vals, default_i_vals, lone_keyword_val, lone_keyword_l_val, lone_keyword_r_val, lone_keyword_c_val, lone_keyword_i_val, lone_keyword_l_vals, lone_keyword_r_vals, lone_keyword_c_vals, lone_keyword_i_vals, enum_c_vals, enum_i_vals, enum, enum_strict, enum_desc, unit_str, citations, deprecation_notice, removed)
creates a keyword object
objects that represent the structure of input sections and the data contained in an input section
subroutine, public section_create(section, location, name, description, n_keywords, n_subsections, repeats, citations)
creates a list of keywords
subroutine, public section_add_keyword(section, keyword)
adds a keyword to the given section
subroutine, public section_add_subsection(section, subsection)
adds a subsection to the given section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
recursive subroutine, public section_release(section)
releases the given keyword list (see doc/ReferenceCounting.html)
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
a wrapper for basic fortran types.
integer, parameter, public real_t
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
Interface to the message passing library MPI.
Provides Cartesian and spherical orbital pointers and indices.
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Utilities for string manipulations.
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
type of a logger, at the moment it contains just a print level starting at which level it should be l...
represent a keyword in the input
represent a section of the input file
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.