Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 cp_log_handling, ONLY: cp_get_default_logger,&
16 : cp_logger_type
17 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
18 : cp_print_key_unit_nr
19 : USE cp_units, ONLY: cp_unit_from_cp2k
20 : USE input_cp2k_poisson, ONLY: create_ewald_section
21 : USE input_enumeration_types, ONLY: enum_i2c,&
22 : enumeration_type
23 : USE input_keyword_types, ONLY: keyword_get,&
24 : keyword_type
25 : USE input_section_types, ONLY: section_get_keyword,&
26 : section_release,&
27 : section_type,&
28 : section_vals_get_subs_vals,&
29 : section_vals_release,&
30 : section_vals_retain,&
31 : section_vals_type,&
32 : section_vals_val_get
33 : USE kinds, ONLY: dp
34 : USE mathconstants, ONLY: twopi
35 : USE mathlib, ONLY: det_3x3
36 : USE message_passing, ONLY: mp_comm_type,&
37 : mp_para_env_release,&
38 : mp_para_env_type
39 : USE pw_grid_info, ONLY: pw_grid_n_for_fft
40 : USE pw_poisson_types, ONLY: do_ewald_ewald,&
41 : do_ewald_none,&
42 : do_ewald_pme,&
43 : do_ewald_spme
44 : #include "./base/base_uses.f90"
45 :
46 : IMPLICIT NONE
47 : PRIVATE
48 :
49 : ! **************************************************************************************************
50 : !> \brief to build arrays of pointers
51 : !> \param ewald_env the pointer to the ewald_env
52 : !> \par History
53 : !> 11/03
54 : !> \author CJM
55 : ! **************************************************************************************************
56 : TYPE ewald_environment_type
57 : PRIVATE
58 : LOGICAL :: do_multipoles = .FALSE. ! Flag for using the multipole code
59 : INTEGER :: do_ipol = -1 ! Solver for induced dipoles
60 : INTEGER :: max_multipole = -1 ! max expansion in the multipoles
61 : INTEGER :: max_ipol_iter = -1 ! max number of interaction for induced dipoles
62 : INTEGER :: ewald_type = -1 ! type of ewald
63 : INTEGER :: gmax(3) = -1 ! max Miller index
64 : INTEGER :: ns_max = -1 ! # grid points for small grid (PME)
65 : INTEGER :: o_spline = -1 ! order of spline (SPME)
66 : REAL(KIND=dp) :: precs = 0.0_dp ! precision achieved when evaluating the real-space part
67 : REAL(KIND=dp) :: alpha = 0.0_dp, rcut = 0.0_dp ! ewald alpha and real-space cutoff
68 : REAL(KIND=dp) :: epsilon = 0.0_dp ! tolerance for small grid (PME)
69 : REAL(KIND=dp) :: eps_pol = 0.0_dp ! tolerance for convergence of induced dipoles
70 : TYPE(mp_para_env_type), POINTER :: para_env => NULL()
71 : TYPE(section_vals_type), POINTER :: poisson_section => NULL()
72 : ! interaction cutoff is required to make the electrostatic interaction
73 : ! continuous at a pair distance equal to rcut. this is ignored by the
74 : ! multipole code and is otherwise only active when SHIFT_CUTOFF is used.
75 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: interaction_cutoffs => NULL()
76 : ! store current cell, used to rebuild lazily.
77 : REAL(KIND=dp), DIMENSION(3, 3) :: cell_hmat = -1.0_dp
78 : END TYPE ewald_environment_type
79 :
80 : ! *** Public data types ***
81 : PUBLIC :: ewald_environment_type
82 :
83 : ! *** Public subroutines ***
84 : PUBLIC :: ewald_env_get, &
85 : ewald_env_set, &
86 : ewald_env_create, &
87 : ewald_env_release, &
88 : read_ewald_section, &
89 : read_ewald_section_tb
90 :
91 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ewald_environment_types'
92 :
93 : CONTAINS
94 :
95 : ! **************************************************************************************************
96 : !> \brief Purpose: Get the EWALD environment.
97 : !> \param ewald_env the pointer to the ewald_env
98 : !> \param ewald_type ...
99 : !> \param alpha ...
100 : !> \param eps_pol ...
101 : !> \param epsilon ...
102 : !> \param gmax ...
103 : !> \param ns_max ...
104 : !> \param o_spline ...
105 : !> \param group ...
106 : !> \param para_env ...
107 : !> \param poisson_section ...
108 : !> \param precs ...
109 : !> \param rcut ...
110 : !> \param do_multipoles ...
111 : !> \param max_multipole ...
112 : !> \param do_ipol ...
113 : !> \param max_ipol_iter ...
114 : !> \param interaction_cutoffs ...
115 : !> \param cell_hmat ...
116 : !> \par History
117 : !> 11/03
118 : !> \author CJM
119 : ! **************************************************************************************************
120 535506 : SUBROUTINE ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, &
121 : gmax, ns_max, o_spline, group, para_env, poisson_section, precs, &
122 : rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, &
123 : interaction_cutoffs, cell_hmat)
124 : TYPE(ewald_environment_type), INTENT(IN) :: ewald_env
125 : INTEGER, OPTIONAL :: ewald_type
126 : REAL(KIND=dp), OPTIONAL :: alpha, eps_pol, epsilon
127 : INTEGER, OPTIONAL :: gmax(3), ns_max, o_spline
128 : TYPE(mp_comm_type), INTENT(OUT), OPTIONAL :: group
129 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
130 : TYPE(section_vals_type), OPTIONAL, POINTER :: poisson_section
131 : REAL(KIND=dp), OPTIONAL :: precs, rcut
132 : LOGICAL, INTENT(OUT), OPTIONAL :: do_multipoles
133 : INTEGER, INTENT(OUT), OPTIONAL :: max_multipole, do_ipol, max_ipol_iter
134 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
135 : POINTER :: interaction_cutoffs
136 : REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL :: cell_hmat
137 :
138 535506 : IF (PRESENT(ewald_type)) ewald_type = ewald_env%ewald_type
139 535506 : IF (PRESENT(do_multipoles)) do_multipoles = ewald_env%do_multipoles
140 535506 : IF (PRESENT(do_ipol)) do_ipol = ewald_env%do_ipol
141 535506 : IF (PRESENT(max_multipole)) max_multipole = ewald_env%max_multipole
142 535506 : IF (PRESENT(max_ipol_iter)) max_ipol_iter = ewald_env%max_ipol_iter
143 535506 : IF (PRESENT(alpha)) alpha = ewald_env%alpha
144 535506 : IF (PRESENT(precs)) precs = ewald_env%precs
145 535506 : IF (PRESENT(rcut)) rcut = ewald_env%rcut
146 535506 : IF (PRESENT(epsilon)) epsilon = ewald_env%epsilon
147 535506 : IF (PRESENT(eps_pol)) eps_pol = ewald_env%eps_pol
148 552622 : IF (PRESENT(gmax)) gmax = ewald_env%gmax
149 535506 : IF (PRESENT(ns_max)) ns_max = ewald_env%ns_max
150 535506 : IF (PRESENT(o_spline)) o_spline = ewald_env%o_spline
151 535506 : IF (PRESENT(group)) group = ewald_env%para_env
152 535506 : IF (PRESENT(para_env)) para_env => ewald_env%para_env
153 535506 : IF (PRESENT(poisson_section)) poisson_section => ewald_env%poisson_section
154 535506 : IF (PRESENT(interaction_cutoffs)) interaction_cutoffs => &
155 79000 : ewald_env%interaction_cutoffs
156 1613765 : IF (PRESENT(cell_hmat)) cell_hmat = ewald_env%cell_hmat
157 535506 : END SUBROUTINE ewald_env_get
158 :
159 : ! **************************************************************************************************
160 : !> \brief Purpose: Set the EWALD environment.
161 : !> \param ewald_env the pointer to the ewald_env
162 : !> \param ewald_type ...
163 : !> \param alpha ...
164 : !> \param epsilon ...
165 : !> \param eps_pol ...
166 : !> \param gmax ...
167 : !> \param ns_max ...
168 : !> \param precs ...
169 : !> \param o_spline ...
170 : !> \param para_env ...
171 : !> \param poisson_section ...
172 : !> \param interaction_cutoffs ...
173 : !> \param cell_hmat ...
174 : !> \par History
175 : !> 11/03
176 : !> \author CJM
177 : ! **************************************************************************************************
178 23020 : SUBROUTINE ewald_env_set(ewald_env, ewald_type, alpha, epsilon, eps_pol, &
179 : gmax, ns_max, precs, o_spline, para_env, poisson_section, &
180 : interaction_cutoffs, cell_hmat)
181 :
182 : TYPE(ewald_environment_type), INTENT(INOUT) :: ewald_env
183 : INTEGER, OPTIONAL :: ewald_type
184 : REAL(KIND=dp), OPTIONAL :: alpha, epsilon, eps_pol
185 : INTEGER, OPTIONAL :: gmax(3), ns_max
186 : REAL(KIND=dp), OPTIONAL :: precs
187 : INTEGER, OPTIONAL :: o_spline
188 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
189 : TYPE(section_vals_type), OPTIONAL, POINTER :: poisson_section
190 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
191 : POINTER :: interaction_cutoffs
192 : REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL :: cell_hmat
193 :
194 23020 : IF (PRESENT(ewald_type)) ewald_env%ewald_type = ewald_type
195 23020 : IF (PRESENT(alpha)) ewald_env%alpha = alpha
196 23020 : IF (PRESENT(precs)) ewald_env%precs = precs
197 23020 : IF (PRESENT(epsilon)) ewald_env%epsilon = epsilon
198 23020 : IF (PRESENT(eps_pol)) ewald_env%eps_pol = eps_pol
199 23020 : IF (PRESENT(gmax)) ewald_env%gmax = gmax
200 23020 : IF (PRESENT(ns_max)) ewald_env%ns_max = ns_max
201 23020 : IF (PRESENT(o_spline)) ewald_env%o_spline = o_spline
202 23020 : IF (PRESENT(para_env)) ewald_env%para_env => para_env
203 23020 : IF (PRESENT(poisson_section)) THEN
204 4271 : CALL section_vals_retain(poisson_section)
205 4271 : CALL section_vals_release(ewald_env%poisson_section)
206 4271 : ewald_env%poisson_section => poisson_section
207 : END IF
208 23020 : IF (PRESENT(interaction_cutoffs)) ewald_env%interaction_cutoffs => &
209 2643 : interaction_cutoffs
210 232398 : IF (PRESENT(cell_hmat)) ewald_env%cell_hmat = cell_hmat
211 23020 : END SUBROUTINE ewald_env_set
212 :
213 : ! **************************************************************************************************
214 : !> \brief allocates and intitializes a ewald_env
215 : !> \param ewald_env the object to create
216 : !> \param para_env ...
217 : !> \par History
218 : !> 12.2002 created [fawzi]
219 : !> \author Fawzi Mohamed
220 : ! **************************************************************************************************
221 72607 : SUBROUTINE ewald_env_create(ewald_env, para_env)
222 : TYPE(ewald_environment_type), INTENT(OUT) :: ewald_env
223 : TYPE(mp_para_env_type), POINTER :: para_env
224 :
225 : NULLIFY (ewald_env%poisson_section)
226 4271 : CALL para_env%retain()
227 4271 : ewald_env%para_env => para_env
228 4271 : NULLIFY (ewald_env%interaction_cutoffs) ! allocated and initialized later
229 4271 : END SUBROUTINE ewald_env_create
230 :
231 : ! **************************************************************************************************
232 : !> \brief releases the given ewald_env (see doc/ReferenceCounting.html)
233 : !> \param ewald_env the object to release
234 : !> \par History
235 : !> 12.2002 created [fawzi]
236 : !> \author Fawzi Mohamed
237 : ! **************************************************************************************************
238 4271 : SUBROUTINE ewald_env_release(ewald_env)
239 : TYPE(ewald_environment_type), INTENT(INOUT) :: ewald_env
240 :
241 4271 : CALL mp_para_env_release(ewald_env%para_env)
242 4271 : CALL section_vals_release(ewald_env%poisson_section)
243 4271 : IF (ASSOCIATED(ewald_env%interaction_cutoffs)) THEN
244 2643 : DEALLOCATE (ewald_env%interaction_cutoffs)
245 : END IF
246 :
247 4271 : END SUBROUTINE ewald_env_release
248 :
249 : ! **************************************************************************************************
250 : !> \brief Purpose: read the EWALD section
251 : !> \param ewald_env the pointer to the ewald_env
252 : !> \param ewald_section ...
253 : !> \author Teodoro Laino [tlaino] -University of Zurich - 2005
254 : ! **************************************************************************************************
255 3835 : SUBROUTINE read_ewald_section(ewald_env, ewald_section)
256 : TYPE(ewald_environment_type), INTENT(INOUT) :: ewald_env
257 : TYPE(section_vals_type), POINTER :: ewald_section
258 :
259 : INTEGER :: iw
260 3835 : INTEGER, DIMENSION(:), POINTER :: gmax_read
261 : LOGICAL :: explicit
262 : REAL(KIND=dp) :: dummy
263 : TYPE(cp_logger_type), POINTER :: logger
264 : TYPE(enumeration_type), POINTER :: enum
265 : TYPE(keyword_type), POINTER :: keyword
266 : TYPE(section_type), POINTER :: section
267 : TYPE(section_vals_type), POINTER :: multipole_section
268 :
269 3835 : NULLIFY (enum, keyword, section, multipole_section)
270 7670 : logger => cp_get_default_logger()
271 3835 : CALL section_vals_val_get(ewald_section, "EWALD_TYPE", i_val=ewald_env%ewald_type)
272 3835 : CALL section_vals_val_get(ewald_section, "ALPHA", r_val=ewald_env%alpha)
273 3835 : CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", r_val=ewald_env%precs)
274 :
275 3835 : IF (ewald_env%ewald_type == do_ewald_none) THEN
276 1050 : ewald_env%rcut = 0.0_dp
277 : ELSE
278 2785 : CALL section_vals_val_get(ewald_section, "RCUT", explicit=explicit)
279 2785 : IF (explicit) THEN
280 8 : CALL section_vals_val_get(ewald_section, "RCUT", r_val=ewald_env%rcut)
281 : ELSE
282 2777 : ewald_env%rcut = find_ewald_optimal_value(ewald_env%precs)/ewald_env%alpha
283 : END IF
284 : END IF
285 : ! we have no defaults for gmax, gmax is only needed for ewald and spme
286 6534 : SELECT CASE (ewald_env%ewald_type)
287 : CASE (do_ewald_ewald, do_ewald_spme)
288 2699 : CALL section_vals_val_get(ewald_section, "GMAX", i_vals=gmax_read)
289 5088 : SELECT CASE (SIZE(gmax_read, 1))
290 : CASE (1)
291 9556 : ewald_env%gmax = gmax_read(1)
292 : CASE (3)
293 1240 : ewald_env%gmax = gmax_read
294 : CASE DEFAULT
295 2699 : CPABORT("")
296 : END SELECT
297 2699 : IF (ewald_env%ewald_type == do_ewald_spme) THEN
298 1878 : CALL section_vals_val_get(ewald_section, "O_SPLINE", i_val=ewald_env%o_spline)
299 : END IF
300 : CASE (do_ewald_pme)
301 86 : CALL section_vals_val_get(ewald_section, "NS_MAX", i_val=ewald_env%ns_max)
302 86 : CALL section_vals_val_get(ewald_section, "EPSILON", r_val=ewald_env%epsilon)
303 : CASE DEFAULT
304 : ! this should not be used for do_ewald_none
305 4200 : ewald_env%gmax = HUGE(0)
306 4885 : ewald_env%ns_max = HUGE(0)
307 : END SELECT
308 :
309 : ! Multipoles
310 3835 : multipole_section => section_vals_get_subs_vals(ewald_section, "MULTIPOLES")
311 3835 : CALL section_vals_val_get(multipole_section, "_SECTION_PARAMETERS_", l_val=ewald_env%do_multipoles)
312 3835 : CALL section_vals_val_get(multipole_section, "POL_SCF", i_val=ewald_env%do_ipol)
313 3835 : CALL section_vals_val_get(multipole_section, "EPS_POL", r_val=ewald_env%eps_pol)
314 3835 : IF (ewald_env%do_multipoles) THEN
315 336 : SELECT CASE (ewald_env%ewald_type)
316 : CASE (do_ewald_ewald)
317 : CALL section_vals_val_get(multipole_section, "MAX_MULTIPOLE_EXPANSION", &
318 168 : i_val=ewald_env%max_multipole)
319 168 : CALL section_vals_val_get(multipole_section, "MAX_IPOL_ITER", i_val=ewald_env%max_ipol_iter)
320 : CASE DEFAULT
321 168 : CPABORT("Multipole code works at the moment only with standard EWALD sums.")
322 : END SELECT
323 : END IF
324 :
325 : iw = cp_print_key_unit_nr(logger, ewald_section, "PRINT%PROGRAM_RUN_INFO", &
326 3835 : extension=".log")
327 3835 : IF (iw > 0) THEN
328 1906 : NULLIFY (keyword, enum)
329 1906 : CALL create_ewald_section(section)
330 1906 : IF (ewald_env%ewald_type /= do_ewald_none) THEN
331 1419 : keyword => section_get_keyword(section, "EWALD_TYPE")
332 1419 : CALL keyword_get(keyword, enum=enum)
333 1419 : WRITE (iw, '(/,T2,"EWALD| ",A,T67,A14 )') 'Summation is done by:', &
334 2838 : ADJUSTR(TRIM(enum_i2c(enum, ewald_env%ewald_type)))
335 1419 : IF (ewald_env%do_multipoles) THEN
336 84 : NULLIFY (keyword, enum)
337 84 : keyword => section_get_keyword(section, "MULTIPOLES%MAX_MULTIPOLE_EXPANSION")
338 84 : CALL keyword_get(keyword, enum=enum)
339 84 : WRITE (iw, '( T2,"EWALD| ",A )') 'Enabled Multipole Method'
340 84 : WRITE (iw, '( T2,"EWALD| ",A,T67,A14 )') 'Max Term in Multipole Expansion :', &
341 168 : ADJUSTR(TRIM(enum_i2c(enum, ewald_env%max_multipole)))
342 84 : WRITE (iw, '( T2,"EWALD| ",A,T67,3I10 )') 'Max number Iterations for IPOL :', &
343 168 : ewald_env%max_ipol_iter
344 : END IF
345 1419 : dummy = cp_unit_from_cp2k(ewald_env%alpha, "angstrom^-1")
346 : WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') &
347 1419 : 'Alpha parameter [', 'ANGSTROM^-1', ']', dummy
348 1419 : dummy = cp_unit_from_cp2k(ewald_env%rcut, "angstrom")
349 : WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') &
350 1419 : 'Real Space Cutoff [', 'ANGSTROM', ']', dummy
351 :
352 1869 : SELECT CASE (ewald_env%ewald_type)
353 : CASE (do_ewald_ewald)
354 : WRITE (iw, '( T2,"EWALD| ",A,T51,3I10 )') &
355 450 : 'G-space max. Miller index', ewald_env%gmax
356 : CASE (do_ewald_pme)
357 : WRITE (iw, '( T2,"EWALD| ",A,T71,I10 )') &
358 43 : 'Max small-grid points (input) ', ewald_env%ns_max
359 : WRITE (iw, '( T2,"EWALD| ",A,T71,E10.4 )') &
360 43 : 'Gaussian tolerance (input) ', ewald_env%epsilon
361 : CASE (do_ewald_spme)
362 : WRITE (iw, '( T2,"EWALD| ",A,T51,3I10 )') &
363 926 : 'G-space max. Miller index', ewald_env%gmax
364 : WRITE (iw, '( T2,"EWALD| ",A,T71,I10 )') &
365 926 : 'Spline interpolation order ', ewald_env%o_spline
366 : CASE DEFAULT
367 1419 : CPABORT("")
368 : END SELECT
369 : ELSE
370 487 : WRITE (iw, '( T2,"EWALD| ",T73, A )') 'not used'
371 : END IF
372 1906 : CALL section_release(section)
373 : END IF
374 : CALL cp_print_key_finished_output(iw, logger, ewald_section, &
375 3835 : "PRINT%PROGRAM_RUN_INFO")
376 :
377 3835 : END SUBROUTINE read_ewald_section
378 :
379 : ! **************************************************************************************************
380 : !> \brief Purpose: read the EWALD section for TB methods
381 : !> \param ewald_env the pointer to the ewald_env
382 : !> \param ewald_section ...
383 : !> \param hmat ...
384 : !> \param silent ...
385 : !> \param pset ...
386 : !> \author JGH
387 : ! **************************************************************************************************
388 3052 : SUBROUTINE read_ewald_section_tb(ewald_env, ewald_section, hmat, silent, pset)
389 : TYPE(ewald_environment_type), INTENT(INOUT) :: ewald_env
390 : TYPE(section_vals_type), POINTER :: ewald_section
391 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat
392 : LOGICAL, INTENT(IN), OPTIONAL :: silent
393 : CHARACTER(LEN=*), OPTIONAL :: pset
394 :
395 : CHARACTER(LEN=5) :: param
396 : INTEGER :: i, iw, n(3)
397 436 : INTEGER, DIMENSION(:), POINTER :: gmax_read
398 : LOGICAL :: do_print, explicit
399 : REAL(KIND=dp) :: alat, cutoff, dummy, omega
400 : TYPE(cp_logger_type), POINTER :: logger
401 :
402 872 : logger => cp_get_default_logger()
403 436 : do_print = .TRUE.
404 436 : IF (PRESENT(silent)) do_print = .NOT. silent
405 436 : param = "none"
406 436 : IF (PRESENT(pset)) param = pset
407 :
408 436 : ewald_env%do_multipoles = .FALSE.
409 436 : ewald_env%do_ipol = 0
410 436 : ewald_env%eps_pol = 1.e-12_dp
411 436 : ewald_env%max_multipole = 0
412 436 : ewald_env%max_ipol_iter = 0
413 436 : ewald_env%epsilon = 1.e-12_dp
414 436 : ewald_env%ns_max = HUGE(0)
415 :
416 436 : CALL section_vals_val_get(ewald_section, "EWALD_TYPE", explicit=explicit)
417 436 : IF (explicit) THEN
418 160 : CALL section_vals_val_get(ewald_section, "EWALD_TYPE", i_val=ewald_env%ewald_type)
419 160 : IF (ewald_env%ewald_type /= do_ewald_spme) THEN
420 0 : CPABORT("TB needs EWALD_TYPE SPME")
421 : END IF
422 : ELSE
423 276 : ewald_env%ewald_type = do_ewald_spme
424 : END IF
425 :
426 436 : CALL section_vals_val_get(ewald_section, "ALPHA", explicit=explicit)
427 436 : IF (explicit) THEN
428 158 : CALL section_vals_val_get(ewald_section, "ALPHA", r_val=ewald_env%alpha)
429 : ELSE
430 198 : SELECT CASE (param)
431 : CASE DEFAULT
432 198 : ewald_env%alpha = 1.0_dp
433 : CASE ("EEQ")
434 80 : omega = ABS(det_3x3(hmat))
435 278 : ewald_env%alpha = SQRT(twopi)/omega**(1./3.)
436 : END SELECT
437 : END IF
438 :
439 436 : CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", explicit=explicit)
440 436 : IF (explicit) THEN
441 0 : CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", r_val=ewald_env%precs)
442 : ELSE
443 436 : CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", r_val=ewald_env%precs)
444 : END IF
445 :
446 436 : CALL section_vals_val_get(ewald_section, "O_SPLINE", explicit=explicit)
447 436 : IF (explicit) THEN
448 108 : CALL section_vals_val_get(ewald_section, "O_SPLINE", i_val=ewald_env%o_spline)
449 : ELSE
450 248 : SELECT CASE (param)
451 : CASE DEFAULT
452 248 : ewald_env%o_spline = 6
453 : CASE ("EEQ")
454 328 : ewald_env%o_spline = 4
455 : END SELECT
456 : END IF
457 :
458 436 : CALL section_vals_val_get(ewald_section, "RCUT", explicit=explicit)
459 436 : IF (explicit) THEN
460 0 : CALL section_vals_val_get(ewald_section, "RCUT", r_val=ewald_env%rcut)
461 : ELSE
462 436 : ewald_env%rcut = find_ewald_optimal_value(ewald_env%precs)/ewald_env%alpha
463 : END IF
464 :
465 436 : CALL section_vals_val_get(ewald_section, "GMAX", explicit=explicit)
466 436 : IF (explicit) THEN
467 156 : CALL section_vals_val_get(ewald_section, "GMAX", i_vals=gmax_read)
468 254 : SELECT CASE (SIZE(gmax_read, 1))
469 : CASE (1)
470 392 : ewald_env%gmax = gmax_read(1)
471 : CASE (3)
472 232 : ewald_env%gmax = gmax_read
473 : CASE DEFAULT
474 156 : CPABORT("")
475 : END SELECT
476 : ELSE
477 200 : SELECT CASE (param)
478 : CASE DEFAULT
479 : ! set GMAX using ECUT=alpha*45 Ry
480 200 : cutoff = 45._dp*ewald_env%alpha
481 : CASE ("EEQ")
482 : ! set GMAX using ECUT=alpha*45 Ry
483 280 : cutoff = 30._dp*ewald_env%alpha
484 : END SELECT
485 1120 : DO i = 1, 3
486 3360 : alat = SUM(hmat(:, i)**2)
487 840 : CPASSERT(alat /= 0._dp)
488 1120 : ewald_env%gmax(i) = 2*FLOOR(SQRT(2.0_dp*cutoff*alat)/twopi) + 1
489 : END DO
490 : END IF
491 1744 : n = ewald_env%gmax
492 1744 : ewald_env%gmax = pw_grid_n_for_fft(n, odd=.TRUE.)
493 :
494 : iw = cp_print_key_unit_nr(logger, ewald_section, "PRINT%PROGRAM_RUN_INFO", &
495 436 : extension=".log")
496 436 : IF (iw > 0 .AND. do_print) THEN
497 193 : WRITE (iw, '(/,T2,"EWALD| ",A,T67,A14 )') 'Summation is done by:', ADJUSTR("SPME")
498 193 : dummy = cp_unit_from_cp2k(ewald_env%alpha, "angstrom^-1")
499 : WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') &
500 193 : 'Alpha parameter [', 'ANGSTROM^-1', ']', dummy
501 193 : dummy = cp_unit_from_cp2k(ewald_env%rcut, "angstrom")
502 : WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') &
503 193 : 'Real Space Cutoff [', 'ANGSTROM', ']', dummy
504 : WRITE (iw, '( T2,"EWALD| ",A,T51,3I10 )') &
505 193 : 'G-space max. Miller index', ewald_env%gmax
506 : WRITE (iw, '( T2,"EWALD| ",A,T71,I10 )') &
507 193 : 'Spline interpolation order ', ewald_env%o_spline
508 : END IF
509 : CALL cp_print_key_finished_output(iw, logger, ewald_section, &
510 436 : "PRINT%PROGRAM_RUN_INFO")
511 :
512 436 : END SUBROUTINE read_ewald_section_tb
513 :
514 : ! **************************************************************************************************
515 : !> \brief triggers (by bisection) the optimal value for EWALD parameter x
516 : !> EXP(-x^2)/x^2 = EWALD_ACCURACY
517 : !> \param precs ...
518 : !> \return ...
519 : !> \author Teodoro Laino [tlaino] - University of Zurich - 12.2007
520 : ! **************************************************************************************************
521 3213 : FUNCTION find_ewald_optimal_value(precs) RESULT(value)
522 : REAL(KIND=dp) :: precs, value
523 :
524 : REAL(KIND=dp) :: func, func1, func2, s, s1, s2
525 :
526 3213 : s = 0.1_dp
527 3213 : func = EXP(-s**2)/s**2 - precs
528 3213 : CPASSERT(func > 0.0_dp)
529 105487 : DO WHILE (func > 0.0_dp)
530 102274 : s = s + 0.1_dp
531 105487 : func = EXP(-s**2)/s**2 - precs
532 : END DO
533 3213 : s2 = s
534 3213 : s1 = s - 0.1_dp
535 : ! Start bisection
536 : DO WHILE (.TRUE.)
537 80353 : func2 = EXP(-s2**2)/s2**2 - precs
538 80353 : func1 = EXP(-s1**2)/s1**2 - precs
539 80353 : CPASSERT(func1 >= 0)
540 80353 : CPASSERT(func2 <= 0)
541 80353 : s = 0.5_dp*(s1 + s2)
542 80353 : func = EXP(-s**2)/s**2 - precs
543 80353 : IF (func > 0.0_dp) THEN
544 : s1 = s
545 29166 : ELSE IF (func < 0.0_dp) THEN
546 29166 : s2 = s
547 : END IF
548 80353 : IF (ABS(func) < 100.0_dp*EPSILON(0.0_dp)) EXIT
549 : END DO
550 3213 : value = s
551 3213 : END FUNCTION find_ewald_optimal_value
552 :
553 0 : END MODULE ewald_environment_types
554 :
|
