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 : !> \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 : ! **************************************************************************************************
14 : MODULE 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 :
101 : CONTAINS
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 557347 : 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 557347 : IF (PRESENT(ewald_type)) ewald_type = ewald_env%ewald_type
147 557347 : IF (PRESENT(do_multipoles)) do_multipoles = ewald_env%do_multipoles
148 557347 : IF (PRESENT(do_ipol)) do_ipol = ewald_env%do_ipol
149 557347 : IF (PRESENT(max_multipole)) max_multipole = ewald_env%max_multipole
150 557347 : IF (PRESENT(max_ipol_iter)) max_ipol_iter = ewald_env%max_ipol_iter
151 557347 : IF (PRESENT(alpha)) alpha = ewald_env%alpha
152 557347 : IF (PRESENT(precs)) precs = ewald_env%precs
153 557347 : IF (PRESENT(rcut)) rcut = ewald_env%rcut
154 557347 : IF (PRESENT(epsilon)) epsilon = ewald_env%epsilon
155 557347 : IF (PRESENT(eps_pol)) eps_pol = ewald_env%eps_pol
156 574663 : IF (PRESENT(gmax)) gmax = ewald_env%gmax
157 557347 : IF (PRESENT(ns_max)) ns_max = ewald_env%ns_max
158 557347 : IF (PRESENT(o_spline)) o_spline = ewald_env%o_spline
159 557347 : IF (PRESENT(group)) group = ewald_env%para_env
160 557347 : IF (PRESENT(para_env)) para_env => ewald_env%para_env
161 557347 : IF (PRESENT(poisson_section)) poisson_section => ewald_env%poisson_section
162 557347 : IF (PRESENT(interaction_cutoffs)) interaction_cutoffs => &
163 81608 : ewald_env%interaction_cutoffs
164 1673410 : IF (PRESENT(cell_hmat)) cell_hmat = ewald_env%cell_hmat
165 557347 : 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 25438 : 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 25438 : IF (PRESENT(ewald_type)) ewald_env%ewald_type = ewald_type
203 25438 : IF (PRESENT(alpha)) ewald_env%alpha = alpha
204 25438 : IF (PRESENT(precs)) ewald_env%precs = precs
205 25438 : IF (PRESENT(epsilon)) ewald_env%epsilon = epsilon
206 25438 : IF (PRESENT(eps_pol)) ewald_env%eps_pol = eps_pol
207 25438 : IF (PRESENT(gmax)) ewald_env%gmax = gmax
208 25438 : IF (PRESENT(ns_max)) ewald_env%ns_max = ns_max
209 25438 : IF (PRESENT(o_spline)) ewald_env%o_spline = o_spline
210 25438 : IF (PRESENT(para_env)) ewald_env%para_env => para_env
211 25438 : IF (PRESENT(poisson_section)) THEN
212 4321 : CALL section_vals_retain(poisson_section)
213 4321 : CALL section_vals_release(ewald_env%poisson_section)
214 4321 : ewald_env%poisson_section => poisson_section
215 : END IF
216 25438 : IF (PRESENT(interaction_cutoffs)) ewald_env%interaction_cutoffs => &
217 2647 : interaction_cutoffs
218 265548 : IF (PRESENT(cell_hmat)) ewald_env%cell_hmat = cell_hmat
219 25438 : 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 73457 : 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 4321 : CALL para_env%retain()
235 4321 : ewald_env%para_env => para_env
236 4321 : NULLIFY (ewald_env%interaction_cutoffs) ! allocated and initialized later
237 4321 : 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 4321 : SUBROUTINE ewald_env_release(ewald_env)
247 : TYPE(ewald_environment_type), INTENT(INOUT) :: ewald_env
248 :
249 4321 : CALL mp_para_env_release(ewald_env%para_env)
250 4321 : CALL section_vals_release(ewald_env%poisson_section)
251 4321 : IF (ASSOCIATED(ewald_env%interaction_cutoffs)) THEN
252 2647 : DEALLOCATE (ewald_env%interaction_cutoffs)
253 : END IF
254 :
255 4321 : 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 3839 : 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 3839 : 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 3839 : NULLIFY (enum, keyword, section, multipole_section)
278 7678 : logger => cp_get_default_logger()
279 3839 : CALL section_vals_val_get(ewald_section, "EWALD_TYPE", i_val=ewald_env%ewald_type)
280 3839 : CALL section_vals_val_get(ewald_section, "ALPHA", r_val=ewald_env%alpha)
281 3839 : CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", r_val=ewald_env%precs)
282 :
283 3839 : IF (ewald_env%ewald_type == do_ewald_none) THEN
284 1046 : ewald_env%rcut = 0.0_dp
285 : ELSE
286 2793 : CALL section_vals_val_get(ewald_section, "RCUT", explicit=explicit)
287 2793 : IF (explicit) THEN
288 8 : CALL section_vals_val_get(ewald_section, "RCUT", r_val=ewald_env%rcut)
289 : ELSE
290 2785 : 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 6546 : SELECT CASE (ewald_env%ewald_type)
295 : CASE (do_ewald_ewald, do_ewald_spme)
296 2707 : CALL section_vals_val_get(ewald_section, "GMAX", i_vals=gmax_read)
297 5098 : SELECT CASE (SIZE(gmax_read, 1))
298 : CASE (1)
299 9564 : ewald_env%gmax = gmax_read(1)
300 : CASE (3)
301 1264 : ewald_env%gmax = gmax_read
302 : CASE DEFAULT
303 2707 : CPABORT("")
304 : END SELECT
305 2707 : IF (ewald_env%ewald_type == do_ewald_spme) THEN
306 1884 : CALL section_vals_val_get(ewald_section, "O_SPLINE", i_val=ewald_env%o_spline)
307 : END IF
308 : CASE (do_ewald_pme)
309 86 : CALL section_vals_val_get(ewald_section, "NS_MAX", i_val=ewald_env%ns_max)
310 86 : 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 4184 : ewald_env%gmax = HUGE(0)
314 4885 : ewald_env%ns_max = HUGE(0)
315 : END SELECT
316 :
317 : ! Multipoles
318 3839 : multipole_section => section_vals_get_subs_vals(ewald_section, "MULTIPOLES")
319 3839 : CALL section_vals_val_get(multipole_section, "_SECTION_PARAMETERS_", l_val=ewald_env%do_multipoles)
320 3839 : CALL section_vals_val_get(multipole_section, "POL_SCF", i_val=ewald_env%do_ipol)
321 3839 : CALL section_vals_val_get(multipole_section, "EPS_POL", r_val=ewald_env%eps_pol)
322 3839 : IF (ewald_env%do_multipoles) THEN
323 336 : SELECT CASE (ewald_env%ewald_type)
324 : CASE (do_ewald_ewald)
325 : CALL section_vals_val_get(multipole_section, "MAX_MULTIPOLE_EXPANSION", &
326 168 : i_val=ewald_env%max_multipole)
327 168 : CALL section_vals_val_get(multipole_section, "MAX_IPOL_ITER", i_val=ewald_env%max_ipol_iter)
328 : CASE DEFAULT
329 168 : 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 3839 : extension=".log")
335 3839 : IF (iw > 0) THEN
336 1908 : NULLIFY (keyword, enum)
337 1908 : CALL create_ewald_section(section)
338 1908 : IF (ewald_env%ewald_type /= do_ewald_none) THEN
339 1423 : keyword => section_get_keyword(section, "EWALD_TYPE")
340 1423 : CALL keyword_get(keyword, enum=enum)
341 1423 : WRITE (iw, '(/,T2,"EWALD| ",A,T67,A14 )') 'Summation is done by:', &
342 2846 : ADJUSTR(TRIM(enum_i2c(enum, ewald_env%ewald_type)))
343 1423 : IF (ewald_env%do_multipoles) THEN
344 84 : NULLIFY (keyword, enum)
345 84 : keyword => section_get_keyword(section, "MULTIPOLES%MAX_MULTIPOLE_EXPANSION")
346 84 : CALL keyword_get(keyword, enum=enum)
347 84 : WRITE (iw, '( T2,"EWALD| ",A )') 'Enabled Multipole Method'
348 84 : WRITE (iw, '( T2,"EWALD| ",A,T67,A14 )') 'Max Term in Multipole Expansion :', &
349 168 : ADJUSTR(TRIM(enum_i2c(enum, ewald_env%max_multipole)))
350 84 : WRITE (iw, '( T2,"EWALD| ",A,T67,3I10 )') 'Max number Iterations for IPOL :', &
351 168 : ewald_env%max_ipol_iter
352 : END IF
353 1423 : dummy = cp_unit_from_cp2k(ewald_env%alpha, "angstrom^-1")
354 : WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') &
355 1423 : 'Alpha parameter [', 'ANGSTROM^-1', ']', dummy
356 1423 : dummy = cp_unit_from_cp2k(ewald_env%rcut, "angstrom")
357 : WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') &
358 1423 : 'Real Space Cutoff [', 'ANGSTROM', ']', dummy
359 :
360 1874 : SELECT CASE (ewald_env%ewald_type)
361 : CASE (do_ewald_ewald)
362 : WRITE (iw, '( T2,"EWALD| ",A,T51,3I10 )') &
363 451 : 'G-space max. Miller index', ewald_env%gmax
364 : CASE (do_ewald_pme)
365 : WRITE (iw, '( T2,"EWALD| ",A,T71,I10 )') &
366 43 : 'Max small-grid points (input) ', ewald_env%ns_max
367 : WRITE (iw, '( T2,"EWALD| ",A,T71,E10.4 )') &
368 43 : 'Gaussian tolerance (input) ', ewald_env%epsilon
369 : CASE (do_ewald_spme)
370 : WRITE (iw, '( T2,"EWALD| ",A,T51,3I10 )') &
371 929 : 'G-space max. Miller index', ewald_env%gmax
372 : WRITE (iw, '( T2,"EWALD| ",A,T71,I10 )') &
373 929 : 'Spline interpolation order ', ewald_env%o_spline
374 : CASE DEFAULT
375 1423 : CPABORT("")
376 : END SELECT
377 : ELSE
378 485 : WRITE (iw, '( T2,"EWALD| ",T73, A )') 'not used'
379 : END IF
380 1908 : CALL section_release(section)
381 : END IF
382 : CALL cp_print_key_finished_output(iw, logger, ewald_section, &
383 3839 : "PRINT%PROGRAM_RUN_INFO")
384 :
385 3839 : 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 3374 : 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 482 : 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 964 : logger => cp_get_default_logger()
415 482 : do_print = .TRUE.
416 482 : IF (PRESENT(silent)) do_print = .NOT. silent
417 482 : param = "none"
418 482 : IF (PRESENT(pset)) param = pset
419 :
420 482 : IF (PRESENT(cell_periodic)) THEN
421 482 : poisson_section => ewald_env%poisson_section
422 482 : IF (ASSOCIATED(poisson_section)) THEN
423 482 : CALL section_vals_val_get(poisson_section, "PERIODIC", i_val=poisson_key)
424 482 : CALL tb_decode_periodicity(poisson_key, poisson_periodic)
425 1910 : 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 6 : " while SUBSYS/CELL/PERIODIC="//TRIM(tb_periodicity_label(cell_periodic))//".")
431 : END IF
432 : END IF
433 : END IF
434 :
435 482 : ewald_env%do_multipoles = .FALSE.
436 482 : ewald_env%do_ipol = 0
437 482 : ewald_env%eps_pol = 1.e-12_dp
438 482 : ewald_env%max_multipole = 0
439 482 : ewald_env%max_ipol_iter = 0
440 482 : ewald_env%epsilon = 1.e-12_dp
441 482 : ewald_env%ns_max = HUGE(0)
442 :
443 482 : CALL section_vals_val_get(ewald_section, "EWALD_TYPE", explicit=explicit)
444 482 : IF (explicit) THEN
445 196 : CALL section_vals_val_get(ewald_section, "EWALD_TYPE", i_val=ewald_env%ewald_type)
446 196 : IF (ewald_env%ewald_type /= do_ewald_spme) THEN
447 0 : CPABORT("TB needs EWALD_TYPE SPME")
448 : END IF
449 : ELSE
450 286 : ewald_env%ewald_type = do_ewald_spme
451 : END IF
452 :
453 482 : CALL section_vals_val_get(ewald_section, "ALPHA", explicit=explicit)
454 482 : IF (explicit) THEN
455 202 : CALL section_vals_val_get(ewald_section, "ALPHA", r_val=ewald_env%alpha)
456 : ELSE
457 194 : SELECT CASE (param)
458 : CASE DEFAULT
459 194 : CALL section_vals_val_get(ewald_section, "ALPHA", r_val=ewald_env%alpha)
460 : CASE ("EEQ")
461 86 : omega = ABS(det_3x3(hmat))
462 280 : ewald_env%alpha = SQRT(twopi)/omega**(1./3.)
463 : END SELECT
464 : END IF
465 :
466 482 : CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", explicit=explicit)
467 482 : IF (explicit) THEN
468 0 : CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", r_val=ewald_env%precs)
469 : ELSE
470 482 : CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", r_val=ewald_env%precs)
471 : END IF
472 :
473 482 : CALL section_vals_val_get(ewald_section, "O_SPLINE", explicit=explicit)
474 482 : IF (explicit) THEN
475 108 : CALL section_vals_val_get(ewald_section, "O_SPLINE", i_val=ewald_env%o_spline)
476 : ELSE
477 278 : SELECT CASE (param)
478 : CASE DEFAULT
479 278 : ewald_env%o_spline = 6
480 : CASE ("EEQ")
481 374 : ewald_env%o_spline = 4
482 : END SELECT
483 : END IF
484 :
485 482 : CALL section_vals_val_get(ewald_section, "RCUT", explicit=explicit)
486 482 : IF (explicit) THEN
487 0 : CALL section_vals_val_get(ewald_section, "RCUT", r_val=ewald_env%rcut)
488 : ELSE
489 482 : ewald_env%rcut = find_ewald_optimal_value(ewald_env%precs)/ewald_env%alpha
490 : END IF
491 :
492 482 : CALL section_vals_val_get(ewald_section, "GMAX", explicit=explicit)
493 482 : IF (explicit) THEN
494 186 : CALL section_vals_val_get(ewald_section, "GMAX", i_vals=gmax_read)
495 314 : SELECT CASE (SIZE(gmax_read, 1))
496 : CASE (1)
497 512 : ewald_env%gmax = gmax_read(1)
498 : CASE (3)
499 232 : ewald_env%gmax = gmax_read
500 : CASE DEFAULT
501 186 : CPABORT("")
502 : END SELECT
503 : ELSE
504 210 : SELECT CASE (param)
505 : CASE DEFAULT
506 : ! set GMAX using ECUT=alpha*45 Ry
507 210 : cutoff = 45._dp*ewald_env%alpha
508 : CASE ("EEQ")
509 : ! set GMAX using ECUT=alpha*45 Ry
510 296 : cutoff = 30._dp*ewald_env%alpha
511 : END SELECT
512 1184 : DO i = 1, 3
513 3552 : alat = SUM(hmat(:, i)**2)
514 888 : CPASSERT(alat /= 0._dp)
515 1184 : ewald_env%gmax(i) = 2*FLOOR(SQRT(2.0_dp*cutoff*alat)/twopi) + 1
516 : END DO
517 : END IF
518 1928 : n = ewald_env%gmax
519 1928 : 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 482 : extension=".log")
523 482 : IF (iw > 0 .AND. do_print) THEN
524 213 : WRITE (iw, '(/,T2,"EWALD| ",A,T67,A14 )') 'Summation is done by:', ADJUSTR("SPME")
525 213 : dummy = cp_unit_from_cp2k(ewald_env%alpha, "angstrom^-1")
526 : WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') &
527 213 : 'Alpha parameter [', 'ANGSTROM^-1', ']', dummy
528 213 : dummy = cp_unit_from_cp2k(ewald_env%rcut, "angstrom")
529 : WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') &
530 213 : 'Real Space Cutoff [', 'ANGSTROM', ']', dummy
531 : WRITE (iw, '( T2,"EWALD| ",A,T51,3I10 )') &
532 213 : 'G-space max. Miller index', ewald_env%gmax
533 : WRITE (iw, '( T2,"EWALD| ",A,T71,I10 )') &
534 213 : 'Spline interpolation order ', ewald_env%o_spline
535 : END IF
536 : CALL cp_print_key_finished_output(iw, logger, ewald_section, &
537 482 : "PRINT%PROGRAM_RUN_INFO")
538 :
539 482 : 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 482 : SUBROUTINE tb_decode_periodicity(periodic_key, periodic)
547 : INTEGER, INTENT(IN) :: periodic_key
548 : INTEGER, DIMENSION(3), INTENT(OUT) :: periodic
549 :
550 482 : SELECT CASE (periodic_key)
551 : CASE (use_perd_x)
552 0 : periodic = [1, 0, 0]
553 : CASE (use_perd_y)
554 0 : periodic = [0, 1, 0]
555 : CASE (use_perd_z)
556 0 : periodic = [0, 0, 1]
557 : CASE (use_perd_xy)
558 0 : periodic = [1, 1, 0]
559 : CASE (use_perd_xz)
560 0 : periodic = [1, 0, 1]
561 : CASE (use_perd_yz)
562 0 : periodic = [0, 1, 1]
563 : CASE (use_perd_xyz)
564 476 : periodic = [1, 1, 1]
565 : CASE (use_perd_none)
566 6 : periodic = [0, 0, 0]
567 : CASE DEFAULT
568 482 : CPABORT("Invalid PERIODIC setting for TB Ewald")
569 : END SELECT
570 :
571 482 : 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 12 : FUNCTION tb_periodicity_label(periodic) RESULT(label)
579 : INTEGER, DIMENSION(3), INTENT(IN) :: periodic
580 : CHARACTER(LEN=4) :: label
581 :
582 12 : label = ""
583 30 : IF (ALL(periodic == 0)) THEN
584 6 : label = "NONE"
585 : ELSE
586 6 : IF (periodic(1) == 1) label = TRIM(label)//"X"
587 6 : IF (periodic(2) == 1) label = TRIM(label)//"Y"
588 6 : IF (periodic(3) == 1) label = TRIM(label)//"Z"
589 : END IF
590 :
591 12 : 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 3267 : 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 3267 : s = 0.1_dp
606 3267 : func = EXP(-s**2)/s**2 - precs
607 3267 : CPASSERT(func > 0.0_dp)
608 107323 : DO WHILE (func > 0.0_dp)
609 104056 : s = s + 0.1_dp
610 107323 : func = EXP(-s**2)/s**2 - precs
611 : END DO
612 3267 : s2 = s
613 3267 : s1 = s - 0.1_dp
614 : ! Start bisection
615 : DO WHILE (.TRUE.)
616 81649 : func2 = EXP(-s2**2)/s2**2 - precs
617 81649 : func1 = EXP(-s1**2)/s1**2 - precs
618 81649 : CPASSERT(func1 >= 0)
619 81649 : CPASSERT(func2 <= 0)
620 81649 : s = 0.5_dp*(s1 + s2)
621 81649 : func = EXP(-s**2)/s**2 - precs
622 81649 : IF (func > 0.0_dp) THEN
623 : s1 = s
624 29598 : ELSE IF (func < 0.0_dp) THEN
625 29598 : s2 = s
626 : END IF
627 81649 : IF (ABS(func) < 100.0_dp*EPSILON(0.0_dp)) EXIT
628 : END DO
629 3267 : value = s
630 3267 : END FUNCTION find_ewald_optimal_value
631 :
632 0 : END MODULE ewald_environment_types
|
