(git:f2099e5)
Loading...
Searching...
No Matches
eri_mme_types.F
Go to the documentation of this file.
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 Types and initialization / release routines for Minimax-Ewald method for electron
10!> repulsion integrals.
11!> \par History
12!> 2015 09 created
13!> \author Patrick Seewald
14! **************************************************************************************************
15
17
25 USE eri_mme_util, ONLY: g_abs_min,&
27 USE kinds, ONLY: dp
28 USE mathlib, ONLY: det_3x3,&
32#include "../base/base_uses.f90"
33
34 IMPLICIT NONE
35
36 PRIVATE
37
38 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .false.
39
40 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eri_mme_types'
41
42 INTEGER, PARAMETER, PUBLIC :: n_minimax_max = 53
43
44 PUBLIC :: eri_mme_param, &
54
55 TYPE minimax_grid
56 REAL(kind=dp) :: cutoff = 0.0_dp
57 INTEGER :: n_minimax = 0
58 REAL(kind=dp), POINTER, &
59 DIMENSION(:) :: minimax_aw => null()
60 REAL(kind=dp) :: error = 0.0_dp
61 END TYPE minimax_grid
62
64 INTEGER :: n_minimax = 0
65 REAL(kind=dp), DIMENSION(3, 3) :: hmat = 0.0_dp, h_inv = 0.0_dp
66 REAL(kind=dp) :: vol = 0.0_dp
67 LOGICAL :: is_ortho = .false.
68 REAL(kind=dp) :: cutoff = 0.0_dp
69 LOGICAL :: do_calib_cutoff = .false.
70 LOGICAL :: do_error_est = .false.
71 LOGICAL :: print_calib = .false.
72 REAL(kind=dp) :: cutoff_min = 0.0_dp, cutoff_max = 0.0_dp, cutoff_delta = 0.0_dp, &
73 cutoff_eps = 0.0_dp
74 REAL(kind=dp) :: err_mm = 0.0_dp, err_c = 0.0_dp
75 REAL(kind=dp) :: mm_delta = 0.0_dp
76 REAL(kind=dp) :: g_min = 0.0_dp, r_min = 0.0_dp
77 LOGICAL :: is_valid = .false.
78 LOGICAL :: debug = .false.
79 REAL(kind=dp) :: debug_delta = 0.0_dp
80 INTEGER :: debug_nsum = 0
81 REAL(kind=dp) :: c_mm = 0.0_dp
82 INTEGER :: unit_nr = -1
83 REAL(kind=dp) :: sum_precision = 0.0_dp
84 INTEGER :: n_grids = 0
85 TYPE(minimax_grid), DIMENSION(:), &
86 ALLOCATABLE :: minimax_grid
87 REAL(kind=dp) :: zet_max = 0.0_dp, zet_min = 0.0_dp
88 INTEGER :: l_mm = -1, l_max_zet = -1
89 INTEGER :: potential = 0
90 REAL(kind=dp) :: pot_par = 0.0_dp
91 END TYPE eri_mme_param
92
93CONTAINS
94
95! **************************************************************************************************
96!> \brief ...
97!> \param param ...
98!> \param n_minimax ...
99!> \param cutoff ...
100!> \param do_calib_cutoff ...
101!> \param do_error_est ...
102!> \param cutoff_min ...
103!> \param cutoff_max ...
104!> \param cutoff_eps ...
105!> \param cutoff_delta ...
106!> \param sum_precision ...
107!> \param debug ...
108!> \param debug_delta ...
109!> \param debug_nsum ...
110!> \param unit_nr ...
111!> \param print_calib ...
112! **************************************************************************************************
113 SUBROUTINE eri_mme_init(param, n_minimax, cutoff, do_calib_cutoff, do_error_est, &
114 cutoff_min, cutoff_max, cutoff_eps, cutoff_delta, sum_precision, &
115 debug, debug_delta, debug_nsum, unit_nr, print_calib)
116 TYPE(eri_mme_param), INTENT(OUT) :: param
117 INTEGER, INTENT(IN) :: n_minimax
118 REAL(kind=dp), INTENT(IN) :: cutoff
119 LOGICAL, INTENT(IN) :: do_calib_cutoff, do_error_est
120 REAL(kind=dp), INTENT(IN) :: cutoff_min, cutoff_max, cutoff_eps, &
121 cutoff_delta, sum_precision
122 LOGICAL, INTENT(IN) :: debug
123 REAL(kind=dp), INTENT(IN) :: debug_delta
124 INTEGER, INTENT(IN) :: debug_nsum, unit_nr
125 LOGICAL, INTENT(IN) :: print_calib
126
127 CHARACTER(len=2) :: string
128
129 WRITE (string, '(I2)') n_minimax_max
130 IF (n_minimax > n_minimax_max) THEN
131 cpabort("The maximum allowed number of minimax points N_MINIMAX is "//trim(string))
132 END IF
133
134 param%n_minimax = n_minimax
135 param%n_grids = 1
136 param%cutoff = cutoff
137 param%do_calib_cutoff = do_calib_cutoff
138 param%do_error_est = do_error_est
139 param%cutoff_min = cutoff_min
140 param%cutoff_max = cutoff_max
141 param%cutoff_eps = cutoff_eps
142 param%cutoff_delta = cutoff_delta
143 param%sum_precision = sum_precision
144 param%debug = debug
145 param%debug_delta = debug_delta
146 param%debug_nsum = debug_nsum
147 param%print_calib = print_calib
148 param%unit_nr = unit_nr
149 param%err_mm = -1.0_dp
150 param%err_c = -1.0_dp
151
152 param%is_valid = .false.
153 ALLOCATE (param%minimax_grid(param%n_grids))
154 END SUBROUTINE eri_mme_init
155
156! **************************************************************************************************
157!> \brief Set parameters for MME method with manual specification of basis parameters.
158!> Takes care of cutoff calibration if requested.
159!> \param param ...
160!> \param hmat ...
161!> \param is_ortho ...
162!> \param zet_min Exponent used to estimate error of minimax approximation.
163!> \param zet_max Exponent used to estimate error of finite cutoff.
164!> \param l_max_zet Total ang. mom. quantum numbers l to be combined with exponents in
165!> zet_max.
166!> \param l_max Maximum total angular momentum quantum number
167!> \param para_env ...
168!> \param potential potential to use. Accepts the following values:
169!> 1: coulomb potential V(r)=1/r
170!> 2: yukawa potential V(r)=e(-a*r)/r
171!> 3: long-range coulomb erf(a*r)/r
172!> \param pot_par potential parameter a for yukawa V(r)=e(-a*r)/r or long-range coulomb V(r)=erf(a*r)/r
173! **************************************************************************************************
174 SUBROUTINE eri_mme_set_params(param, hmat, is_ortho, zet_min, zet_max, l_max_zet, l_max, para_env, &
175 potential, pot_par)
176 TYPE(eri_mme_param), INTENT(INOUT) :: param
177 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN) :: hmat
178 LOGICAL, INTENT(IN) :: is_ortho
179 REAL(kind=dp), INTENT(IN) :: zet_min, zet_max
180 INTEGER, INTENT(IN) :: l_max_zet, l_max
181 TYPE(mp_para_env_type), INTENT(IN), OPTIONAL :: para_env
182 INTEGER, INTENT(IN), OPTIONAL :: potential
183 REAL(kind=dp), INTENT(IN), OPTIONAL :: pot_par
184
185 CHARACTER(LEN=*), PARAMETER :: routinen = 'eri_mme_set_params'
186
187 INTEGER :: handle, l_mm, n_grids
188 LOGICAL :: s_only
189 REAL(kind=dp) :: cutoff
190 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: minimax_aw
191
192 CALL timeset(routinen, handle)
193
194 ! Note: in MP2 default logger hacked and does not use global default print level
195 s_only = l_max == 0
196
197 CALL init_orbital_pointers(3*l_max) ! allow for orbital pointers of combined index
198
199 ! l values for minimax error estimate (l_mm) and for cutoff error estimate (l_max_zet)
200 l_mm = merge(0, 1, s_only)
201
202 ! cell info
203 ! Note: we recompute basic quantities from hmat to avoid dependency on cp2k cell type
204 param%hmat = hmat
205 param%h_inv = inv_3x3(hmat)
206 param%vol = abs(det_3x3(hmat))
207 param%is_ortho = is_ortho
208
209 ! Minimum lattice vectors
210 param%G_min = g_abs_min(param%h_inv)
211 param%R_min = r_abs_min(param%hmat)
212
213 ! Minimum and maximum exponents
214 param%zet_max = zet_max
215 param%zet_min = zet_min
216 param%l_max_zet = l_max_zet
217 param%l_mm = l_mm
218
219 ! cutoff calibration not yet implemented for general cell
220 IF (.NOT. param%is_ortho) THEN
221 param%do_calib_cutoff = .false.
222 param%do_error_est = .false.
223 END IF
224
225 n_grids = param%n_grids
226
227 ! Cutoff calibration and error estimate for orthorhombic cell
228 ! Here we assume Coulomb potential which should give an upper bound error also for the other
229 ! potentials
230 IF (param%do_calib_cutoff) THEN
231 CALL calibrate_cutoff(param%hmat, param%h_inv, param%G_min, param%vol, &
232 zet_min, l_mm, zet_max, l_max_zet, param%n_minimax, &
233 param%cutoff_min, param%cutoff_max, param%cutoff_eps, &
234 param%cutoff_delta, cutoff, param%err_mm, param%err_c, &
235 param%C_mm, para_env, param%print_calib, param%unit_nr)
236
237 param%cutoff = cutoff
238 ELSE IF (param%do_error_est) THEN
239 ALLOCATE (minimax_aw(2*param%n_minimax))
240 CALL cutoff_minimax_error(param%cutoff, param%hmat, param%h_inv, param%vol, param%G_min, &
241 zet_min, l_mm, zet_max, l_max_zet, param%n_minimax, &
242 minimax_aw, param%err_mm, param%err_c, param%C_mm, para_env)
243 DEALLOCATE (minimax_aw)
244 END IF
245
246 param%is_valid = .true.
247
248 CALL eri_mme_set_potential(param, potential=potential, pot_par=pot_par)
249
250 CALL timestop(handle)
251 END SUBROUTINE eri_mme_set_params
252
253! **************************************************************************************************
254!> \brief ...
255!> \param param ...
256!> \param potential potential to use. Accepts the following values:
257!> 1: coulomb potential V(r)=1/r
258!> 2: yukawa potential V(r)=e(-a*r)/r
259!> 3: long-range coulomb erf(a*r)/r
260!> \param pot_par potential parameter a for yukawa V(r)=e(-a*r)/r or long-range coulomb V(r)=erf(a*r)/r
261! **************************************************************************************************
262 SUBROUTINE eri_mme_set_potential(param, potential, pot_par)
263 TYPE(eri_mme_param), INTENT(INOUT) :: param
264 INTEGER, INTENT(IN), OPTIONAL :: potential
265 REAL(kind=dp), INTENT(IN), OPTIONAL :: pot_par
266
267 REAL(kind=dp), PARAMETER :: zet_d = 1.0e-12_dp
268
269 REAL(kind=dp) :: cutoff_max, cutoff_min, cutoff_rel
270 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: minimax_aw
271
272 cpassert(param%is_valid)
273
274 IF (PRESENT(potential)) THEN
275 param%potential = potential
276 ELSE
277 param%potential = eri_mme_coulomb
278 END IF
279
280 IF (PRESENT(pot_par)) THEN
281 param%pot_par = pot_par
282 ELSE
283 param%pot_par = 0.0_dp
284 END IF
285
286 ALLOCATE (minimax_aw(2*param%n_minimax))
287
288 CALL minimax_error(param%cutoff, param%hmat, param%vol, param%G_min, param%zet_min, param%l_mm, &
289 param%n_minimax, minimax_aw, param%err_mm, param%mm_delta, potential=potential, pot_par=pot_par)
290
291 DEALLOCATE (minimax_aw)
292
293 cpassert(param%zet_max + zet_d > param%zet_min)
294 cpassert(param%n_grids >= 1)
295
296 cutoff_max = param%cutoff
297 cutoff_rel = param%cutoff/param%zet_max
298 cutoff_min = param%zet_min*cutoff_rel
299
300 CALL eri_mme_destroy_minimax_grids(param%minimax_grid)
301 ALLOCATE (param%minimax_grid(param%n_grids))
302
303 CALL eri_mme_create_minimax_grids(param%n_grids, param%minimax_grid, param%n_minimax, &
304 cutoff_max, cutoff_min, param%G_min, &
305 param%mm_delta, potential=potential, pot_par=pot_par)
306
307 END SUBROUTINE eri_mme_set_potential
308
309! **************************************************************************************************
310!> \brief ...
311!> \param n_grids ...
312!> \param minimax_grids ...
313!> \param n_minimax ...
314!> \param cutoff_max ...
315!> \param cutoff_min ...
316!> \param G_min ...
317!> \param target_error ...
318!> \param potential ...
319!> \param pot_par ...
320! **************************************************************************************************
321 SUBROUTINE eri_mme_create_minimax_grids(n_grids, minimax_grids, n_minimax, &
322 cutoff_max, cutoff_min, G_min, &
323 target_error, potential, pot_par)
324 INTEGER, INTENT(IN) :: n_grids
325 TYPE(minimax_grid), DIMENSION(n_grids), &
326 INTENT(OUT) :: minimax_grids
327 INTEGER, INTENT(IN) :: n_minimax
328 REAL(kind=dp), INTENT(IN) :: cutoff_max, cutoff_min, g_min, &
329 target_error
330 INTEGER, INTENT(IN), OPTIONAL :: potential
331 REAL(kind=dp), INTENT(IN), OPTIONAL :: pot_par
332
333 REAL(kind=dp), PARAMETER :: err_mm_d = 1.0e-12_dp
334
335 INTEGER :: i_grid, n_mm
336 REAL(kind=dp) :: cutoff, cutoff_delta, err_mm, err_mm_prev
337 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: minimax_aw, minimax_aw_prev
338
339 cutoff_delta = (cutoff_max/cutoff_min)**(1.0_dp/(n_grids))
340 cutoff = cutoff_max
341
342 ALLOCATE (minimax_aw(2*n_minimax))
343 ! for first grid (for max. cutoff) always use default n_minimax
344 CALL get_minimax_coeff_v_gspace(n_minimax, cutoff, g_min, minimax_aw, err_minimax=err_mm, &
345 potential=potential, pot_par=pot_par)
346 cpassert(err_mm < 1.1_dp*target_error + err_mm_d)
347 CALL create_minimax_grid(minimax_grids(n_grids), cutoff, n_minimax, minimax_aw, err_mm)
348 DEALLOCATE (minimax_aw)
349
350 DO i_grid = n_grids - 1, 1, -1
351 DO n_mm = n_minimax, 1, -1
352 ALLOCATE (minimax_aw(2*n_mm))
353 CALL get_minimax_coeff_v_gspace(n_mm, cutoff, g_min, minimax_aw, err_minimax=err_mm, &
354 potential=potential, pot_par=pot_par)
355
356 IF (err_mm > 1.1_dp*target_error) THEN
357 cpassert(n_mm /= n_minimax)
358 CALL create_minimax_grid(minimax_grids(i_grid), cutoff, n_mm + 1, minimax_aw_prev, err_mm_prev)
359
360 DEALLOCATE (minimax_aw)
361 EXIT
362 END IF
363
364 IF (ALLOCATED(minimax_aw_prev)) DEALLOCATE (minimax_aw_prev)
365 ALLOCATE (minimax_aw_prev(2*n_mm))
366 minimax_aw_prev(:) = minimax_aw(:)
367 DEALLOCATE (minimax_aw)
368 err_mm_prev = err_mm
369 END DO
370 cutoff = cutoff/cutoff_delta
371 END DO
372 END SUBROUTINE eri_mme_create_minimax_grids
373
374! **************************************************************************************************
375!> \brief ...
376!> \param minimax_grids ...
377! **************************************************************************************************
378 SUBROUTINE eri_mme_destroy_minimax_grids(minimax_grids)
379 TYPE(minimax_grid), ALLOCATABLE, DIMENSION(:), &
380 INTENT(INOUT) :: minimax_grids
381
382 INTEGER :: igrid
383
384 IF (ALLOCATED(minimax_grids)) THEN
385 DO igrid = 1, SIZE(minimax_grids)
386 IF (ASSOCIATED(minimax_grids(igrid)%minimax_aw)) THEN
387 DEALLOCATE (minimax_grids(igrid)%minimax_aw)
388 END IF
389 END DO
390 DEALLOCATE (minimax_grids)
391 END IF
392 END SUBROUTINE eri_mme_destroy_minimax_grids
393
394! **************************************************************************************************
395!> \brief ...
396!> \param grid ...
397!> \param cutoff ...
398!> \param n_minimax ...
399!> \param minimax_aw ...
400!> \param error ...
401! **************************************************************************************************
402 SUBROUTINE create_minimax_grid(grid, cutoff, n_minimax, minimax_aw, error)
403 TYPE(minimax_grid), INTENT(OUT) :: grid
404 REAL(kind=dp), INTENT(IN) :: cutoff
405 INTEGER, INTENT(IN) :: n_minimax
406 REAL(kind=dp), DIMENSION(2*n_minimax), INTENT(IN) :: minimax_aw
407 REAL(kind=dp), INTENT(IN) :: error
408
409 grid%cutoff = cutoff
410 grid%n_minimax = n_minimax
411 ALLOCATE (grid%minimax_aw(2*n_minimax))
412 grid%minimax_aw(:) = minimax_aw(:)
413 grid%error = error
414
415 END SUBROUTINE create_minimax_grid
416
417! **************************************************************************************************
418!> \brief ...
419!> \param grids ...
420!> \param cutoff ...
421!> \param n_minimax ...
422!> \param minimax_aw ...
423!> \param grid_no ...
424! **************************************************************************************************
425 SUBROUTINE get_minimax_from_cutoff(grids, cutoff, n_minimax, minimax_aw, grid_no)
426 TYPE(minimax_grid), DIMENSION(:), INTENT(IN) :: grids
427 REAL(kind=dp), INTENT(IN) :: cutoff
428 INTEGER, INTENT(OUT) :: n_minimax
429 REAL(kind=dp), DIMENSION(:), INTENT(OUT), POINTER :: minimax_aw
430 INTEGER, INTENT(OUT) :: grid_no
431
432 INTEGER :: igrid
433
434 grid_no = 0
435 DO igrid = 1, SIZE(grids)
436 IF (grids(igrid)%cutoff >= cutoff/2) THEN
437 n_minimax = grids(igrid)%n_minimax
438 minimax_aw => grids(igrid)%minimax_aw
439 grid_no = igrid
440 EXIT
441 END IF
442 END DO
443 IF (grid_no == 0) THEN
444 igrid = SIZE(grids)
445 n_minimax = grids(igrid)%n_minimax
446 minimax_aw => grids(igrid)%minimax_aw
447 grid_no = igrid
448 END IF
449 END SUBROUTINE get_minimax_from_cutoff
450
451! **************************************************************************************************
452!> \brief ...
453!> \param grid ...
454!> \param grid_no ...
455!> \param unit_nr ...
456! **************************************************************************************************
457 SUBROUTINE eri_mme_print_grid_info(grid, grid_no, unit_nr)
458 TYPE(minimax_grid), INTENT(IN) :: grid
459 INTEGER, INTENT(IN) :: grid_no, unit_nr
460
461 IF (unit_nr > 0) THEN
462 WRITE (unit_nr, '(T2, A, 1X, I2)') "ERI_MME | Info for grid no.", grid_no
463 WRITE (unit_nr, '(T2, A, 1X, ES9.2)') "ERI_MME | Cutoff", grid%cutoff
464 WRITE (unit_nr, '(T2, A, 1X, I2)') "ERI_MME | Number of minimax points", grid%n_minimax
465 WRITE (unit_nr, '(T2, A, 1X, 2ES9.2)') "ERI_MME | minimax error", grid%error
466 WRITE (unit_nr, *)
467 END IF
468
469 END SUBROUTINE eri_mme_print_grid_info
470
471! **************************************************************************************************
472!> \brief ...
473!> \param param ...
474! **************************************************************************************************
475 SUBROUTINE eri_mme_release(param)
476 TYPE(eri_mme_param), INTENT(INOUT) :: param
477
478 IF (ALLOCATED(param%minimax_grid)) THEN
479 CALL eri_mme_destroy_minimax_grids(param%minimax_grid)
480 END IF
481 END SUBROUTINE eri_mme_release
482
483END MODULE eri_mme_types
Methods aiming for error estimate and automatic cutoff calibration. integrals.
subroutine, public cutoff_minimax_error(cutoff, hmat, h_inv, vol, g_min, zet_min, l_mm, zet_max, l_max_zet, n_minimax, minimax_aw, err_mm, err_ctff, c_mm, para_env)
Compute upper bounds for the errors of 2-center ERI's (P|P) due to minimax approximation and due to f...
subroutine, public calibrate_cutoff(hmat, h_inv, g_min, vol, zet_min, l_mm, zet_max, l_max_zet, n_minimax, cutoff_l, cutoff_r, tol, delta, cutoff, err_mm, err_c, c_mm, para_env, print_calib, unit_nr)
Find optimal cutoff minimizing errors due to minimax approximation and due to finite cutoff using bis...
subroutine, public minimax_error(cutoff, hmat, vol, g_min, zet_min, l_mm, n_minimax, minimax_aw, err_mm, delta_mm, potential, pot_par)
Minimax error, simple analytical formula Note minimax error may blow up for small exponents....
Methods related to properties of Hermite and Cartesian Gaussian functions.
integer, parameter, public eri_mme_longrange
integer, parameter, public eri_mme_coulomb
subroutine, public get_minimax_coeff_v_gspace(n_minimax, cutoff, g_min, minimax_aw, potential, pot_par, err_minimax)
Get minimax coefficient a_i and w_i for approximating 1/G^2 by sum_i w_i exp(-a_i G^2).
integer, parameter, public eri_mme_yukawa
Types and initialization / release routines for Minimax-Ewald method for electron repulsion integrals...
subroutine, public eri_mme_release(param)
...
subroutine, public get_minimax_from_cutoff(grids, cutoff, n_minimax, minimax_aw, grid_no)
...
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_potential(param, potential, pot_par)
...
integer, parameter, public n_minimax_max
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)
...
Some utility methods used in different contexts.
real(kind=dp) function, public g_abs_min(h_inv)
Find minimum length of G vectors, for a general (not necessarily orthorhombic) cell.
real(kind=dp) function, public r_abs_min(hmat)
Find minimum length of R vectors, for a general (not necessarily orthorhombic) cell.
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
pure real(kind=dp) function, dimension(3, 3), public inv_3x3(a)
Returns the inverse of the 3 x 3 matrix a.
Definition mathlib.F:524
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.
stores all the informations relevant to an mpi environment