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 preps the system for a Monte Carlo run (sets up some environments,
10 : !> calls the routines to read in the MC parameters)...converted
11 : !> from qs_mc.F
12 : !> \par Literature
13 : !> a list of papers for the theory behind various MC moves
14 : !> Books:
15 : !> D. Frenkel, B. Smit: Understanding Molecular Simulation (1996)
16 : !> M.P. Allen, D.J. Tildesley: Computer Simulations of Liquids (1987)
17 : !> \par
18 : !> Aggregation volume bias Monte Carlo (AVBMC):
19 : !> Chen, B.; Siepmann, J.I. J. Phys. Chem. B 2000, 104, 8725.
20 : !> \par
21 : !> Biasing with an inexpensive potential:
22 : !> Iftimie et al. J. Chem. Phys. 2000, 113, 4852.
23 : !> Gelb, L. D. J. Chem. Phys. 2003, 118, 7747.
24 : !> \par
25 : !> Configurational bias Monte Carlo (CBMC):
26 : !> Siepmann, J.I.; Frenkel, D. Mol. Phys. 1992, 75, 59.
27 : !> \par
28 : !> Gibbs ensemble Monte Carlo (GEMC):
29 : !> Panagiotopoulos, A.Z. Mol. Phys. 1987, 61, 813.
30 : !> Panagiotopoulos et al. Mol. Phys. 1988, 63, 527.
31 : !> Smit et al. Mol. Phys. 1989, 68, 931.
32 : !> \par
33 : !> Isobaric-isothermal ensemble:
34 : !> McDonald, I.R. Mol. Phys. 1972, 23, 41.
35 : !> \par
36 : !> Original Monte Carlo paper:
37 : !> Metropolis et al. J. Chem. Phys. 1953, 21, 1087.
38 : !> \author MJM-Oct-15-03
39 : ! **************************************************************************************************
40 : MODULE mc_run
41 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
42 : USE force_env_types, ONLY: force_env_p_type,&
43 : force_env_release,&
44 : force_env_retain,&
45 : force_env_type,&
46 : use_fist_force
47 : USE global_types, ONLY: global_environment_type
48 : USE input_section_types, ONLY: section_type,&
49 : section_vals_get,&
50 : section_vals_get_subs_vals,&
51 : section_vals_type,&
52 : section_vals_val_get
53 : USE kinds, ONLY: dp
54 : USE mc_control, ONLY: mc_create_force_env,&
55 : read_mc_restart
56 : USE mc_ensembles, ONLY: mc_compute_virial,&
57 : mc_run_ensemble
58 : USE mc_environment_types, ONLY: get_mc_env,&
59 : mc_env_create,&
60 : mc_env_release,&
61 : mc_environment_p_type,&
62 : set_mc_env
63 : USE mc_types, ONLY: &
64 : find_mc_rcut, get_mc_molecule_info, get_mc_par, mc_determine_molecule_info, &
65 : mc_input_file_create, mc_input_file_destroy, mc_input_file_type, &
66 : mc_input_parameters_check, mc_molecule_info_destroy, mc_molecule_info_type, &
67 : mc_sim_par_create, mc_sim_par_destroy, mc_simulation_parameters_p_type, read_mc_section, &
68 : set_mc_par
69 : USE message_passing, ONLY: mp_para_env_type
70 : USE parallel_rng_types, ONLY: UNIFORM,&
71 : rng_stream_type
72 : USE physcon, ONLY: angstrom
73 : #include "../../base/base_uses.f90"
74 :
75 : IMPLICIT NONE
76 :
77 : PRIVATE
78 : ! *** Global parameters ***
79 :
80 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_run'
81 :
82 : PUBLIC :: do_mon_car
83 :
84 : !-----------------------------------------------------------------------------!
85 :
86 : CONTAINS
87 :
88 : ! **************************************************************************************************
89 : !> \brief starts the Monte Carlo simulation and determines which ensemble we're
90 : !> running
91 : !> \param force_env_1 the force environment for the simulation, or
92 : !> the force environment for box 1, depending on which
93 : !> ensemble we're running
94 : !> \param globenv the global environment for the simulation
95 : !> \param input_declaration ...
96 : !> \param input_file_name the name of the input file that force_env_1 was
97 : !> created from
98 : !> \author MJM
99 : !> \note
100 : !> Designed for parallel.
101 : ! **************************************************************************************************
102 20 : SUBROUTINE do_mon_car(force_env_1, globenv, input_declaration, input_file_name)
103 :
104 : TYPE(force_env_type), POINTER :: force_env_1
105 : TYPE(global_environment_type), POINTER :: globenv
106 : TYPE(section_type), POINTER :: input_declaration
107 : CHARACTER(LEN=*) :: input_file_name
108 :
109 : CHARACTER(LEN=20) :: ensemble
110 : CHARACTER(LEN=40) :: box2_file, dat_file
111 : INTEGER :: box_number, ibox, isos, iw, nboxes, &
112 : nmol_types
113 20 : INTEGER, DIMENSION(:), POINTER :: nunits_tot
114 : LOGICAL :: lbias, lhmc, lrestart, lskip_box, &
115 : lterminate
116 : REAL(dp) :: rcut
117 20 : REAL(dp), DIMENSION(:, :), POINTER :: empty_coordinates
118 20 : TYPE(force_env_p_type), DIMENSION(:), POINTER :: force_env
119 : TYPE(force_env_type), POINTER :: force_env_2
120 : TYPE(global_environment_type), POINTER :: globenv_2
121 20 : TYPE(mc_environment_p_type), DIMENSION(:), POINTER :: mc_env
122 : TYPE(mc_input_file_type), POINTER :: mc_bias_file, mc_input_file
123 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
124 : TYPE(mc_simulation_parameters_p_type), &
125 20 : DIMENSION(:), POINTER :: mc_par
126 : TYPE(mp_para_env_type), POINTER :: para_env
127 : TYPE(rng_stream_type) :: rng_stream
128 : TYPE(section_vals_type), POINTER :: force_env_section, mc_section, &
129 : root_section
130 :
131 20 : NULLIFY (mc_env, mc_par, force_env_2, force_env_section, &
132 20 : root_section, mc_molecule_info)
133 :
134 20 : CALL force_env_retain(force_env_1)
135 :
136 20 : para_env => force_env_1%para_env
137 20 : iw = cp_logger_get_default_io_unit()
138 20 : force_env_section => force_env_1%force_env_section
139 20 : root_section => force_env_1%root_section
140 20 : CALL section_vals_get(force_env_section, n_repetition=isos)
141 20 : CPASSERT(isos == 1)
142 : ! set some values...will use get_globenv if that ever comes around
143 :
144 : ! initialize the random numbers
145 : rng_stream = rng_stream_type( &
146 : last_rng_stream=force_env_1%globenv%gaussian_rng_stream, &
147 : name="first", &
148 20 : distribution_type=UNIFORM)
149 :
150 : ! need to figure out how many boxes we have, based on the value
151 : ! of mc_par % ensemble
152 20 : NULLIFY (mc_section)
153 : mc_section => section_vals_get_subs_vals(root_section, &
154 20 : "MOTION%MC")
155 : CALL section_vals_val_get(mc_section, "ENSEMBLE", &
156 20 : c_val=ensemble)
157 :
158 : ! now we read in the second force_env, if we have another box
159 14 : SELECT CASE (ensemble)
160 : CASE ("TRADITIONAL", "VIRIAL")
161 14 : nboxes = 1
162 : CASE ("GEMC_NVT", "GEMC_NPT")
163 6 : nboxes = 2
164 : CALL section_vals_val_get(mc_section, "BOX2_FILE_NAME", &
165 6 : c_val=box2_file)
166 : CALL mc_create_force_env(force_env_2, input_declaration, para_env, &
167 26 : box2_file, globenv_new=globenv_2)
168 : END SELECT
169 :
170 : ! now we create the various pointers that contain information for all boxes
171 86 : ALLOCATE (force_env(1:nboxes))
172 14 : SELECT CASE (ensemble)
173 : CASE ("TRADITIONAL", "VIRIAL")
174 14 : force_env(1)%force_env => force_env_1
175 : CASE ("GEMC_NVT", "GEMC_NPT")
176 6 : force_env(1)%force_env => force_env_1
177 20 : force_env(2)%force_env => force_env_2
178 : END SELECT
179 66 : ALLOCATE (mc_par(1:nboxes))
180 66 : ALLOCATE (mc_env(1:nboxes))
181 :
182 : ! now we need the molecule information
183 : ! determine the total number of molecules and atoms
184 : CALL mc_determine_molecule_info(force_env, mc_molecule_info, &
185 20 : coordinates_empty=empty_coordinates)
186 : CALL get_mc_molecule_info(mc_molecule_info, nmol_types=nmol_types, &
187 20 : nunits_tot=nunits_tot)
188 :
189 46 : DO ibox = 1, nboxes
190 :
191 26 : IF (iw > 0) THEN
192 13 : WRITE (iw, *)
193 13 : WRITE (iw, *)
194 13 : WRITE (iw, '(A,I2,A)') '******************************** Begin BOX ', ibox, &
195 26 : ' *******************************'
196 : END IF
197 :
198 : ! allocates an mc_env and sets the variables to zero
199 26 : ALLOCATE (mc_env(ibox)%mc_env)
200 26 : CALL mc_env_create(mc_env(ibox)%mc_env)
201 :
202 : ! now read in the values of the mc_pars
203 : ! creating the mc_par
204 26 : CALL mc_sim_par_create(mc_par(ibox)%mc_par, nmol_types)
205 :
206 : ! attach molecule information to all mc_par structures, so we know what we have to read in
207 26 : CALL set_mc_par(mc_par(ibox)%mc_par, mc_molecule_info=mc_molecule_info)
208 :
209 : ! read the input of the Monte Carlo section
210 26 : force_env_section => force_env(ibox)%force_env%force_env_section
211 26 : root_section => force_env(ibox)%force_env%root_section
212 26 : IF (ibox == 1) THEN
213 : CALL read_mc_section(mc_par(ibox)%mc_par, para_env, globenv, input_file_name, &
214 20 : root_section, force_env_section)
215 : ELSE
216 : CALL read_mc_section(mc_par(ibox)%mc_par, para_env, globenv, box2_file, &
217 6 : root_section, force_env_section)
218 : END IF
219 :
220 : ! get the input file data, in case we need to make a restart...
221 : ! this is also used in the swap move, or anytime we need to make a dat file
222 : ! always judge based on lbias from box 1...in case someone only changes the value
223 : ! for box 1
224 26 : CALL get_mc_par(mc_par(ibox)%mc_par, mc_input_file=mc_input_file, lhmc=lhmc)
225 26 : CALL get_mc_par(mc_par(1)%mc_par, lbias=lbias)
226 26 : IF (ibox == 1) THEN
227 : CALL mc_input_file_create(mc_input_file, &
228 20 : input_file_name, mc_molecule_info, empty_coordinates, lhmc)
229 : ELSE
230 : CALL mc_input_file_create(mc_input_file, &
231 6 : box2_file, mc_molecule_info, empty_coordinates, lhmc)
232 : END IF
233 26 : CALL set_mc_par(mc_par(ibox)%mc_par, mc_input_file=mc_input_file)
234 :
235 26 : IF (lbias) THEN
236 14 : CALL get_mc_par(mc_par(ibox)%mc_par, mc_bias_file=mc_bias_file)
237 : CALL mc_input_file_create(mc_bias_file, &
238 14 : "bias_template.inp", mc_molecule_info, empty_coordinates, lhmc)
239 14 : CALL set_mc_par(mc_par(ibox)%mc_par, mc_bias_file=mc_bias_file)
240 : END IF
241 :
242 : ! check for restart
243 : CALL get_mc_par(mc_par(ibox)%mc_par, lrestart=lrestart, &
244 26 : dat_file=dat_file)
245 72 : IF (lrestart) THEN
246 : CALL read_mc_restart(mc_par(ibox)%mc_par, force_env(ibox)%force_env, &
247 2 : iw, nunits_tot(ibox), rng_stream)
248 : ! release the old force env and make the new one
249 2 : CALL force_env_release(force_env(ibox)%force_env)
250 : CALL mc_create_force_env(force_env(ibox)%force_env, &
251 2 : input_declaration, para_env, dat_file)
252 : END IF
253 :
254 : END DO
255 :
256 : ! figure out if we have an empty box
257 20 : box_number = 0
258 20 : lskip_box = .FALSE.
259 46 : DO ibox = 1, nboxes
260 46 : IF (nunits_tot(ibox) == 0) THEN
261 0 : IF (lskip_box) THEN
262 0 : CPABORT('More than one box has no atoms in it!')
263 : END IF
264 0 : box_number = ibox
265 0 : lskip_box = .TRUE.
266 : END IF
267 : END DO
268 :
269 : ! in case there was a restart, we need to do this again
270 20 : CALL mc_molecule_info_destroy(mc_molecule_info)
271 20 : CALL mc_determine_molecule_info(force_env, mc_molecule_info, box_number=box_number)
272 : CALL get_mc_molecule_info(mc_molecule_info, nmol_types=nmol_types, &
273 20 : nunits_tot=nunits_tot)
274 46 : DO ibox = 1, nboxes
275 46 : CALL set_mc_par(mc_par(ibox)%mc_par, mc_molecule_info=mc_molecule_info)
276 : END DO
277 : ! if we're doing a classical simulation, figure out the largest
278 : ! potential cutoff and write it to the screen
279 20 : IF (force_env(1)%force_env%in_use == use_fist_force) THEN
280 14 : CALL find_mc_rcut(mc_par(1)%mc_par, force_env(1)%force_env, lterminate)
281 14 : CALL get_mc_par(mc_par(1)%mc_par, rcut=rcut)
282 14 : IF (iw > 0) WRITE (iw, '( A,T73,F8.4 )') &
283 7 : ' MC| Interaction cutoff [angstroms]', rcut*angstrom
284 14 : IF (lterminate) THEN
285 0 : CPABORT('Cutoff larger than twice the boxlength')
286 : END IF
287 : END IF
288 :
289 : ! make sure some values are the same between boxes
290 20 : IF (nboxes == 2) THEN
291 6 : CALL equilize_mc_sim_parameters(mc_par, iw)
292 : END IF
293 :
294 : ! now check the input parameters and run the simulation
295 46 : DO ibox = 1, nboxes
296 :
297 26 : CALL mc_input_parameters_check(mc_par(ibox)%mc_par)
298 :
299 : ! attach all the structures to one convientent structure
300 : CALL set_mc_env(mc_env(ibox)%mc_env, mc_par=mc_par(ibox)%mc_par, &
301 46 : force_env=force_env(ibox)%force_env)
302 :
303 : END DO
304 :
305 : ! if we're computing the second virial coefficient, do that instead
306 : ! of running a simulation
307 2 : SELECT CASE (ensemble)
308 : CASE ("VIRIAL")
309 2 : CALL mc_compute_virial(mc_env, rng_stream)
310 : CASE DEFAULT
311 20 : CALL mc_run_ensemble(mc_env, para_env, globenv, input_declaration, nboxes, rng_stream)
312 : END SELECT
313 :
314 : ! get rid of all the MC molecule information
315 20 : CALL get_mc_env(mc_env(1)%mc_env, mc_par=mc_par(1)%mc_par)
316 20 : CALL get_mc_par(mc_par(1)%mc_par, mc_molecule_info=mc_molecule_info)
317 20 : CALL mc_molecule_info_destroy(mc_molecule_info)
318 :
319 46 : DO ibox = 1, nboxes
320 : CALL get_mc_env(mc_env(ibox)%mc_env, &
321 26 : mc_par=mc_par(ibox)%mc_par, force_env=force_env(ibox)%force_env)
322 26 : CALL get_mc_par(mc_par(ibox)%mc_par, mc_input_file=mc_input_file)
323 :
324 26 : CALL mc_input_file_destroy(mc_input_file)
325 26 : IF (lbias) THEN
326 14 : CALL get_mc_par(mc_par(ibox)%mc_par, mc_bias_file=mc_bias_file)
327 14 : CALL mc_input_file_destroy(mc_bias_file)
328 : END IF
329 :
330 26 : CALL mc_sim_par_destroy(mc_par(ibox)%mc_par)
331 26 : CALL mc_env_release(mc_env(ibox)%mc_env)
332 26 : DEALLOCATE (mc_env(ibox)%mc_env)
333 46 : CALL force_env_release(force_env(ibox)%force_env)
334 : END DO
335 :
336 20 : DEALLOCATE (empty_coordinates)
337 20 : DEALLOCATE (mc_par)
338 20 : DEALLOCATE (mc_env)
339 20 : DEALLOCATE (force_env)
340 :
341 560 : END SUBROUTINE do_mon_car
342 :
343 : ! **************************************************************************************************
344 : !> \brief takes some parameters from one set of MC simulation parameters
345 : !> and transfers them to a second set...used so that we're not using
346 : !> different parameters between two simulation boxes, if they should
347 : !> be the same (move probabilities, for instance)
348 : !> \param mc_par the pointer that contains the simulation parameters
349 : !> for both boxes
350 : !> \param iw the unit number that prints to screen
351 : !> \author MJM
352 : !> \note
353 : !> Designed for parallel.
354 : ! **************************************************************************************************
355 12 : SUBROUTINE equilize_mc_sim_parameters(mc_par, iw)
356 : TYPE(mc_simulation_parameters_p_type), &
357 : DIMENSION(:), POINTER :: mc_par
358 : INTEGER, INTENT(IN) :: iw
359 :
360 : CHARACTER(20) :: ensemble
361 : INTEGER :: iprint, iupcltrans, iuptrans, iupvolume, &
362 : nmoves, nstep, nswapmoves
363 6 : INTEGER, DIMENSION(:), POINTER :: avbmc_atom
364 : LOGICAL :: lbias, lrestart, lstop
365 : REAL(dp) :: BETA, pmswap, pmtraion, pmtrans, &
366 : pmvolume, pressure, rcut, temperature
367 6 : REAL(dp), DIMENSION(:), POINTER :: avbmc_rmax, avbmc_rmin, pbias, &
368 6 : pmavbmc_mol, pmrot_mol, pmswap_mol, &
369 6 : pmtraion_mol, pmtrans_mol
370 :
371 6 : IF (iw > 0) THEN
372 3 : WRITE (iw, '( A,A )') 'Ignoring some input for box 2, and ', &
373 6 : 'using the values for box 1 for the following variables:'
374 3 : WRITE (iw, '( A,A )') 'nstep,iupvolume,iuptrans,iupcltrans,nmoves,', &
375 6 : 'nswapmoves,iprint,lbias,lstop,temperature,pressure'
376 3 : WRITE (iw, '( A,A )') 'pmvolume,pmtraion,pmtrans,BETA,rcut,', &
377 6 : 'lrestart'
378 3 : WRITE (iw, '( A,A )') 'pmtraion_mol,pmtrans_mol,pmrot_mol,pmswap_mol,', &
379 6 : 'avbmc_atom'
380 3 : WRITE (iw, '( A,A )') 'avbmc_rmin,avmbc_rmax,pmavbmc_mol,pbias,pmswap'
381 : END IF
382 :
383 : CALL get_mc_par(mc_par(1)%mc_par, nstep=nstep, iupvolume=iupvolume, iupcltrans=iupcltrans, &
384 : iuptrans=iuptrans, nmoves=nmoves, nswapmoves=nswapmoves, &
385 : iprint=iprint, lbias=lbias, lstop=lstop, temperature=temperature, &
386 : pressure=pressure, pmswap=pmswap, pmvolume=pmvolume, &
387 : pmtraion=pmtraion, pmtrans=pmtrans, BETA=BETA, rcut=rcut, &
388 : lrestart=lrestart, pmtraion_mol=pmtraion_mol, pmtrans_mol=pmtrans_mol, &
389 : pmrot_mol=pmrot_mol, pmswap_mol=pmswap_mol, avbmc_atom=avbmc_atom, &
390 : avbmc_rmin=avbmc_rmin, avbmc_rmax=avbmc_rmax, pmavbmc_mol=pmavbmc_mol, &
391 6 : pbias=pbias, ensemble=ensemble)
392 : CALL set_mc_par(mc_par(2)%mc_par, nstep=nstep, iupvolume=iupvolume, iupcltrans=iupcltrans, &
393 : iuptrans=iuptrans, nmoves=nmoves, nswapmoves=nswapmoves, &
394 : iprint=iprint, lbias=lbias, lstop=lstop, temperature=temperature, &
395 : pressure=pressure, pmswap=pmswap, pmvolume=pmvolume, &
396 : pmtraion=pmtraion, pmtrans=pmtrans, BETA=BETA, rcut=rcut, &
397 : lrestart=lrestart, pmtraion_mol=pmtraion_mol, pmtrans_mol=pmtrans_mol, &
398 : pmrot_mol=pmrot_mol, pmswap_mol=pmswap_mol, avbmc_atom=avbmc_atom, &
399 : avbmc_rmin=avbmc_rmin, avbmc_rmax=avbmc_rmax, pmavbmc_mol=pmavbmc_mol, &
400 6 : pbias=pbias, ensemble=ensemble)
401 :
402 6 : END SUBROUTINE equilize_mc_sim_parameters
403 :
404 : END MODULE mc_run
|
