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 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 :
16 : MODULE eri_mme_types
17 :
18 : USE eri_mme_error_control, ONLY: calibrate_cutoff,&
19 : cutoff_minimax_error,&
20 : minimax_error
21 : USE eri_mme_gaussian, ONLY: eri_mme_coulomb,&
22 : eri_mme_longrange,&
23 : eri_mme_yukawa,&
24 : get_minimax_coeff_v_gspace
25 : USE eri_mme_util, ONLY: G_abs_min,&
26 : R_abs_min
27 : USE kinds, ONLY: dp
28 : USE mathlib, ONLY: det_3x3,&
29 : inv_3x3
30 : USE message_passing, ONLY: mp_para_env_type
31 : USE orbital_pointers, ONLY: init_orbital_pointers
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, &
45 : eri_mme_init, &
46 : eri_mme_release, &
47 : eri_mme_set_params, &
48 : eri_mme_print_grid_info, &
49 : get_minimax_from_cutoff, &
50 : eri_mme_coulomb, &
51 : eri_mme_yukawa, &
52 : eri_mme_longrange, &
53 : eri_mme_set_potential
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 :
63 : TYPE eri_mme_param
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 :
93 : CONTAINS
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 1782 : 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 66 : WRITE (string, '(I2)') n_minimax_max
130 66 : IF (n_minimax > n_minimax_max) THEN
131 0 : CPABORT("The maximum allowed number of minimax points N_MINIMAX is "//TRIM(string))
132 : END IF
133 :
134 66 : param%n_minimax = n_minimax
135 66 : param%n_grids = 1
136 66 : param%cutoff = cutoff
137 66 : param%do_calib_cutoff = do_calib_cutoff
138 66 : param%do_error_est = do_error_est
139 66 : param%cutoff_min = cutoff_min
140 66 : param%cutoff_max = cutoff_max
141 66 : param%cutoff_eps = cutoff_eps
142 66 : param%cutoff_delta = cutoff_delta
143 66 : param%sum_precision = sum_precision
144 66 : param%debug = debug
145 66 : param%debug_delta = debug_delta
146 66 : param%debug_nsum = debug_nsum
147 66 : param%print_calib = print_calib
148 66 : param%unit_nr = unit_nr
149 66 : param%err_mm = -1.0_dp
150 66 : param%err_c = -1.0_dp
151 :
152 66 : param%is_valid = .FALSE.
153 132 : ALLOCATE (param%minimax_grid(param%n_grids))
154 66 : 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 158 : 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 158 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: minimax_aw
191 :
192 158 : CALL timeset(routineN, handle)
193 :
194 : ! Note: in MP2 default logger hacked and does not use global default print level
195 158 : s_only = l_max == 0
196 :
197 158 : 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 308 : 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 2054 : param%hmat = hmat
205 2054 : param%h_inv = inv_3x3(hmat)
206 158 : param%vol = ABS(det_3x3(hmat))
207 158 : param%is_ortho = is_ortho
208 :
209 : ! Minimum lattice vectors
210 158 : param%G_min = G_abs_min(param%h_inv)
211 158 : param%R_min = R_abs_min(param%hmat)
212 :
213 : ! Minimum and maximum exponents
214 158 : param%zet_max = zet_max
215 158 : param%zet_min = zet_min
216 158 : param%l_max_zet = l_max_zet
217 158 : param%l_mm = l_mm
218 :
219 : ! cutoff calibration not yet implemented for general cell
220 158 : IF (.NOT. param%is_ortho) THEN
221 36 : param%do_calib_cutoff = .FALSE.
222 36 : param%do_error_est = .FALSE.
223 : END IF
224 :
225 158 : 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 158 : 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 84 : param%C_mm, para_env, param%print_calib, param%unit_nr)
236 :
237 84 : param%cutoff = cutoff
238 74 : ELSE IF (param%do_error_est) THEN
239 114 : 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 38 : minimax_aw, param%err_mm, param%err_c, param%C_mm, para_env)
243 38 : DEALLOCATE (minimax_aw)
244 : END IF
245 :
246 158 : param%is_valid = .TRUE.
247 :
248 158 : CALL eri_mme_set_potential(param, potential=potential, pot_par=pot_par)
249 :
250 158 : CALL timestop(handle)
251 158 : 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 85492 : 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 85492 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: minimax_aw
271 :
272 85492 : CPASSERT(param%is_valid)
273 :
274 85492 : IF (PRESENT(potential)) THEN
275 85342 : param%potential = potential
276 : ELSE
277 150 : param%potential = eri_mme_coulomb
278 : END IF
279 :
280 85492 : IF (PRESENT(pot_par)) THEN
281 85342 : param%pot_par = pot_par
282 : ELSE
283 150 : param%pot_par = 0.0_dp
284 : END IF
285 :
286 256476 : 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 85492 : param%n_minimax, minimax_aw, param%err_mm, param%mm_delta, potential=potential, pot_par=pot_par)
290 :
291 85492 : DEALLOCATE (minimax_aw)
292 :
293 85492 : CPASSERT(param%zet_max + zet_d > param%zet_min)
294 85492 : CPASSERT(param%n_grids >= 1)
295 :
296 85492 : cutoff_max = param%cutoff
297 85492 : cutoff_rel = param%cutoff/param%zet_max
298 85492 : cutoff_min = param%zet_min*cutoff_rel
299 :
300 85492 : CALL eri_mme_destroy_minimax_grids(param%minimax_grid)
301 341968 : 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 85492 : param%mm_delta, potential=potential, pot_par=pot_par)
306 :
307 85492 : 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 170984 : 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 85492 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: minimax_aw, minimax_aw_prev
338 :
339 85492 : cutoff_delta = (cutoff_max/cutoff_min)**(1.0_dp/(n_grids))
340 85492 : cutoff = cutoff_max
341 :
342 256476 : 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 85492 : potential=potential, pot_par=pot_par)
346 85492 : CPASSERT(err_mm < 1.1_dp*target_error + err_mm_d)
347 85492 : CALL create_minimax_grid(minimax_grids(n_grids), cutoff, n_minimax, minimax_aw, err_mm)
348 85492 : DEALLOCATE (minimax_aw)
349 :
350 85492 : DO i_grid = n_grids - 1, 1, -1
351 0 : DO n_mm = n_minimax, 1, -1
352 0 : 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 0 : potential=potential, pot_par=pot_par)
355 :
356 0 : IF (err_mm > 1.1_dp*target_error) THEN
357 0 : CPASSERT(n_mm /= n_minimax)
358 0 : CALL create_minimax_grid(minimax_grids(i_grid), cutoff, n_mm + 1, minimax_aw_prev, err_mm_prev)
359 :
360 0 : DEALLOCATE (minimax_aw)
361 0 : EXIT
362 : END IF
363 :
364 0 : IF (ALLOCATED(minimax_aw_prev)) DEALLOCATE (minimax_aw_prev)
365 0 : ALLOCATE (minimax_aw_prev(2*n_mm))
366 0 : minimax_aw_prev(:) = minimax_aw(:)
367 0 : DEALLOCATE (minimax_aw)
368 0 : err_mm_prev = err_mm
369 : END DO
370 85492 : cutoff = cutoff/cutoff_delta
371 : END DO
372 170984 : END SUBROUTINE eri_mme_create_minimax_grids
373 :
374 : ! **************************************************************************************************
375 : !> \brief ...
376 : !> \param minimax_grids ...
377 : ! **************************************************************************************************
378 85558 : 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 85558 : IF (ALLOCATED(minimax_grids)) THEN
385 171116 : DO igrid = 1, SIZE(minimax_grids)
386 171116 : IF (ASSOCIATED(minimax_grids(igrid)%minimax_aw)) THEN
387 85492 : DEALLOCATE (minimax_grids(igrid)%minimax_aw)
388 : END IF
389 : END DO
390 85558 : DEALLOCATE (minimax_grids)
391 : END IF
392 85558 : 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 85492 : 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 85492 : grid%cutoff = cutoff
410 85492 : grid%n_minimax = n_minimax
411 341968 : ALLOCATE (grid%minimax_aw(2*n_minimax))
412 2745732 : grid%minimax_aw(:) = minimax_aw(:)
413 85492 : grid%error = error
414 :
415 85492 : 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 210488 : 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 210488 : grid_no = 0
435 226830 : DO igrid = 1, SIZE(grids)
436 226830 : IF (grids(igrid)%cutoff >= cutoff/2) THEN
437 194146 : n_minimax = grids(igrid)%n_minimax
438 194146 : minimax_aw => grids(igrid)%minimax_aw
439 194146 : grid_no = igrid
440 194146 : EXIT
441 : END IF
442 : END DO
443 210488 : IF (grid_no == 0) THEN
444 16342 : igrid = SIZE(grids)
445 16342 : n_minimax = grids(igrid)%n_minimax
446 16342 : minimax_aw => grids(igrid)%minimax_aw
447 16342 : grid_no = igrid
448 : END IF
449 210488 : END SUBROUTINE get_minimax_from_cutoff
450 :
451 : ! **************************************************************************************************
452 : !> \brief ...
453 : !> \param grid ...
454 : !> \param grid_no ...
455 : !> \param unit_nr ...
456 : ! **************************************************************************************************
457 0 : 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 0 : IF (unit_nr > 0) THEN
462 0 : WRITE (unit_nr, '(T2, A, 1X, I2)') "ERI_MME | Info for grid no.", grid_no
463 0 : WRITE (unit_nr, '(T2, A, 1X, ES9.2)') "ERI_MME | Cutoff", grid%cutoff
464 0 : WRITE (unit_nr, '(T2, A, 1X, I2)') "ERI_MME | Number of minimax points", grid%n_minimax
465 0 : WRITE (unit_nr, '(T2, A, 1X, 2ES9.2)') "ERI_MME | minimax error", grid%error
466 0 : WRITE (unit_nr, *)
467 : END IF
468 :
469 0 : END SUBROUTINE eri_mme_print_grid_info
470 :
471 : ! **************************************************************************************************
472 : !> \brief ...
473 : !> \param param ...
474 : ! **************************************************************************************************
475 66 : SUBROUTINE eri_mme_release(param)
476 : TYPE(eri_mme_param), INTENT(INOUT) :: param
477 :
478 66 : IF (ALLOCATED(param%minimax_grid)) THEN
479 66 : CALL eri_mme_destroy_minimax_grids(param%minimax_grid)
480 : END IF
481 66 : END SUBROUTINE eri_mme_release
482 :
483 0 : END MODULE eri_mme_types
|
