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