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