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 Setup of regions with different temperature
10 : !> \par History
11 : !> - Added support for Langevin regions (2014/01/08, LT)
12 : !> - Added print subroutine for langevin regions (2014/02/04, LT)
13 : !> - Changed print_thermal_regions to print_thermal_regions_temperature
14 : !> (2014/02/04, LT)
15 : !> \author MI
16 : ! **************************************************************************************************
17 : MODULE thermal_region_utils
18 :
19 : USE bibliography, ONLY: Kantorovich2008,&
20 : Kantorovich2008a,&
21 : cite_reference
22 : USE cp_log_handling, ONLY: cp_get_default_logger,&
23 : cp_logger_get_default_io_unit,&
24 : cp_logger_type,&
25 : cp_to_string
26 : USE cp_output_handling, ONLY: cp_p_file,&
27 : cp_print_key_finished_output,&
28 : cp_print_key_should_output,&
29 : cp_print_key_unit_nr
30 : USE cp_subsys_types, ONLY: cp_subsys_get,&
31 : cp_subsys_type
32 : USE cp_units, ONLY: cp_unit_from_cp2k,&
33 : cp_unit_to_cp2k
34 : USE force_env_types, ONLY: force_env_get,&
35 : force_env_type
36 : USE force_fields_util, ONLY: get_generic_info
37 : USE fparser, ONLY: EvalErrType,&
38 : evalf,&
39 : finalizef,&
40 : initf,&
41 : parsef
42 : USE input_constants, ONLY: do_region_thermal,&
43 : langevin_ensemble,&
44 : npt_f_ensemble,&
45 : npt_i_ensemble,&
46 : npt_ia_ensemble,&
47 : nvt_ensemble
48 : USE input_section_types, ONLY: section_vals_get,&
49 : section_vals_get_subs_vals,&
50 : section_vals_type,&
51 : section_vals_val_get
52 : USE kinds, ONLY: default_string_length,&
53 : dp
54 : USE molecule_kind_list_types, ONLY: molecule_kind_list_type
55 : USE molecule_kind_types, ONLY: molecule_kind_type
56 : USE molecule_list_types, ONLY: molecule_list_type
57 : USE molecule_types, ONLY: get_molecule,&
58 : molecule_type
59 : USE particle_list_types, ONLY: particle_list_type
60 : USE physcon, ONLY: femtoseconds,&
61 : kelvin
62 : USE simpar_types, ONLY: simpar_type
63 : USE string_utilities, ONLY: integer_to_string
64 : USE thermal_region_types, ONLY: allocate_thermal_regions,&
65 : release_thermal_regions,&
66 : thermal_region_type,&
67 : thermal_regions_type
68 : #include "../base/base_uses.f90"
69 :
70 : IMPLICIT NONE
71 :
72 : PRIVATE
73 : PUBLIC :: create_thermal_regions, &
74 : print_thermal_regions_temperature, &
75 : print_thermal_regions_langevin, &
76 : set_t_region_index
77 :
78 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'thermal_region_utils'
79 :
80 : CONTAINS
81 :
82 : ! **************************************************************************************************
83 : !> \brief create thermal_regions
84 : !> \param thermal_regions ...
85 : !> \param md_section ...
86 : !> \param simpar ...
87 : !> \param force_env ...
88 : !> \par History
89 : !> - Added support for Langevin regions (2014/01/08, LT)
90 : !> \author
91 : ! **************************************************************************************************
92 5388 : SUBROUTINE create_thermal_regions(thermal_regions, md_section, simpar, force_env)
93 : TYPE(thermal_regions_type), POINTER :: thermal_regions
94 : TYPE(section_vals_type), POINTER :: md_section
95 : TYPE(simpar_type), POINTER :: simpar
96 : TYPE(force_env_type), POINTER :: force_env
97 :
98 : CHARACTER(LEN=default_string_length) :: my_region
99 : CHARACTER(LEN=default_string_length), &
100 1796 : DIMENSION(:), POINTER :: tmpstringlist
101 : INTEGER :: first_atom, flen, i, ikind, il, imol, &
102 : ipart, ireg, itimes, last_atom, nlist, &
103 : nnpart, nregions, output_unit, region
104 1796 : INTEGER, DIMENSION(:), POINTER :: tmplist
105 : LOGICAL :: apply_thermostat, do_langevin, &
106 : do_langevin_default, do_read_ngr, &
107 : explicit, use_t, use_temp_tol, &
108 : use_tfunc
109 : REAL(KIND=dp) :: temp, temp_tol
110 : TYPE(cp_logger_type), POINTER :: logger
111 : TYPE(cp_subsys_type), POINTER :: subsys
112 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
113 1796 : TYPE(molecule_kind_type), POINTER :: molecule_kind, molecule_kind_set(:)
114 : TYPE(molecule_list_type), POINTER :: molecules
115 1796 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
116 : TYPE(molecule_type), POINTER :: molecule
117 : TYPE(particle_list_type), POINTER :: particles
118 : TYPE(section_vals_type), POINTER :: region_sections, tfunc_section, &
119 : thermal_region_section
120 : TYPE(thermal_region_type), POINTER :: t_region
121 :
122 1796 : NULLIFY (logger, molecule, molecule_kind, molecule_kind_set, molecule_kinds, &
123 1796 : region_sections, t_region, tfunc_section, &
124 1796 : thermal_region_section, particles, subsys, tmplist)
125 3592 : logger => cp_get_default_logger()
126 1796 : output_unit = cp_logger_get_default_io_unit(logger)
127 1796 : ALLOCATE (thermal_regions)
128 1796 : CALL allocate_thermal_regions(thermal_regions)
129 1796 : thermal_region_section => section_vals_get_subs_vals(md_section, "THERMAL_REGION")
130 1796 : CALL section_vals_get(thermal_region_section, explicit=explicit)
131 1796 : IF (explicit) THEN
132 : apply_thermostat = (simpar%ensemble == nvt_ensemble) .OR. &
133 : (simpar%ensemble == npt_f_ensemble) .OR. &
134 : (simpar%ensemble == npt_ia_ensemble) .OR. &
135 14 : (simpar%ensemble == npt_i_ensemble)
136 0 : IF (apply_thermostat) THEN
137 0 : CALL section_vals_val_get(md_section, "THERMOSTAT%REGION", i_val=region)
138 0 : IF (region /= do_region_thermal) THEN
139 : CALL cp_warn(__LOCATION__, &
140 : "An ensemble has been chosen where one or more thermostats are "// &
141 : "used to control temperature, and one or more thermal regions "// &
142 : "are also defined. To ensure consistency between thermostat "// &
143 0 : "regions and thermal regions, set THERMOSTAT/REGION to THERMAL.")
144 : END IF
145 : END IF
146 : CALL section_vals_val_get(thermal_region_section, "FORCE_RESCALING", &
147 14 : l_val=thermal_regions%force_rescaling)
148 : region_sections => section_vals_get_subs_vals(thermal_region_section, &
149 14 : "DEFINE_REGION")
150 14 : CALL section_vals_get(region_sections, n_repetition=nregions)
151 14 : IF (nregions > 0) THEN
152 14 : thermal_regions%nregions = nregions
153 14 : thermal_regions%section => thermal_region_section
154 62 : ALLOCATE (thermal_regions%thermal_region(nregions))
155 14 : simpar%do_thermal_region = .TRUE.
156 14 : CALL force_env_get(force_env, subsys=subsys)
157 : CALL cp_subsys_get(subsys, molecules=molecules, &
158 : molecule_kinds=molecule_kinds, &
159 14 : particles=particles)
160 14 : molecule_kind_set => molecule_kinds%els
161 14 : molecule_set => molecules%els
162 :
163 34 : DO ireg = 1, nregions
164 20 : NULLIFY (t_region)
165 20 : CALL integer_to_string(ireg, my_region)
166 20 : t_region => thermal_regions%thermal_region(ireg)
167 20 : t_region%region_index = ireg
168 20 : NULLIFY (t_region%part_index)
169 :
170 : ! Set up thermal region mapping to particles
171 : CALL section_vals_val_get(region_sections, "LIST", &
172 20 : i_rep_section=ireg, n_rep_val=nlist)
173 20 : IF (nlist > 0) THEN
174 42 : DO il = 1, nlist
175 : CALL section_vals_val_get(region_sections, "LIST", i_rep_section=ireg, &
176 22 : i_rep_val=il, i_vals=tmplist)
177 162 : DO i = 1, SIZE(tmplist)
178 120 : ipart = tmplist(i)
179 142 : CALL set_t_region_index(particles, ipart, ireg)
180 : END DO
181 : END DO
182 : END IF
183 : CALL section_vals_val_get(region_sections, "MOLNAME", &
184 20 : i_rep_section=ireg, n_rep_val=nlist)
185 20 : IF (nlist > 0) THEN
186 0 : DO il = 1, nlist
187 : CALL section_vals_val_get(region_sections, "MOLNAME", i_rep_section=ireg, &
188 0 : i_rep_val=il, c_vals=tmpstringlist)
189 0 : DO i = 1, SIZE(tmpstringlist)
190 0 : DO ikind = 1, SIZE(molecule_kind_set)
191 0 : molecule_kind => molecule_kind_set(ikind)
192 0 : IF (molecule_kind%name == tmpstringlist(i)) THEN
193 0 : DO imol = 1, SIZE(molecule_kind%molecule_list)
194 0 : molecule => molecule_set(molecule_kind%molecule_list(imol))
195 0 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
196 0 : DO ipart = first_atom, last_atom
197 0 : CALL set_t_region_index(particles, ipart, ireg)
198 : END DO
199 : END DO
200 : END IF
201 : END DO
202 : END DO
203 : END DO
204 : END IF
205 : ! TODO: parsing logic of thermal region should be aligned
206 : ! with thermostat region as in SUBROUTINE get_defined_region_info,
207 : ! src/motion/thermostat/thermostat_utils.F, so that
208 : ! MM_SUBSYS and QM_SUBSYS keywords work for thermal regions
209 : CALL section_vals_val_get(region_sections, "QM_SUBSYS", &
210 20 : i_rep_section=ireg, n_rep_val=nlist, explicit=explicit)
211 20 : IF (explicit) THEN
212 0 : CPABORT("QM_SUBSYS not yet implemented for thermal regions")
213 : END IF
214 : CALL section_vals_val_get(region_sections, "MM_SUBSYS", &
215 20 : i_rep_section=ireg, n_rep_val=nlist, explicit=explicit)
216 20 : IF (explicit) THEN
217 0 : CPABORT("MM_SUBSYS not yet implemented for thermal regions")
218 : END IF
219 :
220 : ! Set up t_region%part_index after parsing input
221 20 : nnpart = 0
222 20 : t_region%npart = 0
223 500 : DO ipart = 1, particles%n_els
224 500 : IF (particles%els(ipart)%t_region_index == ireg) nnpart = nnpart + 1
225 : END DO
226 60 : ALLOCATE (t_region%part_index(nnpart))
227 500 : DO ipart = 1, particles%n_els
228 500 : IF (particles%els(ipart)%t_region_index == ireg) THEN
229 120 : t_region%npart = t_region%npart + 1
230 120 : t_region%part_index(t_region%npart) = ipart
231 : END IF
232 : END DO
233 :
234 : ! Set up temperature function, if any
235 20 : tfunc_section => section_vals_get_subs_vals(region_sections, "TFUNC", i_rep_section=ireg)
236 20 : CALL section_vals_get(tfunc_section, explicit=use_tfunc)
237 20 : IF (use_tfunc) THEN
238 : CALL get_generic_info(tfunc_section, "FUNCTION", &
239 : t_region%temperature_function, &
240 : t_region%temperature_function_parameters, &
241 : t_region%temperature_function_values, &
242 0 : input_variables=["S", "T"], i_rep_sec=1)
243 : END IF
244 :
245 : ! Set up t_region%temp_expected
246 : ! Precedence of temperature setting: the first positive value in
247 : ! 1. MD/THERMAL_REGION/DEFINE_REGION/TEMPERATURE
248 : ! 2. The temperature function evaluated at S = STEP_START_VAL and T = 0
249 : ! 3. MD/TEMPERATURE
250 : ! Otherwise it will be set to 0
251 20 : temp = 0.0_dp
252 : CALL section_vals_val_get(region_sections, "TEMPERATURE", &
253 20 : i_rep_section=ireg, explicit=use_t)
254 20 : IF (use_t) THEN
255 : CALL section_vals_val_get(region_sections, "TEMPERATURE", &
256 20 : i_rep_section=ireg, r_val=temp)
257 : END IF
258 20 : IF (temp > 0.0_dp) THEN
259 20 : t_region%temp_expected = temp
260 : ELSE
261 0 : IF (use_tfunc) THEN
262 0 : CALL initf(1)
263 : CALL parsef(1, TRIM(t_region%temperature_function), &
264 0 : t_region%temperature_function_parameters)
265 0 : CALL section_vals_val_get(md_section, "STEP_START_VAL", i_val=itimes)
266 0 : t_region%temperature_function_values(1) = itimes*1.0_dp
267 0 : t_region%temperature_function_values(2) = 0.0_dp
268 0 : temp = cp_unit_to_cp2k(evalf(1, t_region%temperature_function_values), "K")
269 0 : CPASSERT(EvalErrType <= 0)
270 0 : CALL finalizef()
271 : END IF
272 0 : IF (temp > 0.0_dp) THEN
273 0 : t_region%temp_expected = temp
274 : ELSE
275 0 : temp = simpar%temp_ext
276 0 : IF (temp > 0.0_dp) THEN
277 0 : t_region%temp_expected = temp
278 : ELSE
279 : CALL cp_warn(__LOCATION__, &
280 : "For thermal region "//TRIM(my_region)//" , none of the "// &
281 : "input values for the initial temperature, nor the result "// &
282 : "of evaluating temperature function at STEP_START_VAL, is "// &
283 0 : "positive in Kelvin. A value of zero will be used instead.")
284 0 : t_region%temp_expected = 0.0_dp
285 : END IF
286 : END IF
287 : END IF
288 :
289 : ! Set up t_region%temp_tol
290 : ! Precedence of tolerance setting: the first positive value in
291 : ! 1. MD/THERMAL_REGION/DEFINE_REGION/TEMP_TOL
292 : ! 2. MD/TEMP_TOL
293 : ! Otherwise it will be set to 0 (i.e. not used)
294 20 : temp_tol = 0.0_dp
295 : CALL section_vals_val_get(region_sections, "TEMP_TOL", &
296 20 : i_rep_section=ireg, explicit=use_temp_tol)
297 20 : IF (use_temp_tol) THEN
298 : CALL section_vals_val_get(region_sections, "TEMP_TOL", &
299 2 : i_rep_section=ireg, r_val=temp_tol)
300 : END IF
301 174 : IF (temp_tol > 0.0_dp) THEN
302 2 : t_region%temp_tol = temp_tol
303 : ELSE
304 18 : temp_tol = simpar%temp_tol
305 18 : IF (temp_tol > 0.0_dp) THEN
306 0 : t_region%temp_tol = temp_tol
307 : ELSE
308 18 : t_region%temp_tol = 0.0_dp
309 : END IF
310 : END IF
311 : ! End of one thermal region in the loop
312 : END DO
313 :
314 : ! For Langevin ensemble, determine thermal_regions%do_langevin
315 : ! and t_region%noisy_gamma_region
316 14 : IF (simpar%ensemble == langevin_ensemble) THEN
317 12 : CALL cite_reference(Kantorovich2008)
318 12 : CALL cite_reference(Kantorovich2008a)
319 : CALL section_vals_val_get(thermal_region_section, "DO_LANGEVIN_DEFAULT", &
320 12 : l_val=do_langevin_default)
321 36 : ALLOCATE (thermal_regions%do_langevin(particles%n_els))
322 84 : thermal_regions%do_langevin = do_langevin_default
323 40 : DO ireg = 1, nregions
324 16 : NULLIFY (t_region)
325 16 : t_region => thermal_regions%thermal_region(ireg)
326 16 : CALL integer_to_string(ireg, my_region)
327 : CALL section_vals_val_get(region_sections, "DO_LANGEVIN", &
328 16 : i_rep_section=ireg, l_val=do_langevin)
329 112 : DO ipart = 1, particles%n_els
330 112 : IF (particles%els(ipart)%t_region_index == ireg) THEN
331 44 : thermal_regions%do_langevin(ipart) = do_langevin
332 : END IF
333 : END DO
334 16 : CALL section_vals_val_get(region_sections, "NOISY_GAMMA_REGION", i_rep_section=ireg, explicit=do_read_ngr)
335 44 : IF (do_read_ngr) THEN
336 : CALL section_vals_val_get(region_sections, "NOISY_GAMMA_REGION", i_rep_section=ireg, &
337 6 : r_val=t_region%noisy_gamma_region)
338 6 : IF (.NOT. do_langevin) THEN
339 : CALL cp_warn(__LOCATION__, &
340 : "You provided NOISY_GAMMA_REGION but atoms in thermal region "//TRIM(my_region)// &
341 : " will not undergo Langevin MD. "// &
342 0 : "NOISY_GAMMA_REGION will be ignored and its value discarded!")
343 : END IF
344 : ELSE
345 10 : t_region%noisy_gamma_region = simpar%noisy_gamma
346 : END IF
347 : END DO
348 : END IF
349 : ! End of Langevin treatment
350 :
351 : ! Output thermal region temperature and mapping to particles
352 : ! The region indices are assumed to be 0-999
353 14 : IF (output_unit > 0) THEN
354 : WRITE (output_unit, '(/,T2,A)') &
355 7 : "THERMAL_REGION| Temperature value and function for each region [K]"
356 17 : DO ireg = 1, nregions
357 10 : NULLIFY (t_region)
358 10 : t_region => thermal_regions%thermal_region(ireg)
359 : WRITE (output_unit, '(T2,A,T25,I3,T29,A,T34,I6,T41,A,T66,F15.6)') &
360 10 : "THERMAL_REGION| Region", t_region%region_index, &
361 10 : "with", t_region%npart, &
362 10 : "particles initialized at", &
363 20 : cp_unit_from_cp2k(t_region%temp_expected, "K")
364 10 : flen = LEN(TRIM(t_region%temperature_function))
365 17 : IF (flen > 0) THEN
366 0 : DO i = 1, flen, 64
367 : WRITE (output_unit, '(T2,A,T17,A64)') &
368 0 : "THERMAL_REGION|", t_region%temperature_function(i:MIN(i + 64, flen))
369 : END DO
370 : ELSE
371 : WRITE (output_unit, '(T2,A,T66,F15.6)') &
372 10 : "THERMAL_REGION|", cp_unit_from_cp2k(t_region%temp_expected, "K")
373 : END IF
374 : END DO
375 : WRITE (output_unit, '(/,T2,A)') &
376 7 : "THERMAL_REGION| Mapping of thermal region indices to particles"
377 19 : DO ipart = 1, particles%n_els, 16
378 : WRITE (output_unit, '(T2,A,T17,16(" ",I3))') &
379 151 : "THERMAL_REGION|", particles%els(ipart:MIN(ipart + 15, particles%n_els))%t_region_index
380 : END DO
381 : END IF
382 : ELSE
383 0 : CALL release_thermal_regions(thermal_regions)
384 0 : DEALLOCATE (thermal_regions)
385 0 : simpar%do_thermal_region = .FALSE.
386 : END IF
387 : ELSE
388 1782 : CALL release_thermal_regions(thermal_regions)
389 1782 : DEALLOCATE (thermal_regions)
390 1782 : simpar%do_thermal_region = .FALSE.
391 : END IF
392 :
393 1796 : END SUBROUTINE create_thermal_regions
394 :
395 : ! **************************************************************************************************
396 : !> \brief set particles%els(ipart)%t_region_index to ireg
397 : !> \param particles ...
398 : !> \param ipart ...
399 : !> \param ireg ...
400 : !> \par History
401 : !> - Created (04/2026)
402 : !> \author
403 : ! **************************************************************************************************
404 120 : SUBROUTINE set_t_region_index(particles, ipart, ireg)
405 : TYPE(particle_list_type), POINTER :: particles
406 : INTEGER, INTENT(IN) :: ipart, ireg
407 :
408 : CHARACTER(LEN=default_string_length) :: my_part, my_reg, my_region
409 : INTEGER :: iregion
410 :
411 120 : CALL integer_to_string(ipart, my_part)
412 120 : IF ((ipart > 0) .AND. (ipart <= particles%n_els)) THEN
413 120 : CALL integer_to_string(ireg, my_reg)
414 120 : iregion = particles%els(ipart)%t_region_index
415 120 : IF ((iregion == 0) .OR. (iregion == ireg)) THEN
416 : ! No thermal region has been assigned, or
417 : ! the same one has already been assigned
418 120 : particles%els(ipart)%t_region_index = ireg
419 : ELSE
420 0 : CALL integer_to_string(iregion, my_region)
421 : CALL cp_abort(__LOCATION__, &
422 : "Atom "//TRIM(my_part)//" has been "// &
423 : "assigned to two different thermal "// &
424 : "regions "//TRIM(my_region)//" and "// &
425 0 : TRIM(my_reg)//" which is not allowed!")
426 : END IF
427 : ELSE
428 0 : IF (.NOT. ipart > 0) &
429 : CALL cp_abort(__LOCATION__, &
430 0 : "Input atom index "//TRIM(my_part)//" is non-positive!")
431 0 : IF (.NOT. ipart <= particles%n_els) &
432 : CALL cp_abort(__LOCATION__, &
433 0 : "Input atom index "//TRIM(my_part)//" is out of bounds!")
434 : END IF
435 :
436 120 : END SUBROUTINE set_t_region_index
437 :
438 : ! **************************************************************************************************
439 : !> \brief print_thermal_regions_temperature
440 : !> \param thermal_regions : thermal regions type contains information
441 : !> about the regions
442 : !> \param itimes : iteration number of the time step
443 : !> \param time : simulation time of the time step
444 : !> \param pos : file position
445 : !> \param act : file action
446 : !> \par History
447 : !> - added doxygen header and changed subroutine name from
448 : !> print_thermal_regions to print_thermal_regions_temperature
449 : !> (2014/02/04, LT)
450 : !> \author
451 : ! **************************************************************************************************
452 43875 : SUBROUTINE print_thermal_regions_temperature(thermal_regions, itimes, time, pos, act)
453 : TYPE(thermal_regions_type), POINTER :: thermal_regions
454 : INTEGER, INTENT(IN) :: itimes
455 : REAL(KIND=dp), INTENT(IN) :: time
456 : CHARACTER(LEN=default_string_length) :: pos, act
457 :
458 : CHARACTER(LEN=default_string_length) :: fmd
459 : INTEGER :: ireg, nregions, unit
460 : LOGICAL :: new_file
461 43875 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: temp, temp_tfunc
462 : TYPE(cp_logger_type), POINTER :: logger
463 : TYPE(section_vals_type), POINTER :: print_key
464 :
465 43875 : NULLIFY (logger)
466 43875 : logger => cp_get_default_logger()
467 :
468 43875 : IF (ASSOCIATED(thermal_regions)) THEN
469 110 : print_key => section_vals_get_subs_vals(thermal_regions%section, "PRINT%TEMPERATURE")
470 110 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
471 : unit = cp_print_key_unit_nr(logger, thermal_regions%section, "PRINT%TEMPERATURE", &
472 : extension=".tregion", file_position=pos, &
473 42 : file_action=act, is_new_file=new_file)
474 42 : IF (unit > 0) THEN
475 21 : IF (new_file) THEN
476 1 : WRITE (unit, '(A)') "# Temperature per Region"
477 1 : WRITE (unit, '("#",3X,A,2X,A,13X,A)') "Step Nr.", "Time[fs]", "Temp.[K] ...."
478 : END IF
479 21 : nregions = thermal_regions%nregions
480 63 : ALLOCATE (temp(0:nregions))
481 63 : ALLOCATE (temp_tfunc(1:nregions))
482 84 : temp = 0.0_dp
483 21 : temp(0) = thermal_regions%temp_reg0
484 63 : DO ireg = 1, nregions
485 : ! Note the unit conversion here
486 42 : temp(ireg) = thermal_regions%thermal_region(ireg)%temperature
487 63 : temp_tfunc(ireg) = cp_unit_from_cp2k(thermal_regions%thermal_region(ireg)%temp_expected, "K")
488 : END DO
489 21 : fmd = "(I10,F20.3,"//TRIM(ADJUSTL(cp_to_string(2*nregions + 1)))//"F20.10)"
490 : fmd = TRIM(fmd)
491 63 : WRITE (UNIT=unit, FMT=fmd) itimes, time, temp(0), (temp_tfunc(ireg), temp(ireg), ireg=1, nregions)
492 21 : DEALLOCATE (temp)
493 21 : DEALLOCATE (temp_tfunc)
494 : END IF
495 42 : CALL cp_print_key_finished_output(unit, logger, thermal_regions%section, "PRINT%TEMPERATURE")
496 : END IF
497 : END IF
498 43875 : END SUBROUTINE print_thermal_regions_temperature
499 :
500 : ! **************************************************************************************************
501 : !> \brief print out information regarding to langevin regions defined in
502 : !> thermal_regions section
503 : !> \param thermal_regions : thermal regions type containing the relevant
504 : !> langevin regions information
505 : !> \param simpar : wrapper for simulation parameters
506 : !> \param pos : file position
507 : !> \param act : file action
508 : !> \par History
509 : !> - created (2014/02/02, LT)
510 : !> \author Lianheng Tong [LT] (tonglianheng@gmail.com)
511 : ! **************************************************************************************************
512 42 : SUBROUTINE print_thermal_regions_langevin(thermal_regions, simpar, pos, act)
513 : TYPE(thermal_regions_type), POINTER :: thermal_regions
514 : TYPE(simpar_type), POINTER :: simpar
515 : CHARACTER(LEN=default_string_length) :: pos, act
516 :
517 : INTEGER :: ipart, ipart_reg, ireg, natoms, &
518 : print_unit
519 42 : INTEGER, ALLOCATABLE, DIMENSION(:) :: region_id
520 : LOGICAL :: new_file
521 42 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: noisy_gamma_region, temperature
522 : TYPE(cp_logger_type), POINTER :: logger
523 : TYPE(section_vals_type), POINTER :: print_key
524 :
525 42 : NULLIFY (logger)
526 42 : logger => cp_get_default_logger()
527 :
528 42 : IF (ASSOCIATED(thermal_regions)) THEN
529 12 : IF (ASSOCIATED(thermal_regions%do_langevin)) THEN
530 : print_key => section_vals_get_subs_vals(thermal_regions%section, &
531 12 : "PRINT%LANGEVIN_REGIONS")
532 12 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), &
533 : cp_p_file)) THEN
534 : print_unit = cp_print_key_unit_nr(logger, thermal_regions%section, &
535 : "PRINT%LANGEVIN_REGIONS", &
536 : extension=".lgv_regions", &
537 : file_position=pos, file_action=act, &
538 12 : is_new_file=new_file)
539 12 : IF (print_unit > 0) THEN
540 6 : IF (new_file) THEN
541 6 : WRITE (print_unit, '(A)') "# Atoms Undergoing Langevin MD"
542 : WRITE (print_unit, '(A,3X,A,3X,A,3X,A,3X,A,3X,A)') &
543 6 : "#", "Atom_ID", "Region_ID", "Langevin(L)/NVE(N)", "Expected_T[K]", "[NoisyGamma]"
544 : END IF
545 6 : natoms = SIZE(thermal_regions%do_langevin)
546 18 : ALLOCATE (temperature(natoms))
547 18 : ALLOCATE (region_id(natoms))
548 12 : ALLOCATE (noisy_gamma_region(natoms))
549 42 : temperature(:) = simpar%temp_ext
550 42 : region_id(:) = 0
551 42 : noisy_gamma_region(:) = simpar%noisy_gamma
552 14 : DO ireg = 1, thermal_regions%nregions
553 36 : DO ipart_reg = 1, thermal_regions%thermal_region(ireg)%npart
554 22 : ipart = thermal_regions%thermal_region(ireg)%part_index(ipart_reg)
555 22 : temperature(ipart) = thermal_regions%thermal_region(ireg)%temp_expected
556 22 : region_id(ipart) = thermal_regions%thermal_region(ireg)%region_index
557 30 : noisy_gamma_region(ipart) = thermal_regions%thermal_region(ireg)%noisy_gamma_region
558 : END DO
559 : END DO
560 42 : DO ipart = 1, natoms
561 36 : WRITE (print_unit, '(1X,I10,2X)', advance='no') ipart
562 36 : WRITE (print_unit, '(I10,3X)', advance='no') region_id(ipart)
563 42 : IF (thermal_regions%do_langevin(ipart)) THEN
564 24 : WRITE (print_unit, '(A,3X)', advance='no') "L"
565 24 : IF (noisy_gamma_region(ipart) > 0._dp) THEN
566 10 : WRITE (print_unit, '(10X,F20.3,3X,ES9.3)') temperature(ipart)*kelvin, &
567 20 : noisy_gamma_region(ipart)/femtoseconds
568 : ELSE
569 14 : WRITE (print_unit, '(10X,F20.3)') temperature(ipart)*kelvin
570 : END IF
571 : ELSE
572 12 : WRITE (print_unit, '(A,3X)', advance='no') "N"
573 12 : WRITE (print_unit, '(18X,A)') "--"
574 : END IF
575 : END DO
576 6 : DEALLOCATE (region_id)
577 6 : DEALLOCATE (temperature)
578 6 : DEALLOCATE (noisy_gamma_region)
579 : END IF
580 : CALL cp_print_key_finished_output(print_unit, logger, thermal_regions%section, &
581 12 : "PRINT%LANGEVIN_REGIONS")
582 : END IF
583 : END IF
584 : END IF
585 42 : END SUBROUTINE print_thermal_regions_langevin
586 :
587 : END MODULE thermal_region_utils
|
