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 holds all the structure types needed for Monte Carlo, except
10 : !> the mc_environment_type
11 : !> \par History
12 : !> none
13 : !> \author MJM
14 : ! **************************************************************************************************
15 : MODULE mc_types
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind
18 : USE cell_types, ONLY: cell_type,&
19 : get_cell
20 : USE cp_files, ONLY: close_file,&
21 : open_file
22 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
23 : USE cp_subsys_types, ONLY: cp_subsys_get,&
24 : cp_subsys_type
25 : USE cp_units, ONLY: cp_unit_to_cp2k
26 : USE fist_environment_types, ONLY: fist_env_get,&
27 : fist_environment_type
28 : USE fist_nonbond_env_types, ONLY: fist_nonbond_env_get,&
29 : fist_nonbond_env_type
30 : USE force_env_types, ONLY: force_env_get,&
31 : force_env_p_type,&
32 : force_env_type,&
33 : use_fist_force,&
34 : use_qs_force
35 : USE global_types, ONLY: global_environment_type
36 : USE input_constants, ONLY: do_fist,&
37 : do_qs
38 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
39 : section_vals_type,&
40 : section_vals_val_get,&
41 : section_vals_val_set
42 : USE kinds, ONLY: default_path_length,&
43 : default_string_length,&
44 : dp
45 : USE mathconstants, ONLY: pi
46 : USE message_passing, ONLY: mp_comm_type,&
47 : mp_para_env_type
48 : USE molecule_kind_list_types, ONLY: molecule_kind_list_p_type
49 : USE molecule_kind_types, ONLY: atom_type,&
50 : get_molecule_kind,&
51 : molecule_kind_type
52 : USE pair_potential_types, ONLY: pair_potential_pp_type
53 : USE particle_list_types, ONLY: particle_list_p_type
54 : USE physcon, ONLY: boltzmann,&
55 : joule
56 : USE string_utilities, ONLY: uppercase,&
57 : xstring
58 : #include "../../base/base_uses.f90"
59 :
60 : IMPLICIT NONE
61 :
62 : PRIVATE
63 :
64 : ! *** Global parameters ***
65 :
66 : PUBLIC :: mc_simpar_type, &
67 : mc_simulation_parameters_p_type, &
68 : mc_averages_type, mc_averages_p_type, &
69 : mc_moves_type, mc_moves_p_type, accattempt, &
70 : get_mc_par, set_mc_par, read_mc_section, &
71 : find_mc_rcut, mc_molecule_info_type, &
72 : mc_determine_molecule_info, mc_molecule_info_destroy, &
73 : mc_sim_par_create, mc_sim_par_destroy, &
74 : get_mc_molecule_info, &
75 : mc_input_file_type, mc_input_file_create, &
76 : get_mc_input_file, mc_input_file_destroy, &
77 : mc_input_parameters_check, &
78 : mc_ekin_type
79 :
80 : ! **************************************************************************************************
81 : TYPE mc_simpar_type
82 : ! contains all the parameters needed for running Monte Carlo simulations
83 : PRIVATE
84 : INTEGER, DIMENSION(:), POINTER :: avbmc_atom => NULL()
85 : INTEGER :: nstep = 0
86 : INTEGER :: iupvolume = 0
87 : INTEGER :: iuptrans = 0
88 : INTEGER :: iupcltrans = 0
89 : INTEGER :: nvirial = 0
90 : INTEGER :: nbox = 0
91 : INTEGER :: nmoves = 0
92 : INTEGER :: nswapmoves = 0
93 : INTEGER :: rm = 0
94 : INTEGER :: cl = 0
95 : INTEGER :: diff = 0
96 : INTEGER :: nstart = 0
97 : INTEGER :: source = 0
98 : TYPE(mp_comm_type) :: group = mp_comm_type()
99 : INTEGER :: iprint = 0
100 : LOGICAL :: ldiscrete = .FALSE.
101 : LOGICAL :: lbias = .FALSE.
102 : LOGICAL :: ionode = .FALSE.
103 : LOGICAL :: lrestart = .FALSE.
104 : LOGICAL :: lstop = .FALSE.
105 : LOGICAL :: lhmc = .FALSE.
106 : CHARACTER(LEN=20) :: ensemble = ""
107 : CHARACTER(LEN=default_path_length) :: restart_file_name = ""
108 : CHARACTER(LEN=default_path_length) :: molecules_file = ""
109 : CHARACTER(LEN=default_path_length) :: moves_file = ""
110 : CHARACTER(LEN=default_path_length) :: coords_file = ""
111 : CHARACTER(LEN=default_path_length) :: energy_file = ""
112 : CHARACTER(LEN=default_path_length) :: displacement_file = ""
113 : CHARACTER(LEN=default_path_length) :: cell_file = ""
114 : CHARACTER(LEN=default_path_length) :: dat_file = ""
115 : CHARACTER(LEN=default_path_length) :: data_file = ""
116 : CHARACTER(LEN=default_path_length) :: box2_file = ""
117 : CHARACTER(LEN=200) :: fft_lib = ""
118 : CHARACTER(LEN=50) :: PROGRAM = ""
119 : REAL(dp), DIMENSION(:), POINTER :: pmtraion_mol => NULL(), pmtrans_mol => NULL(), &
120 : pmrot_mol => NULL(), pmswap_mol => NULL(), &
121 : pbias => NULL(), pmavbmc_mol => NULL()
122 : REAL(dp) :: discrete_step = 0.0_dp
123 : REAL(dp) :: rmvolume = 0.0_dp, rmcltrans = 0.0_dp
124 : REAL(dp), DIMENSION(:), POINTER :: rmbond => NULL(), rmangle => NULL(), rmdihedral => NULL(), &
125 : rmrot => NULL(), rmtrans => NULL()
126 : REAL(dp), DIMENSION(:), POINTER :: eta => NULL()
127 : REAL(dp) :: temperature = 0.0_dp
128 : REAL(dp) :: pressure = 0.0_dp
129 : REAL(dp) :: rclus = 0.0_dp
130 : REAL(dp) :: pmavbmc = 0.0_dp
131 : REAL(dp) :: pmswap = 0.0_dp
132 : REAL(dp) :: pmvolume = 0.0_dp, pmvol_box = 0.0_dp, pmclus_box = 0.0_dp
133 : REAL(dp) :: pmhmc = 0.0_dp, pmhmc_box = 0.0_dp
134 : REAL(dp) :: pmtraion = 0.0_dp
135 : REAL(dp) :: pmtrans = 0.0_dp
136 : REAL(dp) :: pmcltrans = 0.0_dp
137 : REAL(dp) :: BETA = 0.0_dp
138 : REAL(dp) :: rcut = 0.0_dp
139 : REAL(dp), DIMENSION(:), POINTER :: avbmc_rmin => NULL(), avbmc_rmax => NULL()
140 : REAL(dp), DIMENSION(:), POINTER :: virial_temps => NULL()
141 : TYPE(mc_input_file_type), POINTER :: &
142 : mc_input_file => NULL(), mc_bias_file => NULL()
143 : TYPE(section_vals_type), POINTER :: input_file => NULL()
144 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info => NULL()
145 : REAL(dp) :: exp_min_val = 0.0_dp, exp_max_val = 0.0_dp, min_val = 0.0_dp, max_val = 0.0_dp
146 : INTEGER :: rand2skip = 0
147 : END TYPE mc_simpar_type
148 :
149 : ! **************************************************************************************************
150 : TYPE mc_ekin_type
151 : ! contains the kinetic energy of an MD sequence from hybrid MC
152 : REAL(dp) :: initial_ekin = 0.0_dp, final_ekin = 0.0_dp
153 : END TYPE mc_ekin_type
154 : ! **************************************************************************************************
155 : TYPE mc_input_file_type
156 : ! contains all the text of the input file
157 : PRIVATE
158 : INTEGER :: run_type_row = 0, run_type_column = 0, coord_row_start = 0, coord_row_end = 0, &
159 : cell_row = 0, cell_column = 0, force_eval_row_start = 0, force_eval_row_end = 0, &
160 : global_row_start = 0, global_row_end = 0, in_use = 0, nunits_empty = 0, &
161 : motion_row_end = 0, motion_row_start = 0
162 : INTEGER, DIMENSION(:), POINTER :: mol_set_nmol_row => NULL(), mol_set_nmol_column => NULL()
163 : CHARACTER(LEN=default_string_length), DIMENSION(:), POINTER :: text => NULL()
164 : CHARACTER(LEN=default_string_length), DIMENSION(:), POINTER :: atom_names_empty => NULL()
165 : REAL(dp), DIMENSION(:, :), POINTER :: coordinates_empty => NULL()
166 : END TYPE mc_input_file_type
167 :
168 : ! **************************************************************************************************
169 : TYPE mc_molecule_info_type
170 : ! contains information on molecules...I had to do this because I use multiple force
171 : ! environments, and they know nothing about each other
172 : PRIVATE
173 : INTEGER :: nboxes = 0
174 : CHARACTER(LEN=default_string_length), DIMENSION(:), POINTER :: names => NULL()
175 : CHARACTER(LEN=default_string_length), &
176 : DIMENSION(:, :), POINTER :: atom_names => NULL()
177 : REAL(dp), DIMENSION(:, :), POINTER :: conf_prob => NULL(), mass => NULL()
178 : INTEGER, DIMENSION(:, :), POINTER :: nchains => NULL()
179 : INTEGER :: nmol_types = 0, nchain_total = 0
180 : INTEGER, DIMENSION(:), POINTER :: nunits => NULL(), mol_type => NULL(), nunits_tot => NULL(), in_box => NULL()
181 : END TYPE mc_molecule_info_type
182 :
183 : ! **************************************************************************************************
184 : TYPE mc_simulation_parameters_p_type
185 : ! a pointer for multiple boxes
186 : TYPE(mc_simpar_type), POINTER :: mc_par => NULL()
187 : END TYPE mc_simulation_parameters_p_type
188 :
189 : ! **************************************************************************************************
190 : TYPE mc_averages_type
191 : ! contains some data that is averaged throughout the course of a run
192 : REAL(KIND=dp) :: ave_energy = 0.0_dp
193 : REAL(KIND=dp) :: ave_energy_squared = 0.0_dp
194 : REAL(KIND=dp) :: ave_volume = 0.0_dp
195 : REAL(KIND=dp) :: molecules = 0.0_dp
196 : END TYPE mc_averages_type
197 :
198 : ! **************************************************************************************************
199 : TYPE mc_averages_p_type
200 : TYPE(mc_averages_type), POINTER :: averages => NULL()
201 : END TYPE mc_averages_p_type
202 :
203 : ! **************************************************************************************************
204 : TYPE mc_moves_type
205 : ! contains data for how many moves of a give type we've accepted/rejected
206 : TYPE(accattempt), POINTER :: bias_bond => NULL()
207 : TYPE(accattempt), POINTER :: bias_angle => NULL()
208 : TYPE(accattempt), POINTER :: bias_dihedral => NULL()
209 : TYPE(accattempt), POINTER :: bias_trans => NULL()
210 : TYPE(accattempt), POINTER :: bias_cltrans => NULL()
211 : TYPE(accattempt), POINTER :: bias_rot => NULL()
212 : TYPE(accattempt), POINTER :: bond => NULL()
213 : TYPE(accattempt), POINTER :: angle => NULL()
214 : TYPE(accattempt), POINTER :: dihedral => NULL()
215 : TYPE(accattempt), POINTER :: trans => NULL()
216 : TYPE(accattempt), POINTER :: cltrans => NULL()
217 : TYPE(accattempt), POINTER :: rot => NULL()
218 : TYPE(accattempt), POINTER :: swap => NULL()
219 : TYPE(accattempt), POINTER :: volume => NULL()
220 : TYPE(accattempt), POINTER :: hmc => NULL()
221 : TYPE(accattempt), POINTER :: avbmc_inin => NULL()
222 : TYPE(accattempt), POINTER :: avbmc_outin => NULL()
223 : TYPE(accattempt), POINTER :: avbmc_inout => NULL()
224 : TYPE(accattempt), POINTER :: avbmc_outout => NULL()
225 : TYPE(accattempt), POINTER :: Quickstep => NULL()
226 : REAL(KIND=dp) :: trans_dis = 0.0_dp, qtrans_dis = 0.0_dp
227 : REAL(KIND=dp) :: cltrans_dis = 0.0_dp, qcltrans_dis = 0.0_dp
228 : INTEGER :: empty = 0, grown = 0, empty_conf = 0, empty_avbmc = 0
229 : END TYPE mc_moves_type
230 :
231 : ! **************************************************************************************************
232 : TYPE accattempt
233 : INTEGER :: successes = 0
234 : INTEGER :: qsuccesses = 0
235 : INTEGER :: attempts = 0
236 : END TYPE accattempt
237 :
238 : ! **************************************************************************************************
239 : TYPE mc_moves_p_type
240 : TYPE(mc_moves_type), POINTER :: moves => NULL()
241 : END TYPE mc_moves_p_type
242 :
243 : ! *** Global parameters ***
244 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_types'
245 :
246 : CONTAINS
247 :
248 : ! **************************************************************************************************
249 : !> \brief accesses the private elements of the mc_input_file_type
250 : !> \param mc_input_file the input file you want data for
251 : !>
252 : !> Suitable for parallel.
253 : !> \param run_type_row ...
254 : !> \param run_type_column ...
255 : !> \param coord_row_start ...
256 : !> \param coord_row_end ...
257 : !> \param cell_row ...
258 : !> \param cell_column ...
259 : !> \param force_eval_row_start ...
260 : !> \param force_eval_row_end ...
261 : !> \param global_row_start ...
262 : !> \param global_row_end ...
263 : !> \param mol_set_nmol_row ...
264 : !> \param mol_set_nmol_column ...
265 : !> \param in_use ...
266 : !> \param text ...
267 : !> \param atom_names_empty ...
268 : !> \param nunits_empty ...
269 : !> \param coordinates_empty ...
270 : !> \param motion_row_end ...
271 : !> \param motion_row_start ...
272 : !> \author MJM
273 : ! **************************************************************************************************
274 102 : SUBROUTINE get_mc_input_file(mc_input_file, run_type_row, run_type_column, &
275 : coord_row_start, coord_row_end, cell_row, cell_column, &
276 : force_eval_row_start, force_eval_row_end, global_row_start, global_row_end, &
277 : mol_set_nmol_row, mol_set_nmol_column, in_use, text, atom_names_empty, &
278 : nunits_empty, coordinates_empty, motion_row_end, motion_row_start)
279 :
280 : TYPE(mc_input_file_type), POINTER :: mc_input_file
281 : INTEGER, INTENT(OUT), OPTIONAL :: run_type_row, run_type_column, coord_row_start, &
282 : coord_row_end, cell_row, cell_column, force_eval_row_start, force_eval_row_end, &
283 : global_row_start, global_row_end
284 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: mol_set_nmol_row, mol_set_nmol_column
285 : INTEGER, INTENT(OUT), OPTIONAL :: in_use
286 : CHARACTER(LEN=default_string_length), &
287 : DIMENSION(:), OPTIONAL, POINTER :: text, atom_names_empty
288 : INTEGER, INTENT(OUT), OPTIONAL :: nunits_empty
289 : REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: coordinates_empty
290 : INTEGER, INTENT(OUT), OPTIONAL :: motion_row_end, motion_row_start
291 :
292 102 : IF (PRESENT(in_use)) in_use = mc_input_file%in_use
293 102 : IF (PRESENT(run_type_row)) run_type_row = mc_input_file%run_type_row
294 102 : IF (PRESENT(run_type_column)) run_type_column = mc_input_file%run_type_column
295 102 : IF (PRESENT(coord_row_start)) coord_row_start = mc_input_file%coord_row_start
296 102 : IF (PRESENT(coord_row_end)) coord_row_end = mc_input_file%coord_row_end
297 102 : IF (PRESENT(cell_row)) cell_row = mc_input_file%cell_row
298 102 : IF (PRESENT(cell_column)) cell_column = mc_input_file%cell_column
299 102 : IF (PRESENT(force_eval_row_start)) force_eval_row_start = mc_input_file%force_eval_row_start
300 102 : IF (PRESENT(force_eval_row_end)) force_eval_row_end = mc_input_file%force_eval_row_end
301 102 : IF (PRESENT(global_row_start)) global_row_start = mc_input_file%global_row_start
302 102 : IF (PRESENT(global_row_end)) global_row_end = mc_input_file%global_row_end
303 102 : IF (PRESENT(nunits_empty)) nunits_empty = mc_input_file%nunits_empty
304 102 : IF (PRESENT(motion_row_start)) motion_row_start = mc_input_file%motion_row_start
305 102 : IF (PRESENT(motion_row_end)) motion_row_end = mc_input_file%motion_row_end
306 :
307 102 : IF (PRESENT(mol_set_nmol_row)) mol_set_nmol_row => mc_input_file%mol_set_nmol_row
308 102 : IF (PRESENT(mol_set_nmol_column)) mol_set_nmol_column => mc_input_file%mol_set_nmol_column
309 102 : IF (PRESENT(text)) text => mc_input_file%text
310 102 : IF (PRESENT(atom_names_empty)) atom_names_empty => mc_input_file%atom_names_empty
311 102 : IF (PRESENT(coordinates_empty)) coordinates_empty => mc_input_file%coordinates_empty
312 :
313 102 : END SUBROUTINE GET_MC_INPUT_FILE
314 : ! **************************************************************************************************
315 : !> \brief ...
316 : !> \param mc_par ...
317 : !> \param nstep ...
318 : !> \param nvirial ...
319 : !> \param iuptrans ...
320 : !> \param iupcltrans ...
321 : !> \param iupvolume ...
322 : !> \param nmoves ...
323 : !> \param nswapmoves ...
324 : !> \param rm ...
325 : !> \param cl ...
326 : !> \param diff ...
327 : !> \param nstart ...
328 : !> \param source ...
329 : !> \param group ...
330 : !> \param lbias ...
331 : !> \param ionode ...
332 : !> \param lrestart ...
333 : !> \param lstop ...
334 : !> \param rmvolume ...
335 : !> \param rmcltrans ...
336 : !> \param rmbond ...
337 : !> \param rmangle ...
338 : !> \param rmrot ...
339 : !> \param rmtrans ...
340 : !> \param temperature ...
341 : !> \param pressure ...
342 : !> \param rclus ...
343 : !> \param BETA ...
344 : !> \param pmswap ...
345 : !> \param pmvolume ...
346 : !> \param pmtraion ...
347 : !> \param pmtrans ...
348 : !> \param pmcltrans ...
349 : !> \param ensemble ...
350 : !> \param PROGRAM ...
351 : !> \param restart_file_name ...
352 : !> \param molecules_file ...
353 : !> \param moves_file ...
354 : !> \param coords_file ...
355 : !> \param energy_file ...
356 : !> \param displacement_file ...
357 : !> \param cell_file ...
358 : !> \param dat_file ...
359 : !> \param data_file ...
360 : !> \param box2_file ...
361 : !> \param fft_lib ...
362 : !> \param iprint ...
363 : !> \param rcut ...
364 : !> \param ldiscrete ...
365 : !> \param discrete_step ...
366 : !> \param pmavbmc ...
367 : !> \param pbias ...
368 : !> \param avbmc_atom ...
369 : !> \param avbmc_rmin ...
370 : !> \param avbmc_rmax ...
371 : !> \param rmdihedral ...
372 : !> \param input_file ...
373 : !> \param mc_molecule_info ...
374 : !> \param pmswap_mol ...
375 : !> \param pmavbmc_mol ...
376 : !> \param pmtrans_mol ...
377 : !> \param pmrot_mol ...
378 : !> \param pmtraion_mol ...
379 : !> \param mc_input_file ...
380 : !> \param mc_bias_file ...
381 : !> \param pmvol_box ...
382 : !> \param pmclus_box ...
383 : !> \param virial_temps ...
384 : !> \param exp_min_val ...
385 : !> \param exp_max_val ...
386 : !> \param min_val ...
387 : !> \param max_val ...
388 : !> \param eta ...
389 : !> \param pmhmc ...
390 : !> \param pmhmc_box ...
391 : !> \param lhmc ...
392 : !> \param rand2skip ...
393 : ! **************************************************************************************************
394 2646 : SUBROUTINE get_mc_par(mc_par, nstep, nvirial, iuptrans, iupcltrans, iupvolume, &
395 : nmoves, nswapmoves, rm, cl, diff, nstart, &
396 : source, group, lbias, ionode, lrestart, lstop, rmvolume, rmcltrans, rmbond, rmangle, &
397 : rmrot, rmtrans, temperature, pressure, rclus, BETA, pmswap, pmvolume, pmtraion, pmtrans, &
398 : pmcltrans, ensemble, PROGRAM, restart_file_name, molecules_file, moves_file, coords_file, &
399 : energy_file, displacement_file, cell_file, dat_file, data_file, box2_file, &
400 : fft_lib, iprint, rcut, ldiscrete, discrete_step, &
401 : pmavbmc, pbias, avbmc_atom, avbmc_rmin, avbmc_rmax, rmdihedral, &
402 : input_file, mc_molecule_info, pmswap_mol, pmavbmc_mol, pmtrans_mol, pmrot_mol, &
403 : pmtraion_mol, mc_input_file, mc_bias_file, pmvol_box, pmclus_box, virial_temps, &
404 : exp_min_val, exp_max_val, min_val, max_val, eta, pmhmc, pmhmc_box, lhmc, rand2skip)
405 :
406 : ! **************************************************************************************************
407 : !> \brief accesses the private elements of the mc_parameters_type
408 : !> \param mc_par the structure mc parameters you want
409 : !>
410 : !> Suitable for parallel.
411 : !> \author MJM
412 : TYPE(mc_simpar_type), POINTER :: mc_par
413 : INTEGER, INTENT(OUT), OPTIONAL :: nstep, nvirial, iuptrans, iupcltrans, &
414 : iupvolume, nmoves, nswapmoves, rm, cl, &
415 : diff, nstart, source
416 : TYPE(mp_comm_type), INTENT(OUT), OPTIONAL :: group
417 : LOGICAL, INTENT(OUT), OPTIONAL :: lbias, ionode, lrestart, lstop
418 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: rmvolume, rmcltrans
419 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: rmbond, rmangle, rmrot, rmtrans
420 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: temperature, pressure, rclus, BETA, &
421 : pmswap, pmvolume, pmtraion, pmtrans, &
422 : pmcltrans
423 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: ensemble, PROGRAM, restart_file_name, &
424 : molecules_file, moves_file, coords_file, energy_file, displacement_file, cell_file, &
425 : dat_file, data_file, box2_file, fft_lib
426 : INTEGER, INTENT(OUT), OPTIONAL :: iprint
427 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: rcut
428 : LOGICAL, INTENT(OUT), OPTIONAL :: ldiscrete
429 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: discrete_step, pmavbmc
430 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: pbias
431 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: avbmc_atom
432 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: avbmc_rmin, avbmc_rmax, rmdihedral
433 : TYPE(section_vals_type), OPTIONAL, POINTER :: input_file
434 : TYPE(mc_molecule_info_type), OPTIONAL, POINTER :: mc_molecule_info
435 : REAL(dp), DIMENSION(:), OPTIONAL, POINTER :: pmswap_mol
436 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: pmavbmc_mol, pmtrans_mol, pmrot_mol, &
437 : pmtraion_mol
438 : TYPE(mc_input_file_type), OPTIONAL, POINTER :: mc_input_file, mc_bias_file
439 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: pmvol_box, pmclus_box
440 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: virial_temps
441 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: exp_min_val, exp_max_val, min_val, &
442 : max_val
443 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: eta
444 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: pmhmc, pmhmc_box
445 : LOGICAL, INTENT(OUT), OPTIONAL :: lhmc
446 : INTEGER, INTENT(OUT), OPTIONAL :: rand2skip
447 :
448 3295 : IF (PRESENT(nstep)) nstep = mc_par%nstep
449 3295 : IF (PRESENT(nvirial)) nvirial = mc_par%nvirial
450 3295 : IF (PRESENT(iuptrans)) iuptrans = mc_par%iuptrans
451 3295 : IF (PRESENT(iupcltrans)) iupcltrans = mc_par%iupcltrans
452 3295 : IF (PRESENT(iupvolume)) iupvolume = mc_par%iupvolume
453 3295 : IF (PRESENT(nmoves)) nmoves = mc_par%nmoves
454 3295 : IF (PRESENT(nswapmoves)) nswapmoves = mc_par%nswapmoves
455 3295 : IF (PRESENT(rm)) rm = mc_par%rm
456 3295 : IF (PRESENT(cl)) cl = mc_par%cl
457 3295 : IF (PRESENT(diff)) diff = mc_par%diff
458 3295 : IF (PRESENT(nstart)) nstart = mc_par%nstart
459 3295 : IF (PRESENT(source)) source = mc_par%source
460 3295 : IF (PRESENT(group)) group = mc_par%group
461 3295 : IF (PRESENT(iprint)) iprint = mc_par%iprint
462 :
463 3295 : IF (PRESENT(lbias)) lbias = mc_par%lbias
464 3295 : IF (PRESENT(ionode)) ionode = mc_par%ionode
465 3295 : IF (PRESENT(lrestart)) lrestart = mc_par%lrestart
466 3295 : IF (PRESENT(lstop)) lstop = mc_par%lstop
467 3295 : IF (PRESENT(ldiscrete)) ldiscrete = mc_par%ldiscrete
468 :
469 3295 : IF (PRESENT(virial_temps)) virial_temps => mc_par%virial_temps
470 3295 : IF (PRESENT(avbmc_atom)) avbmc_atom => mc_par%avbmc_atom
471 3295 : IF (PRESENT(avbmc_rmin)) avbmc_rmin => mc_par%avbmc_rmin
472 3295 : IF (PRESENT(avbmc_rmax)) avbmc_rmax => mc_par%avbmc_rmax
473 3295 : IF (PRESENT(rcut)) rcut = mc_par%rcut
474 3295 : IF (PRESENT(discrete_step)) discrete_step = mc_par%discrete_step
475 3295 : IF (PRESENT(rmvolume)) rmvolume = mc_par%rmvolume
476 3295 : IF (PRESENT(rmcltrans)) rmcltrans = mc_par%rmcltrans
477 3295 : IF (PRESENT(temperature)) temperature = mc_par%temperature
478 3295 : IF (PRESENT(pressure)) pressure = mc_par%pressure
479 3295 : IF (PRESENT(rclus)) rclus = mc_par%rclus
480 3295 : IF (PRESENT(BETA)) BETA = mc_par%BETA
481 3295 : IF (PRESENT(exp_max_val)) exp_max_val = mc_par%exp_max_val
482 3295 : IF (PRESENT(exp_min_val)) exp_min_val = mc_par%exp_min_val
483 3295 : IF (PRESENT(max_val)) max_val = mc_par%max_val
484 3295 : IF (PRESENT(min_val)) min_val = mc_par%min_val
485 3295 : IF (PRESENT(pmswap)) pmswap = mc_par%pmswap
486 3295 : IF (PRESENT(pmvolume)) pmvolume = mc_par%pmvolume
487 3295 : IF (PRESENT(lhmc)) lhmc = mc_par%lhmc
488 3295 : IF (PRESENT(pmhmc)) pmhmc = mc_par%pmhmc
489 3295 : IF (PRESENT(pmhmc_box)) pmhmc_box = mc_par%pmhmc_box
490 3295 : IF (PRESENT(pmvol_box)) pmvol_box = mc_par%pmvol_box
491 3295 : IF (PRESENT(pmclus_box)) pmclus_box = mc_par%pmclus_box
492 3295 : IF (PRESENT(pmtraion)) pmtraion = mc_par%pmtraion
493 3295 : IF (PRESENT(pmtrans)) pmtrans = mc_par%pmtrans
494 3295 : IF (PRESENT(pmcltrans)) pmcltrans = mc_par%pmcltrans
495 3295 : IF (PRESENT(pmavbmc)) pmavbmc = mc_par%pmavbmc
496 3295 : IF (PRESENT(pmrot_mol)) pmrot_mol => mc_par%pmrot_mol
497 3295 : IF (PRESENT(pmtrans_mol)) pmtrans_mol => mc_par%pmtrans_mol
498 3295 : IF (PRESENT(pmtraion_mol)) pmtraion_mol => mc_par%pmtraion_mol
499 3295 : IF (PRESENT(pmavbmc_mol)) pmavbmc_mol => mc_par%pmavbmc_mol
500 3295 : IF (PRESENT(pbias)) pbias => mc_par%pbias
501 :
502 3295 : IF (PRESENT(rmbond)) rmbond => mc_par%rmbond
503 3295 : IF (PRESENT(rmangle)) rmangle => mc_par%rmangle
504 3295 : IF (PRESENT(rmdihedral)) rmdihedral => mc_par%rmdihedral
505 3295 : IF (PRESENT(rmrot)) rmrot => mc_par%rmrot
506 3295 : IF (PRESENT(rmtrans)) rmtrans => mc_par%rmtrans
507 :
508 3295 : IF (PRESENT(eta)) eta => mc_par%eta
509 :
510 3295 : IF (PRESENT(ensemble)) ensemble = mc_par%ensemble
511 3295 : IF (PRESENT(PROGRAM)) PROGRAM = mc_par%program
512 3295 : IF (PRESENT(restart_file_name)) restart_file_name = &
513 32 : mc_par%restart_file_name
514 3295 : IF (PRESENT(moves_file)) moves_file = mc_par%moves_file
515 3295 : IF (PRESENT(coords_file)) coords_file = mc_par%coords_file
516 3295 : IF (PRESENT(molecules_file)) molecules_file = mc_par%molecules_file
517 3295 : IF (PRESENT(energy_file)) energy_file = mc_par%energy_file
518 3295 : IF (PRESENT(displacement_file)) displacement_file = &
519 24 : mc_par%displacement_file
520 3295 : IF (PRESENT(cell_file)) cell_file = mc_par%cell_file
521 3295 : IF (PRESENT(dat_file)) dat_file = mc_par%dat_file
522 3295 : IF (PRESENT(data_file)) data_file = mc_par%data_file
523 3295 : IF (PRESENT(box2_file)) box2_file = mc_par%box2_file
524 3295 : IF (PRESENT(fft_lib)) fft_lib = mc_par%fft_lib
525 :
526 3295 : IF (PRESENT(input_file)) input_file => mc_par%input_file
527 3295 : IF (PRESENT(mc_molecule_info)) mc_molecule_info => mc_par%mc_molecule_info
528 3295 : IF (PRESENT(mc_input_file)) mc_input_file => mc_par%mc_input_file
529 3295 : IF (PRESENT(mc_bias_file)) mc_bias_file => mc_par%mc_bias_file
530 :
531 3295 : IF (PRESENT(pmswap_mol)) pmswap_mol => mc_par%pmswap_mol
532 3295 : IF (PRESENT(rand2skip)) rand2skip = mc_par%rand2skip
533 3295 : END SUBROUTINE get_mc_par
534 :
535 : ! **************************************************************************************************
536 : !> \brief ...
537 : !> \param mc_molecule_info ...
538 : !> \param nmol_types ...
539 : !> \param nchain_total ...
540 : !> \param nboxes ...
541 : !> \param names ...
542 : !> \param conf_prob ...
543 : !> \param nchains ...
544 : !> \param nunits ...
545 : !> \param mol_type ...
546 : !> \param nunits_tot ...
547 : !> \param in_box ...
548 : !> \param atom_names ...
549 : !> \param mass ...
550 : ! **************************************************************************************************
551 3616 : SUBROUTINE get_mc_molecule_info(mc_molecule_info, nmol_types, nchain_total, nboxes, &
552 : names, conf_prob, nchains, nunits, mol_type, nunits_tot, in_box, atom_names, &
553 : mass)
554 :
555 : ! **************************************************************************************************
556 : !> \brief accesses the private elements of the mc_molecule_info_type
557 : !> \param mc_molecule_info the structure you want the parameters for
558 : !>
559 : !> Suitable for parallel.
560 : !> \author MJM
561 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
562 : INTEGER, INTENT(OUT), OPTIONAL :: nmol_types, nchain_total, nboxes
563 : CHARACTER(LEN=default_string_length), &
564 : DIMENSION(:), OPTIONAL, POINTER :: names
565 : REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: conf_prob
566 : INTEGER, DIMENSION(:, :), OPTIONAL, POINTER :: nchains
567 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: nunits, mol_type, nunits_tot, in_box
568 : CHARACTER(LEN=default_string_length), &
569 : DIMENSION(:, :), OPTIONAL, POINTER :: atom_names
570 : REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: mass
571 :
572 3616 : IF (PRESENT(nboxes)) nboxes = mc_molecule_info%nboxes
573 3616 : IF (PRESENT(nchain_total)) nchain_total = mc_molecule_info%nchain_total
574 3616 : IF (PRESENT(nmol_types)) nmol_types = mc_molecule_info%nmol_types
575 :
576 3616 : IF (PRESENT(names)) names => mc_molecule_info%names
577 3616 : IF (PRESENT(atom_names)) atom_names => mc_molecule_info%atom_names
578 :
579 3616 : IF (PRESENT(conf_prob)) conf_prob => mc_molecule_info%conf_prob
580 3616 : IF (PRESENT(mass)) mass => mc_molecule_info%mass
581 :
582 3616 : IF (PRESENT(nchains)) nchains => mc_molecule_info%nchains
583 3616 : IF (PRESENT(nunits)) nunits => mc_molecule_info%nunits
584 3616 : IF (PRESENT(mol_type)) mol_type => mc_molecule_info%mol_type
585 3616 : IF (PRESENT(nunits_tot)) nunits_tot => mc_molecule_info%nunits_tot
586 3616 : IF (PRESENT(in_box)) in_box => mc_molecule_info%in_box
587 :
588 3616 : END SUBROUTINE get_mc_molecule_info
589 :
590 : ! **************************************************************************************************
591 : !> \brief changes the private elements of the mc_parameters_type
592 : !> \param mc_par the structure mc parameters you want
593 : !> \param rm ...
594 : !> \param cl ...
595 : !> \param diff ...
596 : !> \param nstart ...
597 : !> \param rmvolume ...
598 : !> \param rmcltrans ...
599 : !> \param rmbond ...
600 : !> \param rmangle ...
601 : !> \param rmdihedral ...
602 : !> \param rmrot ...
603 : !> \param rmtrans ...
604 : !> \param PROGRAM ...
605 : !> \param nmoves ...
606 : !> \param nswapmoves ...
607 : !> \param lstop ...
608 : !> \param temperature ...
609 : !> \param pressure ...
610 : !> \param rclus ...
611 : !> \param iuptrans ...
612 : !> \param iupcltrans ...
613 : !> \param iupvolume ...
614 : !> \param pmswap ...
615 : !> \param pmvolume ...
616 : !> \param pmtraion ...
617 : !> \param pmtrans ...
618 : !> \param pmcltrans ...
619 : !> \param BETA ...
620 : !> \param rcut ...
621 : !> \param iprint ...
622 : !> \param lbias ...
623 : !> \param nstep ...
624 : !> \param lrestart ...
625 : !> \param ldiscrete ...
626 : !> \param discrete_step ...
627 : !> \param pmavbmc ...
628 : !> \param mc_molecule_info ...
629 : !> \param pmavbmc_mol ...
630 : !> \param pmtrans_mol ...
631 : !> \param pmrot_mol ...
632 : !> \param pmtraion_mol ...
633 : !> \param pmswap_mol ...
634 : !> \param avbmc_rmin ...
635 : !> \param avbmc_rmax ...
636 : !> \param avbmc_atom ...
637 : !> \param pbias ...
638 : !> \param ensemble ...
639 : !> \param pmvol_box ...
640 : !> \param pmclus_box ...
641 : !> \param eta ...
642 : !> \param mc_input_file ...
643 : !> \param mc_bias_file ...
644 : !> \param exp_max_val ...
645 : !> \param exp_min_val ...
646 : !> \param min_val ...
647 : !> \param max_val ...
648 : !> \param pmhmc ...
649 : !> \param pmhmc_box ...
650 : !> \param lhmc ...
651 : !> \param ionode ...
652 : !> \param source ...
653 : !> \param group ...
654 : !> \param rand2skip ...
655 : !> Suitable for parallel.
656 : !> \author MJM
657 : ! **************************************************************************************************
658 154 : SUBROUTINE set_mc_par(mc_par, rm, cl, &
659 : diff, nstart, rmvolume, rmcltrans, rmbond, rmangle, rmdihedral, rmrot, rmtrans, PROGRAM, &
660 : nmoves, nswapmoves, lstop, temperature, pressure, rclus, iuptrans, iupcltrans, iupvolume, &
661 : pmswap, pmvolume, pmtraion, pmtrans, pmcltrans, BETA, rcut, iprint, lbias, nstep, &
662 : lrestart, ldiscrete, discrete_step, pmavbmc, mc_molecule_info, &
663 : pmavbmc_mol, pmtrans_mol, pmrot_mol, pmtraion_mol, pmswap_mol, &
664 : avbmc_rmin, avbmc_rmax, avbmc_atom, pbias, ensemble, pmvol_box, pmclus_box, eta, &
665 : mc_input_file, mc_bias_file, exp_max_val, exp_min_val, min_val, max_val, &
666 : pmhmc, pmhmc_box, lhmc, ionode, source, group, rand2skip)
667 :
668 : TYPE(mc_simpar_type), POINTER :: mc_par
669 : INTEGER, INTENT(IN), OPTIONAL :: rm, cl, diff, nstart
670 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: rmvolume, rmcltrans
671 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: rmbond, rmangle, rmdihedral, rmrot, &
672 : rmtrans
673 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: PROGRAM
674 : INTEGER, INTENT(IN), OPTIONAL :: nmoves, nswapmoves
675 : LOGICAL, INTENT(IN), OPTIONAL :: lstop
676 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: temperature, pressure, rclus
677 : INTEGER, INTENT(IN), OPTIONAL :: iuptrans, iupcltrans, iupvolume
678 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: pmswap, pmvolume, pmtraion, pmtrans, &
679 : pmcltrans, BETA, rcut
680 : INTEGER, INTENT(IN), OPTIONAL :: iprint
681 : LOGICAL, INTENT(IN), OPTIONAL :: lbias
682 : INTEGER, INTENT(IN), OPTIONAL :: nstep
683 : LOGICAL, INTENT(IN), OPTIONAL :: lrestart, ldiscrete
684 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: discrete_step, pmavbmc
685 : TYPE(mc_molecule_info_type), OPTIONAL, POINTER :: mc_molecule_info
686 : REAL(dp), DIMENSION(:), OPTIONAL, POINTER :: pmavbmc_mol, pmtrans_mol, pmrot_mol, &
687 : pmtraion_mol, pmswap_mol, avbmc_rmin, &
688 : avbmc_rmax
689 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: avbmc_atom
690 : REAL(dp), DIMENSION(:), OPTIONAL, POINTER :: pbias
691 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: ensemble
692 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: pmvol_box, pmclus_box
693 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: eta
694 : TYPE(mc_input_file_type), OPTIONAL, POINTER :: mc_input_file, mc_bias_file
695 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: exp_max_val, exp_min_val, min_val, &
696 : max_val, pmhmc, pmhmc_box
697 : LOGICAL, INTENT(IN), OPTIONAL :: lhmc, ionode
698 : INTEGER, INTENT(IN), OPTIONAL :: source
699 :
700 : CLASS(mp_comm_type), INTENT(IN), OPTIONAL :: group
701 : INTEGER, INTENT(IN), OPTIONAL :: rand2skip
702 :
703 : INTEGER :: iparm
704 :
705 : ! These are the only values that change during the course of the simulation
706 : ! or are computed outside of this module
707 :
708 8 : IF (PRESENT(nstep)) mc_par%nstep = nstep
709 154 : IF (PRESENT(rm)) mc_par%rm = rm
710 154 : IF (PRESENT(cl)) mc_par%cl = cl
711 154 : IF (PRESENT(diff)) mc_par%diff = diff
712 154 : IF (PRESENT(nstart)) mc_par%nstart = nstart
713 154 : IF (PRESENT(nmoves)) mc_par%nmoves = nmoves
714 154 : IF (PRESENT(nswapmoves)) mc_par%nswapmoves = nswapmoves
715 154 : IF (PRESENT(iprint)) mc_par%iprint = iprint
716 154 : IF (PRESENT(iuptrans)) mc_par%iuptrans = iuptrans
717 154 : IF (PRESENT(iupcltrans)) mc_par%iupcltrans = iupcltrans
718 154 : IF (PRESENT(iupvolume)) mc_par%iupvolume = iupvolume
719 :
720 154 : IF (PRESENT(ldiscrete)) mc_par%ldiscrete = ldiscrete
721 154 : IF (PRESENT(lstop)) mc_par%lstop = lstop
722 154 : IF (PRESENT(lbias)) mc_par%lbias = lbias
723 154 : IF (PRESENT(lrestart)) mc_par%lrestart = lrestart
724 :
725 154 : IF (PRESENT(exp_max_val)) mc_par%exp_max_val = exp_max_val
726 154 : IF (PRESENT(exp_min_val)) mc_par%exp_min_val = exp_min_val
727 154 : IF (PRESENT(max_val)) mc_par%max_val = max_val
728 154 : IF (PRESENT(min_val)) mc_par%min_val = min_val
729 154 : IF (PRESENT(BETA)) mc_par%BETA = BETA
730 154 : IF (PRESENT(temperature)) mc_par%temperature = temperature
731 154 : IF (PRESENT(rcut)) mc_par%rcut = rcut
732 154 : IF (PRESENT(pressure)) mc_par%pressure = pressure
733 154 : IF (PRESENT(rclus)) mc_par%rclus = rclus
734 154 : IF (PRESENT(pmvolume)) mc_par%pmvolume = pmvolume
735 154 : IF (PRESENT(lhmc)) mc_par%lhmc = lhmc
736 154 : IF (PRESENT(pmhmc)) mc_par%pmhmc = pmhmc
737 154 : IF (PRESENT(pmhmc_box)) mc_par%pmhmc_box = pmhmc_box
738 154 : IF (PRESENT(pmvol_box)) mc_par%pmvol_box = pmvol_box
739 154 : IF (PRESENT(pmclus_box)) mc_par%pmclus_box = pmclus_box
740 154 : IF (PRESENT(pmswap)) mc_par%pmswap = pmswap
741 154 : IF (PRESENT(pmtrans)) mc_par%pmtrans = pmtrans
742 154 : IF (PRESENT(pmcltrans)) mc_par%pmcltrans = pmcltrans
743 154 : IF (PRESENT(pmtraion)) mc_par%pmtraion = pmtraion
744 154 : IF (PRESENT(pmavbmc)) mc_par%pmavbmc = pmavbmc
745 :
746 154 : IF (PRESENT(discrete_step)) mc_par%discrete_step = discrete_step
747 154 : IF (PRESENT(rmvolume)) mc_par%rmvolume = rmvolume
748 154 : IF (PRESENT(rmcltrans)) mc_par%rmcltrans = rmcltrans
749 :
750 154 : IF (PRESENT(mc_input_file)) mc_par%mc_input_file => mc_input_file
751 154 : IF (PRESENT(mc_bias_file)) mc_par%mc_bias_file => mc_bias_file
752 :
753 : ! cannot just change the pointers, because then we have memory leaks
754 : ! and the program crashes if we try to release it (in certain cases)
755 154 : IF (PRESENT(eta)) THEN
756 0 : DO iparm = 1, SIZE(eta)
757 0 : mc_par%eta(iparm) = eta(iparm)
758 : END DO
759 : END IF
760 154 : IF (PRESENT(rmbond)) THEN
761 0 : DO iparm = 1, SIZE(rmbond)
762 0 : mc_par%rmbond(iparm) = rmbond(iparm)
763 : END DO
764 : END IF
765 154 : IF (PRESENT(rmangle)) THEN
766 0 : DO iparm = 1, SIZE(rmangle)
767 0 : mc_par%rmangle(iparm) = rmangle(iparm)
768 : END DO
769 : END IF
770 154 : IF (PRESENT(rmdihedral)) THEN
771 0 : DO iparm = 1, SIZE(rmdihedral)
772 0 : mc_par%rmdihedral(iparm) = rmdihedral(iparm)
773 : END DO
774 : END IF
775 154 : IF (PRESENT(rmrot)) THEN
776 0 : DO iparm = 1, SIZE(rmrot)
777 0 : mc_par%rmrot(iparm) = rmrot(iparm)
778 : END DO
779 : END IF
780 154 : IF (PRESENT(rmtrans)) THEN
781 0 : DO iparm = 1, SIZE(rmtrans)
782 0 : mc_par%rmtrans(iparm) = rmtrans(iparm)
783 : END DO
784 : END IF
785 :
786 154 : IF (PRESENT(pmavbmc_mol)) THEN
787 18 : DO iparm = 1, SIZE(pmavbmc_mol)
788 18 : mc_par%pmavbmc_mol(iparm) = pmavbmc_mol(iparm)
789 : END DO
790 : END IF
791 154 : IF (PRESENT(pmtrans_mol)) THEN
792 18 : DO iparm = 1, SIZE(pmtrans_mol)
793 18 : mc_par%pmtrans_mol(iparm) = pmtrans_mol(iparm)
794 : END DO
795 : END IF
796 154 : IF (PRESENT(pmrot_mol)) THEN
797 18 : DO iparm = 1, SIZE(pmrot_mol)
798 18 : mc_par%pmrot_mol(iparm) = pmrot_mol(iparm)
799 : END DO
800 : END IF
801 154 : IF (PRESENT(pmtraion_mol)) THEN
802 18 : DO iparm = 1, SIZE(pmtraion_mol)
803 18 : mc_par%pmtraion_mol(iparm) = pmtraion_mol(iparm)
804 : END DO
805 : END IF
806 154 : IF (PRESENT(pmswap_mol)) THEN
807 18 : DO iparm = 1, SIZE(pmswap_mol)
808 18 : mc_par%pmswap_mol(iparm) = pmswap_mol(iparm)
809 : END DO
810 : END IF
811 :
812 154 : IF (PRESENT(avbmc_atom)) THEN
813 18 : DO iparm = 1, SIZE(avbmc_atom)
814 18 : mc_par%avbmc_atom(iparm) = avbmc_atom(iparm)
815 : END DO
816 : END IF
817 154 : IF (PRESENT(avbmc_rmin)) THEN
818 18 : DO iparm = 1, SIZE(avbmc_rmin)
819 18 : mc_par%avbmc_rmin(iparm) = avbmc_rmin(iparm)
820 : END DO
821 : END IF
822 154 : IF (PRESENT(avbmc_rmax)) THEN
823 18 : DO iparm = 1, SIZE(avbmc_rmax)
824 18 : mc_par%avbmc_rmax(iparm) = avbmc_rmax(iparm)
825 : END DO
826 : END IF
827 154 : IF (PRESENT(pbias)) THEN
828 18 : DO iparm = 1, SIZE(pbias)
829 18 : mc_par%pbias(iparm) = pbias(iparm)
830 : END DO
831 : END IF
832 :
833 154 : IF (PRESENT(PROGRAM)) mc_par%program = PROGRAM
834 154 : IF (PRESENT(ensemble)) mc_par%ensemble = ensemble
835 :
836 154 : IF (PRESENT(mc_molecule_info)) mc_par%mc_molecule_info => mc_molecule_info
837 154 : IF (PRESENT(source)) mc_par%source = source
838 154 : IF (PRESENT(group)) mc_par%group = group
839 154 : IF (PRESENT(ionode)) mc_par%ionode = ionode
840 :
841 154 : IF (PRESENT(rand2skip)) mc_par%rand2skip = rand2skip
842 154 : END SUBROUTINE set_mc_par
843 :
844 : ! **************************************************************************************************
845 : !> \brief creates (allocates) the mc_simulation_parameters type
846 : !> \param mc_par the structure that will be created (allocated)
847 : !> \param nmol_types the number of molecule types in the system
848 : !> \author MJM
849 : ! **************************************************************************************************
850 26 : SUBROUTINE mc_sim_par_create(mc_par, nmol_types)
851 :
852 : TYPE(mc_simpar_type), POINTER :: mc_par
853 : INTEGER, INTENT(IN) :: nmol_types
854 :
855 : NULLIFY (mc_par)
856 :
857 26 : ALLOCATE (mc_par)
858 78 : ALLOCATE (mc_par%pmtraion_mol(1:nmol_types))
859 78 : ALLOCATE (mc_par%pmtrans_mol(1:nmol_types))
860 78 : ALLOCATE (mc_par%pmrot_mol(1:nmol_types))
861 78 : ALLOCATE (mc_par%pmswap_mol(1:nmol_types))
862 :
863 78 : ALLOCATE (mc_par%eta(1:nmol_types))
864 :
865 78 : ALLOCATE (mc_par%rmbond(1:nmol_types))
866 78 : ALLOCATE (mc_par%rmangle(1:nmol_types))
867 78 : ALLOCATE (mc_par%rmdihedral(1:nmol_types))
868 78 : ALLOCATE (mc_par%rmtrans(1:nmol_types))
869 78 : ALLOCATE (mc_par%rmrot(1:nmol_types))
870 :
871 78 : ALLOCATE (mc_par%avbmc_atom(1:nmol_types))
872 78 : ALLOCATE (mc_par%avbmc_rmin(1:nmol_types))
873 78 : ALLOCATE (mc_par%avbmc_rmax(1:nmol_types))
874 78 : ALLOCATE (mc_par%pmavbmc_mol(1:nmol_types))
875 78 : ALLOCATE (mc_par%pbias(1:nmol_types))
876 :
877 26 : ALLOCATE (mc_par%mc_input_file)
878 26 : ALLOCATE (mc_par%mc_bias_file)
879 :
880 26 : END SUBROUTINE mc_sim_par_create
881 :
882 : ! **************************************************************************************************
883 : !> \brief destroys (deallocates) the mc_simulation_parameters type
884 : !> \param mc_par the structure that will be destroyed
885 : !> \author MJM
886 : ! **************************************************************************************************
887 26 : SUBROUTINE mc_sim_par_destroy(mc_par)
888 :
889 : TYPE(mc_simpar_type), POINTER :: mc_par
890 :
891 26 : DEALLOCATE (mc_par%mc_input_file)
892 26 : DEALLOCATE (mc_par%mc_bias_file)
893 :
894 26 : DEALLOCATE (mc_par%pmtraion_mol)
895 26 : DEALLOCATE (mc_par%pmtrans_mol)
896 26 : DEALLOCATE (mc_par%pmrot_mol)
897 26 : DEALLOCATE (mc_par%pmswap_mol)
898 :
899 26 : DEALLOCATE (mc_par%eta)
900 :
901 26 : DEALLOCATE (mc_par%rmbond)
902 26 : DEALLOCATE (mc_par%rmangle)
903 26 : DEALLOCATE (mc_par%rmdihedral)
904 26 : DEALLOCATE (mc_par%rmtrans)
905 26 : DEALLOCATE (mc_par%rmrot)
906 :
907 26 : DEALLOCATE (mc_par%avbmc_atom)
908 26 : DEALLOCATE (mc_par%avbmc_rmin)
909 26 : DEALLOCATE (mc_par%avbmc_rmax)
910 26 : DEALLOCATE (mc_par%pmavbmc_mol)
911 26 : DEALLOCATE (mc_par%pbias)
912 26 : IF (mc_par%ensemble == "VIRIAL") THEN
913 2 : DEALLOCATE (mc_par%virial_temps)
914 : END IF
915 :
916 26 : DEALLOCATE (mc_par)
917 : NULLIFY (mc_par)
918 :
919 26 : END SUBROUTINE mc_sim_par_destroy
920 :
921 : ! **************************************************************************************************
922 : !> \brief creates (allocates) the mc_input_file_type
923 : !> \param mc_input_file the structure that will be created (allocated)
924 : !> \param input_file_name the name of the file to read
925 : !> \param mc_molecule_info ...
926 : !> \param empty_coords ...
927 : !> \param lhmc ...
928 : !> \author MJM
929 : ! **************************************************************************************************
930 40 : SUBROUTINE mc_input_file_create(mc_input_file, input_file_name, &
931 40 : mc_molecule_info, empty_coords, lhmc)
932 :
933 : TYPE(mc_input_file_type), POINTER :: mc_input_file
934 : CHARACTER(LEN=*), INTENT(IN) :: input_file_name
935 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
936 : REAL(dp), DIMENSION(:, :) :: empty_coords
937 : LOGICAL, INTENT(IN) :: lhmc
938 :
939 : CHARACTER(default_string_length) :: cdum, line, method_name
940 : CHARACTER(default_string_length), &
941 40 : DIMENSION(:, :), POINTER :: atom_names
942 : INTEGER :: abc_column, abc_row, cell_column, cell_row, idum, iline, io_stat, irep, iunit, &
943 : iw, line_number, nlines, nmol_types, nstart, row_number, unit
944 40 : INTEGER, DIMENSION(:), POINTER :: nunits
945 :
946 : ! some stuff in case we need to write error messages to the screen
947 :
948 80 : iw = cp_logger_get_default_io_unit()
949 :
950 : ! allocate the array we'll need in case we have an empty box
951 : CALL get_mc_molecule_info(mc_molecule_info, nmol_types=nmol_types, &
952 40 : nunits=nunits, atom_names=atom_names)
953 120 : ALLOCATE (mc_input_file%coordinates_empty(1:3, 1:nunits(1)))
954 120 : ALLOCATE (mc_input_file%atom_names_empty(1:nunits(1)))
955 154 : DO iunit = 1, nunits(1)
956 114 : mc_input_file%atom_names_empty(iunit) = atom_names(iunit, 1)
957 496 : mc_input_file%coordinates_empty(1:3, iunit) = empty_coords(1:3, iunit)
958 : END DO
959 40 : mc_input_file%nunits_empty = nunits(1)
960 :
961 : ! allocate some arrays
962 120 : ALLOCATE (mc_input_file%mol_set_nmol_row(1:nmol_types))
963 80 : ALLOCATE (mc_input_file%mol_set_nmol_column(1:nmol_types))
964 :
965 : ! now read in and store the input file, first finding out how many lines it is
966 : CALL open_file(file_name=input_file_name, &
967 40 : unit_number=unit, file_action='READ', file_status='OLD')
968 :
969 40 : nlines = 0
970 9240 : DO
971 9280 : READ (unit, '(A)', IOSTAT=io_stat) line
972 9280 : IF (io_stat /= 0) EXIT
973 9240 : nlines = nlines + 1
974 : END DO
975 :
976 120 : ALLOCATE (mc_input_file%text(1:nlines))
977 :
978 40 : REWIND (unit)
979 9280 : DO iline = 1, nlines
980 9280 : READ (unit, '(A)') mc_input_file%text(iline)
981 : END DO
982 :
983 : ! close the input file
984 40 : CALL close_file(unit_number=unit)
985 :
986 : ! now we need to find the row and column numbers of a variety of things
987 : ! first, Let's find the first END after GLOBAL
988 : CALL mc_parse_text(mc_input_file%text, 1, nlines, "&GLOBAL", .TRUE., &
989 40 : mc_input_file%global_row_end, idum, start_row_number=mc_input_file%global_row_start)
990 40 : IF (mc_input_file%global_row_end == 0) THEN
991 0 : IF (iw > 0) THEN
992 0 : WRITE (iw, *)
993 0 : WRITE (iw, *) 'File name ', input_file_name
994 : END IF
995 : CALL cp_abort(__LOCATION__, &
996 0 : 'Could not find &END after &GLOBAL (make sure & is in the same column) and last line is blank')
997 : END IF
998 :
999 : ! Let's find the first END after MOTION...we might need this whole section
1000 : ! because of hybrid MD/MC moves, which requires the MD information to always
1001 : ! be attached to the force env
1002 : CALL mc_parse_text(mc_input_file%text, 1, nlines, "&MOTION", .TRUE., &
1003 40 : mc_input_file%motion_row_end, idum, start_row_number=mc_input_file%motion_row_start)
1004 40 : IF (mc_input_file%motion_row_end == 0) THEN
1005 0 : IF (iw > 0) THEN
1006 0 : WRITE (iw, *)
1007 0 : WRITE (iw, *) 'File name ', input_file_name, mc_input_file%motion_row_start
1008 : END IF
1009 : CALL cp_abort(__LOCATION__, &
1010 0 : 'Could not find &END after &MOTION (make sure & is in the same column) and last line is blank')
1011 : END IF
1012 :
1013 : ! find &FORCE_EVAL sections
1014 : CALL mc_parse_text(mc_input_file%text, 1, nlines, "&FORCE_EVAL", .TRUE., &
1015 40 : mc_input_file%force_eval_row_end, idum, start_row_number=mc_input_file%force_eval_row_start)
1016 :
1017 : ! first, let's find the first END after &COORD, and the line of &COORD
1018 : CALL mc_parse_text(mc_input_file%text, 1, nlines, "&COORD", .TRUE., &
1019 40 : mc_input_file%coord_row_end, idum)
1020 : CALL mc_parse_text(mc_input_file%text, 1, nlines, "&COORD", .FALSE., &
1021 40 : mc_input_file%coord_row_start, idum)
1022 :
1023 40 : IF (mc_input_file%coord_row_end == 0) THEN
1024 0 : IF (iw > 0) THEN
1025 0 : WRITE (iw, *)
1026 0 : WRITE (iw, *) 'File name ', input_file_name
1027 : END IF
1028 : CALL cp_abort(__LOCATION__, &
1029 0 : 'Could not find &END after &COORD (make sure & is the first in the same column after &COORD)')
1030 : END IF
1031 :
1032 : ! now the RUN_TYPE
1033 : CALL mc_parse_text(mc_input_file%text, 1, nlines, "RUN_TYPE", .FALSE., &
1034 40 : mc_input_file%run_type_row, mc_input_file%run_type_column)
1035 40 : mc_input_file%run_type_column = mc_input_file%run_type_column + 9
1036 40 : IF (mc_input_file%run_type_row == 0) THEN
1037 0 : IF (iw > 0) THEN
1038 0 : WRITE (iw, *)
1039 0 : WRITE (iw, *) 'File name ', input_file_name
1040 : END IF
1041 : CALL cp_abort(__LOCATION__, &
1042 0 : 'Could not find RUN_TYPE in the input file (should be in the &GLOBAL section).')
1043 : END IF
1044 :
1045 : ! now the CELL stuff..we don't care about REF_CELL since that should
1046 : ! never change if it's there
1047 : CALL mc_parse_text(mc_input_file%text, 1, nlines, "&CELL", .FALSE., &
1048 40 : mc_input_file%cell_row, mc_input_file%cell_column)
1049 : ! now find the ABC input line after CELL
1050 : CALL mc_parse_text(mc_input_file%text, mc_input_file%cell_row + 1, nlines, &
1051 40 : "ABC", .FALSE., abc_row, abc_column)
1052 : ! is there a &CELL inbetween? If so, that ABC will be for the ref_cell
1053 : ! and we need to find the next one
1054 : CALL mc_parse_text(mc_input_file%text, mc_input_file%cell_row + 1, abc_row, &
1055 40 : "&CELL", .FALSE., cell_row, cell_column)
1056 40 : IF (cell_row == 0) THEN ! nothing in between...we found the correct ABC
1057 40 : mc_input_file%cell_row = abc_row
1058 40 : mc_input_file%cell_column = abc_column + 4
1059 : ELSE
1060 : CALL mc_parse_text(mc_input_file%text, abc_row + 1, nlines, &
1061 0 : "ABC", .FALSE., mc_input_file%cell_row, mc_input_file%cell_column)
1062 : END IF
1063 40 : IF (mc_input_file%cell_row == 0) THEN
1064 0 : IF (iw > 0) THEN
1065 0 : WRITE (iw, *)
1066 0 : WRITE (iw, *) 'File name ', input_file_name
1067 : END IF
1068 : CALL cp_abort(__LOCATION__, &
1069 0 : 'Could not find &CELL section in the input file. Or could not find the ABC line after it.')
1070 : END IF
1071 :
1072 : ! now we need all the MOL_SET NMOLS indices
1073 40 : IF (.NOT. lhmc) THEN
1074 38 : irep = 0
1075 38 : nstart = 1
1076 68 : DO
1077 : CALL mc_parse_text(mc_input_file%text, nstart, nlines, "&MOLECULE", &
1078 106 : .FALSE., row_number, idum)
1079 106 : IF (row_number == 0) EXIT
1080 68 : nstart = row_number + 1
1081 68 : irep = irep + 1
1082 : CALL mc_parse_text(mc_input_file%text, nstart, nlines, "NMOL", &
1083 : .FALSE., mc_input_file%mol_set_nmol_row(irep), &
1084 68 : mc_input_file%mol_set_nmol_column(irep))
1085 68 : mc_input_file%mol_set_nmol_column(irep) = mc_input_file%mol_set_nmol_column(irep) + 5
1086 :
1087 : END DO
1088 38 : IF (irep /= nmol_types) THEN
1089 0 : IF (iw > 0) THEN
1090 0 : WRITE (iw, *)
1091 0 : WRITE (iw, *) 'File name ', input_file_name
1092 0 : WRITE (iw, *) 'Number of molecule types ', nmol_types
1093 0 : WRITE (iw, *) '&MOLECULE sections found ', irep
1094 : END IF
1095 : CALL cp_abort( &
1096 : __LOCATION__, &
1097 0 : 'Did not find MOLECULE sections for every molecule in the simulation (make sure both input files have all types)')
1098 : END IF
1099 : END IF
1100 :
1101 : ! now let's find the type of force_env...right now, I'm only concerned with
1102 : ! QS, and Fist, though it should be trivial to add others
1103 : CALL mc_parse_text(mc_input_file%text, 1, nlines, &
1104 40 : "METHOD", .FALSE., line_number, idum)
1105 40 : READ (mc_input_file%text(line_number), *) cdum, method_name
1106 40 : CALL uppercase(method_name)
1107 80 : SELECT CASE (TRIM(ADJUSTL(method_name)))
1108 : CASE ("FIST")
1109 34 : mc_input_file%in_use = use_fist_force
1110 : CASE ("QS", "QUICKSTEP")
1111 6 : mc_input_file%in_use = use_qs_force
1112 : CASE default
1113 80 : CPABORT('Cannot determine the type of force_env we are using (check METHOD)')
1114 : END SELECT
1115 :
1116 40 : END SUBROUTINE mc_input_file_create
1117 :
1118 : ! **************************************************************************************************
1119 : !> \brief destroys (deallocates) things in the mc_input_file_type
1120 : !> \param mc_input_file the structure that will be released (deallocated)
1121 : !> \author MJM
1122 : ! **************************************************************************************************
1123 40 : SUBROUTINE mc_input_file_destroy(mc_input_file)
1124 :
1125 : TYPE(mc_input_file_type), POINTER :: mc_input_file
1126 :
1127 40 : DEALLOCATE (mc_input_file%mol_set_nmol_row)
1128 40 : DEALLOCATE (mc_input_file%mol_set_nmol_column)
1129 40 : DEALLOCATE (mc_input_file%text)
1130 40 : DEALLOCATE (mc_input_file%atom_names_empty)
1131 40 : DEALLOCATE (mc_input_file%coordinates_empty)
1132 :
1133 40 : END SUBROUTINE mc_input_file_destroy
1134 :
1135 : ! **************************************************************************************************
1136 : !> \brief a basic text parser used to find the row and column numbers of various strings
1137 : !> in the input file, to store as indices for when we create a dat_file...
1138 : !> returns 0 for the row if nothing is found
1139 : !> \param text the text to parse
1140 : !> \param nstart the line to start searching from
1141 : !> \param nend the line to end searching
1142 : !> \param string_search the text we're looking for
1143 : !> \param lend if .TRUE., find the &END that comes after string_search...assumes that
1144 : !> the row is the first row with & in the same column as the search string
1145 : !> \param row_number the row the string is first found on
1146 : !> \param column_number the column number that the string first appears on
1147 : !> \param start_row_number ...
1148 : !> \author MJM
1149 : ! **************************************************************************************************
1150 574 : SUBROUTINE mc_parse_text(text, nstart, nend, string_search, lend, &
1151 : row_number, column_number, start_row_number)
1152 :
1153 : CHARACTER(LEN=*), DIMENSION(:), INTENT(IN) :: text
1154 : INTEGER, INTENT(IN) :: nstart, nend
1155 : CHARACTER(LEN=*), INTENT(IN) :: string_search
1156 : LOGICAL, INTENT(IN) :: lend
1157 : INTEGER, INTENT(OUT) :: row_number, column_number
1158 : INTEGER, INTENT(OUT), OPTIONAL :: start_row_number
1159 :
1160 : CHARACTER(default_string_length) :: text_temp
1161 : INTEGER :: iline, index_string, jline, string_size
1162 :
1163 : ! did not see any string utilities in the code to help, so here I go
1164 :
1165 574 : row_number = 0
1166 574 : column_number = 0
1167 :
1168 : ! how long is our string to search for?
1169 574 : string_size = LEN_TRIM(string_search)
1170 29728 : whole_file: DO iline = nstart, nend
1171 :
1172 29650 : index_string = -1
1173 29650 : index_string = INDEX(text(iline), string_search(1:string_size))
1174 :
1175 29728 : IF (index_string > 0) THEN ! we found one
1176 496 : IF (.NOT. lend) THEN
1177 336 : row_number = iline
1178 336 : column_number = index_string
1179 : ELSE
1180 160 : IF (PRESENT(start_row_number)) start_row_number = iline
1181 160 : column_number = index_string
1182 12034 : DO jline = iline + 1, nend
1183 : ! now we find the &END that matches up with this one...
1184 : ! I need proper indentation because I'm not very smart
1185 12034 : text_temp = text(jline)
1186 12034 : IF (text_temp(column_number:column_number) == "&") THEN
1187 160 : row_number = jline
1188 160 : EXIT
1189 : END IF
1190 : END DO
1191 : END IF
1192 : EXIT whole_file
1193 : END IF
1194 : END DO whole_file
1195 :
1196 574 : END SUBROUTINE mc_parse_text
1197 :
1198 : ! **************************************************************************************************
1199 : !> \brief reads in the Monte Carlo simulation parameters from an input file
1200 : !> \param mc_par the structure that will store the parameters
1201 : !> \param para_env ...
1202 : !> \param globenv the global environment for the simulation
1203 : !> \param input_file_name the name of the input_file
1204 : !> \param input_file the structure that contains all the keywords in the input file
1205 : !> \param force_env_section used to grab the type of force_env
1206 : !> \author MJM
1207 : ! **************************************************************************************************
1208 78 : SUBROUTINE read_mc_section(mc_par, para_env, globenv, input_file_name, input_file, force_env_section)
1209 :
1210 : TYPE(mc_simpar_type), POINTER :: mc_par
1211 : TYPE(mp_para_env_type), POINTER :: para_env
1212 : TYPE(global_environment_type), POINTER :: globenv
1213 : CHARACTER(LEN=*), INTENT(IN) :: input_file_name
1214 : TYPE(section_vals_type), POINTER :: input_file, force_env_section
1215 :
1216 : CHARACTER(len=*), PARAMETER :: routineN = 'read_mc_section'
1217 :
1218 : CHARACTER(LEN=5) :: box_string, molecule_string, &
1219 : tab_box_string, tab_string, &
1220 : tab_string_int
1221 : CHARACTER(LEN=default_string_length) :: format_box_string, format_string, &
1222 : format_string_int
1223 : INTEGER :: handle, ia, ie, imol, itype, iw, &
1224 : method_name_id, nmol_types, stop_num
1225 26 : INTEGER, DIMENSION(:), POINTER :: temp_i_array
1226 26 : REAL(dp), DIMENSION(:), POINTER :: temp_r_array
1227 : TYPE(section_vals_type), POINTER :: mc_section
1228 :
1229 : ! begin the timing of the subroutine
1230 :
1231 0 : CPASSERT(ASSOCIATED(input_file))
1232 :
1233 26 : CALL timeset(routineN, handle)
1234 :
1235 26 : NULLIFY (mc_section)
1236 : mc_section => section_vals_get_subs_vals(input_file, &
1237 26 : "MOTION%MC")
1238 :
1239 : ! need the input file sturcutre that we're reading from for when we make
1240 : ! dat files
1241 26 : mc_par%input_file => input_file
1242 :
1243 : ! set the ionode and mepos
1244 26 : mc_par%ionode = para_env%is_source()
1245 26 : mc_par%group = para_env
1246 26 : mc_par%source = para_env%source
1247 :
1248 : ! for convenience
1249 26 : nmol_types = mc_par%mc_molecule_info%nmol_types
1250 :
1251 : !..defaults...most are set in input_cp2k_motion
1252 26 : mc_par%nstart = 0
1253 26 : CALL section_vals_val_get(force_env_section, "METHOD", i_val=method_name_id)
1254 20 : SELECT CASE (method_name_id)
1255 : CASE (do_fist)
1256 20 : mc_par%iprint = 100
1257 : CASE (do_qs)
1258 26 : mc_par%iprint = 1
1259 : END SELECT
1260 :
1261 26 : iw = cp_logger_get_default_io_unit()
1262 26 : IF (iw > 0) WRITE (iw, *)
1263 :
1264 : !..filenames
1265 26 : mc_par%program = input_file_name
1266 26 : CALL xstring(mc_par%program, ia, ie)
1267 26 : mc_par%coords_file = mc_par%program(ia:ie)//'.coordinates'
1268 26 : mc_par%molecules_file = mc_par%program(ia:ie)//'.molecules'
1269 26 : mc_par%moves_file = mc_par%program(ia:ie)//'.moves'
1270 26 : mc_par%energy_file = mc_par%program(ia:ie)//'.energy'
1271 26 : mc_par%cell_file = mc_par%program(ia:ie)//'.cell'
1272 : mc_par%displacement_file = mc_par%program(ia:ie) &
1273 26 : //'.max_displacements'
1274 26 : mc_par%data_file = mc_par%program(ia:ie)//'.data'
1275 26 : stop_num = ie - 3
1276 26 : mc_par%dat_file = mc_par%program(ia:stop_num)//'dat'
1277 :
1278 : ! set them into the input parameter structure as the new defaults
1279 : CALL section_vals_val_set(mc_section, "COORDINATE_FILE_NAME", &
1280 26 : c_val=mc_par%coords_file)
1281 : CALL section_vals_val_set(mc_section, "DATA_FILE_NAME", &
1282 26 : c_val=mc_par%data_file)
1283 : CALL section_vals_val_set(mc_section, "CELL_FILE_NAME", &
1284 26 : c_val=mc_par%cell_file)
1285 : CALL section_vals_val_set(mc_section, "MAX_DISP_FILE_NAME", &
1286 26 : c_val=mc_par%displacement_file)
1287 : CALL section_vals_val_set(mc_section, "MOVES_FILE_NAME", &
1288 26 : c_val=mc_par%moves_file)
1289 : CALL section_vals_val_set(mc_section, "MOLECULES_FILE_NAME", &
1290 26 : c_val=mc_par%molecules_file)
1291 : CALL section_vals_val_set(mc_section, "ENERGY_FILE_NAME", &
1292 26 : c_val=mc_par%energy_file)
1293 :
1294 : ! grab the FFT library name and print level...this is used for writing the dat file
1295 : ! and hopefully will be changed
1296 26 : mc_par%fft_lib = globenv%default_fft_library
1297 :
1298 : ! get all the move probabilities first...if we are not doing certain types of moves, we don't care
1299 : ! if the values for those move parameters are strange
1300 :
1301 : ! find out if we're only doing HMC moves...we can ignore a lot of extra information
1302 : ! then, which would ordinarily be cause for concern
1303 26 : CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMHMC", r_val=mc_par%pmhmc)
1304 26 : CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMSWAP", r_val=mc_par%pmswap)
1305 26 : CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMVOLUME", r_val=mc_par%pmvolume)
1306 26 : CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMAVBMC", r_val=mc_par%pmavbmc)
1307 26 : CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMTRANS", r_val=mc_par%pmtrans)
1308 26 : CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMTRAION", r_val=mc_par%pmtraion)
1309 26 : CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMCLTRANS", r_val=mc_par%pmcltrans)
1310 :
1311 : ! first, grab all the integer values
1312 26 : CALL section_vals_val_get(mc_section, "NSTEP", i_val=mc_par%nstep)
1313 26 : CALL section_vals_val_get(mc_section, "NMOVES", i_val=mc_par%nmoves)
1314 26 : CALL section_vals_val_get(mc_section, "NSWAPMOVES", i_val=mc_par%nswapmoves)
1315 : CALL section_vals_val_get(mc_section, "MOVE_UPDATES%IUPVOLUME", &
1316 26 : i_val=mc_par%iupvolume)
1317 26 : CALL section_vals_val_get(mc_section, "NVIRIAL", i_val=mc_par%nvirial)
1318 : CALL section_vals_val_get(mc_section, "MOVE_UPDATES%IUPTRANS", &
1319 26 : i_val=mc_par%iuptrans)
1320 : CALL section_vals_val_get(mc_section, "MOVE_UPDATES%IUPCLTRANS", &
1321 26 : i_val=mc_par%iupcltrans)
1322 26 : CALL section_vals_val_get(mc_section, "IPRINT", i_val=mc_par%iprint)
1323 : ! now an integer array
1324 26 : CALL section_vals_val_get(mc_section, "AVBMC%AVBMC_ATOM", i_vals=temp_i_array)
1325 :
1326 26 : IF (mc_par%pmhmc - mc_par%pmswap >= 1.0_dp .AND. mc_par%pmhmc == 1.0_dp) THEN
1327 2 : mc_par%lhmc = .TRUE.
1328 : ELSE
1329 24 : mc_par%lhmc = .FALSE.
1330 : END IF
1331 :
1332 26 : IF (.NOT. mc_par%lhmc) THEN
1333 24 : IF (SIZE(temp_i_array) /= nmol_types) THEN
1334 0 : CPABORT('AVBMC_ATOM must have as many values as there are molecule types.')
1335 : END IF
1336 : END IF
1337 68 : DO imol = 1, SIZE(temp_i_array)
1338 68 : mc_par%avbmc_atom(imol) = temp_i_array(imol)
1339 : END DO
1340 :
1341 : ! now the real values
1342 26 : CALL section_vals_val_get(mc_section, "PRESSURE", r_val=mc_par%pressure)
1343 26 : CALL section_vals_val_get(mc_section, "TEMPERATURE", r_val=mc_par%temperature)
1344 26 : CALL section_vals_val_get(mc_section, "RCLUS", r_val=mc_par%rclus)
1345 : CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%BOX_PROBABILITIES%PMCLUS_BOX", &
1346 26 : r_val=mc_par%pmclus_box)
1347 : CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%BOX_PROBABILITIES%PMHMC_BOX", &
1348 26 : r_val=mc_par%pmhmc_box)
1349 : CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%BOX_PROBABILITIES%PMVOL_BOX", &
1350 26 : r_val=mc_par%pmvol_box)
1351 26 : CALL section_vals_val_get(mc_section, "DISCRETE_STEP", r_val=mc_par%discrete_step)
1352 : CALL section_vals_val_get(mc_section, "MAX_DISPLACEMENTS%BOX_DISPLACEMENTS%RMVOLUME", &
1353 26 : r_val=mc_par%rmvolume)
1354 : CALL section_vals_val_get(mc_section, "MAX_DISPLACEMENTS%BOX_DISPLACEMENTS%RMCLTRANS", &
1355 26 : r_val=mc_par%rmcltrans)
1356 :
1357 : ! finally the character values
1358 26 : CALL section_vals_val_get(mc_section, "ENSEMBLE", c_val=mc_par%ensemble)
1359 26 : CALL section_vals_val_get(mc_section, "RESTART_FILE_NAME", c_val=mc_par%restart_file_name)
1360 26 : CALL section_vals_val_get(mc_section, "COORDINATE_FILE_NAME", c_val=mc_par%coords_file)
1361 26 : CALL section_vals_val_get(mc_section, "ENERGY_FILE_NAME", c_val=mc_par%energy_file)
1362 26 : CALL section_vals_val_get(mc_section, "MOVES_FILE_NAME", c_val=mc_par%moves_file)
1363 26 : CALL section_vals_val_get(mc_section, "MOLECULES_FILE_NAME", c_val=mc_par%molecules_file)
1364 26 : CALL section_vals_val_get(mc_section, "CELL_FILE_NAME", c_val=mc_par%cell_file)
1365 26 : CALL section_vals_val_get(mc_section, "DATA_FILE_NAME", c_val=mc_par%data_file)
1366 26 : CALL section_vals_val_get(mc_section, "MAX_DISP_FILE_NAME", c_val=mc_par%displacement_file)
1367 26 : CALL section_vals_val_get(mc_section, "BOX2_FILE_NAME", c_val=mc_par%box2_file)
1368 :
1369 : ! set the values of the arrays...if we just point, we have problems when we start fooling around
1370 : ! with releasing force_envs and wonky values get set (despite that these are private)
1371 26 : IF (mc_par%ensemble == "VIRIAL") THEN
1372 2 : CALL section_vals_val_get(mc_section, "VIRIAL_TEMPS", r_vals=temp_r_array)
1373 : ! yes, I'm allocating here...I cannot find a better place to do it, though
1374 6 : ALLOCATE (mc_par%virial_temps(1:SIZE(temp_r_array)))
1375 6 : DO imol = 1, SIZE(temp_r_array)
1376 6 : mc_par%virial_temps(imol) = temp_r_array(imol)
1377 : END DO
1378 : END IF
1379 : ! all of these arrays should have one value for each type of molecule...so check that
1380 26 : CALL section_vals_val_get(mc_section, "AVBMC%AVBMC_RMIN", r_vals=temp_r_array)
1381 26 : IF (.NOT. mc_par%lhmc) THEN
1382 24 : IF (SIZE(temp_r_array) /= nmol_types) THEN
1383 0 : CPABORT('AVBMC_RMIN must have as many values as there are molecule types.')
1384 : END IF
1385 : END IF
1386 68 : DO imol = 1, SIZE(temp_r_array)
1387 68 : mc_par%avbmc_rmin(imol) = temp_r_array(imol)
1388 : END DO
1389 26 : CALL section_vals_val_get(mc_section, "AVBMC%AVBMC_RMAX", r_vals=temp_r_array)
1390 26 : IF (.NOT. mc_par%lhmc) THEN
1391 24 : IF (SIZE(temp_r_array) /= nmol_types) THEN
1392 0 : CPABORT('AVBMC_RMAXL must have as many values as there are molecule types.')
1393 : END IF
1394 : END IF
1395 68 : DO imol = 1, SIZE(temp_r_array)
1396 68 : mc_par%avbmc_rmax(imol) = temp_r_array(imol)
1397 : END DO
1398 26 : CALL section_vals_val_get(mc_section, "AVBMC%PBIAS", r_vals=temp_r_array)
1399 26 : IF (.NOT. mc_par%lhmc) THEN
1400 24 : IF (SIZE(temp_r_array) /= nmol_types) THEN
1401 0 : CPABORT('PBIAS must have as many values as there are molecule types.')
1402 : END IF
1403 : END IF
1404 68 : DO imol = 1, SIZE(temp_r_array)
1405 68 : mc_par%pbias(imol) = temp_r_array(imol)
1406 : END DO
1407 : CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%MOL_PROBABILITIES%PMAVBMC_MOL", &
1408 26 : r_vals=temp_r_array)
1409 26 : IF (.NOT. mc_par%lhmc) THEN
1410 24 : IF (SIZE(temp_r_array) /= nmol_types) THEN
1411 0 : CPABORT('PMAVBMC_MOL must have as many values as there are molecule types.')
1412 : END IF
1413 : END IF
1414 68 : DO imol = 1, SIZE(temp_r_array)
1415 68 : mc_par%pmavbmc_mol(imol) = temp_r_array(imol)
1416 : END DO
1417 26 : CALL section_vals_val_get(mc_section, "ETA", r_vals=temp_r_array)
1418 26 : IF (.NOT. mc_par%lhmc) THEN
1419 24 : IF (SIZE(temp_r_array) /= nmol_types) THEN
1420 0 : CPABORT('ETA must have as many values as there are molecule types.')
1421 : END IF
1422 : END IF
1423 68 : DO imol = 1, SIZE(temp_r_array)
1424 68 : mc_par%eta(imol) = temp_r_array(imol)
1425 : END DO
1426 : CALL section_vals_val_get(mc_section, "MAX_DISPLACEMENTS%MOL_DISPLACEMENTS%RMBOND", &
1427 26 : r_vals=temp_r_array)
1428 26 : IF (.NOT. mc_par%lhmc) THEN
1429 24 : IF (SIZE(temp_r_array) /= nmol_types) THEN
1430 0 : CPABORT('RMBOND must have as many values as there are molecule types.')
1431 : END IF
1432 : END IF
1433 68 : DO imol = 1, SIZE(temp_r_array)
1434 68 : mc_par%rmbond(imol) = temp_r_array(imol)
1435 : END DO
1436 : CALL section_vals_val_get(mc_section, "MAX_DISPLACEMENTS%MOL_DISPLACEMENTS%RMANGLE", &
1437 26 : r_vals=temp_r_array)
1438 26 : IF (.NOT. mc_par%lhmc) THEN
1439 24 : IF (SIZE(temp_r_array) /= nmol_types) THEN
1440 0 : CPABORT('RMANGLE must have as many values as there are molecule types.')
1441 : END IF
1442 : END IF
1443 68 : DO imol = 1, SIZE(temp_r_array)
1444 68 : mc_par%rmangle(imol) = temp_r_array(imol)
1445 : END DO
1446 : CALL section_vals_val_get(mc_section, "MAX_DISPLACEMENTS%MOL_DISPLACEMENTS%RMDIHEDRAL", &
1447 26 : r_vals=temp_r_array)
1448 26 : IF (.NOT. mc_par%lhmc) THEN
1449 24 : IF (SIZE(temp_r_array) /= nmol_types) THEN
1450 0 : CPABORT('RMDIHEDRAL must have as many values as there are molecule types.')
1451 : END IF
1452 : END IF
1453 68 : DO imol = 1, SIZE(temp_r_array)
1454 68 : mc_par%rmdihedral(imol) = temp_r_array(imol)
1455 : END DO
1456 : CALL section_vals_val_get(mc_section, "MAX_DISPLACEMENTS%MOL_DISPLACEMENTS%RMROT", &
1457 26 : r_vals=temp_r_array)
1458 26 : IF (.NOT. mc_par%lhmc) THEN
1459 24 : IF (SIZE(temp_r_array) /= nmol_types) THEN
1460 0 : CPABORT('RMROT must have as many values as there are molecule types.')
1461 : END IF
1462 : END IF
1463 68 : DO imol = 1, SIZE(temp_r_array)
1464 68 : mc_par%rmrot(imol) = temp_r_array(imol)
1465 : END DO
1466 : CALL section_vals_val_get(mc_section, "MAX_DISPLACEMENTS%MOL_DISPLACEMENTS%RMTRANS", &
1467 26 : r_vals=temp_r_array)
1468 26 : IF (.NOT. mc_par%lhmc) THEN
1469 24 : IF (SIZE(temp_r_array) /= nmol_types) THEN
1470 0 : CPABORT('RMTRANS must have as many values as there are molecule types.')
1471 : END IF
1472 : END IF
1473 68 : DO imol = 1, SIZE(temp_r_array)
1474 68 : mc_par%rmtrans(imol) = temp_r_array(imol)
1475 : END DO
1476 :
1477 : CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%MOL_PROBABILITIES%PMTRAION_MOL", &
1478 26 : r_vals=temp_r_array)
1479 26 : IF (.NOT. mc_par%lhmc) THEN
1480 24 : IF (SIZE(temp_r_array) /= nmol_types) THEN
1481 0 : CPABORT('PMTRAION_MOL must have as many values as there are molecule types.')
1482 : END IF
1483 : END IF
1484 68 : DO imol = 1, SIZE(temp_r_array)
1485 68 : mc_par%pmtraion_mol(imol) = temp_r_array(imol)
1486 : END DO
1487 :
1488 : CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%MOL_PROBABILITIES%PMTRANS_MOL", &
1489 26 : r_vals=temp_r_array)
1490 26 : IF (.NOT. mc_par%lhmc) THEN
1491 24 : IF (SIZE(temp_r_array) /= nmol_types) THEN
1492 0 : CPABORT('PMTRANS_MOL must have as many values as there are molecule types.')
1493 : END IF
1494 : END IF
1495 68 : DO imol = 1, SIZE(temp_r_array)
1496 68 : mc_par%pmtrans_mol(imol) = temp_r_array(imol)
1497 : END DO
1498 :
1499 : CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%MOL_PROBABILITIES%PMROT_MOL", &
1500 26 : r_vals=temp_r_array)
1501 26 : IF (.NOT. mc_par%lhmc) THEN
1502 24 : IF (SIZE(temp_r_array) /= nmol_types) THEN
1503 0 : CPABORT('PMROT_MOL must have as many values as there are molecule types.')
1504 : END IF
1505 : END IF
1506 68 : DO imol = 1, SIZE(temp_r_array)
1507 68 : mc_par%pmrot_mol(imol) = temp_r_array(imol)
1508 : END DO
1509 :
1510 : CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%MOL_PROBABILITIES%PMSWAP_MOL", &
1511 26 : r_vals=temp_r_array)
1512 26 : IF (.NOT. mc_par%lhmc) THEN
1513 24 : IF (SIZE(temp_r_array) /= nmol_types) THEN
1514 0 : CPABORT('PMSWAP_MOL must have as many values as there are molecule types.')
1515 : END IF
1516 : END IF
1517 68 : DO imol = 1, SIZE(temp_r_array)
1518 68 : mc_par%pmswap_mol(imol) = temp_r_array(imol)
1519 : END DO
1520 :
1521 : ! now some logical values
1522 26 : CALL section_vals_val_get(mc_section, "LBIAS", l_val=mc_par%lbias)
1523 26 : CALL section_vals_val_get(mc_section, "LDISCRETE", l_val=mc_par%ldiscrete)
1524 26 : CALL section_vals_val_get(mc_section, "LSTOP", l_val=mc_par%lstop)
1525 26 : CALL section_vals_val_get(mc_section, "RESTART", l_val=mc_par%lrestart)
1526 26 : CALL section_vals_val_get(mc_section, "LBIAS", l_val=mc_par%lbias)
1527 :
1528 : !..end of parsing the input section
1529 :
1530 : !..write some information to output
1531 26 : IF (iw > 0) THEN
1532 13 : WRITE (box_string, '(I2)') mc_par%mc_molecule_info%nboxes
1533 13 : WRITE (molecule_string, '(I2)') mc_par%mc_molecule_info%nmol_types
1534 13 : WRITE (tab_string, '(I4)') 81 - 10*mc_par%mc_molecule_info%nmol_types
1535 13 : WRITE (tab_box_string, '(I4)') 81 - 10*mc_par%mc_molecule_info%nboxes
1536 13 : WRITE (tab_string_int, '(I4)') 81 - 5*mc_par%mc_molecule_info%nmol_types
1537 13 : format_string = "(A,T"//TRIM(ADJUSTL(tab_string))//","//TRIM(ADJUSTL(molecule_string))//"(2X,F8.4))"
1538 13 : format_box_string = "(A,T"//TRIM(ADJUSTL(tab_box_string))//","//TRIM(ADJUSTL(box_string))//"(2X,F8.4))"
1539 13 : format_string_int = "(A,T"//TRIM(ADJUSTL(tab_string))//","//TRIM(ADJUSTL(molecule_string))//"(I3,2X))"
1540 13 : WRITE (iw, '( A )') ' MC| Monte Carlo Protocol '
1541 13 : WRITE (iw, '( A,T71,I10 )') ' MC| total number of steps ', &
1542 26 : mc_par%nstep
1543 13 : WRITE (iw, '( A,T71,F10.3 )') ' MC| pmvolume ', &
1544 26 : mc_par%pmvolume
1545 13 : WRITE (iw, '( A,T71,F10.3 )') ' MC| pmvol_box ', &
1546 26 : mc_par%pmvol_box
1547 13 : WRITE (iw, '( A,T71,F10.3 )') ' MC| pmclus_box ', &
1548 26 : mc_par%pmclus_box
1549 13 : WRITE (iw, '( A,T71,F10.3 )') ' MC| pmhmc ', &
1550 26 : mc_par%pmhmc
1551 13 : WRITE (iw, '( A,T71,F10.3 )') ' MC| pmhmc_box ', &
1552 26 : mc_par%pmhmc_box
1553 13 : WRITE (iw, '( A,T71,F10.3 )') ' MC| pmswap ', &
1554 26 : mc_par%pmswap
1555 13 : WRITE (iw, '( A,T71,F10.3 )') ' MC| pmavbmc ', &
1556 26 : mc_par%pmavbmc
1557 13 : WRITE (iw, '( A,T71,F10.3 )') ' MC| pmtraion ', &
1558 26 : mc_par%pmtraion
1559 13 : WRITE (iw, '( A,T71,F10.3 )') ' MC| pmtrans ', &
1560 26 : mc_par%pmtrans
1561 13 : WRITE (iw, '( A,T71,F10.3 )') ' MC| pmcltrans ', &
1562 26 : mc_par%pmcltrans
1563 13 : WRITE (iw, '( A,T71,I10 )') ' MC| iupvolume ', &
1564 26 : mc_par%iupvolume
1565 13 : WRITE (iw, '( A,T71,I10 )') ' MC| nvirial ', &
1566 26 : mc_par%nvirial
1567 13 : WRITE (iw, '( A,T71,I10 )') ' MC| iuptrans ', &
1568 26 : mc_par%iuptrans
1569 13 : WRITE (iw, '( A,T71,I10 )') ' MC| iupcltrans ', &
1570 26 : mc_par%iupcltrans
1571 13 : WRITE (iw, '( A,T71,I10 )') ' MC| iprint ', &
1572 26 : mc_par%iprint
1573 13 : WRITE (iw, '( A,T61,A20 )') ' MC| ensemble ', &
1574 26 : ADJUSTR(mc_par%ensemble)
1575 : ! format string may not be valid if all the molecules are atomic,
1576 : ! like in HMC
1577 13 : IF (.NOT. mc_par%lhmc) THEN
1578 12 : WRITE (iw, format_string) ' MC| pmswap_mol ', &
1579 44 : mc_par%pmswap_mol
1580 12 : WRITE (iw, format_string) ' MC| pmavbmc_mol ', &
1581 44 : mc_par%pmavbmc_mol
1582 12 : WRITE (iw, format_string) ' MC| pbias ', &
1583 44 : mc_par%pbias
1584 12 : WRITE (iw, format_string) ' MC| pmtraion_mol ', &
1585 44 : mc_par%pmtraion_mol
1586 12 : WRITE (iw, format_string) ' MC| pmtrans_mol ', &
1587 44 : mc_par%pmtrans_mol
1588 12 : WRITE (iw, format_string) ' MC| pmrot_mol ', &
1589 44 : mc_par%pmrot_mol
1590 12 : WRITE (iw, format_string) ' MC| eta [K]', &
1591 44 : mc_par%eta
1592 12 : WRITE (iw, format_string) ' MC| rmbond [angstroms]', &
1593 44 : mc_par%rmbond
1594 12 : WRITE (iw, format_string) ' MC| rmangle [degrees]', &
1595 44 : mc_par%rmangle
1596 12 : WRITE (iw, format_string) ' MC| rmdihedral [degrees]', &
1597 44 : mc_par%rmdihedral
1598 12 : WRITE (iw, format_string) ' MC| rmtrans [angstroms]', &
1599 44 : mc_par%rmtrans
1600 12 : WRITE (iw, format_string) ' MC| rmcltrans [angstroms]', &
1601 24 : mc_par%rmcltrans
1602 12 : WRITE (iw, format_string) ' MC| rmrot [degrees]', &
1603 44 : mc_par%rmrot
1604 12 : WRITE (iw, format_string_int) ' MC| AVBMC target atom ', &
1605 44 : mc_par%avbmc_atom
1606 12 : WRITE (iw, format_string) ' MC| AVBMC inner cutoff [ang]', &
1607 44 : mc_par%avbmc_rmin
1608 12 : WRITE (iw, format_string) ' MC| AVBMC outer cutoff [ang]', &
1609 44 : mc_par%avbmc_rmax
1610 : END IF
1611 13 : IF (mc_par%ensemble == 'GEMC-NVT') THEN
1612 0 : WRITE (iw, '( A,T58,A)') ' MC| Box 2 file', &
1613 0 : TRIM(mc_par%box2_file)
1614 : END IF
1615 13 : WRITE (iw, '( A,T58,A )') ' MC| Name of restart file:', &
1616 26 : TRIM(mc_par%restart_file_name)
1617 13 : WRITE (iw, '( A,A,T44,A )') ' MC| Name of output ', &
1618 13 : 'coordinate file:', &
1619 26 : TRIM(mc_par%coords_file)
1620 13 : WRITE (iw, '( A,T44,A )') ' MC| Name of output data file:', &
1621 26 : TRIM(mc_par%data_file)
1622 13 : WRITE (iw, '( A,A,T44,A )') ' MC| Name of output ', &
1623 13 : 'molecules file:', &
1624 26 : TRIM(mc_par%molecules_file)
1625 13 : WRITE (iw, '( A,T44,A )') ' MC| Name of output moves file:', &
1626 26 : TRIM(mc_par%moves_file)
1627 13 : WRITE (iw, '( A,T44,A )') ' MC| Name of output energy file:', &
1628 26 : TRIM(mc_par%energy_file)
1629 13 : WRITE (iw, '( A,T44,A )') ' MC| Name of output cell file:', &
1630 26 : TRIM(mc_par%cell_file)
1631 13 : WRITE (iw, '( A,A,T44,A )') ' MC| Name of output', &
1632 13 : ' displacement file:', &
1633 26 : TRIM(mc_par%displacement_file)
1634 13 : IF (mc_par%ldiscrete) THEN
1635 0 : WRITE (iw, '(A,A,T71,F10.3)') ' MC| discrete step size', &
1636 0 : '[angstroms]', &
1637 0 : mc_par%discrete_step
1638 : ELSE
1639 13 : WRITE (iw, '( A,A,T71,F10.3 )') ' MC| rmvolume ', &
1640 13 : '[cubic angstroms]', &
1641 26 : mc_par%rmvolume
1642 : END IF
1643 13 : WRITE (iw, '( A,T71,F10.2 )') ' MC| Temperature [K] ', &
1644 26 : mc_par%temperature
1645 13 : WRITE (iw, '( A,T71,F10.5 )') ' MC| Pressure [bar] ', &
1646 26 : mc_par%pressure
1647 13 : WRITE (iw, '( A,T71,F10.5 )') ' MC| Rclus [ang] ', &
1648 26 : mc_par%rclus
1649 13 : IF (mc_par%lrestart) THEN
1650 1 : WRITE (iw, '(A,A)') ' MC| Initial data will be read from a', &
1651 2 : ' restart file.'
1652 : END IF
1653 13 : IF (mc_par%lbias) THEN
1654 7 : WRITE (iw, '(A)') ' MC| The moves will be biased.'
1655 : ELSE
1656 6 : WRITE (iw, '(A)') ' MC| The moves will not be biased,'
1657 : END IF
1658 13 : IF (mc_par%nmoves == 1) THEN
1659 10 : WRITE (iw, '(A,A)') ' MC| A full energy calculation ', &
1660 20 : 'will be done at every step.'
1661 : ELSE
1662 3 : WRITE (iw, '(A,I4,A,A)') ' MC| ', mc_par%nmoves, &
1663 3 : ' moves will be attempted ', &
1664 6 : 'before a Quickstep energy calculation'
1665 3 : WRITE (iw, '(A)') ' MC| takes place.'
1666 : END IF
1667 :
1668 13 : WRITE (iw, '(A,I4,A,A)') ' MC| ', mc_par%nswapmoves, &
1669 13 : ' swap insertions will be attempted ', &
1670 26 : 'per molecular swap move'
1671 :
1672 13 : IF (mc_par%lhmc) THEN
1673 1 : WRITE (iw, '(A)') '********************************************************'
1674 1 : WRITE (iw, '(A)') '************** ONLY DOING HYBRID MONTE CARLO MOVES **************************'
1675 1 : WRITE (iw, '(A)') '********************************************************'
1676 :
1677 : END IF
1678 : END IF
1679 :
1680 : ! figure out what beta (1/kT) is in atomic units (1/Hartree)
1681 26 : mc_par%BETA = 1/mc_par%temperature/boltzmann*joule
1682 :
1683 : ! 0.9_dp is to avoid overflow or underflow
1684 26 : mc_par%exp_max_val = 0.9_dp*LOG(HUGE(0.0_dp))
1685 26 : mc_par%exp_min_val = 0.9_dp*LOG(TINY(0.0_dp))
1686 26 : mc_par%max_val = EXP(mc_par%exp_max_val)
1687 26 : mc_par%min_val = 0.0_dp
1688 :
1689 : ! convert from bar to a.u.
1690 26 : mc_par%pressure = cp_unit_to_cp2k(mc_par%pressure, "bar")
1691 26 : mc_par%rclus = cp_unit_to_cp2k(mc_par%rclus, "angstrom")
1692 : ! convert from angstrom to a.u.
1693 68 : DO itype = 1, mc_par%mc_molecule_info%nmol_types
1694 : ! convert from Kelvin to a.u.
1695 : ! mc_par % eta = mc_par%eta(itype) * boltzmann / joule
1696 : ! convert from degrees to radians
1697 42 : mc_par%rmrot(itype) = mc_par%rmrot(itype)/180.0e0_dp*pi
1698 42 : mc_par%rmangle(itype) = mc_par%rmangle(itype)/180.0e0_dp*pi
1699 42 : mc_par%rmdihedral(itype) = mc_par%rmdihedral(itype)/180.0e0_dp*pi
1700 42 : mc_par%rmtrans(itype) = cp_unit_to_cp2k(mc_par%rmtrans(itype), "angstrom")
1701 42 : mc_par%rmbond(itype) = cp_unit_to_cp2k(mc_par%rmbond(itype), "angstrom")
1702 42 : mc_par%eta(itype) = cp_unit_to_cp2k(mc_par%eta(itype), "K_e")
1703 42 : mc_par%avbmc_rmin(itype) = cp_unit_to_cp2k(mc_par%avbmc_rmin(itype), "angstrom")
1704 68 : mc_par%avbmc_rmax(itype) = cp_unit_to_cp2k(mc_par%avbmc_rmax(itype), "angstrom")
1705 : END DO
1706 26 : mc_par%rmvolume = cp_unit_to_cp2k(mc_par%rmvolume, "angstrom^3")
1707 26 : mc_par%rmcltrans = cp_unit_to_cp2k(mc_par%rmcltrans, "angstrom")
1708 26 : mc_par%discrete_step = cp_unit_to_cp2k(mc_par%discrete_step, "angstrom")
1709 :
1710 : ! end the timing
1711 26 : CALL timestop(handle)
1712 :
1713 26 : END SUBROUTINE read_mc_section
1714 :
1715 : ! **************************************************************************************************
1716 : !> \brief finds the largest interaction cutoff value in a classical simulation
1717 : !> so we know the smallest size we can make the box in a volume move
1718 : !> \param mc_par the structure that will store the parameters
1719 : !> \param force_env the force environment that we'll grab the rcut parameter
1720 : !> out of
1721 : !> \param lterminate set to .TRUE. if one of the sides of the box is
1722 : !> less than twice the cutoff
1723 : !>
1724 : !> Suitable for parallel.
1725 : !> \author MJM
1726 : ! **************************************************************************************************
1727 14 : SUBROUTINE find_mc_rcut(mc_par, force_env, lterminate)
1728 :
1729 : TYPE(mc_simpar_type), INTENT(INOUT) :: mc_par
1730 : TYPE(force_env_type), POINTER :: force_env
1731 : LOGICAL, INTENT(OUT) :: lterminate
1732 :
1733 : INTEGER :: itype, jtype
1734 : REAL(KIND=dp) :: rcutsq_max
1735 : REAL(KIND=dp), DIMENSION(1:3) :: abc
1736 : TYPE(cell_type), POINTER :: cell
1737 : TYPE(fist_environment_type), POINTER :: fist_env
1738 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
1739 : TYPE(pair_potential_pp_type), POINTER :: potparm
1740 :
1741 14 : NULLIFY (cell, potparm, fist_nonbond_env, fist_env)
1742 :
1743 14 : lterminate = .FALSE.
1744 14 : CALL force_env_get(force_env, fist_env=fist_env)
1745 : CALL fist_env_get(fist_env, cell=cell, &
1746 14 : fist_nonbond_env=fist_nonbond_env)
1747 14 : CALL fist_nonbond_env_get(fist_nonbond_env, potparm=potparm)
1748 14 : CALL get_cell(cell, abc=abc)
1749 :
1750 : ! find the largest value of rcutsq
1751 14 : rcutsq_max = 0.0e0_dp
1752 50 : DO itype = 1, SIZE(potparm%pot, 1)
1753 118 : DO jtype = itype, SIZE(potparm%pot, 2)
1754 68 : IF (potparm%pot(itype, jtype)%pot%rcutsq > rcutsq_max) &
1755 50 : rcutsq_max = potparm%pot(itype, jtype)%pot%rcutsq
1756 : END DO
1757 : END DO
1758 :
1759 : ! check to make sure all box dimensions are greater than two times this
1760 : ! value
1761 14 : mc_par%rcut = rcutsq_max**0.5_dp
1762 56 : DO itype = 1, 3
1763 56 : IF (abc(itype) < 2.0_dp*mc_par%rcut) THEN
1764 0 : lterminate = .TRUE.
1765 : END IF
1766 : END DO
1767 :
1768 14 : END SUBROUTINE find_mc_rcut
1769 :
1770 : ! **************************************************************************************************
1771 : !> \brief figures out the number of total molecules, the number of atoms in each
1772 : !> molecule, an array with the molecule types, etc...a lot of information
1773 : !> that we need. I did this because I use multiple force_env (simulation
1774 : !> boxes) for MC, and they don't know about each other.
1775 : !> \param force_env the pointer containing all the force environments in the
1776 : !> simulation
1777 : !> \param mc_molecule_info the structure that will hold all the information
1778 : !> for the molecule types
1779 : !>
1780 : !> Suitable for parallel.
1781 : !> \param box_number ...
1782 : !> \param coordinates_empty ...
1783 : !> \author MJM
1784 : ! **************************************************************************************************
1785 62 : SUBROUTINE mc_determine_molecule_info(force_env, mc_molecule_info, box_number, &
1786 : coordinates_empty)
1787 :
1788 : TYPE(force_env_p_type), DIMENSION(:), POINTER :: force_env
1789 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
1790 : INTEGER, INTENT(IN), OPTIONAL :: box_number
1791 : REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: coordinates_empty
1792 :
1793 : CHARACTER(LEN=default_string_length) :: name
1794 : CHARACTER(LEN=default_string_length), &
1795 62 : ALLOCATABLE, DIMENSION(:) :: names_init
1796 : INTEGER :: first_mol, iatom, ibox, imol, imolecule, itype, iunique, iunit, jtype, last_mol, &
1797 : natom, natoms_large, nbend, nbond, nboxes, nmol_types, nmolecule, ntorsion, ntypes, &
1798 : skip_box, this_molecule, total
1799 62 : INTEGER, DIMENSION(:), POINTER :: molecule_list
1800 : LOGICAL :: lnew_type
1801 62 : TYPE(atom_type), DIMENSION(:), POINTER :: atom_list
1802 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1803 : TYPE(cp_subsys_type), POINTER :: subsys
1804 : TYPE(molecule_kind_list_p_type), DIMENSION(:), &
1805 62 : POINTER :: molecule_kinds
1806 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1807 62 : TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
1808 :
1809 : ! how many simulation boxes do we have?
1810 :
1811 62 : nboxes = SIZE(force_env(:))
1812 :
1813 : ! allocate the pointer
1814 0 : ALLOCATE (mc_molecule_info)
1815 62 : mc_molecule_info%nboxes = nboxes
1816 :
1817 : ! if box_number is present, that box is supposed to be empty,
1818 : ! so we don't count it in anything
1819 62 : IF (PRESENT(box_number)) THEN
1820 20 : skip_box = box_number
1821 : ELSE
1822 : skip_box = 0
1823 : END IF
1824 :
1825 : ! we need to figure out how many different kinds of molecules we have...
1826 : ! the fun stuff comes in when box 1 is missing a molecule type...I'll
1827 : ! distinguish molecules based on names from the .psf files...
1828 : ! this is only really a problem when we have more than one simulation
1829 : ! box
1830 62 : NULLIFY (subsys, molecule_kinds, molecule_kind)
1831 :
1832 282 : ALLOCATE (particles(1:nboxes))
1833 220 : ALLOCATE (molecule_kinds(1:nboxes))
1834 :
1835 : ! find out how many types we have
1836 : ntypes = 0
1837 158 : DO ibox = 1, nboxes
1838 96 : IF (ibox == skip_box) CYCLE
1839 : CALL force_env_get(force_env(ibox)%force_env, &
1840 96 : subsys=subsys)
1841 : CALL cp_subsys_get(subsys, &
1842 : molecule_kinds=molecule_kinds(ibox)%list, &
1843 96 : particles=particles(ibox)%list)
1844 158 : ntypes = ntypes + SIZE(molecule_kinds(ibox)%list%els(:))
1845 : END DO
1846 :
1847 186 : ALLOCATE (names_init(1:ntypes))
1848 :
1849 : ! now get the names of all types so we can figure out how many distinct
1850 : ! types we have
1851 158 : itype = 1
1852 158 : natoms_large = 0
1853 158 : DO ibox = 1, nboxes
1854 96 : IF (ibox == skip_box) CYCLE
1855 330 : DO imolecule = 1, SIZE(molecule_kinds(ibox)%list%els(:))
1856 172 : molecule_kind => molecule_kinds(ibox)%list%els(imolecule)
1857 : CALL get_molecule_kind(molecule_kind, name=names_init(itype), &
1858 172 : natom=natom)
1859 172 : IF (natom > natoms_large) natoms_large = natom
1860 268 : itype = itype + 1
1861 : END DO
1862 : END DO
1863 :
1864 : nmol_types = 0
1865 234 : DO itype = 1, ntypes
1866 : lnew_type = .TRUE.
1867 384 : DO jtype = 1, itype - 1
1868 212 : IF (TRIM(names_init(itype)) == TRIM(names_init(jtype))) &
1869 240 : lnew_type = .FALSE.
1870 : END DO
1871 234 : IF (lnew_type) THEN
1872 104 : nmol_types = nmol_types + 1
1873 : ELSE
1874 68 : names_init(itype) = ''
1875 : END IF
1876 : END DO
1877 :
1878 : ! allocate arrays for the names of the molecules, the number of atoms, and
1879 : ! the conformational probabilities
1880 62 : mc_molecule_info%nmol_types = nmol_types
1881 186 : ALLOCATE (mc_molecule_info%names(1:nmol_types))
1882 186 : ALLOCATE (mc_molecule_info%nunits(1:nmol_types))
1883 186 : ALLOCATE (mc_molecule_info%conf_prob(1:3, 1:nmol_types))
1884 248 : ALLOCATE (mc_molecule_info%nchains(1:nmol_types, 1:nboxes))
1885 186 : ALLOCATE (mc_molecule_info%nunits_tot(1:nboxes))
1886 248 : ALLOCATE (mc_molecule_info%atom_names(1:natoms_large, 1:nmol_types))
1887 248 : ALLOCATE (mc_molecule_info%mass(1:natoms_large, 1:nmol_types))
1888 :
1889 : ! now record all the information for every molecule type
1890 62 : iunique = 0
1891 62 : itype = 0
1892 158 : DO ibox = 1, nboxes
1893 96 : IF (ibox == skip_box) CYCLE
1894 330 : DO imolecule = 1, SIZE(molecule_kinds(ibox)%list%els(:))
1895 172 : itype = itype + 1
1896 172 : IF (names_init(itype) == '') CYCLE
1897 104 : iunique = iunique + 1
1898 104 : mc_molecule_info%names(iunique) = names_init(itype)
1899 104 : molecule_kind => molecule_kinds(ibox)%list%els(imolecule)
1900 : CALL get_molecule_kind(molecule_kind, natom=mc_molecule_info%nunits(iunique), &
1901 104 : nbond=nbond, nbend=nbend, ntorsion=ntorsion, atom_list=atom_list)
1902 :
1903 320 : DO iatom = 1, mc_molecule_info%nunits(iunique)
1904 216 : atomic_kind => atom_list(iatom)%atomic_kind
1905 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1906 : name=mc_molecule_info%atom_names(iatom, iunique), &
1907 320 : mass=mc_molecule_info%mass(iatom, iunique))
1908 : END DO
1909 :
1910 : ! compute the probabilities of doing any particular kind of conformation change
1911 104 : total = nbond + nbend + ntorsion
1912 :
1913 304 : IF (total == 0) THEN
1914 184 : mc_molecule_info%conf_prob(:, iunique) = 0.0e0_dp
1915 : ELSE
1916 58 : mc_molecule_info%conf_prob(1, iunique) = REAL(nbond, dp)/REAL(total, dp)
1917 58 : mc_molecule_info%conf_prob(2, iunique) = REAL(nbend, dp)/REAL(total, dp)
1918 58 : mc_molecule_info%conf_prob(3, iunique) = REAL(ntorsion, dp)/REAL(total, dp)
1919 : END IF
1920 :
1921 : END DO
1922 : END DO
1923 :
1924 : ! figure out the empty coordinates
1925 62 : IF (PRESENT(coordinates_empty)) THEN
1926 60 : ALLOCATE (coordinates_empty(1:3, 1:mc_molecule_info%nunits(1)))
1927 74 : DO iunit = 1, mc_molecule_info%nunits(1)
1928 236 : coordinates_empty(1:3, iunit) = particles(1)%list%els(iunit)%r(1:3)
1929 : END DO
1930 : END IF
1931 :
1932 : ! now we need to figure out the total number of molecules of each kind in
1933 : ! each box, and the total number of interaction sites in each box
1934 330 : mc_molecule_info%nchains(:, :) = 0
1935 166 : DO iunique = 1, nmol_types
1936 338 : DO ibox = 1, nboxes
1937 172 : IF (ibox == skip_box) CYCLE
1938 600 : DO imolecule = 1, SIZE(molecule_kinds(ibox)%list%els(:))
1939 324 : molecule_kind => molecule_kinds(ibox)%list%els(imolecule)
1940 : CALL get_molecule_kind(molecule_kind, nmolecule=nmolecule, &
1941 324 : name=name)
1942 324 : IF (TRIM(name) /= mc_molecule_info%names(iunique)) CYCLE
1943 668 : mc_molecule_info%nchains(iunique, ibox) = mc_molecule_info%nchains(iunique, ibox) + nmolecule
1944 : END DO
1945 : END DO
1946 : END DO
1947 62 : mc_molecule_info%nchain_total = 0
1948 158 : DO ibox = 1, nboxes
1949 96 : mc_molecule_info%nunits_tot(ibox) = 0
1950 96 : IF (ibox == skip_box) CYCLE
1951 268 : DO iunique = 1, nmol_types
1952 : mc_molecule_info%nunits_tot(ibox) = mc_molecule_info%nunits_tot(ibox) + &
1953 268 : mc_molecule_info%nunits(iunique)*mc_molecule_info%nchains(iunique, ibox)
1954 : END DO
1955 330 : mc_molecule_info%nchain_total = mc_molecule_info%nchain_total + SUM(mc_molecule_info%nchains(:, ibox))
1956 : END DO
1957 :
1958 : ! now we need to figure out which type every molecule is,
1959 : ! and which force_env (box) it's in
1960 186 : ALLOCATE (mc_molecule_info%mol_type(1:mc_molecule_info%nchain_total))
1961 186 : ALLOCATE (mc_molecule_info%in_box(1:mc_molecule_info%nchain_total))
1962 :
1963 62 : last_mol = 0
1964 158 : DO ibox = 1, nboxes
1965 96 : IF (ibox == skip_box) CYCLE
1966 96 : first_mol = last_mol + 1
1967 268 : last_mol = first_mol + SUM(mc_molecule_info%nchains(:, ibox)) - 1
1968 5404 : mc_molecule_info%in_box(first_mol:last_mol) = ibox
1969 330 : DO imolecule = 1, SIZE(molecule_kinds(ibox)%list%els(:))
1970 172 : molecule_kind => molecule_kinds(ibox)%list%els(imolecule)
1971 : CALL get_molecule_kind(molecule_kind, nmolecule=nmolecule, &
1972 172 : name=name, molecule_list=molecule_list)
1973 : ! figure out which molecule number this is
1974 172 : this_molecule = 0
1975 496 : DO iunique = 1, nmol_types
1976 496 : IF (TRIM(name) == TRIM(mc_molecule_info%names(iunique))) THEN
1977 172 : this_molecule = iunique
1978 : END IF
1979 : END DO
1980 5748 : DO imol = 1, SIZE(molecule_list(:))
1981 5480 : mc_molecule_info%mol_type(first_mol + molecule_list(imol) - 1) = this_molecule
1982 : END DO
1983 : END DO
1984 : END DO
1985 :
1986 : ! get rid of stuff we don't need
1987 62 : DEALLOCATE (names_init)
1988 62 : DEALLOCATE (molecule_kinds)
1989 62 : DEALLOCATE (particles)
1990 :
1991 62 : END SUBROUTINE MC_DETERMINE_MOLECULE_INFO
1992 :
1993 : ! **************************************************************************************************
1994 : !> \brief deallocates all the arrays in the mc_molecule_info_type
1995 : !> \param mc_molecule_info the structure we wish to deallocate
1996 : !>
1997 : !> Suitable for parallel.
1998 : !> \author MJM
1999 : ! **************************************************************************************************
2000 62 : SUBROUTINE mc_molecule_info_destroy(mc_molecule_info)
2001 :
2002 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
2003 :
2004 62 : DEALLOCATE (mc_molecule_info%nchains)
2005 62 : DEALLOCATE (mc_molecule_info%nunits)
2006 62 : DEALLOCATE (mc_molecule_info%mol_type)
2007 62 : DEALLOCATE (mc_molecule_info%in_box)
2008 62 : DEALLOCATE (mc_molecule_info%names)
2009 62 : DEALLOCATE (mc_molecule_info%atom_names)
2010 62 : DEALLOCATE (mc_molecule_info%conf_prob)
2011 62 : DEALLOCATE (mc_molecule_info%nunits_tot)
2012 62 : DEALLOCATE (mc_molecule_info%mass)
2013 :
2014 62 : DEALLOCATE (mc_molecule_info)
2015 : NULLIFY (mc_molecule_info)
2016 :
2017 62 : END SUBROUTINE mc_molecule_info_destroy
2018 :
2019 : ! **************************************************************************************************
2020 : !> \brief ...
2021 : !> \param mc_par ...
2022 : ! **************************************************************************************************
2023 52 : SUBROUTINE mc_input_parameters_check(mc_par)
2024 :
2025 : ! **************************************************************************************************
2026 : !> \brief accesses the private elements of the mc_molecule_info_type
2027 : !> \param mc_molecule_info the structure you want the parameters for
2028 : !> \param nmol_types the number of molecule types in the simulation
2029 : !>
2030 : !> Suitable for parallel.
2031 : !> \author MJM
2032 : TYPE(mc_simpar_type), POINTER :: mc_par
2033 :
2034 : INTEGER :: imol_type, nboxes, nchain_total, &
2035 : nmol_types
2036 26 : INTEGER, DIMENSION(:), POINTER :: nunits
2037 : LOGICAL :: lhmc
2038 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
2039 :
2040 26 : CALL get_mc_par(mc_par, mc_molecule_info=mc_molecule_info, lhmc=lhmc)
2041 :
2042 : CALL get_mc_molecule_info(mc_molecule_info, nmol_types=nmol_types, &
2043 26 : nboxes=nboxes, nunits=nunits, nchain_total=nchain_total)
2044 :
2045 : ! if we're only doing HMC, we can skip these checks
2046 26 : IF (lhmc) RETURN
2047 :
2048 : ! the final value of these arrays needs to be 1.0, otherwise there is
2049 : ! a chance we won't select a molecule type for a move
2050 24 : IF (mc_par%pmavbmc_mol(nmol_types) < 0.9999E0_dp) THEN
2051 0 : CPABORT('The last value of PMAVBMC_MOL needs to be 1.0')
2052 : END IF
2053 24 : IF (mc_par%pmswap_mol(nmol_types) < 0.9999E0_dp) THEN
2054 0 : CPABORT('The last value of PMSWAP_MOL needs to be 1.0')
2055 : END IF
2056 24 : IF (mc_par%pmtraion_mol(nmol_types) < 0.9999E0_dp) THEN
2057 0 : CPABORT('The last value of PMTRAION_MOL needs to be 1.0')
2058 : END IF
2059 24 : IF (mc_par%pmtrans_mol(nmol_types) < 0.9999E0_dp) THEN
2060 0 : CPABORT('The last value of PMTRANS_MOL needs to be 1.0')
2061 : END IF
2062 24 : IF (mc_par%pmrot_mol(nmol_types) < 0.9999E0_dp) THEN
2063 0 : CPABORT('The last value of PMROT_MOL needs to be 1.0')
2064 : END IF
2065 :
2066 : ! check to make sure we have all the desired values for various
2067 :
2068 : ! some ensembles require multiple boxes or molecule types
2069 28 : SELECT CASE (mc_par%ensemble)
2070 : CASE ("GEMC_NPT")
2071 4 : IF (nmol_types <= 1) &
2072 0 : CPABORT('Cannot have GEMC-NPT simulation with only one molecule type')
2073 4 : IF (nboxes <= 1) &
2074 0 : CPABORT('Cannot have GEMC-NPT simulation with only one box')
2075 : CASE ("GEMC_NVT")
2076 8 : IF (nboxes <= 1) &
2077 0 : CPABORT('Cannot have GEMC-NVT simulation with only one box')
2078 : CASE ("TRADITIONAL")
2079 10 : IF (mc_par%pmswap > 0.0E0_dp) &
2080 0 : CPABORT('You cannot do swap moves in a system with only one box')
2081 : CASE ("VIRIAL")
2082 2 : IF (nchain_total /= 2) &
2083 24 : CPABORT('You need exactly two molecules in the box to compute the second virial.')
2084 : END SELECT
2085 :
2086 : ! can't choose an AVBMC target atom number higher than the number
2087 : ! of units in that molecule
2088 64 : DO imol_type = 1, nmol_types
2089 64 : IF (mc_par%avbmc_atom(imol_type) > nunits(imol_type)) THEN
2090 0 : CPABORT('Cannot have avbmc_atom higher than the number of atoms for that molecule type!')
2091 : END IF
2092 : END DO
2093 :
2094 26 : END SUBROUTINE mc_input_parameters_check
2095 :
2096 0 : END MODULE mc_types
2097 :
|
