(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
18 USE cell_types, ONLY: cell_type,&
20 USE cp_files, ONLY: close_file,&
25 USE cp_units, ONLY: cp_unit_to_cp2k
36 USE input_constants, ONLY: do_fist,&
37 do_qs
42 USE kinds, ONLY: default_path_length,&
44 dp
45 USE mathconstants, ONLY: pi
46 USE message_passing, ONLY: mp_comm_type,&
54 USE physcon, ONLY: boltzmann,&
55 joule
56 USE string_utilities, ONLY: uppercase,&
58#include "../../base/base_uses.f90"
59
60 IMPLICIT NONE
61
62 PRIVATE
63
64! *** Global parameters ***
65
66 PUBLIC :: mc_simpar_type, &
79
80! **************************************************************************************************
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! **************************************************************************************************
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! **************************************************************************************************
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! **************************************************************************************************
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! **************************************************************************************************
185! a pointer for multiple boxes
186 TYPE(mc_simpar_type), POINTER :: mc_par => null()
188
189! **************************************************************************************************
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! **************************************************************************************************
200 TYPE(mc_averages_type), POINTER :: averages => null()
201 END TYPE mc_averages_p_type
202
203! **************************************************************************************************
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! **************************************************************************************************
233 INTEGER :: successes = 0
234 INTEGER :: qsuccesses = 0
235 INTEGER :: attempts = 0
236 END TYPE accattempt
237
238! **************************************************************************************************
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
246CONTAINS
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
2096END 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
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....
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.
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_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
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 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 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
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.
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
represents a system: atoms, molecules, their pos,vel,...
allows for the creation of an array of force_env
wrapper to abstract the force evaluation of the various methods
contains the initially parsed file and the initial parallel environment
stores all the informations relevant to an mpi environment