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_type,&
24 : cp_to_string
25 : USE cp_output_handling, ONLY: cp_p_file,&
26 : cp_print_key_finished_output,&
27 : cp_print_key_should_output,&
28 : cp_print_key_unit_nr
29 : USE cp_subsys_types, ONLY: cp_subsys_get,&
30 : cp_subsys_type
31 : USE cp_units, ONLY: cp_unit_from_cp2k,&
32 : cp_unit_to_cp2k
33 : USE force_env_types, ONLY: force_env_get,&
34 : force_env_type
35 : USE force_fields_util, ONLY: get_generic_info
36 : USE fparser, ONLY: EvalErrType,&
37 : evalf,&
38 : finalizef,&
39 : initf,&
40 : parsef
41 : USE input_constants, ONLY: langevin_ensemble,&
42 : npt_f_ensemble,&
43 : npt_i_ensemble,&
44 : npt_ia_ensemble,&
45 : nvt_ensemble
46 : USE input_section_types, ONLY: section_vals_get,&
47 : section_vals_get_subs_vals,&
48 : section_vals_type,&
49 : section_vals_val_get
50 : USE kinds, ONLY: default_string_length,&
51 : dp
52 : USE memory_utilities, ONLY: reallocate
53 : USE particle_list_types, ONLY: particle_list_type
54 : USE physcon, ONLY: femtoseconds,&
55 : kelvin
56 : USE simpar_types, ONLY: simpar_type
57 : USE string_utilities, ONLY: integer_to_string
58 : USE thermal_region_types, ONLY: allocate_thermal_regions,&
59 : release_thermal_regions,&
60 : thermal_region_type,&
61 : thermal_regions_type
62 : #include "../base/base_uses.f90"
63 :
64 : IMPLICIT NONE
65 :
66 : PRIVATE
67 : PUBLIC :: create_thermal_regions, &
68 : print_thermal_regions_temperature, &
69 : print_thermal_regions_langevin, &
70 : update_thermal_regions_temp_expected
71 :
72 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'thermal_region_utils'
73 :
74 : CONTAINS
75 :
76 : ! **************************************************************************************************
77 : !> \brief create thermal_regions
78 : !> \param thermal_regions ...
79 : !> \param md_section ...
80 : !> \param simpar ...
81 : !> \param force_env ...
82 : !> \par History
83 : !> - Added support for Langevin regions (2014/01/08, LT)
84 : !> \author
85 : ! **************************************************************************************************
86 5364 : SUBROUTINE create_thermal_regions(thermal_regions, md_section, simpar, force_env)
87 : TYPE(thermal_regions_type), POINTER :: thermal_regions
88 : TYPE(section_vals_type), POINTER :: md_section
89 : TYPE(simpar_type), POINTER :: simpar
90 : TYPE(force_env_type), POINTER :: force_env
91 :
92 : CHARACTER(LEN=default_string_length) :: my_region
93 : INTEGER :: i, il, ipart, ireg, itimes, nlist, &
94 : nregions
95 1788 : INTEGER, DIMENSION(:), POINTER :: tmplist
96 : LOGICAL :: apply_thermostat, do_langevin, &
97 : do_langevin_default, do_read_ngr, &
98 : explicit, use_t, use_tfunc
99 : REAL(KIND=dp) :: temp, temp_tol
100 : TYPE(cp_subsys_type), POINTER :: subsys
101 : TYPE(particle_list_type), POINTER :: particles
102 : TYPE(section_vals_type), POINTER :: region_sections, tfunc_section, &
103 : thermal_region_section
104 : TYPE(thermal_region_type), POINTER :: t_region
105 :
106 1788 : NULLIFY (region_sections, t_region, tfunc_section, thermal_region_section, particles, subsys, tmplist)
107 0 : ALLOCATE (thermal_regions)
108 1788 : CALL allocate_thermal_regions(thermal_regions)
109 1788 : thermal_region_section => section_vals_get_subs_vals(md_section, "THERMAL_REGION")
110 1788 : CALL section_vals_get(thermal_region_section, explicit=explicit)
111 1788 : IF (explicit) THEN
112 : apply_thermostat = (simpar%ensemble == nvt_ensemble) .OR. &
113 : (simpar%ensemble == npt_f_ensemble) .OR. &
114 : (simpar%ensemble == npt_ia_ensemble) .OR. &
115 14 : (simpar%ensemble == npt_i_ensemble)
116 : IF (apply_thermostat) THEN
117 : CALL cp_warn(__LOCATION__, &
118 : "With the chosen ensemble the temperature is "// &
119 : "controlled by thermostats. The definition of different thermal "// &
120 0 : "regions might result inconsistent with the presence of thermostats.")
121 : END IF
122 14 : IF (simpar%temp_tol > 0.0_dp) THEN
123 : CALL cp_warn(__LOCATION__, &
124 : "Control of the global temperature by rescaling of the velocity "// &
125 : "is not consistent with the presence of different thermal regions. "// &
126 0 : "The temperature of different regions is rescaled separatedly.")
127 : END IF
128 : CALL section_vals_val_get(thermal_region_section, "FORCE_RESCALING", &
129 14 : l_val=thermal_regions%force_rescaling)
130 : region_sections => section_vals_get_subs_vals(thermal_region_section, &
131 14 : "DEFINE_REGION")
132 14 : CALL section_vals_get(region_sections, n_repetition=nregions)
133 14 : IF (nregions > 0) THEN
134 14 : thermal_regions%nregions = nregions
135 14 : thermal_regions%section => thermal_region_section
136 62 : ALLOCATE (thermal_regions%thermal_region(nregions))
137 14 : CALL force_env_get(force_env, subsys=subsys)
138 14 : CALL cp_subsys_get(subsys, particles=particles)
139 14 : IF (simpar%ensemble == langevin_ensemble) THEN
140 12 : CALL cite_reference(Kantorovich2008)
141 12 : CALL cite_reference(Kantorovich2008a)
142 : CALL section_vals_val_get(thermal_region_section, "DO_LANGEVIN_DEFAULT", &
143 12 : l_val=do_langevin_default)
144 36 : ALLOCATE (thermal_regions%do_langevin(particles%n_els))
145 96 : thermal_regions%do_langevin = do_langevin_default
146 : END IF
147 34 : DO ireg = 1, nregions
148 20 : NULLIFY (t_region)
149 20 : CALL integer_to_string(ireg, my_region)
150 20 : t_region => thermal_regions%thermal_region(ireg)
151 20 : t_region%region_index = ireg
152 :
153 : ! Determine t_region%npart
154 : CALL section_vals_val_get(region_sections, "LIST", &
155 20 : i_rep_section=ireg, n_rep_val=nlist)
156 20 : NULLIFY (t_region%part_index)
157 20 : t_region%npart = 0
158 20 : IF (simpar%ensemble == langevin_ensemble) THEN
159 : CALL section_vals_val_get(region_sections, "DO_LANGEVIN", &
160 16 : i_rep_section=ireg, l_val=do_langevin)
161 : END IF
162 42 : DO il = 1, nlist
163 : CALL section_vals_val_get(region_sections, "LIST", i_rep_section=ireg, &
164 22 : i_rep_val=il, i_vals=tmplist)
165 22 : CALL reallocate(t_region%part_index, 1, t_region%npart + SIZE(tmplist))
166 162 : DO i = 1, SIZE(tmplist)
167 120 : ipart = tmplist(i)
168 120 : CPASSERT(((ipart > 0) .AND. (ipart <= particles%n_els)))
169 120 : t_region%npart = t_region%npart + 1
170 120 : t_region%part_index(t_region%npart) = ipart
171 120 : particles%els(ipart)%t_region_index = ireg
172 142 : IF (simpar%ensemble == langevin_ensemble) THEN
173 44 : thermal_regions%do_langevin(ipart) = do_langevin
174 : END IF
175 : END DO
176 : END DO
177 :
178 : ! Determine t_region%temp_expected
179 20 : tfunc_section => section_vals_get_subs_vals(region_sections, "TFUNC", i_rep_section=ireg)
180 20 : CALL section_vals_val_get(region_sections, "TEMPERATURE", i_rep_section=ireg, explicit=use_t)
181 20 : IF (use_t) CALL section_vals_val_get(region_sections, "TEMPERATURE", i_rep_section=ireg, r_val=temp)
182 20 : CALL section_vals_get(tfunc_section, explicit=use_tfunc)
183 20 : IF (use_tfunc) THEN
184 : CALL get_generic_info(tfunc_section, "FUNCTION", &
185 : t_region%temperature_function, &
186 : t_region%temperature_function_parameters, &
187 : t_region%temperature_function_values, &
188 0 : input_variables=["S", "T"], i_rep_sec=1)
189 0 : CALL initf(1)
190 : CALL parsef(1, TRIM(t_region%temperature_function), &
191 0 : t_region%temperature_function_parameters)
192 0 : IF (.NOT. use_t) THEN ! Evaluate function for initial temperature
193 0 : CALL section_vals_val_get(md_section, "STEP_START_VAL", i_val=itimes)
194 0 : t_region%temperature_function_values(1) = itimes*1.0_dp
195 0 : t_region%temperature_function_values(2) = 0.0_dp
196 0 : temp = cp_unit_to_cp2k(evalf(1, t_region%temperature_function_values), "K")
197 0 : CPASSERT(EvalErrType <= 0)
198 : END IF
199 0 : CALL finalizef()
200 : END IF
201 20 : IF (temp >= 0.0_dp) THEN ! Check for negative temperature in Kelvin
202 20 : t_region%temp_expected = temp
203 : ELSE
204 : CALL cp_abort(__LOCATION__, &
205 : "The initial temperature for region "//TRIM(my_region)//" is "// &
206 : TRIM(ADJUSTL(cp_to_string(cp_unit_from_cp2k(temp, "K"))))// &
207 0 : " K, which is negative and unphysical!")
208 : END IF
209 :
210 : CALL section_vals_val_get(region_sections, "TEMP_TOL", i_rep_section=ireg, &
211 20 : r_val=temp_tol)
212 20 : t_region%temp_tol = temp_tol
213 20 : CALL section_vals_val_get(region_sections, "NOISY_GAMMA_REGION", i_rep_section=ireg, explicit=do_read_ngr)
214 94 : IF (do_read_ngr) THEN
215 : CALL section_vals_val_get(region_sections, "NOISY_GAMMA_REGION", i_rep_section=ireg, &
216 6 : r_val=t_region%noisy_gamma_region)
217 6 : IF (simpar%ensemble == langevin_ensemble) THEN
218 6 : IF (.NOT. do_langevin) THEN
219 : CALL cp_warn(__LOCATION__, &
220 : "You provided NOISY_GAMMA_REGION but atoms in thermal region "//TRIM(my_region)// &
221 : " will not undergo Langevin MD. "// &
222 0 : "NOISY_GAMMA_REGION will be ignored and its value discarded!")
223 : END IF
224 : ELSE
225 : CALL cp_warn(__LOCATION__, &
226 : "You provided NOISY_GAMMA_REGION but the Langevin Ensamble is not selected "// &
227 0 : "NOISY_GAMMA_REGION will be ignored and its value discarded!")
228 : END IF
229 : ELSE
230 14 : t_region%noisy_gamma_region = simpar%noisy_gamma
231 : END IF
232 : END DO
233 14 : simpar%do_thermal_region = .TRUE.
234 : ELSE
235 0 : CALL release_thermal_regions(thermal_regions)
236 0 : DEALLOCATE (thermal_regions)
237 0 : simpar%do_thermal_region = .FALSE.
238 : END IF
239 : ELSE
240 1774 : CALL release_thermal_regions(thermal_regions)
241 1774 : DEALLOCATE (thermal_regions)
242 1774 : simpar%do_thermal_region = .FALSE.
243 : END IF
244 :
245 1788 : END SUBROUTINE create_thermal_regions
246 :
247 : ! **************************************************************************************************
248 : !> \brief update_thermal_regions_temp_expected
249 : !> \param itimes : iteration number of the time step
250 : !> \param thermal_regions : thermal regions type contains information
251 : !> about the regions
252 : !> \date 03.2026
253 : !> \par Note on units: t_region%temperature is updated by
254 : !> SUBROUTINE scale_velocity_region in md_vel_utils.F
255 : !> which is in Kelvin, but t_region%temp_expected is
256 : !> in cp2k internal units. See "Note the unit conversion
257 : !> here" in print_thermal_regions_temperature below.
258 : !> \author HE Zilong
259 : ! **************************************************************************************************
260 96 : SUBROUTINE update_thermal_regions_temp_expected(itimes, thermal_regions)
261 : INTEGER, INTENT(IN) :: itimes
262 : TYPE(thermal_regions_type), POINTER :: thermal_regions
263 :
264 : CHARACTER(LEN=default_string_length) :: my_region
265 : INTEGER :: ireg
266 : TYPE(thermal_region_type), POINTER :: t_region
267 :
268 256 : DO ireg = 1, thermal_regions%nregions
269 160 : NULLIFY (t_region)
270 160 : t_region => thermal_regions%thermal_region(ireg)
271 256 : IF (LEN_TRIM(t_region%temperature_function) > 0) THEN
272 0 : CALL initf(1)
273 0 : CALL parsef(1, TRIM(t_region%temperature_function), t_region%temperature_function_parameters)
274 0 : t_region%temperature_function_values(1) = itimes*1.0_dp
275 0 : t_region%temperature_function_values(2) = t_region%temperature
276 0 : t_region%temp_expected = cp_unit_to_cp2k(evalf(1, t_region%temperature_function_values), "K")
277 0 : CPASSERT(EvalErrType <= 0)
278 0 : CALL finalizef()
279 0 : IF (t_region%temp_expected < 0.0_dp) THEN
280 0 : CALL integer_to_string(ireg, my_region)
281 : CALL cp_abort(__LOCATION__, &
282 : "The temperature evaluated for region "//TRIM(my_region)//" is "// &
283 : TRIM(ADJUSTL(cp_to_string(cp_unit_from_cp2k(t_region%temp_expected, "K"))))// &
284 0 : " K, which is negative and unphysical!")
285 : END IF
286 : END IF
287 : END DO
288 :
289 96 : END SUBROUTINE update_thermal_regions_temp_expected
290 :
291 : ! **************************************************************************************************
292 : !> \brief print_thermal_regions_temperature
293 : !> \param thermal_regions : thermal regions type contains information
294 : !> about the regions
295 : !> \param itimes : iteration number of the time step
296 : !> \param time : simulation time of the time step
297 : !> \param pos : file position
298 : !> \param act : file action
299 : !> \par History
300 : !> - added doxygen header and changed subroutine name from
301 : !> print_thermal_regions to print_thermal_regions_temperature
302 : !> (2014/02/04, LT)
303 : !> \author
304 : ! **************************************************************************************************
305 43819 : SUBROUTINE print_thermal_regions_temperature(thermal_regions, itimes, time, pos, act)
306 : TYPE(thermal_regions_type), POINTER :: thermal_regions
307 : INTEGER, INTENT(IN) :: itimes
308 : REAL(KIND=dp), INTENT(IN) :: time
309 : CHARACTER(LEN=default_string_length) :: pos, act
310 :
311 : CHARACTER(LEN=default_string_length) :: fmd
312 : INTEGER :: ireg, nregions, unit
313 : LOGICAL :: new_file
314 43819 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: temp, temp_tfunc
315 : TYPE(cp_logger_type), POINTER :: logger
316 : TYPE(section_vals_type), POINTER :: print_key
317 :
318 43819 : NULLIFY (logger)
319 43819 : logger => cp_get_default_logger()
320 :
321 43819 : IF (ASSOCIATED(thermal_regions)) THEN
322 110 : print_key => section_vals_get_subs_vals(thermal_regions%section, "PRINT%TEMPERATURE")
323 110 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
324 : unit = cp_print_key_unit_nr(logger, thermal_regions%section, "PRINT%TEMPERATURE", &
325 : extension=".tregion", file_position=pos, &
326 42 : file_action=act, is_new_file=new_file)
327 42 : IF (unit > 0) THEN
328 21 : IF (new_file) THEN
329 1 : WRITE (unit, '(A)') "# Temperature per Region"
330 1 : WRITE (unit, '("#",3X,A,2X,A,13X,A)') "Step Nr.", "Time[fs]", "Temp.[K] ...."
331 : END IF
332 21 : nregions = thermal_regions%nregions
333 63 : ALLOCATE (temp(0:nregions))
334 63 : ALLOCATE (temp_tfunc(1:nregions))
335 84 : temp = 0.0_dp
336 21 : temp(0) = thermal_regions%temp_reg0
337 63 : DO ireg = 1, nregions
338 : ! Note the unit conversion here
339 42 : temp(ireg) = thermal_regions%thermal_region(ireg)%temperature
340 63 : temp_tfunc(ireg) = cp_unit_from_cp2k(thermal_regions%thermal_region(ireg)%temp_expected, "K")
341 : END DO
342 21 : fmd = "(I10,F20.3,"//TRIM(ADJUSTL(cp_to_string(2*nregions + 1)))//"F20.8)"
343 : fmd = TRIM(fmd)
344 63 : WRITE (UNIT=unit, FMT=fmd) itimes, time, temp(0), (temp_tfunc(ireg), temp(ireg), ireg=1, nregions)
345 21 : DEALLOCATE (temp)
346 21 : DEALLOCATE (temp_tfunc)
347 : END IF
348 42 : CALL cp_print_key_finished_output(unit, logger, thermal_regions%section, "PRINT%TEMPERATURE")
349 : END IF
350 : END IF
351 43819 : END SUBROUTINE print_thermal_regions_temperature
352 :
353 : ! **************************************************************************************************
354 : !> \brief print out information regarding to langevin regions defined in
355 : !> thermal_regions section
356 : !> \param thermal_regions : thermal regions type containing the relevant
357 : !> langevin regions information
358 : !> \param simpar : wrapper for simulation parameters
359 : !> \param pos : file position
360 : !> \param act : file action
361 : !> \par History
362 : !> - created (2014/02/02, LT)
363 : !> \author Lianheng Tong [LT] (tonglianheng@gmail.com)
364 : ! **************************************************************************************************
365 42 : SUBROUTINE print_thermal_regions_langevin(thermal_regions, simpar, pos, act)
366 : TYPE(thermal_regions_type), POINTER :: thermal_regions
367 : TYPE(simpar_type), POINTER :: simpar
368 : CHARACTER(LEN=default_string_length) :: pos, act
369 :
370 : INTEGER :: ipart, ipart_reg, ireg, natoms, &
371 : print_unit
372 42 : INTEGER, ALLOCATABLE, DIMENSION(:) :: region_id
373 : LOGICAL :: new_file
374 42 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: noisy_gamma_region, temperature
375 : TYPE(cp_logger_type), POINTER :: logger
376 : TYPE(section_vals_type), POINTER :: print_key
377 :
378 42 : NULLIFY (logger)
379 42 : logger => cp_get_default_logger()
380 :
381 42 : IF (ASSOCIATED(thermal_regions)) THEN
382 12 : IF (ASSOCIATED(thermal_regions%do_langevin)) THEN
383 : print_key => section_vals_get_subs_vals(thermal_regions%section, &
384 12 : "PRINT%LANGEVIN_REGIONS")
385 12 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), &
386 : cp_p_file)) THEN
387 : print_unit = cp_print_key_unit_nr(logger, thermal_regions%section, &
388 : "PRINT%LANGEVIN_REGIONS", &
389 : extension=".lgv_regions", &
390 : file_position=pos, file_action=act, &
391 12 : is_new_file=new_file)
392 12 : IF (print_unit > 0) THEN
393 6 : IF (new_file) THEN
394 6 : WRITE (print_unit, '(A)') "# Atoms Undergoing Langevin MD"
395 : WRITE (print_unit, '(A,3X,A,3X,A,3X,A,3X,A,3X,A)') &
396 6 : "#", "Atom_ID", "Region_ID", "Langevin(L)/NVE(N)", "Expected_T[K]", "[NoisyGamma]"
397 : END IF
398 6 : natoms = SIZE(thermal_regions%do_langevin)
399 18 : ALLOCATE (temperature(natoms))
400 18 : ALLOCATE (region_id(natoms))
401 12 : ALLOCATE (noisy_gamma_region(natoms))
402 42 : temperature(:) = simpar%temp_ext
403 42 : region_id(:) = 0
404 42 : noisy_gamma_region(:) = simpar%noisy_gamma
405 14 : DO ireg = 1, thermal_regions%nregions
406 36 : DO ipart_reg = 1, thermal_regions%thermal_region(ireg)%npart
407 22 : ipart = thermal_regions%thermal_region(ireg)%part_index(ipart_reg)
408 22 : temperature(ipart) = thermal_regions%thermal_region(ireg)%temp_expected
409 22 : region_id(ipart) = thermal_regions%thermal_region(ireg)%region_index
410 30 : noisy_gamma_region(ipart) = thermal_regions%thermal_region(ireg)%noisy_gamma_region
411 : END DO
412 : END DO
413 42 : DO ipart = 1, natoms
414 36 : WRITE (print_unit, '(1X,I10,2X)', advance='no') ipart
415 36 : WRITE (print_unit, '(I10,3X)', advance='no') region_id(ipart)
416 42 : IF (thermal_regions%do_langevin(ipart)) THEN
417 24 : WRITE (print_unit, '(A,3X)', advance='no') "L"
418 24 : IF (noisy_gamma_region(ipart) > 0._dp) THEN
419 10 : WRITE (print_unit, '(10X,F20.3,3X,ES9.3)') temperature(ipart)*kelvin, &
420 20 : noisy_gamma_region(ipart)/femtoseconds
421 : ELSE
422 14 : WRITE (print_unit, '(10X,F20.3)') temperature(ipart)*kelvin
423 : END IF
424 : ELSE
425 12 : WRITE (print_unit, '(A,3X)', advance='no') "N"
426 12 : WRITE (print_unit, '(18X,A)') "--"
427 : END IF
428 : END DO
429 6 : DEALLOCATE (region_id)
430 6 : DEALLOCATE (temperature)
431 6 : DEALLOCATE (noisy_gamma_region)
432 : END IF
433 : CALL cp_print_key_finished_output(print_unit, logger, thermal_regions%section, &
434 12 : "PRINT%LANGEVIN_REGIONS")
435 : END IF
436 : END IF
437 : END IF
438 42 : END SUBROUTINE print_thermal_regions_langevin
439 :
440 : END MODULE thermal_region_utils
|
