(git:ccc2433)
mc_types.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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,&
18  USE cell_types, ONLY: cell_type,&
19  get_cell
20  USE cp_files, ONLY: close_file,&
21  open_file
23  USE cp_subsys_types, ONLY: cp_subsys_get,&
24  cp_subsys_type
25  USE cp_units, ONLY: cp_unit_to_cp2k
27  fist_environment_type
29  fist_nonbond_env_type
30  USE force_env_types, ONLY: force_env_get,&
31  force_env_p_type,&
32  force_env_type,&
35  USE global_types, ONLY: global_environment_type
36  USE input_constants, ONLY: do_fist,&
37  do_qs
39  section_vals_type,&
42  USE kinds, ONLY: default_path_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,&
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, &
71  find_mc_rcut, mc_molecule_info_type, &
75  mc_input_file_type, mc_input_file_create, &
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  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  IF (PRESENT(in_use)) in_use = mc_input_file%in_use
293  IF (PRESENT(run_type_row)) run_type_row = mc_input_file%run_type_row
294  IF (PRESENT(run_type_column)) run_type_column = mc_input_file%run_type_column
295  IF (PRESENT(coord_row_start)) coord_row_start = mc_input_file%coord_row_start
296  IF (PRESENT(coord_row_end)) coord_row_end = mc_input_file%coord_row_end
297  IF (PRESENT(cell_row)) cell_row = mc_input_file%cell_row
298  IF (PRESENT(cell_column)) cell_column = mc_input_file%cell_column
299  IF (PRESENT(force_eval_row_start)) force_eval_row_start = mc_input_file%force_eval_row_start
300  IF (PRESENT(force_eval_row_end)) force_eval_row_end = mc_input_file%force_eval_row_end
301  IF (PRESENT(global_row_start)) global_row_start = mc_input_file%global_row_start
302  IF (PRESENT(global_row_end)) global_row_end = mc_input_file%global_row_end
303  IF (PRESENT(nunits_empty)) nunits_empty = mc_input_file%nunits_empty
304  IF (PRESENT(motion_row_start)) motion_row_start = mc_input_file%motion_row_start
305  IF (PRESENT(motion_row_end)) motion_row_end = mc_input_file%motion_row_end
306 
307  IF (PRESENT(mol_set_nmol_row)) mol_set_nmol_row => mc_input_file%mol_set_nmol_row
308  IF (PRESENT(mol_set_nmol_column)) mol_set_nmol_column => mc_input_file%mol_set_nmol_column
309  IF (PRESENT(text)) text => mc_input_file%text
310  IF (PRESENT(atom_names_empty)) atom_names_empty => mc_input_file%atom_names_empty
311  IF (PRESENT(coordinates_empty)) coordinates_empty => mc_input_file%coordinates_empty
312 
313  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  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  IF (PRESENT(nstep)) nstep = mc_par%nstep
449  IF (PRESENT(nvirial)) nvirial = mc_par%nvirial
450  IF (PRESENT(iuptrans)) iuptrans = mc_par%iuptrans
451  IF (PRESENT(iupcltrans)) iupcltrans = mc_par%iupcltrans
452  IF (PRESENT(iupvolume)) iupvolume = mc_par%iupvolume
453  IF (PRESENT(nmoves)) nmoves = mc_par%nmoves
454  IF (PRESENT(nswapmoves)) nswapmoves = mc_par%nswapmoves
455  IF (PRESENT(rm)) rm = mc_par%rm
456  IF (PRESENT(cl)) cl = mc_par%cl
457  IF (PRESENT(diff)) diff = mc_par%diff
458  IF (PRESENT(nstart)) nstart = mc_par%nstart
459  IF (PRESENT(source)) source = mc_par%source
460  IF (PRESENT(group)) group = mc_par%group
461  IF (PRESENT(iprint)) iprint = mc_par%iprint
462 
463  IF (PRESENT(lbias)) lbias = mc_par%lbias
464  IF (PRESENT(ionode)) ionode = mc_par%ionode
465  IF (PRESENT(lrestart)) lrestart = mc_par%lrestart
466  IF (PRESENT(lstop)) lstop = mc_par%lstop
467  IF (PRESENT(ldiscrete)) ldiscrete = mc_par%ldiscrete
468 
469  IF (PRESENT(virial_temps)) virial_temps => mc_par%virial_temps
470  IF (PRESENT(avbmc_atom)) avbmc_atom => mc_par%avbmc_atom
471  IF (PRESENT(avbmc_rmin)) avbmc_rmin => mc_par%avbmc_rmin
472  IF (PRESENT(avbmc_rmax)) avbmc_rmax => mc_par%avbmc_rmax
473  IF (PRESENT(rcut)) rcut = mc_par%rcut
474  IF (PRESENT(discrete_step)) discrete_step = mc_par%discrete_step
475  IF (PRESENT(rmvolume)) rmvolume = mc_par%rmvolume
476  IF (PRESENT(rmcltrans)) rmcltrans = mc_par%rmcltrans
477  IF (PRESENT(temperature)) temperature = mc_par%temperature
478  IF (PRESENT(pressure)) pressure = mc_par%pressure
479  IF (PRESENT(rclus)) rclus = mc_par%rclus
480  IF (PRESENT(beta)) beta = mc_par%BETA
481  IF (PRESENT(exp_max_val)) exp_max_val = mc_par%exp_max_val
482  IF (PRESENT(exp_min_val)) exp_min_val = mc_par%exp_min_val
483  IF (PRESENT(max_val)) max_val = mc_par%max_val
484  IF (PRESENT(min_val)) min_val = mc_par%min_val
485  IF (PRESENT(pmswap)) pmswap = mc_par%pmswap
486  IF (PRESENT(pmvolume)) pmvolume = mc_par%pmvolume
487  IF (PRESENT(lhmc)) lhmc = mc_par%lhmc
488  IF (PRESENT(pmhmc)) pmhmc = mc_par%pmhmc
489  IF (PRESENT(pmhmc_box)) pmhmc_box = mc_par%pmhmc_box
490  IF (PRESENT(pmvol_box)) pmvol_box = mc_par%pmvol_box
491  IF (PRESENT(pmclus_box)) pmclus_box = mc_par%pmclus_box
492  IF (PRESENT(pmtraion)) pmtraion = mc_par%pmtraion
493  IF (PRESENT(pmtrans)) pmtrans = mc_par%pmtrans
494  IF (PRESENT(pmcltrans)) pmcltrans = mc_par%pmcltrans
495  IF (PRESENT(pmavbmc)) pmavbmc = mc_par%pmavbmc
496  IF (PRESENT(pmrot_mol)) pmrot_mol => mc_par%pmrot_mol
497  IF (PRESENT(pmtrans_mol)) pmtrans_mol => mc_par%pmtrans_mol
498  IF (PRESENT(pmtraion_mol)) pmtraion_mol => mc_par%pmtraion_mol
499  IF (PRESENT(pmavbmc_mol)) pmavbmc_mol => mc_par%pmavbmc_mol
500  IF (PRESENT(pbias)) pbias => mc_par%pbias
501 
502  IF (PRESENT(rmbond)) rmbond => mc_par%rmbond
503  IF (PRESENT(rmangle)) rmangle => mc_par%rmangle
504  IF (PRESENT(rmdihedral)) rmdihedral => mc_par%rmdihedral
505  IF (PRESENT(rmrot)) rmrot => mc_par%rmrot
506  IF (PRESENT(rmtrans)) rmtrans => mc_par%rmtrans
507 
508  IF (PRESENT(eta)) eta => mc_par%eta
509 
510  IF (PRESENT(ensemble)) ensemble = mc_par%ensemble
511  IF (PRESENT(program)) PROGRAM = mc_par%program
512  IF (PRESENT(restart_file_name)) restart_file_name = &
513  mc_par%restart_file_name
514  IF (PRESENT(moves_file)) moves_file = mc_par%moves_file
515  IF (PRESENT(coords_file)) coords_file = mc_par%coords_file
516  IF (PRESENT(molecules_file)) molecules_file = mc_par%molecules_file
517  IF (PRESENT(energy_file)) energy_file = mc_par%energy_file
518  IF (PRESENT(displacement_file)) displacement_file = &
519  mc_par%displacement_file
520  IF (PRESENT(cell_file)) cell_file = mc_par%cell_file
521  IF (PRESENT(dat_file)) dat_file = mc_par%dat_file
522  IF (PRESENT(data_file)) data_file = mc_par%data_file
523  IF (PRESENT(box2_file)) box2_file = mc_par%box2_file
524  IF (PRESENT(fft_lib)) fft_lib = mc_par%fft_lib
525 
526  IF (PRESENT(input_file)) input_file => mc_par%input_file
527  IF (PRESENT(mc_molecule_info)) mc_molecule_info => mc_par%mc_molecule_info
528  IF (PRESENT(mc_input_file)) mc_input_file => mc_par%mc_input_file
529  IF (PRESENT(mc_bias_file)) mc_bias_file => mc_par%mc_bias_file
530 
531  IF (PRESENT(pmswap_mol)) pmswap_mol => mc_par%pmswap_mol
532  IF (PRESENT(rand2skip)) rand2skip = mc_par%rand2skip
533  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  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  IF (PRESENT(nboxes)) nboxes = mc_molecule_info%nboxes
573  IF (PRESENT(nchain_total)) nchain_total = mc_molecule_info%nchain_total
574  IF (PRESENT(nmol_types)) nmol_types = mc_molecule_info%nmol_types
575 
576  IF (PRESENT(names)) names => mc_molecule_info%names
577  IF (PRESENT(atom_names)) atom_names => mc_molecule_info%atom_names
578 
579  IF (PRESENT(conf_prob)) conf_prob => mc_molecule_info%conf_prob
580  IF (PRESENT(mass)) mass => mc_molecule_info%mass
581 
582  IF (PRESENT(nchains)) nchains => mc_molecule_info%nchains
583  IF (PRESENT(nunits)) nunits => mc_molecule_info%nunits
584  IF (PRESENT(mol_type)) mol_type => mc_molecule_info%mol_type
585  IF (PRESENT(nunits_tot)) nunits_tot => mc_molecule_info%nunits_tot
586  IF (PRESENT(in_box)) in_box => mc_molecule_info%in_box
587 
588  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  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  IF (PRESENT(nstep)) mc_par%nstep = nstep
709  IF (PRESENT(rm)) mc_par%rm = rm
710  IF (PRESENT(cl)) mc_par%cl = cl
711  IF (PRESENT(diff)) mc_par%diff = diff
712  IF (PRESENT(nstart)) mc_par%nstart = nstart
713  IF (PRESENT(nmoves)) mc_par%nmoves = nmoves
714  IF (PRESENT(nswapmoves)) mc_par%nswapmoves = nswapmoves
715  IF (PRESENT(iprint)) mc_par%iprint = iprint
716  IF (PRESENT(iuptrans)) mc_par%iuptrans = iuptrans
717  IF (PRESENT(iupcltrans)) mc_par%iupcltrans = iupcltrans
718  IF (PRESENT(iupvolume)) mc_par%iupvolume = iupvolume
719 
720  IF (PRESENT(ldiscrete)) mc_par%ldiscrete = ldiscrete
721  IF (PRESENT(lstop)) mc_par%lstop = lstop
722  IF (PRESENT(lbias)) mc_par%lbias = lbias
723  IF (PRESENT(lrestart)) mc_par%lrestart = lrestart
724 
725  IF (PRESENT(exp_max_val)) mc_par%exp_max_val = exp_max_val
726  IF (PRESENT(exp_min_val)) mc_par%exp_min_val = exp_min_val
727  IF (PRESENT(max_val)) mc_par%max_val = max_val
728  IF (PRESENT(min_val)) mc_par%min_val = min_val
729  IF (PRESENT(beta)) mc_par%BETA = beta
730  IF (PRESENT(temperature)) mc_par%temperature = temperature
731  IF (PRESENT(rcut)) mc_par%rcut = rcut
732  IF (PRESENT(pressure)) mc_par%pressure = pressure
733  IF (PRESENT(rclus)) mc_par%rclus = rclus
734  IF (PRESENT(pmvolume)) mc_par%pmvolume = pmvolume
735  IF (PRESENT(lhmc)) mc_par%lhmc = lhmc
736  IF (PRESENT(pmhmc)) mc_par%pmhmc = pmhmc
737  IF (PRESENT(pmhmc_box)) mc_par%pmhmc_box = pmhmc_box
738  IF (PRESENT(pmvol_box)) mc_par%pmvol_box = pmvol_box
739  IF (PRESENT(pmclus_box)) mc_par%pmclus_box = pmclus_box
740  IF (PRESENT(pmswap)) mc_par%pmswap = pmswap
741  IF (PRESENT(pmtrans)) mc_par%pmtrans = pmtrans
742  IF (PRESENT(pmcltrans)) mc_par%pmcltrans = pmcltrans
743  IF (PRESENT(pmtraion)) mc_par%pmtraion = pmtraion
744  IF (PRESENT(pmavbmc)) mc_par%pmavbmc = pmavbmc
745 
746  IF (PRESENT(discrete_step)) mc_par%discrete_step = discrete_step
747  IF (PRESENT(rmvolume)) mc_par%rmvolume = rmvolume
748  IF (PRESENT(rmcltrans)) mc_par%rmcltrans = rmcltrans
749 
750  IF (PRESENT(mc_input_file)) mc_par%mc_input_file => mc_input_file
751  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  IF (PRESENT(eta)) THEN
756  DO iparm = 1, SIZE(eta)
757  mc_par%eta(iparm) = eta(iparm)
758  END DO
759  END IF
760  IF (PRESENT(rmbond)) THEN
761  DO iparm = 1, SIZE(rmbond)
762  mc_par%rmbond(iparm) = rmbond(iparm)
763  END DO
764  END IF
765  IF (PRESENT(rmangle)) THEN
766  DO iparm = 1, SIZE(rmangle)
767  mc_par%rmangle(iparm) = rmangle(iparm)
768  END DO
769  END IF
770  IF (PRESENT(rmdihedral)) THEN
771  DO iparm = 1, SIZE(rmdihedral)
772  mc_par%rmdihedral(iparm) = rmdihedral(iparm)
773  END DO
774  END IF
775  IF (PRESENT(rmrot)) THEN
776  DO iparm = 1, SIZE(rmrot)
777  mc_par%rmrot(iparm) = rmrot(iparm)
778  END DO
779  END IF
780  IF (PRESENT(rmtrans)) THEN
781  DO iparm = 1, SIZE(rmtrans)
782  mc_par%rmtrans(iparm) = rmtrans(iparm)
783  END DO
784  END IF
785 
786  IF (PRESENT(pmavbmc_mol)) THEN
787  DO iparm = 1, SIZE(pmavbmc_mol)
788  mc_par%pmavbmc_mol(iparm) = pmavbmc_mol(iparm)
789  END DO
790  END IF
791  IF (PRESENT(pmtrans_mol)) THEN
792  DO iparm = 1, SIZE(pmtrans_mol)
793  mc_par%pmtrans_mol(iparm) = pmtrans_mol(iparm)
794  END DO
795  END IF
796  IF (PRESENT(pmrot_mol)) THEN
797  DO iparm = 1, SIZE(pmrot_mol)
798  mc_par%pmrot_mol(iparm) = pmrot_mol(iparm)
799  END DO
800  END IF
801  IF (PRESENT(pmtraion_mol)) THEN
802  DO iparm = 1, SIZE(pmtraion_mol)
803  mc_par%pmtraion_mol(iparm) = pmtraion_mol(iparm)
804  END DO
805  END IF
806  IF (PRESENT(pmswap_mol)) THEN
807  DO iparm = 1, SIZE(pmswap_mol)
808  mc_par%pmswap_mol(iparm) = pmswap_mol(iparm)
809  END DO
810  END IF
811 
812  IF (PRESENT(avbmc_atom)) THEN
813  DO iparm = 1, SIZE(avbmc_atom)
814  mc_par%avbmc_atom(iparm) = avbmc_atom(iparm)
815  END DO
816  END IF
817  IF (PRESENT(avbmc_rmin)) THEN
818  DO iparm = 1, SIZE(avbmc_rmin)
819  mc_par%avbmc_rmin(iparm) = avbmc_rmin(iparm)
820  END DO
821  END IF
822  IF (PRESENT(avbmc_rmax)) THEN
823  DO iparm = 1, SIZE(avbmc_rmax)
824  mc_par%avbmc_rmax(iparm) = avbmc_rmax(iparm)
825  END DO
826  END IF
827  IF (PRESENT(pbias)) THEN
828  DO iparm = 1, SIZE(pbias)
829  mc_par%pbias(iparm) = pbias(iparm)
830  END DO
831  END IF
832 
833  IF (PRESENT(program)) mc_par%program = PROGRAM
834  IF (PRESENT(ensemble)) mc_par%ensemble = ensemble
835 
836  IF (PRESENT(mc_molecule_info)) mc_par%mc_molecule_info => mc_molecule_info
837  IF (PRESENT(source)) mc_par%source = source
838  IF (PRESENT(group)) mc_par%group = group
839  IF (PRESENT(ionode)) mc_par%ionode = ionode
840 
841  IF (PRESENT(rand2skip)) mc_par%rand2skip = rand2skip
842  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  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  ALLOCATE (mc_par)
858  ALLOCATE (mc_par%pmtraion_mol(1:nmol_types))
859  ALLOCATE (mc_par%pmtrans_mol(1:nmol_types))
860  ALLOCATE (mc_par%pmrot_mol(1:nmol_types))
861  ALLOCATE (mc_par%pmswap_mol(1:nmol_types))
862 
863  ALLOCATE (mc_par%eta(1:nmol_types))
864 
865  ALLOCATE (mc_par%rmbond(1:nmol_types))
866  ALLOCATE (mc_par%rmangle(1:nmol_types))
867  ALLOCATE (mc_par%rmdihedral(1:nmol_types))
868  ALLOCATE (mc_par%rmtrans(1:nmol_types))
869  ALLOCATE (mc_par%rmrot(1:nmol_types))
870 
871  ALLOCATE (mc_par%avbmc_atom(1:nmol_types))
872  ALLOCATE (mc_par%avbmc_rmin(1:nmol_types))
873  ALLOCATE (mc_par%avbmc_rmax(1:nmol_types))
874  ALLOCATE (mc_par%pmavbmc_mol(1:nmol_types))
875  ALLOCATE (mc_par%pbias(1:nmol_types))
876 
877  ALLOCATE (mc_par%mc_input_file)
878  ALLOCATE (mc_par%mc_bias_file)
879 
880  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  SUBROUTINE mc_sim_par_destroy(mc_par)
888 
889  TYPE(mc_simpar_type), POINTER :: mc_par
890 
891  DEALLOCATE (mc_par%mc_input_file)
892  DEALLOCATE (mc_par%mc_bias_file)
893 
894  DEALLOCATE (mc_par%pmtraion_mol)
895  DEALLOCATE (mc_par%pmtrans_mol)
896  DEALLOCATE (mc_par%pmrot_mol)
897  DEALLOCATE (mc_par%pmswap_mol)
898 
899  DEALLOCATE (mc_par%eta)
900 
901  DEALLOCATE (mc_par%rmbond)
902  DEALLOCATE (mc_par%rmangle)
903  DEALLOCATE (mc_par%rmdihedral)
904  DEALLOCATE (mc_par%rmtrans)
905  DEALLOCATE (mc_par%rmrot)
906 
907  DEALLOCATE (mc_par%avbmc_atom)
908  DEALLOCATE (mc_par%avbmc_rmin)
909  DEALLOCATE (mc_par%avbmc_rmax)
910  DEALLOCATE (mc_par%pmavbmc_mol)
911  DEALLOCATE (mc_par%pbias)
912  IF (mc_par%ensemble == "VIRIAL") THEN
913  DEALLOCATE (mc_par%virial_temps)
914  END IF
915 
916  DEALLOCATE (mc_par)
917  NULLIFY (mc_par)
918 
919  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  SUBROUTINE mc_input_file_create(mc_input_file, input_file_name, &
931  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  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  INTEGER, DIMENSION(:), POINTER :: nunits
945 
946 ! some stuff in case we need to write error messages to the screen
947 
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  nunits=nunits, atom_names=atom_names)
953  ALLOCATE (mc_input_file%coordinates_empty(1:3, 1:nunits(1)))
954  ALLOCATE (mc_input_file%atom_names_empty(1:nunits(1)))
955  DO iunit = 1, nunits(1)
956  mc_input_file%atom_names_empty(iunit) = atom_names(iunit, 1)
957  mc_input_file%coordinates_empty(1:3, iunit) = empty_coords(1:3, iunit)
958  END DO
959  mc_input_file%nunits_empty = nunits(1)
960 
961 ! allocate some arrays
962  ALLOCATE (mc_input_file%mol_set_nmol_row(1:nmol_types))
963  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  unit_number=unit, file_action='READ', file_status='OLD')
968 
969  nlines = 0
970  DO
971  READ (unit, '(A)', iostat=io_stat) line
972  IF (io_stat .NE. 0) EXIT
973  nlines = nlines + 1
974  END DO
975 
976  ALLOCATE (mc_input_file%text(1:nlines))
977 
978  rewind(unit)
979  DO iline = 1, nlines
980  READ (unit, '(A)') mc_input_file%text(iline)
981  END DO
982 
983 ! close the input file
984  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  mc_input_file%global_row_end, idum, start_row_number=mc_input_file%global_row_start)
990  IF (mc_input_file%global_row_end == 0) THEN
991  IF (iw > 0) THEN
992  WRITE (iw, *)
993  WRITE (iw, *) 'File name ', input_file_name
994  END IF
995  CALL cp_abort(__location__, &
996  '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  mc_input_file%motion_row_end, idum, start_row_number=mc_input_file%motion_row_start)
1004  IF (mc_input_file%motion_row_end == 0) THEN
1005  IF (iw > 0) THEN
1006  WRITE (iw, *)
1007  WRITE (iw, *) 'File name ', input_file_name, mc_input_file%motion_row_start
1008  END IF
1009  CALL cp_abort(__location__, &
1010  '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  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  mc_input_file%coord_row_end, idum)
1020  CALL mc_parse_text(mc_input_file%text, 1, nlines, "&COORD", .false., &
1021  mc_input_file%coord_row_start, idum)
1022 
1023  IF (mc_input_file%coord_row_end == 0) THEN
1024  IF (iw > 0) THEN
1025  WRITE (iw, *)
1026  WRITE (iw, *) 'File name ', input_file_name
1027  END IF
1028  CALL cp_abort(__location__, &
1029  '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  mc_input_file%run_type_row, mc_input_file%run_type_column)
1035  mc_input_file%run_type_column = mc_input_file%run_type_column + 9
1036  IF (mc_input_file%run_type_row == 0) THEN
1037  IF (iw > 0) THEN
1038  WRITE (iw, *)
1039  WRITE (iw, *) 'File name ', input_file_name
1040  END IF
1041  CALL cp_abort(__location__, &
1042  '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  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  "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  "&CELL", .false., cell_row, cell_column)
1056  IF (cell_row == 0) THEN ! nothing in between...we found the correct ABC
1057  mc_input_file%cell_row = abc_row
1058  mc_input_file%cell_column = abc_column + 4
1059  ELSE
1060  CALL mc_parse_text(mc_input_file%text, abc_row + 1, nlines, &
1061  "ABC", .false., mc_input_file%cell_row, mc_input_file%cell_column)
1062  END IF
1063  IF (mc_input_file%cell_row == 0) THEN
1064  IF (iw > 0) THEN
1065  WRITE (iw, *)
1066  WRITE (iw, *) 'File name ', input_file_name
1067  END IF
1068  CALL cp_abort(__location__, &
1069  '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  IF (.NOT. lhmc) THEN
1074  irep = 0
1075  nstart = 1
1076  DO
1077  CALL mc_parse_text(mc_input_file%text, nstart, nlines, "&MOLECULE", &
1078  .false., row_number, idum)
1079  IF (row_number == 0) EXIT
1080  nstart = row_number + 1
1081  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  mc_input_file%mol_set_nmol_column(irep))
1085  mc_input_file%mol_set_nmol_column(irep) = mc_input_file%mol_set_nmol_column(irep) + 5
1086 
1087  END DO
1088  IF (irep .NE. nmol_types) THEN
1089  IF (iw > 0) THEN
1090  WRITE (iw, *)
1091  WRITE (iw, *) 'File name ', input_file_name
1092  WRITE (iw, *) 'Number of molecule types ', nmol_types
1093  WRITE (iw, *) '&MOLECULE sections found ', irep
1094  END IF
1095  CALL cp_abort( &
1096  __location__, &
1097  '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  "METHOD", .false., line_number, idum)
1105  READ (mc_input_file%text(line_number), *) cdum, method_name
1106  CALL uppercase(method_name)
1107  SELECT CASE (trim(adjustl(method_name)))
1108  CASE ("FIST")
1109  mc_input_file%in_use = use_fist_force
1110  CASE ("QS", "QUICKSTEP")
1111  mc_input_file%in_use = use_qs_force
1112  CASE default
1113  cpabort('Cannot determine the type of force_env we are using (check METHOD)')
1114  END SELECT
1115 
1116  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  SUBROUTINE mc_input_file_destroy(mc_input_file)
1124 
1125  TYPE(mc_input_file_type), POINTER :: mc_input_file
1126 
1127  DEALLOCATE (mc_input_file%mol_set_nmol_row)
1128  DEALLOCATE (mc_input_file%mol_set_nmol_column)
1129  DEALLOCATE (mc_input_file%text)
1130  DEALLOCATE (mc_input_file%atom_names_empty)
1131  DEALLOCATE (mc_input_file%coordinates_empty)
1132 
1133  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  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  row_number = 0
1166  column_number = 0
1167 
1168 ! how long is our string to search for?
1169  string_size = len_trim(string_search)
1170  whole_file: DO iline = nstart, nend
1171 
1172  index_string = -1
1173  index_string = index(text(iline), string_search(1:string_size))
1174 
1175  IF (index_string .GT. 0) THEN ! we found one
1176  IF (.NOT. lend) THEN
1177  row_number = iline
1178  column_number = index_string
1179  ELSE
1180  IF (PRESENT(start_row_number)) start_row_number = iline
1181  column_number = index_string
1182  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  text_temp = text(jline)
1186  IF (text_temp(column_number:column_number) .EQ. "&") THEN
1187  row_number = jline
1188  EXIT
1189  END IF
1190  END DO
1191  END IF
1192  EXIT whole_file
1193  END IF
1194  END DO whole_file
1195 
1196  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  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  INTEGER, DIMENSION(:), POINTER :: temp_i_array
1226  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  cpassert(ASSOCIATED(input_file))
1232 
1233  CALL timeset(routinen, handle)
1234 
1235  NULLIFY (mc_section)
1236  mc_section => section_vals_get_subs_vals(input_file, &
1237  "MOTION%MC")
1238 
1239 ! need the input file sturcutre that we're reading from for when we make
1240 ! dat files
1241  mc_par%input_file => input_file
1242 
1243 ! set the ionode and mepos
1244  mc_par%ionode = para_env%is_source()
1245  mc_par%group = para_env
1246  mc_par%source = para_env%source
1247 
1248 ! for convenience
1249  nmol_types = mc_par%mc_molecule_info%nmol_types
1250 
1251 !..defaults...most are set in input_cp2k_motion
1252  mc_par%nstart = 0
1253  CALL section_vals_val_get(force_env_section, "METHOD", i_val=method_name_id)
1254  SELECT CASE (method_name_id)
1255  CASE (do_fist)
1256  mc_par%iprint = 100
1257  CASE (do_qs)
1258  mc_par%iprint = 1
1259  END SELECT
1260 
1262  IF (iw > 0) WRITE (iw, *)
1263 
1264 !..filenames
1265  mc_par%program = input_file_name
1266  CALL xstring(mc_par%program, ia, ie)
1267  mc_par%coords_file = mc_par%program(ia:ie)//'.coordinates'
1268  mc_par%molecules_file = mc_par%program(ia:ie)//'.molecules'
1269  mc_par%moves_file = mc_par%program(ia:ie)//'.moves'
1270  mc_par%energy_file = mc_par%program(ia:ie)//'.energy'
1271  mc_par%cell_file = mc_par%program(ia:ie)//'.cell'
1272  mc_par%displacement_file = mc_par%program(ia:ie) &
1273  //'.max_displacements'
1274  mc_par%data_file = mc_par%program(ia:ie)//'.data'
1275  stop_num = ie - 3
1276  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  c_val=mc_par%coords_file)
1281  CALL section_vals_val_set(mc_section, "DATA_FILE_NAME", &
1282  c_val=mc_par%data_file)
1283  CALL section_vals_val_set(mc_section, "CELL_FILE_NAME", &
1284  c_val=mc_par%cell_file)
1285  CALL section_vals_val_set(mc_section, "MAX_DISP_FILE_NAME", &
1286  c_val=mc_par%displacement_file)
1287  CALL section_vals_val_set(mc_section, "MOVES_FILE_NAME", &
1288  c_val=mc_par%moves_file)
1289  CALL section_vals_val_set(mc_section, "MOLECULES_FILE_NAME", &
1290  c_val=mc_par%molecules_file)
1291  CALL section_vals_val_set(mc_section, "ENERGY_FILE_NAME", &
1292  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  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  CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMHMC", r_val=mc_par%pmhmc)
1304  CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMSWAP", r_val=mc_par%pmswap)
1305  CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMVOLUME", r_val=mc_par%pmvolume)
1306  CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMAVBMC", r_val=mc_par%pmavbmc)
1307  CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMTRANS", r_val=mc_par%pmtrans)
1308  CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMTRAION", r_val=mc_par%pmtraion)
1309  CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMCLTRANS", r_val=mc_par%pmcltrans)
1310 
1311  ! first, grab all the integer values
1312  CALL section_vals_val_get(mc_section, "NSTEP", i_val=mc_par%nstep)
1313  CALL section_vals_val_get(mc_section, "NMOVES", i_val=mc_par%nmoves)
1314  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  i_val=mc_par%iupvolume)
1317  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  i_val=mc_par%iuptrans)
1320  CALL section_vals_val_get(mc_section, "MOVE_UPDATES%IUPCLTRANS", &
1321  i_val=mc_par%iupcltrans)
1322  CALL section_vals_val_get(mc_section, "IPRINT", i_val=mc_par%iprint)
1323 ! now an integer array
1324  CALL section_vals_val_get(mc_section, "AVBMC%AVBMC_ATOM", i_vals=temp_i_array)
1325 
1326  IF (mc_par%pmhmc - mc_par%pmswap >= 1.0_dp .AND. mc_par%pmhmc == 1.0_dp) THEN
1327  mc_par%lhmc = .true.
1328  ELSE
1329  mc_par%lhmc = .false.
1330  END IF
1331 
1332  IF (.NOT. mc_par%lhmc) THEN
1333  IF (SIZE(temp_i_array) .NE. nmol_types) THEN
1334  cpabort('AVBMC_ATOM must have as many values as there are molecule types.')
1335  END IF
1336  END IF
1337  DO imol = 1, SIZE(temp_i_array)
1338  mc_par%avbmc_atom(imol) = temp_i_array(imol)
1339  END DO
1340 
1341 ! now the real values
1342  CALL section_vals_val_get(mc_section, "PRESSURE", r_val=mc_par%pressure)
1343  CALL section_vals_val_get(mc_section, "TEMPERATURE", r_val=mc_par%temperature)
1344  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  r_val=mc_par%pmclus_box)
1347  CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%BOX_PROBABILITIES%PMHMC_BOX", &
1348  r_val=mc_par%pmhmc_box)
1349  CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%BOX_PROBABILITIES%PMVOL_BOX", &
1350  r_val=mc_par%pmvol_box)
1351  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  r_val=mc_par%rmvolume)
1354  CALL section_vals_val_get(mc_section, "MAX_DISPLACEMENTS%BOX_DISPLACEMENTS%RMCLTRANS", &
1355  r_val=mc_par%rmcltrans)
1356 
1357 ! finally the character values
1358  CALL section_vals_val_get(mc_section, "ENSEMBLE", c_val=mc_par%ensemble)
1359  CALL section_vals_val_get(mc_section, "RESTART_FILE_NAME", c_val=mc_par%restart_file_name)
1360  CALL section_vals_val_get(mc_section, "COORDINATE_FILE_NAME", c_val=mc_par%coords_file)
1361  CALL section_vals_val_get(mc_section, "ENERGY_FILE_NAME", c_val=mc_par%energy_file)
1362  CALL section_vals_val_get(mc_section, "MOVES_FILE_NAME", c_val=mc_par%moves_file)
1363  CALL section_vals_val_get(mc_section, "MOLECULES_FILE_NAME", c_val=mc_par%molecules_file)
1364  CALL section_vals_val_get(mc_section, "CELL_FILE_NAME", c_val=mc_par%cell_file)
1365  CALL section_vals_val_get(mc_section, "DATA_FILE_NAME", c_val=mc_par%data_file)
1366  CALL section_vals_val_get(mc_section, "MAX_DISP_FILE_NAME", c_val=mc_par%displacement_file)
1367  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  IF (mc_par%ensemble == "VIRIAL") THEN
1372  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  ALLOCATE (mc_par%virial_temps(1:SIZE(temp_r_array)))
1375  DO imol = 1, SIZE(temp_r_array)
1376  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  CALL section_vals_val_get(mc_section, "AVBMC%AVBMC_RMIN", r_vals=temp_r_array)
1381  IF (.NOT. mc_par%lhmc) THEN
1382  IF (SIZE(temp_r_array) .NE. nmol_types) THEN
1383  cpabort('AVBMC_RMIN must have as many values as there are molecule types.')
1384  END IF
1385  END IF
1386  DO imol = 1, SIZE(temp_r_array)
1387  mc_par%avbmc_rmin(imol) = temp_r_array(imol)
1388  END DO
1389  CALL section_vals_val_get(mc_section, "AVBMC%AVBMC_RMAX", r_vals=temp_r_array)
1390  IF (.NOT. mc_par%lhmc) THEN
1391  IF (SIZE(temp_r_array) .NE. nmol_types) THEN
1392  cpabort('AVBMC_RMAXL must have as many values as there are molecule types.')
1393  END IF
1394  END IF
1395  DO imol = 1, SIZE(temp_r_array)
1396  mc_par%avbmc_rmax(imol) = temp_r_array(imol)
1397  END DO
1398  CALL section_vals_val_get(mc_section, "AVBMC%PBIAS", r_vals=temp_r_array)
1399  IF (.NOT. mc_par%lhmc) THEN
1400  IF (SIZE(temp_r_array) .NE. nmol_types) THEN
1401  cpabort('PBIAS must have as many values as there are molecule types.')
1402  END IF
1403  END IF
1404  DO imol = 1, SIZE(temp_r_array)
1405  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  r_vals=temp_r_array)
1409  IF (.NOT. mc_par%lhmc) THEN
1410  IF (SIZE(temp_r_array) .NE. nmol_types) THEN
1411  cpabort('PMAVBMC_MOL must have as many values as there are molecule types.')
1412  END IF
1413  END IF
1414  DO imol = 1, SIZE(temp_r_array)
1415  mc_par%pmavbmc_mol(imol) = temp_r_array(imol)
1416  END DO
1417  CALL section_vals_val_get(mc_section, "ETA", r_vals=temp_r_array)
1418  IF (.NOT. mc_par%lhmc) THEN
1419  IF (SIZE(temp_r_array) .NE. nmol_types) THEN
1420  cpabort('ETA must have as many values as there are molecule types.')
1421  END IF
1422  END IF
1423  DO imol = 1, SIZE(temp_r_array)
1424  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  r_vals=temp_r_array)
1428  IF (.NOT. mc_par%lhmc) THEN
1429  IF (SIZE(temp_r_array) .NE. nmol_types) THEN
1430  cpabort('RMBOND must have as many values as there are molecule types.')
1431  END IF
1432  END IF
1433  DO imol = 1, SIZE(temp_r_array)
1434  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  r_vals=temp_r_array)
1438  IF (.NOT. mc_par%lhmc) THEN
1439  IF (SIZE(temp_r_array) .NE. nmol_types) THEN
1440  cpabort('RMANGLE must have as many values as there are molecule types.')
1441  END IF
1442  END IF
1443  DO imol = 1, SIZE(temp_r_array)
1444  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  r_vals=temp_r_array)
1448  IF (.NOT. mc_par%lhmc) THEN
1449  IF (SIZE(temp_r_array) .NE. nmol_types) THEN
1450  cpabort('RMDIHEDRAL must have as many values as there are molecule types.')
1451  END IF
1452  END IF
1453  DO imol = 1, SIZE(temp_r_array)
1454  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  r_vals=temp_r_array)
1458  IF (.NOT. mc_par%lhmc) THEN
1459  IF (SIZE(temp_r_array) .NE. nmol_types) THEN
1460  cpabort('RMROT must have as many values as there are molecule types.')
1461  END IF
1462  END IF
1463  DO imol = 1, SIZE(temp_r_array)
1464  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  r_vals=temp_r_array)
1468  IF (.NOT. mc_par%lhmc) THEN
1469  IF (SIZE(temp_r_array) .NE. nmol_types) THEN
1470  cpabort('RMTRANS must have as many values as there are molecule types.')
1471  END IF
1472  END IF
1473  DO imol = 1, SIZE(temp_r_array)
1474  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  r_vals=temp_r_array)
1479  IF (.NOT. mc_par%lhmc) THEN
1480  IF (SIZE(temp_r_array) .NE. nmol_types) THEN
1481  cpabort('PMTRAION_MOL must have as many values as there are molecule types.')
1482  END IF
1483  END IF
1484  DO imol = 1, SIZE(temp_r_array)
1485  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  r_vals=temp_r_array)
1490  IF (.NOT. mc_par%lhmc) THEN
1491  IF (SIZE(temp_r_array) .NE. nmol_types) THEN
1492  cpabort('PMTRANS_MOL must have as many values as there are molecule types.')
1493  END IF
1494  END IF
1495  DO imol = 1, SIZE(temp_r_array)
1496  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  r_vals=temp_r_array)
1501  IF (.NOT. mc_par%lhmc) THEN
1502  IF (SIZE(temp_r_array) .NE. nmol_types) THEN
1503  cpabort('PMROT_MOL must have as many values as there are molecule types.')
1504  END IF
1505  END IF
1506  DO imol = 1, SIZE(temp_r_array)
1507  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  r_vals=temp_r_array)
1512  IF (.NOT. mc_par%lhmc) THEN
1513  IF (SIZE(temp_r_array) .NE. nmol_types) THEN
1514  cpabort('PMSWAP_MOL must have as many values as there are molecule types.')
1515  END IF
1516  END IF
1517  DO imol = 1, SIZE(temp_r_array)
1518  mc_par%pmswap_mol(imol) = temp_r_array(imol)
1519  END DO
1520 
1521 ! now some logical values
1522  CALL section_vals_val_get(mc_section, "LBIAS", l_val=mc_par%lbias)
1523  CALL section_vals_val_get(mc_section, "LDISCRETE", l_val=mc_par%ldiscrete)
1524  CALL section_vals_val_get(mc_section, "LSTOP", l_val=mc_par%lstop)
1525  CALL section_vals_val_get(mc_section, "RESTART", l_val=mc_par%lrestart)
1526  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  IF (iw > 0) THEN
1532  WRITE (box_string, '(I2)') mc_par%mc_molecule_info%nboxes
1533  WRITE (molecule_string, '(I2)') mc_par%mc_molecule_info%nmol_types
1534  WRITE (tab_string, '(I4)') 81 - 10*mc_par%mc_molecule_info%nmol_types
1535  WRITE (tab_box_string, '(I4)') 81 - 10*mc_par%mc_molecule_info%nboxes
1536  WRITE (tab_string_int, '(I4)') 81 - 5*mc_par%mc_molecule_info%nmol_types
1537  format_string = "(A,T"//trim(adjustl(tab_string))//","//trim(adjustl(molecule_string))//"(2X,F8.4))"
1538  format_box_string = "(A,T"//trim(adjustl(tab_box_string))//","//trim(adjustl(box_string))//"(2X,F8.4))"
1539  format_string_int = "(A,T"//trim(adjustl(tab_string))//","//trim(adjustl(molecule_string))//"(I3,2X))"
1540  WRITE (iw, '( A )') ' MC| Monte Carlo Protocol '
1541  WRITE (iw, '( A,T71,I10 )') ' MC| total number of steps ', &
1542  mc_par%nstep
1543  WRITE (iw, '( A,T71,F10.3 )') ' MC| pmvolume ', &
1544  mc_par%pmvolume
1545  WRITE (iw, '( A,T71,F10.3 )') ' MC| pmvol_box ', &
1546  mc_par%pmvol_box
1547  WRITE (iw, '( A,T71,F10.3 )') ' MC| pmclus_box ', &
1548  mc_par%pmclus_box
1549  WRITE (iw, '( A,T71,F10.3 )') ' MC| pmhmc ', &
1550  mc_par%pmhmc
1551  WRITE (iw, '( A,T71,F10.3 )') ' MC| pmhmc_box ', &
1552  mc_par%pmhmc_box
1553  WRITE (iw, '( A,T71,F10.3 )') ' MC| pmswap ', &
1554  mc_par%pmswap
1555  WRITE (iw, '( A,T71,F10.3 )') ' MC| pmavbmc ', &
1556  mc_par%pmavbmc
1557  WRITE (iw, '( A,T71,F10.3 )') ' MC| pmtraion ', &
1558  mc_par%pmtraion
1559  WRITE (iw, '( A,T71,F10.3 )') ' MC| pmtrans ', &
1560  mc_par%pmtrans
1561  WRITE (iw, '( A,T71,F10.3 )') ' MC| pmcltrans ', &
1562  mc_par%pmcltrans
1563  WRITE (iw, '( A,T71,I10 )') ' MC| iupvolume ', &
1564  mc_par%iupvolume
1565  WRITE (iw, '( A,T71,I10 )') ' MC| nvirial ', &
1566  mc_par%nvirial
1567  WRITE (iw, '( A,T71,I10 )') ' MC| iuptrans ', &
1568  mc_par%iuptrans
1569  WRITE (iw, '( A,T71,I10 )') ' MC| iupcltrans ', &
1570  mc_par%iupcltrans
1571  WRITE (iw, '( A,T71,I10 )') ' MC| iprint ', &
1572  mc_par%iprint
1573  WRITE (iw, '( A,T61,A20 )') ' MC| ensemble ', &
1574  adjustr(mc_par%ensemble)
1575 ! format string may not be valid if all the molecules are atomic,
1576 ! like in HMC
1577  IF (.NOT. mc_par%lhmc) THEN
1578  WRITE (iw, format_string) ' MC| pmswap_mol ', &
1579  mc_par%pmswap_mol
1580  WRITE (iw, format_string) ' MC| pmavbmc_mol ', &
1581  mc_par%pmavbmc_mol
1582  WRITE (iw, format_string) ' MC| pbias ', &
1583  mc_par%pbias
1584  WRITE (iw, format_string) ' MC| pmtraion_mol ', &
1585  mc_par%pmtraion_mol
1586  WRITE (iw, format_string) ' MC| pmtrans_mol ', &
1587  mc_par%pmtrans_mol
1588  WRITE (iw, format_string) ' MC| pmrot_mol ', &
1589  mc_par%pmrot_mol
1590  WRITE (iw, format_string) ' MC| eta [K]', &
1591  mc_par%eta
1592  WRITE (iw, format_string) ' MC| rmbond [angstroms]', &
1593  mc_par%rmbond
1594  WRITE (iw, format_string) ' MC| rmangle [degrees]', &
1595  mc_par%rmangle
1596  WRITE (iw, format_string) ' MC| rmdihedral [degrees]', &
1597  mc_par%rmdihedral
1598  WRITE (iw, format_string) ' MC| rmtrans [angstroms]', &
1599  mc_par%rmtrans
1600  WRITE (iw, format_string) ' MC| rmcltrans [angstroms]', &
1601  mc_par%rmcltrans
1602  WRITE (iw, format_string) ' MC| rmrot [degrees]', &
1603  mc_par%rmrot
1604  WRITE (iw, format_string_int) ' MC| AVBMC target atom ', &
1605  mc_par%avbmc_atom
1606  WRITE (iw, format_string) ' MC| AVBMC inner cutoff [ang]', &
1607  mc_par%avbmc_rmin
1608  WRITE (iw, format_string) ' MC| AVBMC outer cutoff [ang]', &
1609  mc_par%avbmc_rmax
1610  END IF
1611  IF (mc_par%ensemble .EQ. 'GEMC-NVT') THEN
1612  WRITE (iw, '( A,T58,A)') ' MC| Box 2 file', &
1613  trim(mc_par%box2_file)
1614  END IF
1615  WRITE (iw, '( A,T58,A )') ' MC| Name of restart file:', &
1616  trim(mc_par%restart_file_name)
1617  WRITE (iw, '( A,A,T44,A )') ' MC| Name of output ', &
1618  'coordinate file:', &
1619  trim(mc_par%coords_file)
1620  WRITE (iw, '( A,T44,A )') ' MC| Name of output data file:', &
1621  trim(mc_par%data_file)
1622  WRITE (iw, '( A,A,T44,A )') ' MC| Name of output ', &
1623  'molecules file:', &
1624  trim(mc_par%molecules_file)
1625  WRITE (iw, '( A,T44,A )') ' MC| Name of output moves file:', &
1626  trim(mc_par%moves_file)
1627  WRITE (iw, '( A,T44,A )') ' MC| Name of output energy file:', &
1628  trim(mc_par%energy_file)
1629  WRITE (iw, '( A,T44,A )') ' MC| Name of output cell file:', &
1630  trim(mc_par%cell_file)
1631  WRITE (iw, '( A,A,T44,A )') ' MC| Name of output', &
1632  ' displacement file:', &
1633  trim(mc_par%displacement_file)
1634  IF (mc_par%ldiscrete) THEN
1635  WRITE (iw, '(A,A,T71,F10.3)') ' MC| discrete step size', &
1636  '[angstroms]', &
1637  mc_par%discrete_step
1638  ELSE
1639  WRITE (iw, '( A,A,T71,F10.3 )') ' MC| rmvolume ', &
1640  '[cubic angstroms]', &
1641  mc_par%rmvolume
1642  END IF
1643  WRITE (iw, '( A,T71,F10.2 )') ' MC| Temperature [K] ', &
1644  mc_par%temperature
1645  WRITE (iw, '( A,T71,F10.5 )') ' MC| Pressure [bar] ', &
1646  mc_par%pressure
1647  WRITE (iw, '( A,T71,F10.5 )') ' MC| Rclus [ang] ', &
1648  mc_par%rclus
1649  IF (mc_par%lrestart) THEN
1650  WRITE (iw, '(A,A)') ' MC| Initial data will be read from a', &
1651  ' restart file.'
1652  END IF
1653  IF (mc_par%lbias) THEN
1654  WRITE (iw, '(A)') ' MC| The moves will be biased.'
1655  ELSE
1656  WRITE (iw, '(A)') ' MC| The moves will not be biased,'
1657  END IF
1658  IF (mc_par%nmoves .EQ. 1) THEN
1659  WRITE (iw, '(A,A)') ' MC| A full energy calculation ', &
1660  'will be done at every step.'
1661  ELSE
1662  WRITE (iw, '(A,I4,A,A)') ' MC| ', mc_par%nmoves, &
1663  ' moves will be attempted ', &
1664  'before a Quickstep energy calculation'
1665  WRITE (iw, '(A)') ' MC| takes place.'
1666  END IF
1667 
1668  WRITE (iw, '(A,I4,A,A)') ' MC| ', mc_par%nswapmoves, &
1669  ' swap insertions will be attempted ', &
1670  'per molecular swap move'
1671 
1672  IF (mc_par%lhmc) THEN
1673  WRITE (iw, '(A)') '********************************************************'
1674  WRITE (iw, '(A)') '************** ONLY DOING HYBRID MONTE CARLO MOVES **************************'
1675  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  mc_par%BETA = 1/mc_par%temperature/boltzmann*joule
1682 
1683  ! 0.9_dp is to avoid overflow or underflow
1684  mc_par%exp_max_val = 0.9_dp*log(huge(0.0_dp))
1685  mc_par%exp_min_val = 0.9_dp*log(tiny(0.0_dp))
1686  mc_par%max_val = exp(mc_par%exp_max_val)
1687  mc_par%min_val = 0.0_dp
1688 
1689 ! convert from bar to a.u.
1690  mc_par%pressure = cp_unit_to_cp2k(mc_par%pressure, "bar")
1691  mc_par%rclus = cp_unit_to_cp2k(mc_par%rclus, "angstrom")
1692 ! convert from angstrom to a.u.
1693  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  mc_par%rmrot(itype) = mc_par%rmrot(itype)/180.0e0_dp*pi
1698  mc_par%rmangle(itype) = mc_par%rmangle(itype)/180.0e0_dp*pi
1699  mc_par%rmdihedral(itype) = mc_par%rmdihedral(itype)/180.0e0_dp*pi
1700  mc_par%rmtrans(itype) = cp_unit_to_cp2k(mc_par%rmtrans(itype), "angstrom")
1701  mc_par%rmbond(itype) = cp_unit_to_cp2k(mc_par%rmbond(itype), "angstrom")
1702  mc_par%eta(itype) = cp_unit_to_cp2k(mc_par%eta(itype), "K_e")
1703  mc_par%avbmc_rmin(itype) = cp_unit_to_cp2k(mc_par%avbmc_rmin(itype), "angstrom")
1704  mc_par%avbmc_rmax(itype) = cp_unit_to_cp2k(mc_par%avbmc_rmax(itype), "angstrom")
1705  END DO
1706  mc_par%rmvolume = cp_unit_to_cp2k(mc_par%rmvolume, "angstrom^3")
1707  mc_par%rmcltrans = cp_unit_to_cp2k(mc_par%rmcltrans, "angstrom")
1708  mc_par%discrete_step = cp_unit_to_cp2k(mc_par%discrete_step, "angstrom")
1709 
1710 ! end the timing
1711  CALL timestop(handle)
1712 
1713  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  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  NULLIFY (cell, potparm, fist_nonbond_env, fist_env)
1742 
1743  lterminate = .false.
1744  CALL force_env_get(force_env, fist_env=fist_env)
1745  CALL fist_env_get(fist_env, cell=cell, &
1746  fist_nonbond_env=fist_nonbond_env)
1747  CALL fist_nonbond_env_get(fist_nonbond_env, potparm=potparm)
1748  CALL get_cell(cell, abc=abc)
1749 
1750 ! find the largest value of rcutsq
1751  rcutsq_max = 0.0e0_dp
1752  DO itype = 1, SIZE(potparm%pot, 1)
1753  DO jtype = itype, SIZE(potparm%pot, 2)
1754  IF (potparm%pot(itype, jtype)%pot%rcutsq .GT. rcutsq_max) &
1755  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  mc_par%rcut = rcutsq_max**0.5_dp
1762  DO itype = 1, 3
1763  IF (abc(itype) .LT. 2.0_dp*mc_par%rcut) THEN
1764  lterminate = .true.
1765  END IF
1766  END DO
1767 
1768  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  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  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  INTEGER, DIMENSION(:), POINTER :: molecule_list
1800  LOGICAL :: lnew_type
1801  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  POINTER :: molecule_kinds
1806  TYPE(molecule_kind_type), POINTER :: molecule_kind
1807  TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
1808 
1809 ! how many simulation boxes do we have?
1810 
1811  nboxes = SIZE(force_env(:))
1812 
1813 ! allocate the pointer
1814  ALLOCATE (mc_molecule_info)
1815  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  IF (PRESENT(box_number)) THEN
1820  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  NULLIFY (subsys, molecule_kinds, molecule_kind)
1831 
1832  ALLOCATE (particles(1:nboxes))
1833  ALLOCATE (molecule_kinds(1:nboxes))
1834 
1835  ! find out how many types we have
1836  ntypes = 0
1837  DO ibox = 1, nboxes
1838  IF (ibox == skip_box) cycle
1839  CALL force_env_get(force_env(ibox)%force_env, &
1840  subsys=subsys)
1841  CALL cp_subsys_get(subsys, &
1842  molecule_kinds=molecule_kinds(ibox)%list, &
1843  particles=particles(ibox)%list)
1844  ntypes = ntypes + SIZE(molecule_kinds(ibox)%list%els(:))
1845  END DO
1846 
1847  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  itype = 1
1852  natoms_large = 0
1853  DO ibox = 1, nboxes
1854  IF (ibox == skip_box) cycle
1855  DO imolecule = 1, SIZE(molecule_kinds(ibox)%list%els(:))
1856  molecule_kind => molecule_kinds(ibox)%list%els(imolecule)
1857  CALL get_molecule_kind(molecule_kind, name=names_init(itype), &
1858  natom=natom)
1859  IF (natom .GT. natoms_large) natoms_large = natom
1860  itype = itype + 1
1861  END DO
1862  END DO
1863 
1864  nmol_types = 0
1865  DO itype = 1, ntypes
1866  lnew_type = .true.
1867  DO jtype = 1, itype - 1
1868  IF (trim(names_init(itype)) .EQ. trim(names_init(jtype))) &
1869  lnew_type = .false.
1870  END DO
1871  IF (lnew_type) THEN
1872  nmol_types = nmol_types + 1
1873  ELSE
1874  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  mc_molecule_info%nmol_types = nmol_types
1881  ALLOCATE (mc_molecule_info%names(1:nmol_types))
1882  ALLOCATE (mc_molecule_info%nunits(1:nmol_types))
1883  ALLOCATE (mc_molecule_info%conf_prob(1:3, 1:nmol_types))
1884  ALLOCATE (mc_molecule_info%nchains(1:nmol_types, 1:nboxes))
1885  ALLOCATE (mc_molecule_info%nunits_tot(1:nboxes))
1886  ALLOCATE (mc_molecule_info%atom_names(1:natoms_large, 1:nmol_types))
1887  ALLOCATE (mc_molecule_info%mass(1:natoms_large, 1:nmol_types))
1888 
1889 ! now record all the information for every molecule type
1890  iunique = 0
1891  itype = 0
1892  DO ibox = 1, nboxes
1893  IF (ibox == skip_box) cycle
1894  DO imolecule = 1, SIZE(molecule_kinds(ibox)%list%els(:))
1895  itype = itype + 1
1896  IF (names_init(itype) .EQ. '') cycle
1897  iunique = iunique + 1
1898  mc_molecule_info%names(iunique) = names_init(itype)
1899  molecule_kind => molecule_kinds(ibox)%list%els(imolecule)
1900  CALL get_molecule_kind(molecule_kind, natom=mc_molecule_info%nunits(iunique), &
1901  nbond=nbond, nbend=nbend, ntorsion=ntorsion, atom_list=atom_list)
1902 
1903  DO iatom = 1, mc_molecule_info%nunits(iunique)
1904  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  mass=mc_molecule_info%mass(iatom, iunique))
1908  END DO
1909 
1910 ! compute the probabilities of doing any particular kind of conformation change
1911  total = nbond + nbend + ntorsion
1912 
1913  IF (total == 0) THEN
1914  mc_molecule_info%conf_prob(:, iunique) = 0.0e0_dp
1915  ELSE
1916  mc_molecule_info%conf_prob(1, iunique) = real(nbond, dp)/real(total, dp)
1917  mc_molecule_info%conf_prob(2, iunique) = real(nbend, dp)/real(total, dp)
1918  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  IF (PRESENT(coordinates_empty)) THEN
1926  ALLOCATE (coordinates_empty(1:3, 1:mc_molecule_info%nunits(1)))
1927  DO iunit = 1, mc_molecule_info%nunits(1)
1928  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  mc_molecule_info%nchains(:, :) = 0
1935  DO iunique = 1, nmol_types
1936  DO ibox = 1, nboxes
1937  IF (ibox == skip_box) cycle
1938  DO imolecule = 1, SIZE(molecule_kinds(ibox)%list%els(:))
1939  molecule_kind => molecule_kinds(ibox)%list%els(imolecule)
1940  CALL get_molecule_kind(molecule_kind, nmolecule=nmolecule, &
1941  name=name)
1942  IF (trim(name) .NE. mc_molecule_info%names(iunique)) cycle
1943  mc_molecule_info%nchains(iunique, ibox) = mc_molecule_info%nchains(iunique, ibox) + nmolecule
1944  END DO
1945  END DO
1946  END DO
1947  mc_molecule_info%nchain_total = 0
1948  DO ibox = 1, nboxes
1949  mc_molecule_info%nunits_tot(ibox) = 0
1950  IF (ibox == skip_box) cycle
1951  DO iunique = 1, nmol_types
1952  mc_molecule_info%nunits_tot(ibox) = mc_molecule_info%nunits_tot(ibox) + &
1953  mc_molecule_info%nunits(iunique)*mc_molecule_info%nchains(iunique, ibox)
1954  END DO
1955  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  ALLOCATE (mc_molecule_info%mol_type(1:mc_molecule_info%nchain_total))
1961  ALLOCATE (mc_molecule_info%in_box(1:mc_molecule_info%nchain_total))
1962 
1963  last_mol = 0
1964  DO ibox = 1, nboxes
1965  IF (ibox == skip_box) cycle
1966  first_mol = last_mol + 1
1967  last_mol = first_mol + sum(mc_molecule_info%nchains(:, ibox)) - 1
1968  mc_molecule_info%in_box(first_mol:last_mol) = ibox
1969  DO imolecule = 1, SIZE(molecule_kinds(ibox)%list%els(:))
1970  molecule_kind => molecule_kinds(ibox)%list%els(imolecule)
1971  CALL get_molecule_kind(molecule_kind, nmolecule=nmolecule, &
1972  name=name, molecule_list=molecule_list)
1973 ! figure out which molecule number this is
1974  this_molecule = 0
1975  DO iunique = 1, nmol_types
1976  IF (trim(name) .EQ. trim(mc_molecule_info%names(iunique))) THEN
1977  this_molecule = iunique
1978  END IF
1979  END DO
1980  DO imol = 1, SIZE(molecule_list(:))
1981  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  DEALLOCATE (names_init)
1988  DEALLOCATE (molecule_kinds)
1989  DEALLOCATE (particles)
1990 
1991  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  SUBROUTINE mc_molecule_info_destroy(mc_molecule_info)
2001 
2002  TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
2003 
2004  DEALLOCATE (mc_molecule_info%nchains)
2005  DEALLOCATE (mc_molecule_info%nunits)
2006  DEALLOCATE (mc_molecule_info%mol_type)
2007  DEALLOCATE (mc_molecule_info%in_box)
2008  DEALLOCATE (mc_molecule_info%names)
2009  DEALLOCATE (mc_molecule_info%atom_names)
2010  DEALLOCATE (mc_molecule_info%conf_prob)
2011  DEALLOCATE (mc_molecule_info%nunits_tot)
2012  DEALLOCATE (mc_molecule_info%mass)
2013 
2014  DEALLOCATE (mc_molecule_info)
2015  NULLIFY (mc_molecule_info)
2016 
2017  END SUBROUTINE mc_molecule_info_destroy
2018 
2019 ! **************************************************************************************************
2020 !> \brief ...
2021 !> \param mc_par ...
2022 ! **************************************************************************************************
2023  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  INTEGER, DIMENSION(:), POINTER :: nunits
2037  LOGICAL :: lhmc
2038  TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
2039 
2040  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  nboxes=nboxes, nunits=nunits, nchain_total=nchain_total)
2044 
2045 ! if we're only doing HMC, we can skip these checks
2046  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  IF (mc_par%pmavbmc_mol(nmol_types) .LT. 0.9999e0_dp) THEN
2051  cpabort('The last value of PMAVBMC_MOL needs to be 1.0')
2052  END IF
2053  IF (mc_par%pmswap_mol(nmol_types) .LT. 0.9999e0_dp) THEN
2054  cpabort('The last value of PMSWAP_MOL needs to be 1.0')
2055  END IF
2056  IF (mc_par%pmtraion_mol(nmol_types) .LT. 0.9999e0_dp) THEN
2057  cpabort('The last value of PMTRAION_MOL needs to be 1.0')
2058  END IF
2059  IF (mc_par%pmtrans_mol(nmol_types) .LT. 0.9999e0_dp) THEN
2060  cpabort('The last value of PMTRANS_MOL needs to be 1.0')
2061  END IF
2062  IF (mc_par%pmrot_mol(nmol_types) .LT. 0.9999e0_dp) THEN
2063  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  SELECT CASE (mc_par%ensemble)
2070  CASE ("GEMC_NPT")
2071  IF (nmol_types .LE. 1) &
2072  cpabort('Cannot have GEMC-NPT simulation with only one molecule type')
2073  IF (nboxes .LE. 1) &
2074  cpabort('Cannot have GEMC-NPT simulation with only one box')
2075  CASE ("GEMC_NVT")
2076  IF (nboxes .LE. 1) &
2077  cpabort('Cannot have GEMC-NVT simulation with only one box')
2078  CASE ("TRADITIONAL")
2079  IF (mc_par%pmswap .GT. 0.0e0_dp) &
2080  cpabort('You cannot do swap moves in a system with only one box')
2081  CASE ("VIRIAL")
2082  IF (nchain_total .NE. 2) &
2083  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  DO imol_type = 1, nmol_types
2089  IF (mc_par%avbmc_atom(imol_type) .GT. nunits(imol_type)) THEN
2090  cpabort('Cannot have avbmc_atom higher than the number of atoms for that molecule type!')
2091  END IF
2092  END DO
2093 
2094  END SUBROUTINE mc_input_parameters_check
2095 
2096 END MODULE mc_types
2097 
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Handles all functions related to the CELL.
Definition: cell_types.F:15
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition: cell_types.F:195
Utility routines to open and close files. Tracking of preconnections.
Definition: cp_files.F:16
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition: cp_files.F:308
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition: cp_files.F:119
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
types that represent a subsys, i.e. a part of the system
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
unit conversion facility
Definition: cp_units.F:30
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Definition: cp_units.F:1150
Definition: fft_lib.F:7
subroutine, public fist_env_get(fist_env, atomic_kind_set, particle_set, ewald_pw, local_particles, local_molecules, molecule_kind_set, molecule_set, cell, cell_ref, ewald_env, fist_nonbond_env, thermo, para_env, subsys, qmmm, qmmm_env, input, shell_model, shell_model_ad, shell_particle_set, core_particle_set, multipoles, results, exclusions, efield)
Purpose: Get the FIST environment.
subroutine, public fist_nonbond_env_get(fist_nonbond_env, potparm14, potparm, nonbonded, rlist_cut, rlist_lowsq, aup, lup, ei_scale14, vdw_scale14, shift_cutoff, do_electrostatics, r_last_update, r_last_update_pbc, rshell_last_update_pbc, rcore_last_update_pbc, cell_last_update, num_update, last_update, counter, natom_types, long_range_correction, ij_kind_full_fac, eam_data, quip_data, nequip_data, allegro_data, deepmd_data, charges)
sets a fist_nonbond_env
Interface for the force calculations.
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env)
returns various attributes about the force environment
integer, parameter, public use_qs_force
integer, parameter, public use_fist_force
Define type storing the global information of a run. Keep the amount of stored data small....
Definition: global_types.F:21
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_fist
integer, parameter, public do_qs
objects that represent the structure of input sections and the data contained in an input section
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
integer, parameter, public default_path_length
Definition: kinds.F:58
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public pi
holds all the structure types needed for Monte Carlo, except the mc_environment_type
Definition: mc_types.F:15
subroutine, public get_mc_molecule_info(mc_molecule_info, nmol_types, nchain_total, nboxes, names, conf_prob, nchains, nunits, mol_type, nunits_tot, in_box, atom_names, mass)
...
Definition: mc_types.F:554
subroutine, public read_mc_section(mc_par, para_env, globenv, input_file_name, input_file, force_env_section)
reads in the Monte Carlo simulation parameters from an input file
Definition: mc_types.F:1209
subroutine, public mc_determine_molecule_info(force_env, mc_molecule_info, box_number, coordinates_empty)
figures out the number of total molecules, the number of atoms in each molecule, an array with the mo...
Definition: mc_types.F:1787
subroutine, public mc_molecule_info_destroy(mc_molecule_info)
deallocates all the arrays in the mc_molecule_info_type
Definition: mc_types.F:2001
subroutine, public mc_sim_par_destroy(mc_par)
destroys (deallocates) the mc_simulation_parameters type
Definition: mc_types.F:888
subroutine, public set_mc_par(mc_par, rm, cl, diff, nstart, rmvolume, rmcltrans, rmbond, rmangle, rmdihedral, rmrot, rmtrans, PROGRAM, nmoves, nswapmoves, lstop, temperature, pressure, rclus, iuptrans, iupcltrans, iupvolume, pmswap, pmvolume, pmtraion, pmtrans, pmcltrans, BETA, rcut, iprint, lbias, nstep, lrestart, ldiscrete, discrete_step, pmavbmc, mc_molecule_info, pmavbmc_mol, pmtrans_mol, pmrot_mol, pmtraion_mol, pmswap_mol, avbmc_rmin, avbmc_rmax, avbmc_atom, pbias, ensemble, pmvol_box, pmclus_box, eta, mc_input_file, mc_bias_file, exp_max_val, exp_min_val, min_val, max_val, pmhmc, pmhmc_box, lhmc, ionode, source, group, rand2skip)
changes the private elements of the mc_parameters_type
Definition: mc_types.F:667
subroutine, public mc_input_parameters_check(mc_par)
...
Definition: mc_types.F:2024
subroutine, public get_mc_input_file(mc_input_file, run_type_row, run_type_column, coord_row_start, coord_row_end, cell_row, cell_column, force_eval_row_start, force_eval_row_end, global_row_start, global_row_end, mol_set_nmol_row, mol_set_nmol_column, in_use, text, atom_names_empty, nunits_empty, coordinates_empty, motion_row_end, motion_row_start)
accesses the private elements of the mc_input_file_type
Definition: mc_types.F:279
subroutine, public mc_input_file_destroy(mc_input_file)
destroys (deallocates) things in the mc_input_file_type
Definition: mc_types.F:1124
subroutine, public find_mc_rcut(mc_par, force_env, lterminate)
finds the largest interaction cutoff value in a classical simulation so we know the smallest size we ...
Definition: mc_types.F:1728
subroutine, public mc_sim_par_create(mc_par, nmol_types)
creates (allocates) the mc_simulation_parameters type
Definition: mc_types.F:851
subroutine, public mc_input_file_create(mc_input_file, input_file_name, mc_molecule_info, empty_coords, lhmc)
creates (allocates) the mc_input_file_type
Definition: mc_types.F:932
subroutine, public get_mc_par(mc_par, nstep, nvirial, iuptrans, iupcltrans, iupvolume, nmoves, nswapmoves, rm, cl, diff, nstart, source, group, lbias, ionode, lrestart, lstop, rmvolume, rmcltrans, rmbond, rmangle, rmrot, rmtrans, temperature, pressure, rclus, BETA, pmswap, pmvolume, pmtraion, pmtrans, pmcltrans, ensemble, PROGRAM, restart_file_name, molecules_file, moves_file, coords_file, energy_file, displacement_file, cell_file, dat_file, data_file, box2_file, fft_lib, iprint, rcut, ldiscrete, discrete_step, pmavbmc, pbias, avbmc_atom, avbmc_rmin, avbmc_rmax, rmdihedral, input_file, mc_molecule_info, pmswap_mol, pmavbmc_mol, pmtrans_mol, pmrot_mol, pmtraion_mol, mc_input_file, mc_bias_file, pmvol_box, pmclus_box, virial_temps, exp_min_val, exp_max_val, min_val, max_val, eta, pmhmc, pmhmc_box, lhmc, rand2skip)
...
Definition: mc_types.F:405
Interface to the message passing library MPI.
represent a simple array based list of the given type
Define the molecule kind structure types and the corresponding functionality.
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
represent a simple array based list of the given type
Definition of physical constants:
Definition: physcon.F:68
real(kind=dp), parameter, public boltzmann
Definition: physcon.F:129
real(kind=dp), parameter, public joule
Definition: physcon.F:159
Utilities for string manipulations.
elemental subroutine, public xstring(string, ia, ib)
...
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.