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 : !> \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_thermo_al, do_thermo_csvr, do_thermo_gle, do_thermo_nose, &
51 : do_thermo_same_as_part, npe_f_ensemble, npe_i_ensemble, npt_f_ensemble, npt_i_ensemble, &
52 : 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 15921 : 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, save_mem, shell_adiabatic, &
126 : 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, thermo_baro_section, thermo_part_section, thermo_shell_section, &
139 : work_section
140 :
141 1769 : NULLIFY (qmmm_env, cell)
142 1769 : ALLOCATE (thermostats)
143 1769 : CALL allocate_thermostats(thermostats)
144 1769 : adiabatic_fast_section => section_vals_get_subs_vals(md_section, "ADIABATIC_DYNAMICS%THERMOSTAT_FAST")
145 1769 : adiabatic_slow_section => section_vals_get_subs_vals(md_section, "ADIABATIC_DYNAMICS%THERMOSTAT_SLOW")
146 1769 : thermo_part_section => section_vals_get_subs_vals(md_section, "THERMOSTAT")
147 1769 : thermo_shell_section => section_vals_get_subs_vals(md_section, "SHELL%THERMOSTAT")
148 1769 : thermo_baro_section => section_vals_get_subs_vals(md_section, "BAROSTAT%THERMOSTAT")
149 1769 : barostat_section => section_vals_get_subs_vals(md_section, "BAROSTAT")
150 1769 : print_section => section_vals_get_subs_vals(md_section, "PRINT")
151 :
152 1769 : CALL force_env_get(force_env, qmmm_env=qmmm_env, subsys=subsys, cell=cell)
153 1769 : CALL section_vals_get(barostat_section, explicit=explicit_barostat_section)
154 1769 : CALL section_vals_val_get(global_section, "SAVE_MEM", l_val=save_mem)
155 1769 : CALL section_vals_get(thermo_part_section, explicit=explicit_part)
156 1769 : CALL section_vals_get(thermo_shell_section, explicit=explicit_shell)
157 1769 : CALL section_vals_get(thermo_baro_section, explicit=explicit_baro)
158 1769 : CALL section_vals_get(adiabatic_fast_section, explicit=explicit_adiabatic_fast)
159 1769 : CALL section_vals_get(adiabatic_slow_section, explicit=explicit_adiabatic_slow)
160 :
161 1769 : apply_thermo_adiabatic = (simpar%ensemble == nvt_adiabatic_ensemble)
162 :
163 : apply_thermo_baro = (simpar%ensemble == npt_f_ensemble) .OR. &
164 : (simpar%ensemble == npt_ia_ensemble) .OR. &
165 : (simpar%ensemble == npt_i_ensemble) .AND. &
166 1769 : (.NOT. apply_thermo_adiabatic)
167 :
168 : apply_general_thermo = apply_thermo_baro .OR. (simpar%ensemble == nvt_ensemble) .AND. &
169 1617 : (.NOT. apply_thermo_adiabatic)
170 :
171 : apply_thermo_shell = (simpar%ensemble == nve_ensemble) .OR. &
172 : (simpar%ensemble == nvt_ensemble) .OR. &
173 : (simpar%ensemble == npt_f_ensemble) .OR. &
174 : (simpar%ensemble == npt_i_ensemble) .OR. &
175 : (simpar%ensemble == npt_ia_ensemble) .OR. &
176 : (simpar%ensemble == npe_i_ensemble) .OR. &
177 : (simpar%ensemble == npe_f_ensemble) .AND. &
178 1769 : (.NOT. apply_thermo_adiabatic)
179 :
180 1769 : binary_restart_file_name = ""
181 : CALL section_vals_val_get(force_env%root_section, "EXT_RESTART%BINARY_RESTART_FILE_NAME", &
182 1769 : c_val=binary_restart_file_name)
183 :
184 : ! Compute Degrees of Freedom
185 1769 : IF (simpar%ensemble == nvt_adiabatic_ensemble) THEN
186 0 : CALL cite_reference(VandeVondele2002)
187 0 : region = do_region_global
188 0 : region_section_fast => section_vals_get_subs_vals(adiabatic_fast_section, "DEFINE_REGION")
189 0 : region_section_slow => section_vals_get_subs_vals(adiabatic_slow_section, "DEFINE_REGION")
190 0 : IF (explicit_adiabatic_fast) CALL section_vals_val_get(adiabatic_fast_section, "REGION", i_val=region)
191 0 : IF (explicit_adiabatic_slow) CALL section_vals_val_get(adiabatic_slow_section, "REGION", i_val=region)
192 : CALL cp_subsys_get(subsys, molecule_kinds=molecule_kinds, local_molecules=local_molecules, &
193 0 : molecules=molecules, gci=gci, particles=particles)
194 : CALL compute_nfree(cell, simpar, molecule_kinds%els, &
195 0 : print_section, particles, gci)
196 0 : IF (explicit_adiabatic_fast .AND. explicit_adiabatic_slow) THEN
197 0 : IF (apply_thermo_adiabatic) THEN
198 0 : ALLOCATE (thermostats%thermostat_fast)
199 : CALL create_thermostat_type(thermostats%thermostat_fast, simpar, adiabatic_fast_section, &
200 0 : label="FAST")
201 0 : ALLOCATE (thermostats%thermostat_slow)
202 : CALL create_thermostat_type(thermostats%thermostat_slow, simpar, adiabatic_slow_section, &
203 0 : label="SLOW")
204 : CALL setup_adiabatic_thermostat_info(thermostats%thermostat_info_fast, &
205 : molecule_kinds%els, local_molecules, molecules, particles, &
206 : region, simpar%ensemble, region_sections=region_section_fast, &
207 0 : qmmm_env=qmmm_env)
208 :
209 : CALL setup_adiabatic_thermostat_info(thermostats%thermostat_info_slow, &
210 : molecule_kinds%els, local_molecules, molecules, particles, &
211 : region, simpar%ensemble, region_sections=region_section_slow, &
212 0 : qmmm_env=qmmm_env)
213 :
214 : ! Initialize or possibly restart Nose on Particles
215 0 : work_section => section_vals_get_subs_vals(adiabatic_fast_section, "NOSE")
216 : CALL initialize_nhc_fast(thermostats%thermostat_info_fast, simpar, local_molecules, &
217 : molecules%els, molecule_kinds%els, para_env, globenv, &
218 : thermostats%thermostat_fast%nhc, nose_section=work_section, gci=gci, &
219 0 : save_mem=save_mem)
220 0 : work_section => section_vals_get_subs_vals(adiabatic_slow_section, "NOSE")
221 : CALL initialize_nhc_slow(thermostats%thermostat_info_slow, simpar, local_molecules, &
222 : molecules%els, molecule_kinds%els, para_env, globenv, &
223 : thermostats%thermostat_slow%nhc, nose_section=work_section, gci=gci, &
224 0 : save_mem=save_mem)
225 : END IF
226 : ELSE
227 : CALL cp_warn(__LOCATION__, &
228 : "Adiabatic Thermostat has been defined but the ensemble provided "// &
229 0 : "does not support thermostat for Particles! Ignoring thermostat input.")
230 : END IF
231 0 : CALL release_thermostat_info(thermostats%thermostat_info_fast)
232 0 : DEALLOCATE (thermostats%thermostat_info_fast)
233 0 : CALL release_thermostat_info(thermostats%thermostat_info_slow)
234 0 : DEALLOCATE (thermostats%thermostat_info_fast)
235 : ELSE
236 1769 : region = do_region_global
237 1769 : region_sections => section_vals_get_subs_vals(thermo_part_section, "DEFINE_REGION")
238 1769 : IF (explicit_part) CALL section_vals_val_get(thermo_part_section, "REGION", i_val=region)
239 : CALL cp_subsys_get(subsys, molecule_kinds=molecule_kinds, local_molecules=local_molecules, &
240 1769 : molecules=molecules, gci=gci, particles=particles)
241 : CALL compute_degrees_of_freedom(thermostats, cell, simpar, molecule_kinds%els, &
242 : local_molecules, molecules, particles, print_section, region_sections, gci, &
243 1769 : region, qmmm_env)
244 :
245 : ! Particles
246 : ! For constant temperature ensembles the thermostat is activated by default
247 1769 : IF (explicit_part) THEN
248 546 : IF (apply_general_thermo) THEN
249 528 : ALLOCATE (thermostats%thermostat_part)
250 : CALL create_thermostat_type(thermostats%thermostat_part, simpar, thermo_part_section, &
251 528 : label="PARTICLES")
252 : ! Initialize thermostat
253 528 : IF (thermostats%thermostat_part%type_of_thermostat == do_thermo_nose) THEN
254 : ! Initialize or possibly restart Nose on Particles
255 396 : work_section => section_vals_get_subs_vals(thermo_part_section, "NOSE")
256 : CALL initialize_nhc_part(thermostats%thermostat_info_part, simpar, local_molecules, &
257 : molecules%els, molecule_kinds%els, para_env, globenv, &
258 : thermostats%thermostat_part%nhc, nose_section=work_section, gci=gci, &
259 396 : save_mem=save_mem, binary_restart_file_name=binary_restart_file_name)
260 132 : ELSE IF (thermostats%thermostat_part%type_of_thermostat == do_thermo_csvr) THEN
261 : ! Initialize or possibly restart CSVR thermostat on Particles
262 124 : work_section => section_vals_get_subs_vals(thermo_part_section, "CSVR")
263 : CALL initialize_csvr_part(thermostats%thermostat_info_part, simpar, local_molecules, &
264 : molecules%els, molecule_kinds%els, para_env, &
265 : thermostats%thermostat_part%csvr, csvr_section=work_section, &
266 124 : gci=gci)
267 8 : ELSE IF (thermostats%thermostat_part%type_of_thermostat == do_thermo_al) THEN
268 : ! Initialize or possibly restart Ad-Langevin thermostat on Particles
269 4 : work_section => section_vals_get_subs_vals(thermo_part_section, "AD_LANGEVIN")
270 : CALL initialize_al_part(thermostats%thermostat_info_part, simpar, local_molecules, &
271 : molecules%els, molecule_kinds%els, para_env, &
272 : thermostats%thermostat_part%al, al_section=work_section, &
273 4 : gci=gci)
274 4 : ELSE IF (thermostats%thermostat_part%type_of_thermostat == do_thermo_gle) THEN
275 : ! Initialize or possibly restart GLE thermostat on Particles
276 4 : work_section => section_vals_get_subs_vals(thermo_part_section, "GLE")
277 : CALL initialize_gle_part(thermostats%thermostat_info_part, simpar, local_molecules, &
278 : molecules%els, molecule_kinds%els, para_env, &
279 : thermostats%thermostat_part%gle, gle_section=work_section, &
280 4 : gci=gci, save_mem=save_mem)
281 : END IF
282 : CALL thermostat_info(thermostats%thermostat_part, "PARTICLES", thermo_part_section, &
283 528 : simpar, para_env)
284 : ELSE
285 : CALL cp_warn(__LOCATION__, &
286 : "Thermostat for Particles has been defined but the ensemble provided "// &
287 18 : "does not support thermostat for Particles! Ignoring thermostat input.")
288 : END IF
289 1223 : ELSE IF (apply_general_thermo) THEN
290 : CALL cp_abort(__LOCATION__, &
291 : "One constant temperature ensemble has been required, but no thermostat for the "// &
292 : "particles has been defined. You may want to change your input and add a "// &
293 0 : "THERMOSTAT section in the MD section.")
294 : END IF
295 :
296 : ! Core-Shell Model
297 1769 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
298 1769 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, shell_present=shell_present)
299 1769 : IF (shell_present) THEN
300 132 : IF (explicit_shell) THEN
301 : ! The thermostat is activated only if explicitely required
302 : ! It can be used to thermalize the shell-core motion when the temperature is not constant (nve, npe)
303 46 : IF (apply_thermo_shell) THEN
304 46 : ALLOCATE (thermostats%thermostat_shell)
305 : CALL create_thermostat_type(thermostats%thermostat_shell, simpar, thermo_shell_section, &
306 46 : label="SHELL")
307 46 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, shell_adiabatic=shell_adiabatic)
308 46 : region_sections => section_vals_get_subs_vals(thermo_shell_section, "DEFINE_REGION")
309 46 : CALL section_vals_val_get(thermo_shell_section, "REGION", i_val=region)
310 : CALL setup_thermostat_info( &
311 : thermostats%thermostat_info_shell, molecule_kinds%els, &
312 : local_molecules, molecules, particles, region, simpar%ensemble, shell=shell_adiabatic, &
313 46 : region_sections=region_sections, qmmm_env=qmmm_env)
314 46 : IF (shell_adiabatic) THEN
315 : ! Initialize thermostat
316 46 : IF (thermostats%thermostat_shell%type_of_thermostat == do_thermo_nose) THEN
317 : ! Initialize or possibly restart Nose on Shells
318 40 : work_section => section_vals_get_subs_vals(thermo_shell_section, "NOSE")
319 : CALL initialize_nhc_shell(thermostats%thermostat_info_shell, simpar, local_molecules, &
320 : molecules%els, molecule_kinds%els, para_env, globenv, &
321 : thermostats%thermostat_shell%nhc, nose_section=work_section, gci=gci, &
322 40 : save_mem=save_mem, binary_restart_file_name=binary_restart_file_name)
323 6 : ELSE IF (thermostats%thermostat_shell%type_of_thermostat == do_thermo_csvr) THEN
324 : ! Initialize or possibly restart CSVR thermostat on Shells
325 6 : work_section => section_vals_get_subs_vals(thermo_shell_section, "CSVR")
326 : CALL initialize_csvr_shell(thermostats%thermostat_info_shell, simpar, local_molecules, &
327 : molecules%els, molecule_kinds%els, para_env, &
328 6 : thermostats%thermostat_shell%csvr, csvr_section=work_section, gci=gci)
329 : END IF
330 : CALL thermostat_info(thermostats%thermostat_shell, "CORE-SHELL", thermo_shell_section, &
331 46 : simpar, para_env)
332 : ELSE
333 : CALL cp_warn(__LOCATION__, &
334 : "Thermostat for Core-Shell motion only with adiabatic shell-model. "// &
335 : "Continuing calculation ignoring the thermostat info! No Thermostat "// &
336 0 : "applied to Shells!")
337 0 : CALL release_thermostat_type(thermostats%thermostat_shell)
338 0 : DEALLOCATE (thermostats%thermostat_shell)
339 0 : CALL release_thermostat_info(thermostats%thermostat_info_shell)
340 0 : DEALLOCATE (thermostats%thermostat_info_shell)
341 : END IF
342 : ELSE
343 : CALL cp_warn(__LOCATION__, &
344 : "Thermostat for Shells has been defined but for the selected ensemble the adiabatic "// &
345 0 : "shell model has not been implemented! Ignoring thermostat input.")
346 : END IF
347 : END IF
348 1637 : ELSE IF (explicit_shell) THEN
349 : CALL cp_warn(__LOCATION__, &
350 : "Thermostat for Shells has been defined but the system provided "// &
351 0 : "does not contain any Shells! Ignoring thermostat input.")
352 : END IF
353 :
354 : ! Barostat Temperature (not necessarily to be controlled by a thermostat)
355 1769 : IF (explicit_barostat_section) THEN
356 174 : simpar%temp_baro_ext = simpar%temp_ext
357 174 : CALL section_vals_val_get(md_section, "BAROSTAT%TEMPERATURE", n_rep_val=n_rep)
358 174 : IF (n_rep /= 0) THEN
359 2 : CALL section_vals_val_get(md_section, "BAROSTAT%TEMPERATURE", r_val=simpar%temp_baro_ext)
360 2 : CPASSERT(simpar%temp_baro_ext >= 0.0_dp)
361 : END IF
362 :
363 : ! Setup Barostat Thermostat
364 174 : IF (apply_thermo_baro) THEN
365 : ! Check if we use the same thermostat as particles
366 152 : CALL section_vals_val_get(thermo_baro_section, "TYPE", i_val=thermostat_type)
367 152 : work_section => thermo_baro_section
368 152 : IF (thermostat_type == do_thermo_same_as_part) work_section => thermo_part_section
369 :
370 152 : ALLOCATE (thermostats%thermostat_baro)
371 : CALL create_thermostat_type(thermostats%thermostat_baro, simpar, work_section, skip_region=.TRUE., &
372 152 : label="BAROSTAT")
373 : ! Initialize thermostat
374 152 : IF (thermostats%thermostat_baro%type_of_thermostat == do_thermo_nose) THEN
375 : ! Initialize or possibly restart Nose on Barostat
376 120 : work_section => section_vals_get_subs_vals(thermo_baro_section, "NOSE")
377 : CALL initialize_nhc_baro(simpar, para_env, globenv, thermostats%thermostat_baro%nhc, &
378 120 : nose_section=work_section, save_mem=save_mem)
379 32 : ELSE IF (thermostats%thermostat_baro%type_of_thermostat == do_thermo_csvr) THEN
380 : ! Initialize or possibly restart CSVR thermostat on Barostat
381 32 : work_section => section_vals_get_subs_vals(thermo_baro_section, "CSVR")
382 : CALL initialize_csvr_baro(simpar, thermostats%thermostat_baro%csvr, &
383 32 : csvr_section=work_section)
384 : END IF
385 : CALL thermostat_info(thermostats%thermostat_baro, "BAROSTAT", thermo_baro_section, &
386 152 : simpar, para_env)
387 : ! If thermostat for barostat uses a diffent kind than the one of the particles
388 : ! let's update infos in the input structure..
389 152 : IF (thermostat_type == do_thermo_same_as_part) THEN
390 148 : CALL update_thermo_baro_section(thermostats%thermostat_baro, thermo_baro_section)
391 : END IF
392 : ELSE
393 22 : IF (explicit_baro) THEN
394 : CALL cp_warn(__LOCATION__, &
395 : "Thermostat for Barostat has been defined but the ensemble provided "// &
396 0 : "does not support thermostat for Barostat! Ignoring thermostat input.")
397 : END IF
398 : ! Let's remove the section
399 22 : CALL section_vals_remove_values(thermo_baro_section)
400 : END IF
401 : END IF
402 :
403 : ! Release the thermostats info..
404 1769 : CALL release_thermostat_info(thermostats%thermostat_info_part)
405 1769 : DEALLOCATE (thermostats%thermostat_info_part)
406 1769 : CALL release_thermostat_info(thermostats%thermostat_info_shell)
407 1769 : DEALLOCATE (thermostats%thermostat_info_shell)
408 :
409 : END IF ! Adiabitic_NVT screening
410 : ! If no thermostats have been allocated deallocate the full structure
411 : IF ((.NOT. ASSOCIATED(thermostats%thermostat_part)) .AND. &
412 : (.NOT. ASSOCIATED(thermostats%thermostat_shell)) .AND. &
413 : (.NOT. ASSOCIATED(thermostats%thermostat_baro)) .AND. &
414 1769 : (.NOT. ASSOCIATED(thermostats%thermostat_fast)) .AND. &
415 : (.NOT. ASSOCIATED(thermostats%thermostat_slow))) THEN
416 1225 : CALL release_thermostats(thermostats)
417 1225 : DEALLOCATE (thermostats)
418 : END IF
419 :
420 1769 : END SUBROUTINE create_thermostats
421 :
422 : ! **************************************************************************************************
423 : !> \brief ...
424 : !> \param thermostat ...
425 : !> \param section ...
426 : !> \par History
427 : !> 10.2007 created [tlaino]
428 : !> \author Teodoro Laino
429 : ! **************************************************************************************************
430 148 : SUBROUTINE update_thermo_baro_section(thermostat, section)
431 : TYPE(thermostat_type), POINTER :: thermostat
432 : TYPE(section_vals_type), POINTER :: section
433 :
434 : TYPE(section_vals_type), POINTER :: work_section
435 :
436 148 : CALL section_vals_val_set(section, "TYPE", i_val=thermostat%type_of_thermostat)
437 264 : SELECT CASE (thermostat%type_of_thermostat)
438 : CASE (do_thermo_nose)
439 116 : work_section => section_vals_get_subs_vals(section, "NOSE")
440 116 : CALL section_vals_val_set(work_section, "LENGTH", i_val=thermostat%nhc%nhc_len)
441 116 : CALL section_vals_val_set(work_section, "YOSHIDA", i_val=thermostat%nhc%nyosh)
442 116 : CALL section_vals_val_set(work_section, "TIMECON", r_val=thermostat%nhc%tau_nhc)
443 116 : CALL section_vals_val_set(work_section, "MTS", i_val=thermostat%nhc%nc)
444 : CASE (do_thermo_csvr)
445 32 : work_section => section_vals_get_subs_vals(section, "CSVR")
446 32 : CALL section_vals_val_set(work_section, "TIMECON", r_val=thermostat%csvr%tau_csvr)
447 : CASE (do_thermo_al)
448 0 : work_section => section_vals_get_subs_vals(section, "AD_LANGEVIN")
449 0 : CALL section_vals_val_set(work_section, "TIMECON_NH", r_val=thermostat%al%tau_nh)
450 148 : CALL section_vals_val_set(work_section, "TIMECON_LANGEVIN", r_val=thermostat%al%tau_langevin)
451 : END SELECT
452 148 : END SUBROUTINE update_thermo_baro_section
453 :
454 : ! **************************************************************************************************
455 : !> \brief ...
456 : !> \param thermostat ...
457 : !> \param label ...
458 : !> \param section ...
459 : !> \param simpar ...
460 : !> \param para_env ...
461 : !> \par History
462 : !> 10.2007 created [tlaino]
463 : !> \author Teodoro Laino
464 : ! **************************************************************************************************
465 1452 : SUBROUTINE thermostat_info(thermostat, label, section, simpar, para_env)
466 : TYPE(thermostat_type), POINTER :: thermostat
467 : CHARACTER(LEN=*), INTENT(IN) :: label
468 : TYPE(section_vals_type), POINTER :: section
469 : TYPE(simpar_type), POINTER :: simpar
470 : TYPE(mp_para_env_type), POINTER :: para_env
471 :
472 : INTEGER :: iw
473 : REAL(KIND=dp) :: kin_energy, pot_energy, tmp
474 : TYPE(cp_logger_type), POINTER :: logger
475 :
476 726 : NULLIFY (logger)
477 726 : logger => cp_get_default_logger()
478 726 : iw = cp_print_key_unit_nr(logger, section, "PRINT%THERMOSTAT_INFO", extension=".log")
479 : ! Total Tehrmostat Energy
480 726 : CALL get_thermostat_energies(thermostat, pot_energy, kin_energy, para_env)
481 726 : IF (iw > 0) THEN
482 : WRITE (iw, '(/,T2,A)') &
483 360 : 'THERMOSTAT| Thermostat information for '//TRIM(label)
484 635 : SELECT CASE (thermostat%type_of_thermostat)
485 : CASE (do_thermo_nose)
486 : WRITE (iw, '(T2,A,T63,A)') &
487 275 : 'THERMOSTAT| Type of thermostat', 'Nose-Hoover-Chains'
488 : WRITE (iw, '(T2,A,T71,I10)') &
489 275 : 'THERMOSTAT| Nose-Hoover-Chain length', thermostat%nhc%nhc_len
490 275 : tmp = cp_unit_from_cp2k(thermostat%nhc%tau_nhc, 'fs')
491 : WRITE (iw, '(T2,A,T61,F20.6)') &
492 275 : 'THERMOSTAT| Nose-Hoover-Chain time constant [fs]', tmp
493 : WRITE (iw, '(T2,A,T71,I10)') &
494 275 : 'THERMOSTAT| Order of Yoshida integrator', thermostat%nhc%nyosh
495 : WRITE (iw, '(T2,A,T71,I10)') &
496 275 : 'THERMOSTAT| Number of multiple time steps', thermostat%nhc%nc
497 : WRITE (iw, '(T2,A,T61,E20.12)') &
498 275 : 'THERMOSTAT| Initial potential energy', pot_energy, &
499 550 : 'THERMOSTAT| Initial kinetic energy', kin_energy
500 : CASE (do_thermo_csvr)
501 : WRITE (iw, '(T2,A,T44,A)') &
502 81 : 'THERMOSTAT| Type of thermostat', 'Canonical Sampling/Velocity Rescaling'
503 81 : tmp = cp_unit_from_cp2k(thermostat%csvr%tau_csvr, 'fs')*0.5_dp*simpar%dt
504 : WRITE (iw, '(T2,A,T61,F20.6)') &
505 81 : 'THERMOSTAT| CSVR time constant [fs]', tmp
506 : WRITE (iw, '(T2,A,T61,E20.12)') &
507 81 : 'THERMOSTAT| Initial kinetic energy', kin_energy
508 : CASE (do_thermo_al)
509 : WRITE (iw, '(T2,A,T44,A)') &
510 2 : 'THERMOSTAT| Type of thermostat', 'Adaptive Langevin'
511 2 : tmp = cp_unit_from_cp2k(thermostat%al%tau_nh, 'fs')
512 : WRITE (iw, '(T2,A,T61,F20.6)') &
513 2 : 'THERMOSTAT| Time constant of Nose-Hoover part [fs]', tmp
514 2 : tmp = cp_unit_from_cp2k(thermostat%al%tau_langevin, 'fs')
515 : WRITE (iw, '(T2,A,T61,F20.6)') &
516 362 : 'THERMOSTAT| Time constant of Langevin part [fs]', tmp
517 : END SELECT
518 : WRITE (iw, '(T2,A)') &
519 360 : 'THERMOSTAT| End of thermostat information for '//TRIM(label)
520 : END IF
521 726 : CALL cp_print_key_finished_output(iw, logger, section, "PRINT%THERMOSTAT_INFO")
522 :
523 726 : END SUBROUTINE thermostat_info
524 :
525 : ! **************************************************************************************************
526 : !> \brief ...
527 : !> \param thermostat ...
528 : !> \param npt ...
529 : !> \param group ...
530 : !> \par History
531 : !> 10.2007 created [tlaino]
532 : !> \author Teodoro Laino
533 : ! **************************************************************************************************
534 4920 : SUBROUTINE apply_thermostat_baro(thermostat, npt, group)
535 : TYPE(thermostat_type), POINTER :: thermostat
536 : TYPE(npt_info_type), DIMENSION(:, :), POINTER :: npt
537 :
538 : CLASS(mp_comm_type), INTENT(IN) :: group
539 :
540 4920 : IF (ASSOCIATED(thermostat)) THEN
541 4360 : IF (thermostat%type_of_thermostat == do_thermo_nose) THEN
542 : ! Apply Nose-Hoover Thermostat
543 3276 : CPASSERT(ASSOCIATED(thermostat%nhc))
544 3276 : CALL lnhc_barostat(thermostat%nhc, npt, group)
545 1084 : ELSE IF (thermostat%type_of_thermostat == do_thermo_csvr) THEN
546 : ! Apply CSVR Thermostat
547 1084 : CPASSERT(ASSOCIATED(thermostat%csvr))
548 1084 : CALL csvr_barostat(thermostat%csvr, npt, group)
549 : END IF
550 : END IF
551 4920 : END SUBROUTINE apply_thermostat_baro
552 :
553 : ! **************************************************************************************************
554 : !> \brief ...
555 : !> \param thermostat ...
556 : !> \param force_env ...
557 : !> \param molecule_kind_set ...
558 : !> \param molecule_set ...
559 : !> \param particle_set ...
560 : !> \param local_molecules ...
561 : !> \param local_particles ...
562 : !> \param group ...
563 : !> \param shell_adiabatic ...
564 : !> \param shell_particle_set ...
565 : !> \param core_particle_set ...
566 : !> \param vel ...
567 : !> \param shell_vel ...
568 : !> \param core_vel ...
569 : !> \par History
570 : !> 10.2007 created [tlaino]
571 : !> \author Teodoro Laino
572 : ! **************************************************************************************************
573 19480 : SUBROUTINE apply_thermostat_particles(thermostat, force_env, molecule_kind_set, molecule_set, &
574 : particle_set, local_molecules, local_particles, &
575 : group, shell_adiabatic, shell_particle_set, &
576 19480 : core_particle_set, vel, shell_vel, core_vel)
577 :
578 : TYPE(thermostat_type), POINTER :: thermostat
579 : TYPE(force_env_type), POINTER :: force_env
580 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
581 : TYPE(molecule_type), POINTER :: molecule_set(:)
582 : TYPE(particle_type), POINTER :: particle_set(:)
583 : TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
584 :
585 : CLASS(mp_comm_type), INTENT(IN) :: group
586 : LOGICAL, INTENT(IN), OPTIONAL :: shell_adiabatic
587 : TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set(:), &
588 : core_particle_set(:)
589 : REAL(KIND=dp), INTENT(INOUT), OPTIONAL :: vel(:, :), shell_vel(:, :), &
590 : core_vel(:, :)
591 :
592 19480 : IF (ASSOCIATED(thermostat)) THEN
593 19480 : IF (thermostat%type_of_thermostat == do_thermo_nose) THEN
594 : ! Apply Nose-Hoover Thermostat
595 13268 : CPASSERT(ASSOCIATED(thermostat%nhc))
596 : CALL lnhc_particles(thermostat%nhc, molecule_kind_set, molecule_set, &
597 : particle_set, local_molecules, group, shell_adiabatic, shell_particle_set, &
598 43778 : core_particle_set, vel, shell_vel, core_vel)
599 6212 : ELSE IF (thermostat%type_of_thermostat == do_thermo_csvr) THEN
600 : ! Apply CSVR Thermostat
601 5396 : CPASSERT(ASSOCIATED(thermostat%csvr))
602 : CALL csvr_particles(thermostat%csvr, molecule_kind_set, molecule_set, &
603 : particle_set, local_molecules, group, shell_adiabatic, shell_particle_set, &
604 18726 : core_particle_set, vel, shell_vel, core_vel)
605 816 : ELSE IF (thermostat%type_of_thermostat == do_thermo_al) THEN
606 : ! Apply AD_LANGEVIN Thermostat
607 16 : CPASSERT(ASSOCIATED(thermostat%al))
608 : CALL al_particles(thermostat%al, force_env, molecule_kind_set, molecule_set, &
609 24 : particle_set, local_molecules, local_particles, group, vel)
610 800 : ELSE IF (thermostat%type_of_thermostat == do_thermo_gle) THEN
611 : ! Apply GLE Thermostat
612 800 : CPASSERT(ASSOCIATED(thermostat%gle))
613 : CALL gle_particles(thermostat%gle, molecule_kind_set, molecule_set, &
614 : particle_set, local_molecules, group, shell_adiabatic, shell_particle_set, &
615 2800 : core_particle_set, vel, shell_vel, core_vel)
616 : END IF
617 : END IF
618 19480 : END SUBROUTINE apply_thermostat_particles
619 :
620 : ! **************************************************************************************************
621 : !> \brief ...
622 : !> \param thermostat ...
623 : !> \param atomic_kind_set ...
624 : !> \param particle_set ...
625 : !> \param local_particles ...
626 : !> \param group ...
627 : !> \param shell_particle_set ...
628 : !> \param core_particle_set ...
629 : !> \param vel ...
630 : !> \param shell_vel ...
631 : !> \param core_vel ...
632 : !> \par History
633 : !> 10.2007 created [tlaino]
634 : !> \author Teodoro Laino
635 : ! **************************************************************************************************
636 5860 : SUBROUTINE apply_thermostat_shells(thermostat, atomic_kind_set, particle_set, &
637 5860 : local_particles, group, shell_particle_set, core_particle_set, vel, shell_vel, &
638 5860 : core_vel)
639 :
640 : TYPE(thermostat_type), POINTER :: thermostat
641 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
642 : TYPE(particle_type), POINTER :: particle_set(:)
643 : TYPE(distribution_1d_type), POINTER :: local_particles
644 :
645 : CLASS(mp_comm_type), INTENT(IN) :: group
646 : TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set(:), &
647 : core_particle_set(:)
648 : REAL(KIND=dp), INTENT(INOUT), OPTIONAL :: vel(:, :), shell_vel(:, :), &
649 : core_vel(:, :)
650 :
651 5860 : IF (ASSOCIATED(thermostat)) THEN
652 1600 : IF (thermostat%type_of_thermostat == do_thermo_nose) THEN
653 : ! Apply Nose-Hoover Thermostat
654 1360 : CPASSERT(ASSOCIATED(thermostat%nhc))
655 : CALL lnhc_shells(thermostat%nhc, atomic_kind_set, particle_set, local_particles, &
656 3400 : group, shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
657 240 : ELSE IF (thermostat%type_of_thermostat == do_thermo_csvr) THEN
658 : ! Apply CSVR Thermostat
659 240 : CPASSERT(ASSOCIATED(thermostat%csvr))
660 : CALL csvr_shells(thermostat%csvr, atomic_kind_set, particle_set, local_particles, &
661 600 : group, shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
662 : END IF
663 : END IF
664 5860 : END SUBROUTINE apply_thermostat_shells
665 :
666 : END MODULE thermostat_methods
|
