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 contains some general routines for dealing with the restart
10 : !> files and creating force_env for MC use
11 : !> \par History
12 : !> none
13 : !> \author MJM
14 : ! **************************************************************************************************
15 : MODULE mc_control
16 : USE atomic_kind_types, ONLY: get_atomic_kind
17 : USE cell_types, ONLY: cell_type,&
18 : get_cell
19 : USE cp_files, ONLY: close_file,&
20 : open_file
21 : USE cp_log_handling, ONLY: cp_logger_get_default_unit_nr
22 : USE cp_subsys_types, ONLY: cp_subsys_get,&
23 : cp_subsys_type
24 : USE f77_interface, ONLY: create_force_env,&
25 : destroy_force_env,&
26 : f_env_add_defaults,&
27 : f_env_rm_defaults,&
28 : f_env_type
29 : USE force_env_types, ONLY: force_env_get,&
30 : force_env_retain,&
31 : force_env_type
32 : USE global_types, ONLY: global_environment_type
33 : USE input_section_types, ONLY: section_type
34 : USE kinds, ONLY: default_path_length,&
35 : default_string_length,&
36 : dp
37 : USE mc_misc, ONLY: mc_make_dat_file_new
38 : USE mc_types, ONLY: get_mc_molecule_info,&
39 : get_mc_par,&
40 : mc_input_file_type,&
41 : mc_molecule_info_type,&
42 : mc_simpar_type,&
43 : set_mc_par
44 : USE message_passing, ONLY: mp_comm_type,&
45 : mp_para_env_type
46 : USE molecule_kind_list_types, ONLY: molecule_kind_list_type
47 : USE molecule_kind_types, ONLY: atom_type,&
48 : get_molecule_kind,&
49 : molecule_kind_type
50 : USE parallel_rng_types, ONLY: rng_stream_type
51 : USE particle_list_types, ONLY: particle_list_type
52 : USE physcon, ONLY: angstrom
53 : #include "../../base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 :
57 : PRIVATE
58 : ! *** Global parameters ***
59 :
60 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_control'
61 :
62 : PUBLIC :: write_mc_restart, read_mc_restart, mc_create_force_env, &
63 : mc_create_bias_force_env
64 :
65 : CONTAINS
66 :
67 : ! **************************************************************************************************
68 : !> \brief writes the coordinates of the current step to a file that can
69 : !> be read in at the start of the next simulation
70 : !> \param nnstep how many steps the simulation has run
71 : !>
72 : !> Only use in serial.
73 : !> \param mc_par the mc parameters for the force env
74 : !> \param nchains ...
75 : !> \param force_env the force environment to write the coords from
76 : !> \author MJM
77 : ! **************************************************************************************************
78 30 : SUBROUTINE write_mc_restart(nnstep, mc_par, nchains, force_env)
79 :
80 : INTEGER, INTENT(IN) :: nnstep
81 : TYPE(mc_simpar_type), POINTER :: mc_par
82 : INTEGER, DIMENSION(:), INTENT(IN) :: nchains
83 : TYPE(force_env_type), POINTER :: force_env
84 :
85 : CHARACTER(len=*), PARAMETER :: routineN = 'write_mc_restart'
86 :
87 : CHARACTER(LEN=20) :: ensemble
88 : CHARACTER(LEN=default_path_length) :: restart_file_name
89 : CHARACTER(LEN=default_string_length) :: name
90 : INTEGER :: handle, ichain, imol_type, iparticle, &
91 : iunit, natom, nmol_types, nmolecule, &
92 : nunits_tot, unit
93 : REAL(KIND=dp) :: temperature
94 : REAL(KIND=dp), DIMENSION(1:3) :: abc
95 30 : TYPE(atom_type), DIMENSION(:), POINTER :: atom_list
96 : TYPE(cell_type), POINTER :: cell
97 : TYPE(cp_subsys_type), POINTER :: subsys
98 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
99 : TYPE(molecule_kind_type), POINTER :: molecule_kind
100 : TYPE(particle_list_type), POINTER :: particles
101 :
102 30 : CALL timeset(routineN, handle)
103 :
104 : ! get some data from mc_par
105 : CALL get_mc_par(mc_par, restart_file_name=restart_file_name, temperature=temperature, &
106 30 : ensemble=ensemble)
107 :
108 : ! open the file and write some simulation parameters
109 : CALL open_file(file_name=restart_file_name, unit_number=unit, &
110 30 : file_action='WRITE', file_status='REPLACE')
111 :
112 : ! get the cell length and coordinates
113 30 : CALL force_env_get(force_env, cell=cell, subsys=subsys)
114 30 : CALL get_cell(cell, abc=abc)
115 : CALL cp_subsys_get(subsys, &
116 : molecule_kinds=molecule_kinds, &
117 30 : particles=particles)
118 :
119 30 : nunits_tot = SIZE(particles%els(:))
120 79 : IF (SUM(nchains(:)) == 0) nunits_tot = 0
121 30 : WRITE (unit, *) nnstep
122 30 : WRITE (unit, *) temperature, nunits_tot
123 30 : WRITE (unit, *) ensemble
124 30 : WRITE (unit, *) nchains(:)
125 120 : WRITE (unit, '(3(F10.6,3X))') abc(1:3)*angstrom ! in angstroms
126 30 : WRITE (unit, *)
127 :
128 : ! can't do a simple particles%els%atomic_kind%element_symbol because
129 : ! of the classical force_env
130 30 : IF (nunits_tot > 0) THEN
131 30 : nmol_types = SIZE(molecule_kinds%els(:))
132 30 : iparticle = 1
133 79 : DO imol_type = 1, nmol_types
134 49 : molecule_kind => molecule_kinds%els(imol_type)
135 : CALL get_molecule_kind(molecule_kind, atom_list=atom_list, &
136 49 : nmolecule=nmolecule, natom=natom)
137 : ! write the coordinates out
138 3602 : DO ichain = 1, nmolecule
139 7391 : DO iunit = 1, natom
140 3819 : CALL get_atomic_kind(atom_list(iunit)%atomic_kind, name=name)
141 : WRITE (unit, '(1X,A,1X,3(F15.10,3X))') &
142 3819 : TRIM(ADJUSTL(name)), &
143 19095 : particles%els(iparticle)%r(1:3)*angstrom
144 7342 : iparticle = iparticle + 1
145 : END DO
146 : END DO
147 : END DO
148 : END IF
149 :
150 30 : CALL close_file(unit_number=unit)
151 :
152 : ! end the timing
153 30 : CALL timestop(handle)
154 :
155 30 : END SUBROUTINE write_mc_restart
156 :
157 : ! **************************************************************************************************
158 : !> \brief reads the input coordinates of the simulation from a file written above
159 : !> \param mc_par the mc parameters for the force env
160 : !> \param force_env the force environment to write the coords from
161 : !> \param iw the unit to write an error message to, in case current
162 : !> simulation parameters don't match what's in the restart file
163 : !> \param mc_nunits_tot ...
164 : !> \param rng_stream the stream we pull random numbers from
165 : !>
166 : !> Used in parallel.
167 : !> \author MJM
168 : ! **************************************************************************************************
169 2 : SUBROUTINE read_mc_restart(mc_par, force_env, iw, mc_nunits_tot, rng_stream)
170 :
171 : TYPE(mc_simpar_type), POINTER :: mc_par
172 : TYPE(force_env_type), POINTER :: force_env
173 : INTEGER, INTENT(IN) :: iw
174 : INTEGER, INTENT(INOUT) :: mc_nunits_tot
175 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
176 :
177 : CHARACTER(len=*), PARAMETER :: routineN = 'read_mc_restart'
178 :
179 : CHARACTER(default_string_length), &
180 2 : DIMENSION(:, :), POINTER :: atom_names
181 : CHARACTER(LEN=20) :: ensemble, mc_ensemble
182 2 : CHARACTER(LEN=5), ALLOCATABLE, DIMENSION(:) :: atom_symbols
183 : CHARACTER(LEN=default_path_length) :: dat_file, restart_file_name
184 : INTEGER :: handle, i, ipart, iunit, nmol_types, &
185 : nstart, nunits_tot, print_level, &
186 : source, unit
187 2 : INTEGER, DIMENSION(:), POINTER :: nchains, nunits
188 : LOGICAL :: ionode
189 : REAL(KIND=dp) :: mc_temp, rand, temperature
190 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r
191 : REAL(KIND=dp), DIMENSION(1:3) :: abc, box_length
192 : TYPE(cell_type), POINTER :: cell
193 : TYPE(cp_subsys_type), POINTER :: subsys
194 : TYPE(mc_input_file_type), POINTER :: mc_input_file
195 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
196 : TYPE(mp_comm_type) :: group
197 : TYPE(particle_list_type), POINTER :: particles
198 :
199 2 : CALL timeset(routineN, handle)
200 :
201 : ! get some stuff from the mc_par
202 : CALL get_mc_par(mc_par, restart_file_name=restart_file_name, temperature=mc_temp, &
203 : ensemble=mc_ensemble, mc_molecule_info=mc_molecule_info, &
204 : ionode=ionode, dat_file=dat_file, &
205 2 : group=group, source=source, mc_input_file=mc_input_file)
206 : CALL get_mc_molecule_info(mc_molecule_info, nunits=nunits, &
207 2 : nmol_types=nmol_types, atom_names=atom_names)
208 :
209 6 : ALLOCATE (nchains(1:nmol_types))
210 :
211 : ! currently a hack, printlevel should be intern to the print_keys
212 2 : print_level = 1
213 :
214 2 : IF (ionode) THEN
215 : ! open the file and read some simulation parameters
216 : CALL open_file(file_name=restart_file_name, unit_number=unit, &
217 1 : file_action='READ', file_status='OLD')
218 :
219 1 : READ (unit, *) nstart
220 1 : READ (unit, *) temperature, nunits_tot
221 1 : READ (unit, *) ensemble
222 1 : READ (unit, *) nchains(1:nmol_types)
223 : END IF
224 2 : CALL group%bcast(nstart, source)
225 2 : CALL group%bcast(temperature, source)
226 2 : CALL group%bcast(nunits_tot, source)
227 2 : CALL group%bcast(ensemble, source)
228 6 : CALL group%bcast(nchains, source)
229 :
230 : ! do some checking
231 2 : IF (ABS(temperature - mc_temp) > 0.01E0_dp) THEN
232 0 : IF (ionode) THEN
233 0 : WRITE (iw, *) 'The temperature in the restart file is ', &
234 0 : 'not the same as the input file.'
235 0 : WRITE (iw, *) 'Input file temperature =', mc_temp
236 0 : WRITE (iw, *) 'Restart file temperature =', temperature
237 : END IF
238 0 : CPABORT("Temperature difference between restart and input")
239 : END IF
240 2 : IF (nunits_tot /= mc_nunits_tot) THEN
241 0 : IF (ionode) THEN
242 0 : WRITE (iw, *) 'The total number of units in the restart ', &
243 0 : 'file is not the same as the input file.'
244 0 : WRITE (iw, *) 'Input file units =', mc_nunits_tot
245 0 : WRITE (iw, *) 'Restart file units =', nunits_tot
246 : END IF
247 0 : mc_nunits_tot = nunits_tot
248 : END IF
249 2 : IF (ensemble /= mc_ensemble) THEN
250 0 : IF (ionode) THEN
251 0 : WRITE (iw, *) 'The ensemble in the restart file is ', &
252 0 : 'not the same as the input file.'
253 0 : WRITE (iw, *) 'Input file ensemble =', mc_ensemble
254 0 : WRITE (iw, *) 'Restart file ensemble =', ensemble
255 : END IF
256 0 : CPABORT("Ensembles different between restart and input")
257 : END IF
258 :
259 : ! get the cell length and coordinates
260 2 : CALL force_env_get(force_env, cell=cell, subsys=subsys)
261 2 : CALL get_cell(cell, abc=abc)
262 : CALL cp_subsys_get(subsys, &
263 2 : particles=particles)
264 :
265 2 : IF (ionode) THEN
266 1 : READ (unit, *) box_length(1:3) ! in angstroms
267 1 : READ (unit, *)
268 4 : box_length(1:3) = box_length(1:3)/angstrom ! convert to a.u.
269 : END IF
270 2 : CALL group%bcast(box_length, source)
271 : IF (ABS(box_length(1) - abc(1)) > 0.0001E0_dp .OR. &
272 2 : ABS(box_length(2) - abc(2)) > 0.0001E0_dp .OR. &
273 : ABS(box_length(3) - abc(3)) > 0.0001E0_dp) THEN
274 2 : IF (ionode) THEN
275 1 : WRITE (iw, *) 'The cell length in the restart file is ', &
276 2 : 'not the same as the input file.'
277 4 : WRITE (iw, *) 'Input file cell length =', abc(1:3)*angstrom
278 4 : WRITE (iw, *) 'Restart file cell length =', box_length(1:3)*angstrom
279 : END IF
280 : END IF
281 :
282 : ! allocate the array holding the coordinates, and read in the coordinates,
283 : ! and write the dat file so we can make a new force_env
284 4 : IF (SUM(nchains(:)) == 0) THEN
285 0 : ALLOCATE (r(3, nunits(1)))
286 0 : ALLOCATE (atom_symbols(nunits(1)))
287 :
288 0 : DO iunit = 1, nunits(1)
289 0 : r(1:3, iunit) = [REAL(iunit, dp), REAL(iunit, dp), REAL(iunit, dp)]
290 0 : IF (LEN_TRIM(atom_names(iunit, 1)) > LEN(atom_symbols(iunit))) THEN
291 0 : CPWARN("The atom symbol will be truncated.")
292 : END IF
293 0 : atom_symbols(iunit) = TRIM(atom_names(iunit, 1))
294 : END DO
295 :
296 0 : IF (ionode) THEN
297 : CALL mc_make_dat_file_new(r(:, :), atom_symbols, 0, &
298 0 : box_length(:), dat_file, nchains(:), mc_input_file)
299 0 : CALL close_file(unit_number=unit)
300 : END IF
301 : ELSE
302 6 : ALLOCATE (r(3, nunits_tot))
303 6 : ALLOCATE (atom_symbols(nunits_tot))
304 :
305 2 : IF (ionode) THEN
306 10 : DO ipart = 1, nunits_tot
307 9 : READ (unit, *) atom_symbols(ipart), r(1:3, ipart)
308 37 : r(1:3, ipart) = r(1:3, ipart)/angstrom
309 : END DO
310 :
311 1 : CALL close_file(unit_number=unit)
312 :
313 : CALL mc_make_dat_file_new(r(:, :), atom_symbols, nunits_tot, &
314 1 : box_length(:), dat_file, nchains(:), mc_input_file)
315 :
316 : END IF
317 : END IF
318 :
319 2 : CALL set_mc_par(mc_par, nstart=nstart)
320 :
321 : ! advance the random number sequence based on the restart step
322 2 : IF (ionode) THEN
323 5 : DO i = 1, nstart + 1
324 5 : rand = rng_stream%next()
325 : END DO
326 : END IF
327 :
328 : ! end the timing
329 2 : CALL timestop(handle)
330 :
331 : ! deallcoate
332 2 : DEALLOCATE (nchains)
333 2 : DEALLOCATE (r)
334 2 : DEALLOCATE (atom_symbols)
335 :
336 2 : END SUBROUTINE read_mc_restart
337 :
338 : ! **************************************************************************************************
339 : !> \brief creates a force environment for any of the different kinds of
340 : !> MC simulations we can do (FIST, QS)
341 : !> \param force_env the force environment to create
342 : !> \param input_declaration ...
343 : !> \param para_env ...
344 : !> \param input_file_name ...
345 : !> \param globenv_new the global environment parameters
346 : !> \author MJM
347 : !> \note Suitable for parallel.
348 : ! **************************************************************************************************
349 210 : SUBROUTINE mc_create_force_env(force_env, input_declaration, para_env, input_file_name, &
350 : globenv_new)
351 :
352 : TYPE(force_env_type), POINTER :: force_env
353 : TYPE(section_type), POINTER :: input_declaration
354 : TYPE(mp_para_env_type), POINTER :: para_env
355 : CHARACTER(LEN=*), INTENT(IN) :: input_file_name
356 : TYPE(global_environment_type), OPTIONAL, POINTER :: globenv_new
357 :
358 : INTEGER :: f_env_id, ierr, output_unit
359 : TYPE(f_env_type), POINTER :: f_env
360 :
361 210 : output_unit = cp_logger_get_default_unit_nr()
362 : CALL create_force_env(f_env_id, &
363 : input_declaration=input_declaration, &
364 : input_path=input_file_name, &
365 : output_unit=output_unit, &
366 210 : mpi_comm=para_env)
367 :
368 210 : CALL f_env_add_defaults(f_env_id, f_env)
369 210 : force_env => f_env%force_env
370 210 : CALL force_env_retain(force_env)
371 210 : CALL f_env_rm_defaults(f_env)
372 210 : CALL destroy_force_env(f_env_id, ierr, .FALSE.)
373 210 : IF (ierr /= 0) CPABORT("mc_create_force_env: destroy_force_env failed")
374 :
375 210 : IF (PRESENT(globenv_new)) &
376 6 : CALL force_env_get(force_env, globenv=globenv_new)
377 :
378 210 : END SUBROUTINE mc_create_force_env
379 :
380 : ! **************************************************************************************************
381 : !> \brief essentially copies the cell size and coordinates of one force env
382 : !> to another that we will use to bias some moves with
383 : !> \param bias_env the force environment to create
384 : !> \param r ...
385 : !> \param atom_symbols ...
386 : !> \param nunits_tot ...
387 : !> \param para_env ...
388 : !> \param box_length ...
389 : !> \param nchains ...
390 : !> \param input_declaration ...
391 : !> \param mc_input_file ...
392 : !> \param ionode ...
393 : !> \author MJM
394 : !> \note Suitable for parallel.
395 : ! **************************************************************************************************
396 114 : SUBROUTINE mc_create_bias_force_env(bias_env, r, atom_symbols, nunits_tot, &
397 : para_env, box_length, nchains, input_declaration, mc_input_file, ionode)
398 :
399 : TYPE(force_env_type), POINTER :: bias_env
400 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: r
401 : CHARACTER(default_string_length), DIMENSION(:), &
402 : INTENT(IN) :: atom_symbols
403 : INTEGER, INTENT(IN) :: nunits_tot
404 : TYPE(mp_para_env_type), POINTER :: para_env
405 : REAL(KIND=dp), DIMENSION(1:3), INTENT(IN) :: box_length
406 : INTEGER, DIMENSION(:), POINTER :: nchains
407 : TYPE(section_type), POINTER :: input_declaration
408 : TYPE(mc_input_file_type), POINTER :: mc_input_file
409 : LOGICAL, INTENT(IN) :: ionode
410 :
411 114 : IF (ionode) &
412 : CALL mc_make_dat_file_new(r(:, :), atom_symbols, nunits_tot, &
413 57 : box_length(:), 'bias_temp.dat', nchains(:), mc_input_file)
414 :
415 114 : CALL mc_create_force_env(bias_env, input_declaration, para_env, 'bias_temp.dat')
416 :
417 114 : END SUBROUTINE mc_create_bias_force_env
418 :
419 : END MODULE mc_control
|
