(git:4cf809f)
Loading...
Searching...
No Matches
ewald_environment_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!> \par History
10!> JGH FEB-13-2007 : Distributed/replicated realspace grids
11!> Teodoro Laino [tlaino] - University of Zurich - 12.2007
12!> \author CJM NOV-30-2003
13! **************************************************************************************************
14MODULE ewald_environment_types
15 USE cell_types, ONLY: use_perd_none,&
16 use_perd_x,&
17 use_perd_xy,&
18 use_perd_xyz,&
19 use_perd_xz,&
20 use_perd_y,&
21 use_perd_yz,&
22 use_perd_z
23 USE cp_log_handling, ONLY: cp_get_default_logger,&
24 cp_logger_type
25 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
26 cp_print_key_unit_nr
27 USE cp_units, ONLY: cp_unit_from_cp2k
28 USE input_cp2k_poisson, ONLY: create_ewald_section
29 USE input_enumeration_types, ONLY: enum_i2c,&
30 enumeration_type
31 USE input_keyword_types, ONLY: keyword_get,&
32 keyword_type
33 USE input_section_types, ONLY: section_get_keyword,&
34 section_release,&
35 section_type,&
36 section_vals_get_subs_vals,&
37 section_vals_release,&
38 section_vals_retain,&
39 section_vals_type,&
40 section_vals_val_get
41 USE kinds, ONLY: dp
42 USE mathconstants, ONLY: twopi
43 USE mathlib, ONLY: det_3x3
44 USE message_passing, ONLY: mp_comm_type,&
45 mp_para_env_release,&
46 mp_para_env_type
47 USE pw_grid_info, ONLY: pw_grid_n_for_fft
48 USE pw_poisson_types, ONLY: do_ewald_ewald,&
49 do_ewald_none,&
50 do_ewald_pme,&
51 do_ewald_spme
52#include "./base/base_uses.f90"
53
54 IMPLICIT NONE
55 PRIVATE
56
57! **************************************************************************************************
58!> \brief to build arrays of pointers
59!> \param ewald_env the pointer to the ewald_env
60!> \par History
61!> 11/03
62!> \author CJM
63! **************************************************************************************************
64 TYPE ewald_environment_type
65 PRIVATE
66 LOGICAL :: do_multipoles = .false. ! Flag for using the multipole code
67 INTEGER :: do_ipol = -1 ! Solver for induced dipoles
68 INTEGER :: max_multipole = -1 ! max expansion in the multipoles
69 INTEGER :: max_ipol_iter = -1 ! max number of interaction for induced dipoles
70 INTEGER :: ewald_type = -1 ! type of ewald
71 INTEGER :: gmax(3) = -1 ! max Miller index
72 INTEGER :: ns_max = -1 ! # grid points for small grid (PME)
73 INTEGER :: o_spline = -1 ! order of spline (SPME)
74 REAL(kind=dp) :: precs = 0.0_dp ! precision achieved when evaluating the real-space part
75 REAL(kind=dp) :: alpha = 0.0_dp, rcut = 0.0_dp ! ewald alpha and real-space cutoff
76 REAL(kind=dp) :: epsilon = 0.0_dp ! tolerance for small grid (PME)
77 REAL(kind=dp) :: eps_pol = 0.0_dp ! tolerance for convergence of induced dipoles
78 TYPE(mp_para_env_type), POINTER :: para_env => null()
79 TYPE(section_vals_type), POINTER :: poisson_section => null()
80 ! interaction cutoff is required to make the electrostatic interaction
81 ! continuous at a pair distance equal to rcut. this is ignored by the
82 ! multipole code and is otherwise only active when SHIFT_CUTOFF is used.
83 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: interaction_cutoffs => null()
84 ! store current cell, used to rebuild lazily.
85 REAL(kind=dp), DIMENSION(3, 3) :: cell_hmat = -1.0_dp
86 END TYPE ewald_environment_type
87
88! *** Public data types ***
89 PUBLIC :: ewald_environment_type
90
91! *** Public subroutines ***
92 PUBLIC :: ewald_env_get, &
93 ewald_env_set, &
94 ewald_env_create, &
95 ewald_env_release, &
96 read_ewald_section, &
97 read_ewald_section_tb
98
99 CHARACTER(len=*), PARAMETER, PRIVATE :: modulen = 'ewald_environment_types'
100
101CONTAINS
102
103! **************************************************************************************************
104!> \brief Purpose: Get the EWALD environment.
105!> \param ewald_env the pointer to the ewald_env
106!> \param ewald_type ...
107!> \param alpha ...
108!> \param eps_pol ...
109!> \param epsilon ...
110!> \param gmax ...
111!> \param ns_max ...
112!> \param o_spline ...
113!> \param group ...
114!> \param para_env ...
115!> \param poisson_section ...
116!> \param precs ...
117!> \param rcut ...
118!> \param do_multipoles ...
119!> \param max_multipole ...
120!> \param do_ipol ...
121!> \param max_ipol_iter ...
122!> \param interaction_cutoffs ...
123!> \param cell_hmat ...
124!> \par History
125!> 11/03
126!> \author CJM
127! **************************************************************************************************
128 SUBROUTINE ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, &
129 gmax, ns_max, o_spline, group, para_env, poisson_section, precs, &
130 rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, &
131 interaction_cutoffs, cell_hmat)
132 TYPE(ewald_environment_type), INTENT(IN) :: ewald_env
133 INTEGER, OPTIONAL :: ewald_type
134 REAL(kind=dp), OPTIONAL :: alpha, eps_pol, epsilon
135 INTEGER, OPTIONAL :: gmax(3), ns_max, o_spline
136 TYPE(mp_comm_type), INTENT(OUT), OPTIONAL :: group
137 TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
138 TYPE(section_vals_type), OPTIONAL, POINTER :: poisson_section
139 REAL(kind=dp), OPTIONAL :: precs, rcut
140 LOGICAL, INTENT(OUT), OPTIONAL :: do_multipoles
141 INTEGER, INTENT(OUT), OPTIONAL :: max_multipole, do_ipol, max_ipol_iter
142 REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL, &
143 POINTER :: interaction_cutoffs
144 REAL(kind=dp), DIMENSION(3, 3), OPTIONAL :: cell_hmat
145
146 IF (PRESENT(ewald_type)) ewald_type = ewald_env%ewald_type
147 IF (PRESENT(do_multipoles)) do_multipoles = ewald_env%do_multipoles
148 IF (PRESENT(do_ipol)) do_ipol = ewald_env%do_ipol
149 IF (PRESENT(max_multipole)) max_multipole = ewald_env%max_multipole
150 IF (PRESENT(max_ipol_iter)) max_ipol_iter = ewald_env%max_ipol_iter
151 IF (PRESENT(alpha)) alpha = ewald_env%alpha
152 IF (PRESENT(precs)) precs = ewald_env%precs
153 IF (PRESENT(rcut)) rcut = ewald_env%rcut
154 IF (PRESENT(epsilon)) epsilon = ewald_env%epsilon
155 IF (PRESENT(eps_pol)) eps_pol = ewald_env%eps_pol
156 IF (PRESENT(gmax)) gmax = ewald_env%gmax
157 IF (PRESENT(ns_max)) ns_max = ewald_env%ns_max
158 IF (PRESENT(o_spline)) o_spline = ewald_env%o_spline
159 IF (PRESENT(group)) group = ewald_env%para_env
160 IF (PRESENT(para_env)) para_env => ewald_env%para_env
161 IF (PRESENT(poisson_section)) poisson_section => ewald_env%poisson_section
162 IF (PRESENT(interaction_cutoffs)) interaction_cutoffs => &
163 ewald_env%interaction_cutoffs
164 IF (PRESENT(cell_hmat)) cell_hmat = ewald_env%cell_hmat
165 END SUBROUTINE ewald_env_get
166
167! **************************************************************************************************
168!> \brief Purpose: Set the EWALD environment.
169!> \param ewald_env the pointer to the ewald_env
170!> \param ewald_type ...
171!> \param alpha ...
172!> \param epsilon ...
173!> \param eps_pol ...
174!> \param gmax ...
175!> \param ns_max ...
176!> \param precs ...
177!> \param o_spline ...
178!> \param para_env ...
179!> \param poisson_section ...
180!> \param interaction_cutoffs ...
181!> \param cell_hmat ...
182!> \par History
183!> 11/03
184!> \author CJM
185! **************************************************************************************************
186 SUBROUTINE ewald_env_set(ewald_env, ewald_type, alpha, epsilon, eps_pol, &
187 gmax, ns_max, precs, o_spline, para_env, poisson_section, &
188 interaction_cutoffs, cell_hmat)
189
190 TYPE(ewald_environment_type), INTENT(INOUT) :: ewald_env
191 INTEGER, OPTIONAL :: ewald_type
192 REAL(kind=dp), OPTIONAL :: alpha, epsilon, eps_pol
193 INTEGER, OPTIONAL :: gmax(3), ns_max
194 REAL(kind=dp), OPTIONAL :: precs
195 INTEGER, OPTIONAL :: o_spline
196 TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
197 TYPE(section_vals_type), OPTIONAL, POINTER :: poisson_section
198 REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL, &
199 POINTER :: interaction_cutoffs
200 REAL(kind=dp), DIMENSION(3, 3), OPTIONAL :: cell_hmat
201
202 IF (PRESENT(ewald_type)) ewald_env%ewald_type = ewald_type
203 IF (PRESENT(alpha)) ewald_env%alpha = alpha
204 IF (PRESENT(precs)) ewald_env%precs = precs
205 IF (PRESENT(epsilon)) ewald_env%epsilon = epsilon
206 IF (PRESENT(eps_pol)) ewald_env%eps_pol = eps_pol
207 IF (PRESENT(gmax)) ewald_env%gmax = gmax
208 IF (PRESENT(ns_max)) ewald_env%ns_max = ns_max
209 IF (PRESENT(o_spline)) ewald_env%o_spline = o_spline
210 IF (PRESENT(para_env)) ewald_env%para_env => para_env
211 IF (PRESENT(poisson_section)) THEN
212 CALL section_vals_retain(poisson_section)
213 CALL section_vals_release(ewald_env%poisson_section)
214 ewald_env%poisson_section => poisson_section
215 END IF
216 IF (PRESENT(interaction_cutoffs)) ewald_env%interaction_cutoffs => &
217 interaction_cutoffs
218 IF (PRESENT(cell_hmat)) ewald_env%cell_hmat = cell_hmat
219 END SUBROUTINE ewald_env_set
220
221! **************************************************************************************************
222!> \brief allocates and intitializes a ewald_env
223!> \param ewald_env the object to create
224!> \param para_env ...
225!> \par History
226!> 12.2002 created [fawzi]
227!> \author Fawzi Mohamed
228! **************************************************************************************************
229 SUBROUTINE ewald_env_create(ewald_env, para_env)
230 TYPE(ewald_environment_type), INTENT(OUT) :: ewald_env
231 TYPE(mp_para_env_type), POINTER :: para_env
232
233 NULLIFY (ewald_env%poisson_section)
234 CALL para_env%retain()
235 ewald_env%para_env => para_env
236 NULLIFY (ewald_env%interaction_cutoffs) ! allocated and initialized later
237 END SUBROUTINE ewald_env_create
238
239! **************************************************************************************************
240!> \brief releases the given ewald_env (see doc/ReferenceCounting.html)
241!> \param ewald_env the object to release
242!> \par History
243!> 12.2002 created [fawzi]
244!> \author Fawzi Mohamed
245! **************************************************************************************************
246 SUBROUTINE ewald_env_release(ewald_env)
247 TYPE(ewald_environment_type), INTENT(INOUT) :: ewald_env
248
249 CALL mp_para_env_release(ewald_env%para_env)
250 CALL section_vals_release(ewald_env%poisson_section)
251 IF (ASSOCIATED(ewald_env%interaction_cutoffs)) THEN
252 DEALLOCATE (ewald_env%interaction_cutoffs)
253 END IF
254
255 END SUBROUTINE ewald_env_release
256
257! **************************************************************************************************
258!> \brief Purpose: read the EWALD section
259!> \param ewald_env the pointer to the ewald_env
260!> \param ewald_section ...
261!> \author Teodoro Laino [tlaino] -University of Zurich - 2005
262! **************************************************************************************************
263 SUBROUTINE read_ewald_section(ewald_env, ewald_section)
264 TYPE(ewald_environment_type), INTENT(INOUT) :: ewald_env
265 TYPE(section_vals_type), POINTER :: ewald_section
266
267 INTEGER :: iw
268 INTEGER, DIMENSION(:), POINTER :: gmax_read
269 LOGICAL :: explicit
270 REAL(kind=dp) :: dummy
271 TYPE(cp_logger_type), POINTER :: logger
272 TYPE(enumeration_type), POINTER :: enum
273 TYPE(keyword_type), POINTER :: keyword
274 TYPE(section_type), POINTER :: section
275 TYPE(section_vals_type), POINTER :: multipole_section
276
277 NULLIFY (enum, keyword, section, multipole_section)
278 logger => cp_get_default_logger()
279 CALL section_vals_val_get(ewald_section, "EWALD_TYPE", i_val=ewald_env%ewald_type)
280 CALL section_vals_val_get(ewald_section, "ALPHA", r_val=ewald_env%alpha)
281 CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", r_val=ewald_env%precs)
282
283 IF (ewald_env%ewald_type == do_ewald_none) THEN
284 ewald_env%rcut = 0.0_dp
285 ELSE
286 CALL section_vals_val_get(ewald_section, "RCUT", explicit=explicit)
287 IF (explicit) THEN
288 CALL section_vals_val_get(ewald_section, "RCUT", r_val=ewald_env%rcut)
289 ELSE
290 ewald_env%rcut = find_ewald_optimal_value(ewald_env%precs)/ewald_env%alpha
291 END IF
292 END IF
293 ! we have no defaults for gmax, gmax is only needed for ewald and spme
294 SELECT CASE (ewald_env%ewald_type)
295 CASE (do_ewald_ewald, do_ewald_spme)
296 CALL section_vals_val_get(ewald_section, "GMAX", i_vals=gmax_read)
297 SELECT CASE (SIZE(gmax_read, 1))
298 CASE (1)
299 ewald_env%gmax = gmax_read(1)
300 CASE (3)
301 ewald_env%gmax = gmax_read
302 CASE DEFAULT
303 cpabort("")
304 END SELECT
305 IF (ewald_env%ewald_type == do_ewald_spme) THEN
306 CALL section_vals_val_get(ewald_section, "O_SPLINE", i_val=ewald_env%o_spline)
307 END IF
308 CASE (do_ewald_pme)
309 CALL section_vals_val_get(ewald_section, "NS_MAX", i_val=ewald_env%ns_max)
310 CALL section_vals_val_get(ewald_section, "EPSILON", r_val=ewald_env%epsilon)
311 CASE DEFAULT
312 ! this should not be used for do_ewald_none
313 ewald_env%gmax = huge(0)
314 ewald_env%ns_max = huge(0)
315 END SELECT
316
317 ! Multipoles
318 multipole_section => section_vals_get_subs_vals(ewald_section, "MULTIPOLES")
319 CALL section_vals_val_get(multipole_section, "_SECTION_PARAMETERS_", l_val=ewald_env%do_multipoles)
320 CALL section_vals_val_get(multipole_section, "POL_SCF", i_val=ewald_env%do_ipol)
321 CALL section_vals_val_get(multipole_section, "EPS_POL", r_val=ewald_env%eps_pol)
322 IF (ewald_env%do_multipoles) THEN
323 SELECT CASE (ewald_env%ewald_type)
324 CASE (do_ewald_ewald)
325 CALL section_vals_val_get(multipole_section, "MAX_MULTIPOLE_EXPANSION", &
326 i_val=ewald_env%max_multipole)
327 CALL section_vals_val_get(multipole_section, "MAX_IPOL_ITER", i_val=ewald_env%max_ipol_iter)
328 CASE DEFAULT
329 cpabort("Multipole code works at the moment only with standard EWALD sums.")
330 END SELECT
331 END IF
332
333 iw = cp_print_key_unit_nr(logger, ewald_section, "PRINT%PROGRAM_RUN_INFO", &
334 extension=".log")
335 IF (iw > 0) THEN
336 NULLIFY (keyword, enum)
337 CALL create_ewald_section(section)
338 IF (ewald_env%ewald_type /= do_ewald_none) THEN
339 keyword => section_get_keyword(section, "EWALD_TYPE")
340 CALL keyword_get(keyword, enum=enum)
341 WRITE (iw, '(/,T2,"EWALD| ",A,T67,A14 )') 'Summation is done by:', &
342 adjustr(trim(enum_i2c(enum, ewald_env%ewald_type)))
343 IF (ewald_env%do_multipoles) THEN
344 NULLIFY (keyword, enum)
345 keyword => section_get_keyword(section, "MULTIPOLES%MAX_MULTIPOLE_EXPANSION")
346 CALL keyword_get(keyword, enum=enum)
347 WRITE (iw, '( T2,"EWALD| ",A )') 'Enabled Multipole Method'
348 WRITE (iw, '( T2,"EWALD| ",A,T67,A14 )') 'Max Term in Multipole Expansion :', &
349 adjustr(trim(enum_i2c(enum, ewald_env%max_multipole)))
350 WRITE (iw, '( T2,"EWALD| ",A,T67,3I10 )') 'Max number Iterations for IPOL :', &
351 ewald_env%max_ipol_iter
352 END IF
353 dummy = cp_unit_from_cp2k(ewald_env%alpha, "angstrom^-1")
354 WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') &
355 'Alpha parameter [', 'ANGSTROM^-1', ']', dummy
356 dummy = cp_unit_from_cp2k(ewald_env%rcut, "angstrom")
357 WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') &
358 'Real Space Cutoff [', 'ANGSTROM', ']', dummy
359
360 SELECT CASE (ewald_env%ewald_type)
361 CASE (do_ewald_ewald)
362 WRITE (iw, '( T2,"EWALD| ",A,T51,3I10 )') &
363 'G-space max. Miller index', ewald_env%gmax
364 CASE (do_ewald_pme)
365 WRITE (iw, '( T2,"EWALD| ",A,T71,I10 )') &
366 'Max small-grid points (input) ', ewald_env%ns_max
367 WRITE (iw, '( T2,"EWALD| ",A,T71,E10.4 )') &
368 'Gaussian tolerance (input) ', ewald_env%epsilon
369 CASE (do_ewald_spme)
370 WRITE (iw, '( T2,"EWALD| ",A,T51,3I10 )') &
371 'G-space max. Miller index', ewald_env%gmax
372 WRITE (iw, '( T2,"EWALD| ",A,T71,I10 )') &
373 'Spline interpolation order ', ewald_env%o_spline
374 CASE DEFAULT
375 cpabort("")
376 END SELECT
377 ELSE
378 WRITE (iw, '( T2,"EWALD| ",T73, A )') 'not used'
379 END IF
380 CALL section_release(section)
381 END IF
382 CALL cp_print_key_finished_output(iw, logger, ewald_section, &
383 "PRINT%PROGRAM_RUN_INFO")
384
385 END SUBROUTINE read_ewald_section
386
387! **************************************************************************************************
388!> \brief Purpose: read the EWALD section for TB methods
389!> \param ewald_env the pointer to the ewald_env
390!> \param ewald_section ...
391!> \param hmat ...
392!> \param silent ...
393!> \param pset ...
394!> \param cell_periodic ...
395!> \author JGH
396! **************************************************************************************************
397 SUBROUTINE read_ewald_section_tb(ewald_env, ewald_section, hmat, silent, pset, cell_periodic)
398 TYPE(ewald_environment_type), INTENT(INOUT) :: ewald_env
399 TYPE(section_vals_type), POINTER :: ewald_section
400 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN) :: hmat
401 LOGICAL, INTENT(IN), OPTIONAL :: silent
402 CHARACTER(LEN=*), OPTIONAL :: pset
403 INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: cell_periodic
404
405 CHARACTER(LEN=5) :: param
406 INTEGER :: i, iw, n(3), poisson_key
407 INTEGER, DIMENSION(3) :: poisson_periodic
408 INTEGER, DIMENSION(:), POINTER :: gmax_read
409 LOGICAL :: do_print, explicit
410 REAL(kind=dp) :: alat, cutoff, dummy, omega
411 TYPE(cp_logger_type), POINTER :: logger
412 TYPE(section_vals_type), POINTER :: poisson_section
413
414 logger => cp_get_default_logger()
415 do_print = .true.
416 IF (PRESENT(silent)) do_print = .NOT. silent
417 param = "none"
418 IF (PRESENT(pset)) param = pset
419
420 IF (PRESENT(cell_periodic)) THEN
421 poisson_section => ewald_env%poisson_section
422 IF (ASSOCIATED(poisson_section)) THEN
423 CALL section_vals_val_get(poisson_section, "PERIODIC", i_val=poisson_key)
424 CALL tb_decode_periodicity(poisson_key, poisson_periodic)
425 IF (any(cell_periodic /= poisson_periodic)) THEN
426 CALL cp_warn(__location__, &
427 "The selected periodicities in SUBSYS/CELL and DFT/POISSON do not match. "// &
428 "The TB Ewald electrostatics will use DFT/POISSON/PERIODIC="// &
429 trim(tb_periodicity_label(poisson_periodic))// &
430 " while SUBSYS/CELL/PERIODIC="//trim(tb_periodicity_label(cell_periodic))//".")
431 END IF
432 END IF
433 END IF
434
435 ewald_env%do_multipoles = .false.
436 ewald_env%do_ipol = 0
437 ewald_env%eps_pol = 1.e-12_dp
438 ewald_env%max_multipole = 0
439 ewald_env%max_ipol_iter = 0
440 ewald_env%epsilon = 1.e-12_dp
441 ewald_env%ns_max = huge(0)
442
443 CALL section_vals_val_get(ewald_section, "EWALD_TYPE", explicit=explicit)
444 IF (explicit) THEN
445 CALL section_vals_val_get(ewald_section, "EWALD_TYPE", i_val=ewald_env%ewald_type)
446 IF (ewald_env%ewald_type /= do_ewald_spme) THEN
447 cpabort("TB needs EWALD_TYPE SPME")
448 END IF
449 ELSE
450 ewald_env%ewald_type = do_ewald_spme
451 END IF
452
453 CALL section_vals_val_get(ewald_section, "ALPHA", explicit=explicit)
454 IF (explicit) THEN
455 CALL section_vals_val_get(ewald_section, "ALPHA", r_val=ewald_env%alpha)
456 ELSE
457 SELECT CASE (param)
458 CASE DEFAULT
459 CALL section_vals_val_get(ewald_section, "ALPHA", r_val=ewald_env%alpha)
460 CASE ("EEQ")
461 omega = abs(det_3x3(hmat))
462 ewald_env%alpha = sqrt(twopi)/omega**(1./3.)
463 END SELECT
464 END IF
465
466 CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", explicit=explicit)
467 IF (explicit) THEN
468 CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", r_val=ewald_env%precs)
469 ELSE
470 CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", r_val=ewald_env%precs)
471 END IF
472
473 CALL section_vals_val_get(ewald_section, "O_SPLINE", explicit=explicit)
474 IF (explicit) THEN
475 CALL section_vals_val_get(ewald_section, "O_SPLINE", i_val=ewald_env%o_spline)
476 ELSE
477 SELECT CASE (param)
478 CASE DEFAULT
479 ewald_env%o_spline = 6
480 CASE ("EEQ")
481 ewald_env%o_spline = 4
482 END SELECT
483 END IF
484
485 CALL section_vals_val_get(ewald_section, "RCUT", explicit=explicit)
486 IF (explicit) THEN
487 CALL section_vals_val_get(ewald_section, "RCUT", r_val=ewald_env%rcut)
488 ELSE
489 ewald_env%rcut = find_ewald_optimal_value(ewald_env%precs)/ewald_env%alpha
490 END IF
491
492 CALL section_vals_val_get(ewald_section, "GMAX", explicit=explicit)
493 IF (explicit) THEN
494 CALL section_vals_val_get(ewald_section, "GMAX", i_vals=gmax_read)
495 SELECT CASE (SIZE(gmax_read, 1))
496 CASE (1)
497 ewald_env%gmax = gmax_read(1)
498 CASE (3)
499 ewald_env%gmax = gmax_read
500 CASE DEFAULT
501 cpabort("")
502 END SELECT
503 ELSE
504 SELECT CASE (param)
505 CASE DEFAULT
506 ! set GMAX using ECUT=alpha*45 Ry
507 cutoff = 45._dp*ewald_env%alpha
508 CASE ("EEQ")
509 ! set GMAX using ECUT=alpha*45 Ry
510 cutoff = 30._dp*ewald_env%alpha
511 END SELECT
512 DO i = 1, 3
513 alat = sum(hmat(:, i)**2)
514 cpassert(alat /= 0._dp)
515 ewald_env%gmax(i) = 2*floor(sqrt(2.0_dp*cutoff*alat)/twopi) + 1
516 END DO
517 END IF
518 n = ewald_env%gmax
519 ewald_env%gmax = pw_grid_n_for_fft(n, odd=.true.)
520
521 iw = cp_print_key_unit_nr(logger, ewald_section, "PRINT%PROGRAM_RUN_INFO", &
522 extension=".log")
523 IF (iw > 0 .AND. do_print) THEN
524 WRITE (iw, '(/,T2,"EWALD| ",A,T67,A14 )') 'Summation is done by:', adjustr("SPME")
525 dummy = cp_unit_from_cp2k(ewald_env%alpha, "angstrom^-1")
526 WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') &
527 'Alpha parameter [', 'ANGSTROM^-1', ']', dummy
528 dummy = cp_unit_from_cp2k(ewald_env%rcut, "angstrom")
529 WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') &
530 'Real Space Cutoff [', 'ANGSTROM', ']', dummy
531 WRITE (iw, '( T2,"EWALD| ",A,T51,3I10 )') &
532 'G-space max. Miller index', ewald_env%gmax
533 WRITE (iw, '( T2,"EWALD| ",A,T71,I10 )') &
534 'Spline interpolation order ', ewald_env%o_spline
535 END IF
536 CALL cp_print_key_finished_output(iw, logger, ewald_section, &
537 "PRINT%PROGRAM_RUN_INFO")
538
539 END SUBROUTINE read_ewald_section_tb
540
541! **************************************************************************************************
542!> \brief Decode CP2K PERIODIC enum into the 3-component periodicity mask.
543!> \param periodic_key input enum value
544!> \param periodic periodicity mask
545! **************************************************************************************************
546 SUBROUTINE tb_decode_periodicity(periodic_key, periodic)
547 INTEGER, INTENT(IN) :: periodic_key
548 INTEGER, DIMENSION(3), INTENT(OUT) :: periodic
549
550 SELECT CASE (periodic_key)
551 CASE (use_perd_x)
552 periodic = [1, 0, 0]
553 CASE (use_perd_y)
554 periodic = [0, 1, 0]
555 CASE (use_perd_z)
556 periodic = [0, 0, 1]
557 CASE (use_perd_xy)
558 periodic = [1, 1, 0]
559 CASE (use_perd_xz)
560 periodic = [1, 0, 1]
561 CASE (use_perd_yz)
562 periodic = [0, 1, 1]
563 CASE (use_perd_xyz)
564 periodic = [1, 1, 1]
565 CASE (use_perd_none)
566 periodic = [0, 0, 0]
567 CASE DEFAULT
568 cpabort("Invalid PERIODIC setting for TB Ewald")
569 END SELECT
570
571 END SUBROUTINE tb_decode_periodicity
572
573! **************************************************************************************************
574!> \brief Return a compact label for a 3-component periodicity mask.
575!> \param periodic periodicity mask
576!> \return label
577! **************************************************************************************************
578 FUNCTION tb_periodicity_label(periodic) RESULT(label)
579 INTEGER, DIMENSION(3), INTENT(IN) :: periodic
580 CHARACTER(LEN=4) :: label
581
582 label = ""
583 IF (all(periodic == 0)) THEN
584 label = "NONE"
585 ELSE
586 IF (periodic(1) == 1) label = trim(label)//"X"
587 IF (periodic(2) == 1) label = trim(label)//"Y"
588 IF (periodic(3) == 1) label = trim(label)//"Z"
589 END IF
590
591 END FUNCTION tb_periodicity_label
592
593! **************************************************************************************************
594!> \brief triggers (by bisection) the optimal value for EWALD parameter x
595!> EXP(-x^2)/x^2 = EWALD_ACCURACY
596!> \param precs ...
597!> \return ...
598!> \author Teodoro Laino [tlaino] - University of Zurich - 12.2007
599! **************************************************************************************************
600 FUNCTION find_ewald_optimal_value(precs) RESULT(value)
601 REAL(kind=dp) :: precs, value
602
603 REAL(kind=dp) :: func, func1, func2, s, s1, s2
604
605 s = 0.1_dp
606 func = exp(-s**2)/s**2 - precs
607 cpassert(func > 0.0_dp)
608 DO WHILE (func > 0.0_dp)
609 s = s + 0.1_dp
610 func = exp(-s**2)/s**2 - precs
611 END DO
612 s2 = s
613 s1 = s - 0.1_dp
614 ! Start bisection
615 DO WHILE (.true.)
616 func2 = exp(-s2**2)/s2**2 - precs
617 func1 = exp(-s1**2)/s1**2 - precs
618 cpassert(func1 >= 0)
619 cpassert(func2 <= 0)
620 s = 0.5_dp*(s1 + s2)
621 func = exp(-s**2)/s**2 - precs
622 IF (func > 0.0_dp) THEN
623 s1 = s
624 ELSE IF (func < 0.0_dp) THEN
625 s2 = s
626 END IF
627 IF (abs(func) < 100.0_dp*epsilon(0.0_dp)) EXIT
628 END DO
629 value = s
630 END FUNCTION find_ewald_optimal_value
631
632END MODULE ewald_environment_types
mathlib::det_3x3
Definition mathlib.F:69
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::use_perd_xyz
integer, parameter, public use_perd_xyz
Definition cell_types.F:42
cell_types::use_perd_y
integer, parameter, public use_perd_y
Definition cell_types.F:42
cell_types::use_perd_xz
integer, parameter, public use_perd_xz
Definition cell_types.F:42
cell_types::use_perd_x
integer, parameter, public use_perd_x
Definition cell_types.F:42
cell_types::use_perd_z
integer, parameter, public use_perd_z
Definition cell_types.F:42
cell_types::use_perd_yz
integer, parameter, public use_perd_yz
Definition cell_types.F:42
cell_types::use_perd_none
integer, parameter, public use_perd_none
Definition cell_types.F:42
cell_types::use_perd_xy
integer, parameter, public use_perd_xy
Definition cell_types.F:42
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
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)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
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,...
Definition cp_output_handling.F:1067
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1239
ewald_environment_types
Definition ewald_environment_types.F:14
ewald_environment_types::ewald_env_set
subroutine, public ewald_env_set(ewald_env, ewald_type, alpha, epsilon, eps_pol, gmax, ns_max, precs, o_spline, para_env, poisson_section, interaction_cutoffs, cell_hmat)
Purpose: Set the EWALD environment.
Definition ewald_environment_types.F:189
ewald_environment_types::ewald_env_create
subroutine, public ewald_env_create(ewald_env, para_env)
allocates and intitializes a ewald_env
Definition ewald_environment_types.F:230
ewald_environment_types::read_ewald_section
subroutine, public read_ewald_section(ewald_env, ewald_section)
Purpose: read the EWALD section.
Definition ewald_environment_types.F:264
ewald_environment_types::read_ewald_section_tb
subroutine, public read_ewald_section_tb(ewald_env, ewald_section, hmat, silent, pset, cell_periodic)
Purpose: read the EWALD section for TB methods.
Definition ewald_environment_types.F:398
ewald_environment_types::ewald_env_release
subroutine, public ewald_env_release(ewald_env)
releases the given ewald_env (see doc/ReferenceCounting.html)
Definition ewald_environment_types.F:247
ewald_environment_types::ewald_env_get
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
Definition ewald_environment_types.F:132
input_cp2k_poisson
function that build the poisson section of the input
Definition input_cp2k_poisson.F:14
input_cp2k_poisson::create_ewald_section
subroutine, public create_ewald_section(section)
...
Definition input_cp2k_poisson.F:283
input_enumeration_types
represents an enumeration, i.e. a mapping between integers and strings
Definition input_enumeration_types.F:14
input_enumeration_types::enum_i2c
character(len=default_string_length) function, public enum_i2c(enum, i)
maps an integer to a string
Definition input_enumeration_types.F:152
input_keyword_types
represents keywords in an input
Definition input_keyword_types.F:14
input_keyword_types::keyword_get
subroutine, public keyword_get(keyword, names, usage, description, type_of_var, n_var, default_value, lone_keyword_value, repeats, enum, citations)
...
Definition input_keyword_types.F:481
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_retain
subroutine, public section_vals_retain(section_vals)
retains the given section values (see doc/ReferenceCounting.html)
Definition input_section_types.F:638
input_section_types::section_vals_get_subs_vals
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
Definition input_section_types.F:731
input_section_types::section_release
recursive subroutine, public section_release(section)
releases the given keyword list (see doc/ReferenceCounting.html)
Definition input_section_types.F:238
input_section_types::section_get_keyword
recursive type(keyword_type) function, pointer, public section_get_keyword(section, keyword_name)
returns the requested keyword
Definition input_section_types.F:478
input_section_types::section_vals_val_get
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
Definition input_section_types.F:1047
input_section_types::section_vals_release
recursive subroutine, public section_vals_release(section_vals)
releases the given object
Definition input_section_types.F:653
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
message_passing::mp_para_env_release
subroutine, public mp_para_env_release(para_env)
releases the para object (to be called when you don't want anymore the shared copy of this object)
Definition message_passing.F:2011
pw_grid_info
This module returns additional info on PW grids.
Definition pw_grid_info.F:14
pw_grid_info::pw_grid_n_for_fft
integer function, dimension(3), public pw_grid_n_for_fft(n, odd)
returns the closest number of points >= n, on which you can perform ffts
Definition pw_grid_info.F:187
pw_poisson_types
functions related to the poisson solver on regular grids
Definition pw_poisson_types.F:22
pw_poisson_types::do_ewald_pme
integer, parameter, public do_ewald_pme
Definition pw_poisson_types.F:76
pw_poisson_types::do_ewald_ewald
integer, parameter, public do_ewald_ewald
Definition pw_poisson_types.F:76
pw_poisson_types::do_ewald_none
integer, parameter, public do_ewald_none
Definition pw_poisson_types.F:76
pw_poisson_types::do_ewald_spme
integer, parameter, public do_ewald_spme
Definition pw_poisson_types.F:76
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
ewald_environment_types::ewald_environment_type
to build arrays of pointers
Definition ewald_environment_types.F:64
input_enumeration_types::enumeration_type
Definition input_enumeration_types.F:43
input_keyword_types::keyword_type
represent a keyword in the input
Definition input_keyword_types.F:89
input_section_types::section_type
represent a section of the input file
Definition input_section_types.F:102
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_comm_type
Definition message_passing.F:155
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743