Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Methods for Thermostats
10 : !> \author teo [tlaino] - University of Zurich - 10.2007
11 : ! **************************************************************************************************
12 : MODULE thermostat_methods
13 : USE al_system_dynamics, ONLY: al_particles
14 : USE al_system_init, ONLY: initialize_al_part
15 : USE atomic_kind_list_types, ONLY: atomic_kind_list_type
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind_set
18 : USE bibliography, ONLY: VandeVondele2002,&
19 : cite_reference
20 : USE cell_types, ONLY: cell_type
21 : USE cp_log_handling, ONLY: cp_get_default_logger,&
22 : cp_logger_type
23 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
24 : cp_print_key_unit_nr
25 : USE cp_subsys_types, ONLY: cp_subsys_get,&
26 : cp_subsys_type
27 : USE cp_units, ONLY: cp_unit_from_cp2k
28 : USE csvr_system_dynamics, ONLY: csvr_barostat,&
29 : csvr_particles,&
30 : csvr_shells
31 : USE csvr_system_init, ONLY: initialize_csvr_baro,&
32 : initialize_csvr_part,&
33 : initialize_csvr_shell
34 : USE distribution_1d_types, ONLY: distribution_1d_type
35 : USE extended_system_dynamics, ONLY: lnhc_barostat,&
36 : lnhc_particles,&
37 : lnhc_shells
38 : USE extended_system_init, ONLY: initialize_nhc_baro,&
39 : initialize_nhc_fast,&
40 : initialize_nhc_part,&
41 : initialize_nhc_shell,&
42 : initialize_nhc_slow
43 : USE extended_system_types, ONLY: npt_info_type
44 : USE force_env_types, ONLY: force_env_get,&
45 : force_env_type
46 : USE gle_system_dynamics, ONLY: gle_particles,&
47 : initialize_gle_part
48 : USE global_types, ONLY: global_environment_type
49 : USE input_constants, ONLY: &
50 : do_region_global, do_region_thermal, do_thermo_al, do_thermo_csvr, do_thermo_gle, &
51 : do_thermo_nose, do_thermo_same_as_part, npe_f_ensemble, npe_i_ensemble, npt_f_ensemble, &
52 : npt_i_ensemble, npt_ia_ensemble, nve_ensemble, nvt_adiabatic_ensemble, nvt_ensemble
53 : USE input_section_types, ONLY: section_vals_get,&
54 : section_vals_get_subs_vals,&
55 : section_vals_remove_values,&
56 : section_vals_type,&
57 : section_vals_val_get,&
58 : section_vals_val_set
59 : USE kinds, ONLY: default_path_length,&
60 : dp
61 : USE message_passing, ONLY: mp_comm_type,&
62 : mp_para_env_type
63 : USE molecule_kind_list_types, ONLY: molecule_kind_list_type
64 : USE molecule_kind_types, ONLY: molecule_kind_type
65 : USE molecule_list_types, ONLY: molecule_list_type
66 : USE molecule_types, ONLY: global_constraint_type,&
67 : molecule_type
68 : USE particle_list_types, ONLY: particle_list_type
69 : USE particle_types, ONLY: particle_type
70 : USE qmmm_types, ONLY: qmmm_env_type
71 : USE simpar_types, ONLY: simpar_type
72 : USE thermostat_types, ONLY: allocate_thermostats,&
73 : create_thermostat_type,&
74 : release_thermostat_info,&
75 : release_thermostat_type,&
76 : release_thermostats,&
77 : thermostat_type,&
78 : thermostats_type
79 : USE thermostat_utils, ONLY: compute_degrees_of_freedom,&
80 : compute_nfree,&
81 : get_thermostat_energies,&
82 : setup_adiabatic_thermostat_info,&
83 : setup_thermostat_info
84 : #include "../../base/base_uses.f90"
85 :
86 : IMPLICIT NONE
87 :
88 : PRIVATE
89 : PUBLIC :: create_thermostats, &
90 : apply_thermostat_baro, &
91 : apply_thermostat_particles, &
92 : apply_thermostat_shells
93 :
94 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'thermostat_methods'
95 :
96 : CONTAINS
97 :
98 : ! **************************************************************************************************
99 : !> \brief ...
100 : !> \param thermostats ...
101 : !> \param md_section ...
102 : !> \param force_env ...
103 : !> \param simpar ...
104 : !> \param para_env ...
105 : !> \param globenv ...
106 : !> \param global_section ...
107 : !> \par History
108 : !> 10.2007 created [tlaino]
109 : !> \author Teodoro Laino
110 : ! **************************************************************************************************
111 17960 : SUBROUTINE create_thermostats(thermostats, md_section, force_env, simpar, &
112 : para_env, globenv, global_section)
113 : TYPE(thermostats_type), POINTER :: thermostats
114 : TYPE(section_vals_type), POINTER :: md_section
115 : TYPE(force_env_type), POINTER :: force_env
116 : TYPE(simpar_type), POINTER :: simpar
117 : TYPE(mp_para_env_type), POINTER :: para_env
118 : TYPE(global_environment_type), POINTER :: globenv
119 : TYPE(section_vals_type), POINTER :: global_section
120 :
121 : CHARACTER(LEN=default_path_length) :: binary_restart_file_name
122 : INTEGER :: n_rep, region, thermostat_type
123 : LOGICAL :: apply_general_thermo, apply_thermo_adiabatic, apply_thermo_baro, &
124 : apply_thermo_shell, explicit_adiabatic_fast, explicit_adiabatic_slow, explicit_baro, &
125 : explicit_barostat_section, explicit_part, explicit_shell, explicit_thermal, save_mem, &
126 : shell_adiabatic, shell_present
127 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
128 : TYPE(cell_type), POINTER :: cell
129 : TYPE(cp_subsys_type), POINTER :: subsys
130 : TYPE(distribution_1d_type), POINTER :: local_molecules
131 : TYPE(global_constraint_type), POINTER :: gci
132 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
133 : TYPE(molecule_list_type), POINTER :: molecules
134 : TYPE(particle_list_type), POINTER :: particles
135 : TYPE(qmmm_env_type), POINTER :: qmmm_env
136 : TYPE(section_vals_type), POINTER :: adiabatic_fast_section, adiabatic_slow_section, &
137 : barostat_section, print_section, region_section_fast, region_section_slow, &
138 : region_sections, thermal_region_section, thermo_baro_section, thermo_part_section, &
139 : thermo_shell_section, work_section
140 :
141 1796 : NULLIFY (qmmm_env, cell)
142 1796 : ALLOCATE (thermostats)
143 1796 : CALL allocate_thermostats(thermostats)
144 1796 : adiabatic_fast_section => section_vals_get_subs_vals(md_section, "ADIABATIC_DYNAMICS%THERMOSTAT_FAST")
145 1796 : adiabatic_slow_section => section_vals_get_subs_vals(md_section, "ADIABATIC_DYNAMICS%THERMOSTAT_SLOW")
146 1796 : thermal_region_section => section_vals_get_subs_vals(md_section, "THERMAL_REGION")
147 1796 : thermo_part_section => section_vals_get_subs_vals(md_section, "THERMOSTAT")
148 1796 : thermo_shell_section => section_vals_get_subs_vals(md_section, "SHELL%THERMOSTAT")
149 1796 : thermo_baro_section => section_vals_get_subs_vals(md_section, "BAROSTAT%THERMOSTAT")
150 1796 : barostat_section => section_vals_get_subs_vals(md_section, "BAROSTAT")
151 1796 : print_section => section_vals_get_subs_vals(md_section, "PRINT")
152 :
153 1796 : CALL force_env_get(force_env, qmmm_env=qmmm_env, subsys=subsys, cell=cell)
154 1796 : CALL section_vals_get(barostat_section, explicit=explicit_barostat_section)
155 1796 : CALL section_vals_val_get(global_section, "SAVE_MEM", l_val=save_mem)
156 1796 : CALL section_vals_get(thermal_region_section, explicit=explicit_thermal)
157 1796 : CALL section_vals_get(thermo_part_section, explicit=explicit_part)
158 1796 : CALL section_vals_get(thermo_shell_section, explicit=explicit_shell)
159 1796 : CALL section_vals_get(thermo_baro_section, explicit=explicit_baro)
160 1796 : CALL section_vals_get(adiabatic_fast_section, explicit=explicit_adiabatic_fast)
161 1796 : CALL section_vals_get(adiabatic_slow_section, explicit=explicit_adiabatic_slow)
162 :
163 1796 : apply_thermo_adiabatic = (simpar%ensemble == nvt_adiabatic_ensemble)
164 :
165 : apply_thermo_baro = (simpar%ensemble == npt_f_ensemble) .OR. &
166 : (simpar%ensemble == npt_ia_ensemble) .OR. &
167 : (simpar%ensemble == npt_i_ensemble) .AND. &
168 1796 : (.NOT. apply_thermo_adiabatic)
169 :
170 : apply_general_thermo = apply_thermo_baro .OR. (simpar%ensemble == nvt_ensemble) .AND. &
171 1644 : (.NOT. apply_thermo_adiabatic)
172 :
173 : apply_thermo_shell = (simpar%ensemble == nve_ensemble) .OR. &
174 : (simpar%ensemble == nvt_ensemble) .OR. &
175 : (simpar%ensemble == npt_f_ensemble) .OR. &
176 : (simpar%ensemble == npt_i_ensemble) .OR. &
177 : (simpar%ensemble == npt_ia_ensemble) .OR. &
178 : (simpar%ensemble == npe_i_ensemble) .OR. &
179 : (simpar%ensemble == npe_f_ensemble) .AND. &
180 1796 : (.NOT. apply_thermo_adiabatic)
181 :
182 1796 : binary_restart_file_name = ""
183 : CALL section_vals_val_get(force_env%root_section, "EXT_RESTART%BINARY_RESTART_FILE_NAME", &
184 1796 : c_val=binary_restart_file_name)
185 :
186 : ! Compute Degrees of Freedom
187 1796 : IF (simpar%ensemble == nvt_adiabatic_ensemble) THEN
188 0 : CALL cite_reference(VandeVondele2002)
189 0 : region = do_region_global
190 0 : region_section_fast => section_vals_get_subs_vals(adiabatic_fast_section, "DEFINE_REGION")
191 0 : region_section_slow => section_vals_get_subs_vals(adiabatic_slow_section, "DEFINE_REGION")
192 0 : IF (explicit_adiabatic_fast) CALL section_vals_val_get(adiabatic_fast_section, "REGION", i_val=region)
193 0 : IF (explicit_adiabatic_slow) CALL section_vals_val_get(adiabatic_slow_section, "REGION", i_val=region)
194 : CALL cp_subsys_get(subsys, molecule_kinds=molecule_kinds, local_molecules=local_molecules, &
195 0 : molecules=molecules, gci=gci, particles=particles)
196 : CALL compute_nfree(cell, simpar, molecule_kinds%els, &
197 0 : print_section, particles, gci)
198 0 : IF (explicit_adiabatic_fast .AND. explicit_adiabatic_slow) THEN
199 0 : IF (apply_thermo_adiabatic) THEN
200 0 : ALLOCATE (thermostats%thermostat_fast)
201 : CALL create_thermostat_type(thermostats%thermostat_fast, simpar, adiabatic_fast_section, &
202 0 : label="FAST")
203 0 : ALLOCATE (thermostats%thermostat_slow)
204 : CALL create_thermostat_type(thermostats%thermostat_slow, simpar, adiabatic_slow_section, &
205 0 : label="SLOW")
206 : CALL setup_adiabatic_thermostat_info(thermostats%thermostat_info_fast, &
207 : molecule_kinds%els, local_molecules, molecules, particles, &
208 : region, simpar%ensemble, region_sections=region_section_fast, &
209 0 : qmmm_env=qmmm_env)
210 :
211 : CALL setup_adiabatic_thermostat_info(thermostats%thermostat_info_slow, &
212 : molecule_kinds%els, local_molecules, molecules, particles, &
213 : region, simpar%ensemble, region_sections=region_section_slow, &
214 0 : qmmm_env=qmmm_env)
215 :
216 : ! Initialize or possibly restart Nose on Particles
217 0 : work_section => section_vals_get_subs_vals(adiabatic_fast_section, "NOSE")
218 : CALL initialize_nhc_fast(thermostats%thermostat_info_fast, simpar, local_molecules, &
219 : molecules%els, molecule_kinds%els, para_env, globenv, &
220 : thermostats%thermostat_fast%nhc, nose_section=work_section, gci=gci, &
221 0 : save_mem=save_mem)
222 0 : work_section => section_vals_get_subs_vals(adiabatic_slow_section, "NOSE")
223 : CALL initialize_nhc_slow(thermostats%thermostat_info_slow, simpar, local_molecules, &
224 : molecules%els, molecule_kinds%els, para_env, globenv, &
225 : thermostats%thermostat_slow%nhc, nose_section=work_section, gci=gci, &
226 0 : save_mem=save_mem)
227 : END IF
228 : ELSE
229 : CALL cp_warn(__LOCATION__, &
230 : "Adiabatic Thermostat has been defined but the ensemble provided "// &
231 0 : "does not support thermostat for Particles! Ignoring thermostat input.")
232 : END IF
233 0 : CALL release_thermostat_info(thermostats%thermostat_info_fast)
234 0 : DEALLOCATE (thermostats%thermostat_info_fast)
235 0 : CALL release_thermostat_info(thermostats%thermostat_info_slow)
236 0 : DEALLOCATE (thermostats%thermostat_info_fast)
237 : ELSE
238 1796 : IF (explicit_part) THEN
239 544 : CALL section_vals_val_get(thermo_part_section, "REGION", i_val=region)
240 : ELSE
241 1252 : region = do_region_global
242 : END IF
243 : ! Pass region and region_sections to compute_degrees_of_freedom()
244 : ! which calls setup_thermostat_info() which calls get_defined_region_info()
245 : ! in thermostat_utils.F
246 1796 : IF (region == do_region_thermal) THEN
247 0 : IF (explicit_thermal) THEN
248 0 : region_sections => section_vals_get_subs_vals(thermal_region_section, "DEFINE_REGION")
249 : ELSE
250 : CALL cp_abort(__LOCATION__, &
251 : "Thermostat region THERMAL must be used with explicitly defined "// &
252 0 : "thermal regions in MD/THERMAL_REGION section, but none is found!")
253 : END IF
254 : ELSE
255 1796 : region_sections => section_vals_get_subs_vals(thermo_part_section, "DEFINE_REGION")
256 : END IF
257 : CALL cp_subsys_get(subsys, molecule_kinds=molecule_kinds, local_molecules=local_molecules, &
258 1796 : molecules=molecules, gci=gci, particles=particles)
259 : CALL compute_degrees_of_freedom(thermostats, cell, simpar, molecule_kinds%els, &
260 : local_molecules, molecules, particles, print_section, region_sections, gci, &
261 1796 : region, qmmm_env)
262 :
263 : ! Particles
264 : ! For constant temperature ensembles the thermostat is activated by default
265 1796 : IF (explicit_part) THEN
266 544 : IF (apply_general_thermo) THEN
267 526 : ALLOCATE (thermostats%thermostat_part)
268 : CALL create_thermostat_type(thermostats%thermostat_part, simpar, thermo_part_section, &
269 526 : label="PARTICLES")
270 : ! Initialize thermostat
271 526 : IF (thermostats%thermostat_part%type_of_thermostat == do_thermo_nose) THEN
272 : ! Initialize or possibly restart Nose on Particles
273 396 : work_section => section_vals_get_subs_vals(thermo_part_section, "NOSE")
274 : CALL initialize_nhc_part(thermostats%thermostat_info_part, simpar, local_molecules, &
275 : molecules%els, molecule_kinds%els, para_env, globenv, &
276 : thermostats%thermostat_part%nhc, nose_section=work_section, gci=gci, &
277 396 : save_mem=save_mem, binary_restart_file_name=binary_restart_file_name)
278 130 : ELSE IF (thermostats%thermostat_part%type_of_thermostat == do_thermo_csvr) THEN
279 : ! Initialize or possibly restart CSVR thermostat on Particles
280 122 : work_section => section_vals_get_subs_vals(thermo_part_section, "CSVR")
281 : CALL initialize_csvr_part(thermostats%thermostat_info_part, simpar, local_molecules, &
282 : molecules%els, molecule_kinds%els, para_env, &
283 : thermostats%thermostat_part%csvr, csvr_section=work_section, &
284 122 : gci=gci)
285 8 : ELSE IF (thermostats%thermostat_part%type_of_thermostat == do_thermo_al) THEN
286 : ! Initialize or possibly restart Ad-Langevin thermostat on Particles
287 4 : work_section => section_vals_get_subs_vals(thermo_part_section, "AD_LANGEVIN")
288 : CALL initialize_al_part(thermostats%thermostat_info_part, simpar, local_molecules, &
289 : molecules%els, molecule_kinds%els, para_env, &
290 : thermostats%thermostat_part%al, al_section=work_section, &
291 4 : gci=gci)
292 4 : ELSE IF (thermostats%thermostat_part%type_of_thermostat == do_thermo_gle) THEN
293 : ! Initialize or possibly restart GLE thermostat on Particles
294 4 : work_section => section_vals_get_subs_vals(thermo_part_section, "GLE")
295 : CALL initialize_gle_part(thermostats%thermostat_info_part, simpar, local_molecules, &
296 : molecules%els, molecule_kinds%els, para_env, &
297 : thermostats%thermostat_part%gle, gle_section=work_section, &
298 4 : gci=gci, save_mem=save_mem)
299 : END IF
300 : CALL thermostat_info(thermostats%thermostat_part, "PARTICLES", thermo_part_section, &
301 526 : simpar, para_env)
302 : ELSE
303 : CALL cp_warn(__LOCATION__, &
304 : "Thermostat for Particles has been defined but the ensemble provided "// &
305 18 : "does not support thermostat for Particles! Ignoring thermostat input.")
306 : END IF
307 1252 : ELSE IF (apply_general_thermo) THEN
308 : CALL cp_abort(__LOCATION__, &
309 : "One constant temperature ensemble has been required, but no thermostat for the "// &
310 : "particles has been defined. You may want to change your input and add a "// &
311 0 : "THERMOSTAT section in the MD section.")
312 : END IF
313 :
314 : ! Core-Shell Model
315 1796 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
316 1796 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, shell_present=shell_present)
317 1796 : IF (shell_present) THEN
318 132 : IF (explicit_shell) THEN
319 : ! The thermostat is activated only if explicitely required
320 : ! It can be used to thermalize the shell-core motion when the temperature is not constant (nve, npe)
321 46 : IF (apply_thermo_shell) THEN
322 46 : ALLOCATE (thermostats%thermostat_shell)
323 : CALL create_thermostat_type(thermostats%thermostat_shell, simpar, thermo_shell_section, &
324 46 : label="SHELL")
325 46 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, shell_adiabatic=shell_adiabatic)
326 46 : region_sections => section_vals_get_subs_vals(thermo_shell_section, "DEFINE_REGION")
327 46 : CALL section_vals_val_get(thermo_shell_section, "REGION", i_val=region)
328 : CALL setup_thermostat_info( &
329 : thermostats%thermostat_info_shell, molecule_kinds%els, &
330 : local_molecules, molecules, particles, region, simpar%ensemble, shell=shell_adiabatic, &
331 46 : region_sections=region_sections, qmmm_env=qmmm_env)
332 46 : IF (shell_adiabatic) THEN
333 : ! Initialize thermostat
334 46 : IF (thermostats%thermostat_shell%type_of_thermostat == do_thermo_nose) THEN
335 : ! Initialize or possibly restart Nose on Shells
336 40 : work_section => section_vals_get_subs_vals(thermo_shell_section, "NOSE")
337 : CALL initialize_nhc_shell(thermostats%thermostat_info_shell, simpar, local_molecules, &
338 : molecules%els, molecule_kinds%els, para_env, globenv, &
339 : thermostats%thermostat_shell%nhc, nose_section=work_section, gci=gci, &
340 40 : save_mem=save_mem, binary_restart_file_name=binary_restart_file_name)
341 6 : ELSE IF (thermostats%thermostat_shell%type_of_thermostat == do_thermo_csvr) THEN
342 : ! Initialize or possibly restart CSVR thermostat on Shells
343 6 : work_section => section_vals_get_subs_vals(thermo_shell_section, "CSVR")
344 : CALL initialize_csvr_shell(thermostats%thermostat_info_shell, simpar, local_molecules, &
345 : molecules%els, molecule_kinds%els, para_env, &
346 6 : thermostats%thermostat_shell%csvr, csvr_section=work_section, gci=gci)
347 : END IF
348 : CALL thermostat_info(thermostats%thermostat_shell, "CORE-SHELL", thermo_shell_section, &
349 46 : simpar, para_env)
350 : ELSE
351 : CALL cp_warn(__LOCATION__, &
352 : "Thermostat for Core-Shell motion only with adiabatic shell-model. "// &
353 : "Continuing calculation ignoring the thermostat info! No Thermostat "// &
354 0 : "applied to Shells!")
355 0 : CALL release_thermostat_type(thermostats%thermostat_shell)
356 0 : DEALLOCATE (thermostats%thermostat_shell)
357 0 : CALL release_thermostat_info(thermostats%thermostat_info_shell)
358 0 : DEALLOCATE (thermostats%thermostat_info_shell)
359 : END IF
360 : ELSE
361 : CALL cp_warn(__LOCATION__, &
362 : "Thermostat for Shells has been defined but for the selected ensemble the adiabatic "// &
363 0 : "shell model has not been implemented! Ignoring thermostat input.")
364 : END IF
365 : END IF
366 1664 : ELSE IF (explicit_shell) THEN
367 : CALL cp_warn(__LOCATION__, &
368 : "Thermostat for Shells has been defined but the system provided "// &
369 0 : "does not contain any Shells! Ignoring thermostat input.")
370 : END IF
371 :
372 : ! Barostat Temperature (not necessarily to be controlled by a thermostat)
373 1796 : IF (explicit_barostat_section) THEN
374 174 : simpar%temp_baro_ext = simpar%temp_ext
375 174 : CALL section_vals_val_get(md_section, "BAROSTAT%TEMPERATURE", n_rep_val=n_rep)
376 174 : IF (n_rep /= 0) THEN
377 2 : CALL section_vals_val_get(md_section, "BAROSTAT%TEMPERATURE", r_val=simpar%temp_baro_ext)
378 2 : CPASSERT(simpar%temp_baro_ext >= 0.0_dp)
379 : END IF
380 :
381 : ! Setup Barostat Thermostat
382 174 : IF (apply_thermo_baro) THEN
383 : ! Check if we use the same thermostat as particles
384 152 : CALL section_vals_val_get(thermo_baro_section, "TYPE", i_val=thermostat_type)
385 152 : work_section => thermo_baro_section
386 152 : IF (thermostat_type == do_thermo_same_as_part) work_section => thermo_part_section
387 :
388 152 : ALLOCATE (thermostats%thermostat_baro)
389 : CALL create_thermostat_type(thermostats%thermostat_baro, simpar, work_section, skip_region=.TRUE., &
390 152 : label="BAROSTAT")
391 : ! Initialize thermostat
392 152 : IF (thermostats%thermostat_baro%type_of_thermostat == do_thermo_nose) THEN
393 : ! Initialize or possibly restart Nose on Barostat
394 120 : work_section => section_vals_get_subs_vals(thermo_baro_section, "NOSE")
395 : CALL initialize_nhc_baro(simpar, para_env, globenv, thermostats%thermostat_baro%nhc, &
396 120 : nose_section=work_section, save_mem=save_mem)
397 32 : ELSE IF (thermostats%thermostat_baro%type_of_thermostat == do_thermo_csvr) THEN
398 : ! Initialize or possibly restart CSVR thermostat on Barostat
399 32 : work_section => section_vals_get_subs_vals(thermo_baro_section, "CSVR")
400 : CALL initialize_csvr_baro(simpar, thermostats%thermostat_baro%csvr, &
401 32 : csvr_section=work_section)
402 : END IF
403 : CALL thermostat_info(thermostats%thermostat_baro, "BAROSTAT", thermo_baro_section, &
404 152 : simpar, para_env)
405 : ! If thermostat for barostat uses a diffent kind than the one of the particles
406 : ! let's update infos in the input structure..
407 152 : IF (thermostat_type == do_thermo_same_as_part) THEN
408 148 : CALL update_thermo_baro_section(thermostats%thermostat_baro, thermo_baro_section)
409 : END IF
410 : ELSE
411 22 : IF (explicit_baro) THEN
412 : CALL cp_warn(__LOCATION__, &
413 : "Thermostat for Barostat has been defined but the ensemble provided "// &
414 0 : "does not support thermostat for Barostat! Ignoring thermostat input.")
415 : END IF
416 : ! Let's remove the section
417 22 : CALL section_vals_remove_values(thermo_baro_section)
418 : END IF
419 : END IF
420 :
421 : ! Release the thermostats info..
422 1796 : CALL release_thermostat_info(thermostats%thermostat_info_part)
423 1796 : DEALLOCATE (thermostats%thermostat_info_part)
424 1796 : CALL release_thermostat_info(thermostats%thermostat_info_shell)
425 1796 : DEALLOCATE (thermostats%thermostat_info_shell)
426 :
427 : END IF ! Adiabitic_NVT screening
428 : ! If no thermostats have been allocated deallocate the full structure
429 : IF ((.NOT. ASSOCIATED(thermostats%thermostat_part)) .AND. &
430 : (.NOT. ASSOCIATED(thermostats%thermostat_shell)) .AND. &
431 : (.NOT. ASSOCIATED(thermostats%thermostat_baro)) .AND. &
432 1796 : (.NOT. ASSOCIATED(thermostats%thermostat_fast)) .AND. &
433 : (.NOT. ASSOCIATED(thermostats%thermostat_slow))) THEN
434 1254 : CALL release_thermostats(thermostats)
435 1254 : DEALLOCATE (thermostats)
436 : END IF
437 :
438 1796 : END SUBROUTINE create_thermostats
439 :
440 : ! **************************************************************************************************
441 : !> \brief ...
442 : !> \param thermostat ...
443 : !> \param section ...
444 : !> \par History
445 : !> 10.2007 created [tlaino]
446 : !> \author Teodoro Laino
447 : ! **************************************************************************************************
448 148 : SUBROUTINE update_thermo_baro_section(thermostat, section)
449 : TYPE(thermostat_type), POINTER :: thermostat
450 : TYPE(section_vals_type), POINTER :: section
451 :
452 : TYPE(section_vals_type), POINTER :: work_section
453 :
454 148 : CALL section_vals_val_set(section, "TYPE", i_val=thermostat%type_of_thermostat)
455 264 : SELECT CASE (thermostat%type_of_thermostat)
456 : CASE (do_thermo_nose)
457 116 : work_section => section_vals_get_subs_vals(section, "NOSE")
458 116 : CALL section_vals_val_set(work_section, "LENGTH", i_val=thermostat%nhc%nhc_len)
459 116 : CALL section_vals_val_set(work_section, "YOSHIDA", i_val=thermostat%nhc%nyosh)
460 116 : CALL section_vals_val_set(work_section, "TIMECON", r_val=thermostat%nhc%tau_nhc)
461 116 : CALL section_vals_val_set(work_section, "MTS", i_val=thermostat%nhc%nc)
462 : CASE (do_thermo_csvr)
463 32 : work_section => section_vals_get_subs_vals(section, "CSVR")
464 32 : CALL section_vals_val_set(work_section, "TIMECON", r_val=thermostat%csvr%tau_csvr)
465 : CASE (do_thermo_al)
466 0 : work_section => section_vals_get_subs_vals(section, "AD_LANGEVIN")
467 0 : CALL section_vals_val_set(work_section, "TIMECON_NH", r_val=thermostat%al%tau_nh)
468 148 : CALL section_vals_val_set(work_section, "TIMECON_LANGEVIN", r_val=thermostat%al%tau_langevin)
469 : END SELECT
470 148 : END SUBROUTINE update_thermo_baro_section
471 :
472 : ! **************************************************************************************************
473 : !> \brief ...
474 : !> \param thermostat ...
475 : !> \param label ...
476 : !> \param section ...
477 : !> \param simpar ...
478 : !> \param para_env ...
479 : !> \par History
480 : !> 10.2007 created [tlaino]
481 : !> \author Teodoro Laino
482 : ! **************************************************************************************************
483 1448 : SUBROUTINE thermostat_info(thermostat, label, section, simpar, para_env)
484 : TYPE(thermostat_type), POINTER :: thermostat
485 : CHARACTER(LEN=*), INTENT(IN) :: label
486 : TYPE(section_vals_type), POINTER :: section
487 : TYPE(simpar_type), POINTER :: simpar
488 : TYPE(mp_para_env_type), POINTER :: para_env
489 :
490 : INTEGER :: iw
491 : REAL(KIND=dp) :: kin_energy, pot_energy, tmp
492 : TYPE(cp_logger_type), POINTER :: logger
493 :
494 724 : NULLIFY (logger)
495 724 : logger => cp_get_default_logger()
496 724 : iw = cp_print_key_unit_nr(logger, section, "PRINT%THERMOSTAT_INFO", extension=".log")
497 : ! Total Tehrmostat Energy
498 724 : CALL get_thermostat_energies(thermostat, pot_energy, kin_energy, para_env)
499 724 : IF (iw > 0) THEN
500 : WRITE (iw, '(/,T2,A)') &
501 359 : 'THERMOSTAT| Thermostat information for '//TRIM(label)
502 634 : SELECT CASE (thermostat%type_of_thermostat)
503 : CASE (do_thermo_nose)
504 : WRITE (iw, '(T2,A,T63,A)') &
505 275 : 'THERMOSTAT| Type of thermostat', 'Nose-Hoover-Chains'
506 : WRITE (iw, '(T2,A,T71,I10)') &
507 275 : 'THERMOSTAT| Nose-Hoover-Chain length', thermostat%nhc%nhc_len
508 275 : tmp = cp_unit_from_cp2k(thermostat%nhc%tau_nhc, 'fs')
509 : WRITE (iw, '(T2,A,T61,F20.6)') &
510 275 : 'THERMOSTAT| Nose-Hoover-Chain time constant [fs]', tmp
511 : WRITE (iw, '(T2,A,T71,I10)') &
512 275 : 'THERMOSTAT| Order of Yoshida integrator', thermostat%nhc%nyosh
513 : WRITE (iw, '(T2,A,T71,I10)') &
514 275 : 'THERMOSTAT| Number of multiple time steps', thermostat%nhc%nc
515 : WRITE (iw, '(T2,A,T61,E20.12)') &
516 275 : 'THERMOSTAT| Initial potential energy', pot_energy, &
517 550 : 'THERMOSTAT| Initial kinetic energy', kin_energy
518 : CASE (do_thermo_csvr)
519 : WRITE (iw, '(T2,A,T44,A)') &
520 80 : 'THERMOSTAT| Type of thermostat', 'Canonical Sampling/Velocity Rescaling'
521 80 : tmp = cp_unit_from_cp2k(thermostat%csvr%tau_csvr, 'fs')*0.5_dp*simpar%dt
522 : WRITE (iw, '(T2,A,T61,F20.6)') &
523 80 : 'THERMOSTAT| CSVR time constant [fs]', tmp
524 : WRITE (iw, '(T2,A,T61,E20.12)') &
525 80 : 'THERMOSTAT| Initial kinetic energy', kin_energy
526 : CASE (do_thermo_al)
527 : WRITE (iw, '(T2,A,T44,A)') &
528 2 : 'THERMOSTAT| Type of thermostat', 'Adaptive Langevin'
529 2 : tmp = cp_unit_from_cp2k(thermostat%al%tau_nh, 'fs')
530 : WRITE (iw, '(T2,A,T61,F20.6)') &
531 2 : 'THERMOSTAT| Time constant of Nose-Hoover part [fs]', tmp
532 2 : tmp = cp_unit_from_cp2k(thermostat%al%tau_langevin, 'fs')
533 : WRITE (iw, '(T2,A,T61,F20.6)') &
534 361 : 'THERMOSTAT| Time constant of Langevin part [fs]', tmp
535 : END SELECT
536 : WRITE (iw, '(T2,A)') &
537 359 : 'THERMOSTAT| End of thermostat information for '//TRIM(label)
538 : END IF
539 724 : CALL cp_print_key_finished_output(iw, logger, section, "PRINT%THERMOSTAT_INFO")
540 :
541 724 : END SUBROUTINE thermostat_info
542 :
543 : ! **************************************************************************************************
544 : !> \brief ...
545 : !> \param thermostat ...
546 : !> \param npt ...
547 : !> \param group ...
548 : !> \par History
549 : !> 10.2007 created [tlaino]
550 : !> \author Teodoro Laino
551 : ! **************************************************************************************************
552 4920 : SUBROUTINE apply_thermostat_baro(thermostat, npt, group)
553 : TYPE(thermostat_type), POINTER :: thermostat
554 : TYPE(npt_info_type), DIMENSION(:, :), POINTER :: npt
555 :
556 : CLASS(mp_comm_type), INTENT(IN) :: group
557 :
558 4920 : IF (ASSOCIATED(thermostat)) THEN
559 4360 : IF (thermostat%type_of_thermostat == do_thermo_nose) THEN
560 : ! Apply Nose-Hoover Thermostat
561 3276 : CPASSERT(ASSOCIATED(thermostat%nhc))
562 3276 : CALL lnhc_barostat(thermostat%nhc, npt, group)
563 1084 : ELSE IF (thermostat%type_of_thermostat == do_thermo_csvr) THEN
564 : ! Apply CSVR Thermostat
565 1084 : CPASSERT(ASSOCIATED(thermostat%csvr))
566 1084 : CALL csvr_barostat(thermostat%csvr, npt, group)
567 : END IF
568 : END IF
569 4920 : END SUBROUTINE apply_thermostat_baro
570 :
571 : ! **************************************************************************************************
572 : !> \brief ...
573 : !> \param thermostat ...
574 : !> \param force_env ...
575 : !> \param molecule_kind_set ...
576 : !> \param molecule_set ...
577 : !> \param particle_set ...
578 : !> \param local_molecules ...
579 : !> \param local_particles ...
580 : !> \param group ...
581 : !> \param shell_adiabatic ...
582 : !> \param shell_particle_set ...
583 : !> \param core_particle_set ...
584 : !> \param vel ...
585 : !> \param shell_vel ...
586 : !> \param core_vel ...
587 : !> \par History
588 : !> 10.2007 created [tlaino]
589 : !> \author Teodoro Laino
590 : ! **************************************************************************************************
591 19480 : SUBROUTINE apply_thermostat_particles(thermostat, force_env, molecule_kind_set, molecule_set, &
592 : particle_set, local_molecules, local_particles, &
593 : group, shell_adiabatic, shell_particle_set, &
594 19480 : core_particle_set, vel, shell_vel, core_vel)
595 :
596 : TYPE(thermostat_type), POINTER :: thermostat
597 : TYPE(force_env_type), POINTER :: force_env
598 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
599 : TYPE(molecule_type), POINTER :: molecule_set(:)
600 : TYPE(particle_type), POINTER :: particle_set(:)
601 : TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
602 :
603 : CLASS(mp_comm_type), INTENT(IN) :: group
604 : LOGICAL, INTENT(IN), OPTIONAL :: shell_adiabatic
605 : TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set(:), &
606 : core_particle_set(:)
607 : REAL(KIND=dp), INTENT(INOUT), OPTIONAL :: vel(:, :), shell_vel(:, :), &
608 : core_vel(:, :)
609 :
610 19480 : IF (ASSOCIATED(thermostat)) THEN
611 19480 : IF (thermostat%type_of_thermostat == do_thermo_nose) THEN
612 : ! Apply Nose-Hoover Thermostat
613 13268 : CPASSERT(ASSOCIATED(thermostat%nhc))
614 : CALL lnhc_particles(thermostat%nhc, molecule_kind_set, molecule_set, &
615 : particle_set, local_molecules, group, shell_adiabatic, shell_particle_set, &
616 43778 : core_particle_set, vel, shell_vel, core_vel)
617 6212 : ELSE IF (thermostat%type_of_thermostat == do_thermo_csvr) THEN
618 : ! Apply CSVR Thermostat
619 5396 : CPASSERT(ASSOCIATED(thermostat%csvr))
620 : CALL csvr_particles(thermostat%csvr, molecule_kind_set, molecule_set, &
621 : particle_set, local_molecules, group, shell_adiabatic, shell_particle_set, &
622 18726 : core_particle_set, vel, shell_vel, core_vel)
623 816 : ELSE IF (thermostat%type_of_thermostat == do_thermo_al) THEN
624 : ! Apply AD_LANGEVIN Thermostat
625 16 : CPASSERT(ASSOCIATED(thermostat%al))
626 : CALL al_particles(thermostat%al, force_env, molecule_kind_set, molecule_set, &
627 24 : particle_set, local_molecules, local_particles, group, vel)
628 800 : ELSE IF (thermostat%type_of_thermostat == do_thermo_gle) THEN
629 : ! Apply GLE Thermostat
630 800 : CPASSERT(ASSOCIATED(thermostat%gle))
631 : CALL gle_particles(thermostat%gle, molecule_kind_set, molecule_set, &
632 : particle_set, local_molecules, group, shell_adiabatic, shell_particle_set, &
633 2800 : core_particle_set, vel, shell_vel, core_vel)
634 : END IF
635 : END IF
636 19480 : END SUBROUTINE apply_thermostat_particles
637 :
638 : ! **************************************************************************************************
639 : !> \brief ...
640 : !> \param thermostat ...
641 : !> \param atomic_kind_set ...
642 : !> \param particle_set ...
643 : !> \param local_particles ...
644 : !> \param group ...
645 : !> \param shell_particle_set ...
646 : !> \param core_particle_set ...
647 : !> \param vel ...
648 : !> \param shell_vel ...
649 : !> \param core_vel ...
650 : !> \par History
651 : !> 10.2007 created [tlaino]
652 : !> \author Teodoro Laino
653 : ! **************************************************************************************************
654 5860 : SUBROUTINE apply_thermostat_shells(thermostat, atomic_kind_set, particle_set, &
655 5860 : local_particles, group, shell_particle_set, core_particle_set, vel, shell_vel, &
656 5860 : core_vel)
657 :
658 : TYPE(thermostat_type), POINTER :: thermostat
659 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
660 : TYPE(particle_type), POINTER :: particle_set(:)
661 : TYPE(distribution_1d_type), POINTER :: local_particles
662 :
663 : CLASS(mp_comm_type), INTENT(IN) :: group
664 : TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set(:), &
665 : core_particle_set(:)
666 : REAL(KIND=dp), INTENT(INOUT), OPTIONAL :: vel(:, :), shell_vel(:, :), &
667 : core_vel(:, :)
668 :
669 5860 : IF (ASSOCIATED(thermostat)) THEN
670 1600 : IF (thermostat%type_of_thermostat == do_thermo_nose) THEN
671 : ! Apply Nose-Hoover Thermostat
672 1360 : CPASSERT(ASSOCIATED(thermostat%nhc))
673 : CALL lnhc_shells(thermostat%nhc, atomic_kind_set, particle_set, local_particles, &
674 3400 : group, shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
675 240 : ELSE IF (thermostat%type_of_thermostat == do_thermo_csvr) THEN
676 : ! Apply CSVR Thermostat
677 240 : CPASSERT(ASSOCIATED(thermostat%csvr))
678 : CALL csvr_shells(thermostat%csvr, atomic_kind_set, particle_set, local_particles, &
679 600 : group, shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
680 : END IF
681 : END IF
682 5860 : END SUBROUTINE apply_thermostat_shells
683 :
684 : END MODULE thermostat_methods
|
