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 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 force_env_types, ONLY: force_env_get,&
32 : force_env_type
33 : USE input_constants, ONLY: langevin_ensemble,&
34 : npt_f_ensemble,&
35 : npt_i_ensemble,&
36 : npt_ia_ensemble,&
37 : nvt_ensemble
38 : USE input_section_types, ONLY: section_vals_get,&
39 : section_vals_get_subs_vals,&
40 : section_vals_type,&
41 : section_vals_val_get
42 : USE kinds, ONLY: default_string_length,&
43 : dp
44 : USE memory_utilities, ONLY: reallocate
45 : USE particle_list_types, ONLY: particle_list_type
46 : USE physcon, ONLY: femtoseconds,&
47 : kelvin
48 : USE simpar_types, ONLY: simpar_type
49 : USE string_utilities, ONLY: integer_to_string
50 : USE thermal_region_types, ONLY: allocate_thermal_regions,&
51 : release_thermal_regions,&
52 : thermal_region_type,&
53 : thermal_regions_type
54 : #include "../base/base_uses.f90"
55 :
56 : IMPLICIT NONE
57 :
58 : PRIVATE
59 : PUBLIC :: create_thermal_regions, &
60 : print_thermal_regions_temperature, &
61 : print_thermal_regions_langevin
62 :
63 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'thermal_region_utils'
64 :
65 : CONTAINS
66 :
67 : ! **************************************************************************************************
68 : !> \brief create thermal_regions
69 : !> \param thermal_regions ...
70 : !> \param md_section ...
71 : !> \param simpar ...
72 : !> \param force_env ...
73 : !> \par History
74 : !> - Added support for Langevin regions (2014/01/08, LT)
75 : !> \author
76 : ! **************************************************************************************************
77 5307 : SUBROUTINE create_thermal_regions(thermal_regions, md_section, simpar, force_env)
78 : TYPE(thermal_regions_type), POINTER :: thermal_regions
79 : TYPE(section_vals_type), POINTER :: md_section
80 : TYPE(simpar_type), POINTER :: simpar
81 : TYPE(force_env_type), POINTER :: force_env
82 :
83 : CHARACTER(LEN=default_string_length) :: my_region
84 : INTEGER :: i, il, ipart, ireg, nlist, nregions
85 1769 : INTEGER, DIMENSION(:), POINTER :: tmplist
86 : LOGICAL :: apply_thermostat, do_langevin, &
87 : do_langevin_default, do_read_ngr, &
88 : explicit
89 : REAL(KIND=dp) :: temp, temp_tol
90 : TYPE(cp_subsys_type), POINTER :: subsys
91 : TYPE(particle_list_type), POINTER :: particles
92 : TYPE(section_vals_type), POINTER :: region_sections, thermal_region_section
93 : TYPE(thermal_region_type), POINTER :: t_region
94 :
95 1769 : NULLIFY (region_sections, t_region, thermal_region_section, particles, subsys, tmplist)
96 0 : ALLOCATE (thermal_regions)
97 1769 : CALL allocate_thermal_regions(thermal_regions)
98 1769 : thermal_region_section => section_vals_get_subs_vals(md_section, "THERMAL_REGION")
99 1769 : CALL section_vals_get(thermal_region_section, explicit=explicit)
100 1769 : IF (explicit) THEN
101 : apply_thermostat = (simpar%ensemble == nvt_ensemble) .OR. &
102 : (simpar%ensemble == npt_f_ensemble) .OR. &
103 : (simpar%ensemble == npt_ia_ensemble) .OR. &
104 14 : (simpar%ensemble == npt_i_ensemble)
105 : IF (apply_thermostat) THEN
106 : CALL cp_warn(__LOCATION__, &
107 : "With the chosen ensemble the temperature is "// &
108 : "controlled by thermostats. The definition of different thermal "// &
109 0 : "regions might result inconsistent with the presence of thermostats.")
110 : END IF
111 14 : IF (simpar%temp_tol > 0.0_dp) THEN
112 : CALL cp_warn(__LOCATION__, &
113 : "Control of the global temperature by rescaling of the velocity "// &
114 : "is not consistent with the presence of different thermal regions. "// &
115 0 : "The temperature of different regions is rescaled separatedly.")
116 : END IF
117 : CALL section_vals_val_get(thermal_region_section, "FORCE_RESCALING", &
118 14 : l_val=thermal_regions%force_rescaling)
119 : region_sections => section_vals_get_subs_vals(thermal_region_section, &
120 14 : "DEFINE_REGION")
121 14 : CALL section_vals_get(region_sections, n_repetition=nregions)
122 14 : IF (nregions > 0) THEN
123 14 : thermal_regions%nregions = nregions
124 14 : thermal_regions%section => thermal_region_section
125 62 : ALLOCATE (thermal_regions%thermal_region(nregions))
126 14 : CALL force_env_get(force_env, subsys=subsys)
127 14 : CALL cp_subsys_get(subsys, particles=particles)
128 14 : IF (simpar%ensemble == langevin_ensemble) THEN
129 12 : CALL cite_reference(Kantorovich2008)
130 12 : CALL cite_reference(Kantorovich2008a)
131 : CALL section_vals_val_get(thermal_region_section, "DO_LANGEVIN_DEFAULT", &
132 12 : l_val=do_langevin_default)
133 36 : ALLOCATE (thermal_regions%do_langevin(particles%n_els))
134 96 : thermal_regions%do_langevin = do_langevin_default
135 : END IF
136 34 : DO ireg = 1, nregions
137 20 : NULLIFY (t_region)
138 20 : t_region => thermal_regions%thermal_region(ireg)
139 20 : t_region%region_index = ireg
140 : CALL section_vals_val_get(region_sections, "LIST", &
141 20 : i_rep_section=ireg, n_rep_val=nlist)
142 20 : NULLIFY (t_region%part_index)
143 20 : t_region%npart = 0
144 20 : IF (simpar%ensemble == langevin_ensemble) THEN
145 : CALL section_vals_val_get(region_sections, "DO_LANGEVIN", &
146 16 : i_rep_section=ireg, l_val=do_langevin)
147 : END IF
148 42 : DO il = 1, nlist
149 : CALL section_vals_val_get(region_sections, "LIST", i_rep_section=ireg, &
150 22 : i_rep_val=il, i_vals=tmplist)
151 22 : CALL reallocate(t_region%part_index, 1, t_region%npart + SIZE(tmplist))
152 162 : DO i = 1, SIZE(tmplist)
153 120 : ipart = tmplist(i)
154 120 : CPASSERT(((ipart > 0) .AND. (ipart <= particles%n_els)))
155 120 : t_region%npart = t_region%npart + 1
156 120 : t_region%part_index(t_region%npart) = ipart
157 120 : particles%els(ipart)%t_region_index = ireg
158 142 : IF (simpar%ensemble == langevin_ensemble) THEN
159 44 : thermal_regions%do_langevin(ipart) = do_langevin
160 : END IF
161 : END DO
162 : END DO
163 : CALL section_vals_val_get(region_sections, "TEMPERATURE", i_rep_section=ireg, &
164 20 : r_val=temp)
165 20 : t_region%temp_expected = temp
166 : CALL section_vals_val_get(region_sections, "TEMP_TOL", i_rep_section=ireg, &
167 20 : r_val=temp_tol)
168 20 : t_region%temp_tol = temp_tol
169 20 : CALL section_vals_val_get(region_sections, "NOISY_GAMMA_REGION", i_rep_section=ireg, explicit=do_read_ngr)
170 54 : IF (do_read_ngr) THEN
171 : CALL section_vals_val_get(region_sections, "NOISY_GAMMA_REGION", i_rep_section=ireg, &
172 6 : r_val=t_region%noisy_gamma_region)
173 6 : IF (simpar%ensemble == langevin_ensemble) THEN
174 6 : IF (.NOT. do_langevin) THEN
175 0 : CALL integer_to_string(ireg, my_region)
176 : CALL cp_warn(__LOCATION__, &
177 : "You provided NOISY_GAMMA_REGION but atoms in thermal region "//TRIM(my_region)// &
178 : " will not undergo Langevin MD. "// &
179 0 : "NOISY_GAMMA_REGION will be ignored and its value discarded!")
180 : END IF
181 : ELSE
182 : CALL cp_warn(__LOCATION__, &
183 : "You provided NOISY_GAMMA_REGION but the Langevin Ensamble is not selected "// &
184 0 : "NOISY_GAMMA_REGION will be ignored and its value discarded!")
185 : END IF
186 : ELSE
187 14 : t_region%noisy_gamma_region = simpar%noisy_gamma
188 : END IF
189 : END DO
190 14 : simpar%do_thermal_region = .TRUE.
191 : ELSE
192 0 : CALL release_thermal_regions(thermal_regions)
193 0 : DEALLOCATE (thermal_regions)
194 0 : simpar%do_thermal_region = .FALSE.
195 : END IF
196 : ELSE
197 1755 : CALL release_thermal_regions(thermal_regions)
198 1755 : DEALLOCATE (thermal_regions)
199 1755 : simpar%do_thermal_region = .FALSE.
200 : END IF
201 :
202 1769 : END SUBROUTINE create_thermal_regions
203 :
204 : ! **************************************************************************************************
205 : !> \brief print_thermal_regions_temperature
206 : !> \param thermal_regions : thermal regions type contains information
207 : !> about the regions
208 : !> \param itimes : iteration number of the time step
209 : !> \param time : simulation time of the time step
210 : !> \param pos : file position
211 : !> \param act : file action
212 : !> \par History
213 : !> - added doxygen header and changed subroutine name from
214 : !> print_thermal_regions to print_thermal_regions_temperature
215 : !> (2014/02/04, LT)
216 : !> \author
217 : ! **************************************************************************************************
218 42149 : SUBROUTINE print_thermal_regions_temperature(thermal_regions, itimes, time, pos, act)
219 : TYPE(thermal_regions_type), POINTER :: thermal_regions
220 : INTEGER, INTENT(IN) :: itimes
221 : REAL(KIND=dp), INTENT(IN) :: time
222 : CHARACTER(LEN=default_string_length) :: pos, act
223 :
224 : CHARACTER(LEN=default_string_length) :: fmd
225 : INTEGER :: ireg, nregions, unit
226 : LOGICAL :: new_file
227 42149 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: temp
228 : TYPE(cp_logger_type), POINTER :: logger
229 : TYPE(section_vals_type), POINTER :: print_key
230 :
231 42149 : NULLIFY (logger)
232 42149 : logger => cp_get_default_logger()
233 :
234 42149 : IF (ASSOCIATED(thermal_regions)) THEN
235 110 : print_key => section_vals_get_subs_vals(thermal_regions%section, "PRINT%TEMPERATURE")
236 110 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
237 : unit = cp_print_key_unit_nr(logger, thermal_regions%section, "PRINT%TEMPERATURE", &
238 : extension=".tregion", file_position=pos, &
239 42 : file_action=act, is_new_file=new_file)
240 42 : IF (unit > 0) THEN
241 21 : IF (new_file) THEN
242 1 : WRITE (unit, '(A)') "# Temperature per Region"
243 1 : WRITE (unit, '("#",3X,A,2X,A,13X,A)') "Step Nr.", "Time[fs]", "Temp.[K] ...."
244 : END IF
245 21 : nregions = thermal_regions%nregions
246 63 : ALLOCATE (temp(0:nregions))
247 84 : temp = 0.0_dp
248 21 : temp(0) = thermal_regions%temp_reg0
249 63 : DO ireg = 1, nregions
250 63 : temp(ireg) = thermal_regions%thermal_region(ireg)%temperature
251 : END DO
252 21 : fmd = "(I10,F20.3,"//TRIM(ADJUSTL(cp_to_string(nregions + 1)))//"F20.6)"
253 : fmd = TRIM(fmd)
254 21 : WRITE (UNIT=unit, FMT=fmd) itimes, time, temp(0:nregions)
255 21 : DEALLOCATE (temp)
256 : END IF
257 42 : CALL cp_print_key_finished_output(unit, logger, thermal_regions%section, "PRINT%TEMPERATURE")
258 : END IF
259 : END IF
260 42149 : END SUBROUTINE print_thermal_regions_temperature
261 :
262 : ! **************************************************************************************************
263 : !> \brief print out information regarding to langevin regions defined in
264 : !> thermal_regions section
265 : !> \param thermal_regions : thermal regions type containing the relevant
266 : !> langevin regions information
267 : !> \param simpar : wrapper for simulation parameters
268 : !> \param pos : file position
269 : !> \param act : file action
270 : !> \par History
271 : !> - created (2014/02/02, LT)
272 : !> \author Lianheng Tong [LT] (tonglianheng@gmail.com)
273 : ! **************************************************************************************************
274 42 : SUBROUTINE print_thermal_regions_langevin(thermal_regions, simpar, pos, act)
275 : TYPE(thermal_regions_type), POINTER :: thermal_regions
276 : TYPE(simpar_type), POINTER :: simpar
277 : CHARACTER(LEN=default_string_length) :: pos, act
278 :
279 : INTEGER :: ipart, ipart_reg, ireg, natoms, &
280 : print_unit
281 42 : INTEGER, ALLOCATABLE, DIMENSION(:) :: region_id
282 : LOGICAL :: new_file
283 42 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: noisy_gamma_region, temperature
284 : TYPE(cp_logger_type), POINTER :: logger
285 : TYPE(section_vals_type), POINTER :: print_key
286 :
287 42 : NULLIFY (logger)
288 42 : logger => cp_get_default_logger()
289 :
290 42 : IF (ASSOCIATED(thermal_regions)) THEN
291 12 : IF (ASSOCIATED(thermal_regions%do_langevin)) THEN
292 : print_key => section_vals_get_subs_vals(thermal_regions%section, &
293 12 : "PRINT%LANGEVIN_REGIONS")
294 12 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), &
295 : cp_p_file)) THEN
296 : print_unit = cp_print_key_unit_nr(logger, thermal_regions%section, &
297 : "PRINT%LANGEVIN_REGIONS", &
298 : extension=".lgv_regions", &
299 : file_position=pos, file_action=act, &
300 12 : is_new_file=new_file)
301 12 : IF (print_unit > 0) THEN
302 6 : IF (new_file) THEN
303 6 : WRITE (print_unit, '(A)') "# Atoms Undergoing Langevin MD"
304 : WRITE (print_unit, '(A,3X,A,3X,A,3X,A,3X,A,3X,A)') &
305 6 : "#", "Atom_ID", "Region_ID", "Langevin(L)/NVE(N)", "Expected_T[K]", "[NoisyGamma]"
306 : END IF
307 6 : natoms = SIZE(thermal_regions%do_langevin)
308 18 : ALLOCATE (temperature(natoms))
309 18 : ALLOCATE (region_id(natoms))
310 12 : ALLOCATE (noisy_gamma_region(natoms))
311 42 : temperature(:) = simpar%temp_ext
312 42 : region_id(:) = 0
313 42 : noisy_gamma_region(:) = simpar%noisy_gamma
314 14 : DO ireg = 1, thermal_regions%nregions
315 36 : DO ipart_reg = 1, thermal_regions%thermal_region(ireg)%npart
316 22 : ipart = thermal_regions%thermal_region(ireg)%part_index(ipart_reg)
317 22 : temperature(ipart) = thermal_regions%thermal_region(ireg)%temp_expected
318 22 : region_id(ipart) = thermal_regions%thermal_region(ireg)%region_index
319 30 : noisy_gamma_region(ipart) = thermal_regions%thermal_region(ireg)%noisy_gamma_region
320 : END DO
321 : END DO
322 42 : DO ipart = 1, natoms
323 36 : WRITE (print_unit, '(1X,I10,2X)', advance='no') ipart
324 36 : WRITE (print_unit, '(I10,3X)', advance='no') region_id(ipart)
325 42 : IF (thermal_regions%do_langevin(ipart)) THEN
326 24 : WRITE (print_unit, '(A,3X)', advance='no') "L"
327 24 : IF (noisy_gamma_region(ipart) > 0._dp) THEN
328 10 : WRITE (print_unit, '(10X,F20.3,3X,ES9.3)') temperature(ipart)*kelvin, &
329 20 : noisy_gamma_region(ipart)/femtoseconds
330 : ELSE
331 14 : WRITE (print_unit, '(10X,F20.3)') temperature(ipart)*kelvin
332 : END IF
333 : ELSE
334 12 : WRITE (print_unit, '(A,3X)', advance='no') "N"
335 12 : WRITE (print_unit, '(18X,A)') "--"
336 : END IF
337 : END DO
338 6 : DEALLOCATE (region_id)
339 6 : DEALLOCATE (temperature)
340 6 : DEALLOCATE (noisy_gamma_region)
341 : END IF
342 : CALL cp_print_key_finished_output(print_unit, logger, thermal_regions%section, &
343 12 : "PRINT%LANGEVIN_REGIONS")
344 : END IF
345 : END IF
346 : END IF
347 42 : END SUBROUTINE print_thermal_regions_langevin
348 :
349 : END MODULE thermal_region_utils
|
