(git:b279b6b)
mc_ensembles.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 Used to run the bulk of the MC simulation, doing things like
10 !> choosing move types and writing data to files
11 !> \author Matthew J. McGrath (09.26.2003)
12 !>
13 !> REVISIONS
14 !> 09.10.05 MJM combined the two subroutines in this module into one
15 ! **************************************************************************************************
16 MODULE mc_ensembles
17  USE cell_types, ONLY: cell_p_type,&
18  get_cell
19  USE cp_external_control, ONLY: external_control
20  USE cp_files, ONLY: close_file,&
21  open_file
22  USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
23  USE cp_subsys_types, ONLY: cp_subsys_get,&
24  cp_subsys_p_type,&
25  cp_subsys_type
26  USE cp_units, ONLY: cp_unit_from_cp2k
27  USE force_env_methods, ONLY: force_env_calc_energy_force
28  USE force_env_types, ONLY: force_env_get,&
29  force_env_p_type,&
30  force_env_release
31  USE global_types, ONLY: global_environment_type
32  USE input_constants, ONLY: dump_xmol
33  USE input_section_types, ONLY: section_type,&
34  section_vals_type,&
35  section_vals_val_get
36  USE kinds, ONLY: default_string_length,&
37  dp
38  USE machine, ONLY: m_flush
39  USE mathconstants, ONLY: pi
40  USE mc_control, ONLY: mc_create_bias_force_env,&
41  write_mc_restart
42  USE mc_coordinates, ONLY: check_for_overlap,&
43  create_discrete_array,&
44  find_mc_test_molecule,&
45  get_center_of_mass,&
46  mc_coordinate_fold,&
47  rotate_molecule
48  USE mc_environment_types, ONLY: get_mc_env,&
49  mc_environment_p_type,&
50  set_mc_env
51  USE mc_ge_moves, ONLY: mc_ge_swap_move,&
52  mc_ge_volume_move,&
53  mc_quickstep_move
54  USE mc_misc, ONLY: final_mc_write,&
55  mc_averages_create,&
56  mc_averages_release
57  USE mc_move_control, ONLY: init_mc_moves,&
58  mc_move_update,&
59  mc_moves_release,&
60  write_move_stats
61  USE mc_moves, ONLY: mc_avbmc_move,&
62  mc_cluster_translation,&
63  mc_conformation_change,&
64  mc_hmc_move,&
65  mc_molecule_rotation,&
66  mc_molecule_translation,&
67  mc_volume_move
68  USE mc_types, ONLY: get_mc_molecule_info,&
69  get_mc_par,&
70  mc_averages_p_type,&
71  mc_input_file_type,&
72  mc_molecule_info_type,&
73  mc_moves_p_type,&
74  mc_simulation_parameters_p_type,&
75  set_mc_par
76  USE message_passing, ONLY: mp_comm_type,&
77  mp_para_env_type
78  USE parallel_rng_types, ONLY: rng_stream_type
79  USE particle_list_types, ONLY: particle_list_p_type,&
80  particle_list_type
81  USE particle_methods, ONLY: write_particle_coordinates
82  USE physcon, ONLY: angstrom,&
83  boltzmann,&
84  joule,&
85  n_avogadro
86 #include "../../base/base_uses.f90"
87 
88  IMPLICIT NONE
89 
90  PRIVATE
91 
92 ! *** Global parameters ***
93 
94  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_ensembles'
95  LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
96 
97  PUBLIC :: mc_run_ensemble, mc_compute_virial
98 
99 CONTAINS
100 
101 ! **************************************************************************************************
102 !> \brief directs the program in running one or two box MC simulations
103 !> \param mc_env a pointer that contains all mc_env for all the simulation
104 !> boxes
105 !> \param para_env ...
106 !> \param globenv the global environment for the simulation
107 !> \param input_declaration ...
108 !> \param nboxes the number of simulation boxes
109 !> \param rng_stream the stream we pull random numbers from
110 !>
111 !> Suitable for parallel.
112 !> \author MJM
113 ! **************************************************************************************************
114  SUBROUTINE mc_run_ensemble(mc_env, para_env, globenv, input_declaration, nboxes, rng_stream)
115 
116  TYPE(mc_environment_p_type), DIMENSION(:), POINTER :: mc_env
117  TYPE(mp_para_env_type), POINTER :: para_env
118  TYPE(global_environment_type), POINTER :: globenv
119  TYPE(section_type), POINTER :: input_declaration
120  INTEGER, INTENT(IN) :: nboxes
121  TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
122 
123  CHARACTER(len=*), PARAMETER :: routinen = 'mc_run_ensemble'
124 
125  CHARACTER(default_string_length), ALLOCATABLE, &
126  DIMENSION(:) :: atom_names_box
127  CHARACTER(default_string_length), &
128  DIMENSION(:, :), POINTER :: atom_names
129  CHARACTER(LEN=20) :: ensemble
130  CHARACTER(LEN=40) :: cbox, cstep, fft_lib, move_type, &
131  move_type_avbmc
132  INTEGER, DIMENSION(:, :), POINTER :: nchains
133  INTEGER, DIMENSION(:), POINTER :: avbmc_atom, mol_type, nchains_box, &
134  nunits, nunits_tot
135  INTEGER, DIMENSION(1:nboxes) :: box_flag, cl, data_unit, diff, istep, &
136  move_unit, rm
137  INTEGER, DIMENSION(1:3, 1:2) :: discrete_array
138  INTEGER :: atom_number, box_number, cell_unit, com_crd, com_ene, com_mol, end_mol, handle, &
139  ibox, idum, imol_type, imolecule, imove, iparticle, iprint, itype, iunit, iuptrans, &
140  iupvolume, iw, jbox, jdum, molecule_type, molecule_type_swap, molecule_type_target, &
141  nchain_total, nmol_types, nmoves, nnstep, nstart, nstep, source, start_atom, &
142  start_atom_swap, start_atom_target, start_mol
143  CHARACTER(LEN=default_string_length) :: unit_str
144  CHARACTER(LEN=40), DIMENSION(1:nboxes) :: cell_file, coords_file, data_file, &
145  displacement_file, energy_file, &
146  molecules_file, moves_file
147  LOGICAL :: ionode, lbias, ldiscrete, lhmc, &
148  lnew_bias_env, loverlap, lreject, &
149  lstop, print_kind, should_stop
150  REAL(dp), DIMENSION(:), POINTER :: pbias, pmavbmc_mol, pmclus_box, &
151  pmhmc_box, pmrot_mol, pmtraion_mol, &
152  pmtrans_mol, pmvol_box
153  REAL(dp), DIMENSION(:, :), POINTER :: conf_prob, mass
154  REAL(kind=dp) :: discrete_step, pmavbmc, pmcltrans, &
155  pmhmc, pmswap, pmtraion, pmtrans, &
156  pmvolume, rand, test_energy, unit_conv
157  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: r_temp
158  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r_old
159  REAL(kind=dp), DIMENSION(1:3, 1:nboxes) :: abc
160  REAL(kind=dp), DIMENSION(1:nboxes) :: bias_energy, energy_check, final_energy, &
161  initial_energy, last_bias_energy, &
162  old_energy
163  TYPE(cell_p_type), DIMENSION(:), POINTER :: cell
164  TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: oldsys
165  TYPE(cp_subsys_type), POINTER :: biassys
166  TYPE(force_env_p_type), DIMENSION(:), POINTER :: bias_env, force_env
167  TYPE(mc_averages_p_type), DIMENSION(:), POINTER :: averages
168  TYPE(mc_input_file_type), POINTER :: mc_bias_file
169  TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
170  TYPE(mc_moves_p_type), DIMENSION(:), POINTER :: test_moves
171  TYPE(mc_moves_p_type), DIMENSION(:, :), POINTER :: move_updates, moves
172  TYPE(mc_simulation_parameters_p_type), &
173  DIMENSION(:), POINTER :: mc_par
174  TYPE(mp_comm_type) :: group
175  TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles_old
176  TYPE(particle_list_type), POINTER :: particles_bias
177  TYPE(section_vals_type), POINTER :: root_section
178 
179  CALL timeset(routinen, handle)
180 
181  ! nullify some pointers
182  NULLIFY (moves, move_updates, test_moves, root_section)
183 
184  ! allocate a whole bunch of stuff based on how many boxes we have
185  ALLOCATE (force_env(1:nboxes))
186  ALLOCATE (bias_env(1:nboxes))
187  ALLOCATE (cell(1:nboxes))
188  ALLOCATE (particles_old(1:nboxes))
189  ALLOCATE (oldsys(1:nboxes))
190  ALLOCATE (averages(1:nboxes))
191  ALLOCATE (mc_par(1:nboxes))
192  ALLOCATE (pmvol_box(1:nboxes))
193  ALLOCATE (pmclus_box(1:nboxes))
194  ALLOCATE (pmhmc_box(1:nboxes))
195 
196  DO ibox = 1, nboxes
197  CALL get_mc_env(mc_env(ibox)%mc_env, &
198  mc_par=mc_par(ibox)%mc_par, &
199  force_env=force_env(ibox)%force_env)
200  END DO
201 
202  ! Gather units of measure for output (if available)
203  root_section => force_env(1)%force_env%root_section
204  CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%UNIT", &
205  c_val=unit_str)
206  unit_conv = cp_unit_from_cp2k(1.0_dp, trim(unit_str))
207  CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%PRINT_ATOM_KIND", &
208  l_val=print_kind)
209 
210  ! get some data out of mc_par
211  CALL get_mc_par(mc_par(1)%mc_par, &
212  ionode=ionode, source=source, group=group, &
213  data_file=data_file(1), moves_file=moves_file(1), &
214  cell_file=cell_file(1), coords_file=coords_file(1), &
215  energy_file=energy_file(1), displacement_file=displacement_file(1), &
216  lstop=lstop, nstep=nstep, nstart=nstart, pmvolume=pmvolume, pmhmc=pmhmc, &
217  molecules_file=molecules_file(1), pmswap=pmswap, nmoves=nmoves, &
218  pmtraion=pmtraion, pmtrans=pmtrans, pmcltrans=pmcltrans, iuptrans=iuptrans, &
219  iupvolume=iupvolume, ldiscrete=ldiscrete, pmtraion_mol=pmtraion_mol, &
220  lbias=lbias, iprint=iprint, pmavbmc_mol=pmavbmc_mol, &
221  discrete_step=discrete_step, fft_lib=fft_lib, avbmc_atom=avbmc_atom, &
222  pmavbmc=pmavbmc, pbias=pbias, mc_molecule_info=mc_molecule_info, &
223  pmrot_mol=pmrot_mol, pmtrans_mol=pmtrans_mol, pmvol_box=pmvol_box(1), &
224  pmclus_box=pmclus_box(1), ensemble=ensemble, pmhmc_box=pmhmc_box(1), lhmc=lhmc)
225 
226  ! get some data from the molecule types
227  CALL get_mc_molecule_info(mc_molecule_info, conf_prob=conf_prob, &
228  nchains=nchains, nmol_types=nmol_types, nunits_tot=nunits_tot, &
229  mol_type=mol_type, nchain_total=nchain_total, nunits=nunits, &
230  atom_names=atom_names, mass=mass)
231 
232  ! allocate some stuff based on the number of molecule types we have
233  ALLOCATE (moves(1:nmol_types, 1:nboxes))
234  ALLOCATE (move_updates(1:nmol_types, 1:nboxes))
235 
236  IF (nboxes .GT. 1) THEN
237  DO ibox = 2, nboxes
238  CALL get_mc_par(mc_par(ibox)%mc_par, &
239  data_file=data_file(ibox), &
240  moves_file=moves_file(ibox), &
241  cell_file=cell_file(ibox), coords_file=coords_file(ibox), &
242  energy_file=energy_file(ibox), &
243  displacement_file=displacement_file(ibox), &
244  molecules_file=molecules_file(ibox), pmvol_box=pmvol_box(ibox), &
245  pmclus_box=pmclus_box(ibox), pmhmc_box=pmhmc_box(ibox))
246  END DO
247  END IF
248 
249  ! this is a check we can't do in the input checking
250  IF (pmvol_box(nboxes) .LT. 1.0e0_dp) THEN
251  cpabort('The last value of PMVOL_BOX needs to be 1.0')
252  END IF
253  IF (pmclus_box(nboxes) .LT. 1.0e0_dp) THEN
254  cpabort('The last value of PMVOL_BOX needs to be 1.0')
255  END IF
256  IF (pmhmc_box(nboxes) .LT. 1.0e0_dp) THEN
257  cpabort('The last value of PMHMC_BOX needs to be 1.0')
258  END IF
259 
260  ! allocate the particle positions array for broadcasting
261  ALLOCATE (r_old(3, sum(nunits_tot), 1:nboxes))
262 
263  ! figure out what the default write unit is
264  iw = cp_logger_get_default_io_unit()
265 
266  IF (iw > 0) THEN
267  WRITE (iw, *)
268  WRITE (iw, *)
269  WRITE (iw, *) 'Beginning the Monte Carlo calculation.'
270  WRITE (iw, *)
271  WRITE (iw, *)
272  END IF
273 
274  ! initialize running average variables
275  energy_check(:) = 0.0e0_dp
276  box_flag(:) = 0
277  istep(:) = 0
278 
279  DO ibox = 1, nboxes
280  ! initialize the moves array, the arrays for updating maximum move
281  ! displacements, and the averages array
282  DO itype = 1, nmol_types
283  CALL init_mc_moves(moves(itype, ibox)%moves)
284  CALL init_mc_moves(move_updates(itype, ibox)%moves)
285  END DO
286  CALL mc_averages_create(averages(ibox)%averages)
287 
288  ! find the energy of the initial configuration
289  IF (sum(nchains(:, ibox)) .NE. 0) THEN
290  CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
291  calc_force=.false.)
292  CALL force_env_get(force_env(ibox)%force_env, &
293  potential_energy=old_energy(ibox))
294  ELSE
295  old_energy(ibox) = 0.0e0_dp
296  END IF
297  initial_energy(ibox) = old_energy(ibox)
298 
299 ! don't care about overlaps if we're only doing HMC
300 
301  IF (.NOT. lhmc) THEN
302  ! check for overlaps
303  start_mol = 1
304  DO jbox = 1, ibox - 1
305  start_mol = start_mol + sum(nchains(:, jbox))
306  END DO
307  end_mol = start_mol + sum(nchains(:, ibox)) - 1
308  CALL check_for_overlap(force_env(ibox)%force_env, nchains(:, ibox), &
309  nunits, loverlap, mol_type(start_mol:end_mol))
310  IF (loverlap) cpabort("overlap in an initial configuration")
311  END IF
312 
313  ! get the subsystems and the cell information
314  CALL force_env_get(force_env(ibox)%force_env, &
315  subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
316  CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
317  CALL cp_subsys_get(oldsys(ibox)%subsys, &
318  particles=particles_old(ibox)%list)
319  ! record the old coordinates, in case a move is rejected
320  DO iparticle = 1, nunits_tot(ibox)
321  r_old(1:3, iparticle, ibox) = &
322  particles_old(ibox)%list%els(iparticle)%r(1:3)
323  END DO
324 
325  ! find the bias energy of the initial run
326  IF (lbias) THEN
327  ! determine the atom names of every particle
328  ALLOCATE (atom_names_box(1:nunits_tot(ibox)))
329 
330  atom_number = 1
331  DO imolecule = 1, sum(nchains(:, ibox))
332  DO iunit = 1, nunits(mol_type(imolecule + start_mol - 1))
333  atom_names_box(atom_number) = &
334  atom_names(iunit, mol_type(imolecule + start_mol - 1))
335  atom_number = atom_number + 1
336  END DO
337  END DO
338 
339  CALL get_mc_par(mc_par(ibox)%mc_par, mc_bias_file=mc_bias_file)
340  nchains_box => nchains(:, ibox)
341  CALL mc_create_bias_force_env(bias_env(ibox)%force_env, &
342  r_old(:, :, ibox), atom_names_box(:), nunits_tot(ibox), &
343  para_env, abc(:, ibox), nchains_box, input_declaration, mc_bias_file, &
344  ionode)
345  IF (sum(nchains(:, ibox)) .NE. 0) THEN
346  CALL force_env_calc_energy_force(bias_env(ibox)%force_env, &
347  calc_force=.false.)
348  CALL force_env_get(bias_env(ibox)%force_env, &
349  potential_energy=last_bias_energy(ibox))
350 
351  ELSE
352  last_bias_energy(ibox) = 0.0e0_dp
353  END IF
354  bias_energy(ibox) = last_bias_energy(ibox)
355  DEALLOCATE (atom_names_box)
356  END IF
357  lnew_bias_env = .false.
358 
359  END DO
360 
361  ! back to seriel for a bunch of I/O stuff
362  IF (ionode) THEN
363 
364  ! record the combined energies,coordinates, and cell lengths
365  CALL open_file(file_name='mc_cell_length', &
366  unit_number=cell_unit, file_position='APPEND', &
367  file_action='WRITE', file_status='UNKNOWN')
368  CALL open_file(file_name='mc_energies', &
369  unit_number=com_ene, file_position='APPEND', &
370  file_action='WRITE', file_status='UNKNOWN')
371  CALL open_file(file_name='mc_coordinates', &
372  unit_number=com_crd, file_position='APPEND', &
373  file_action='WRITE', file_status='UNKNOWN')
374  CALL open_file(file_name='mc_molecules', &
375  unit_number=com_mol, file_position='APPEND', &
376  file_action='WRITE', file_status='UNKNOWN')
377  WRITE (com_ene, *) 'Initial Energies: ', &
378  old_energy(1:nboxes)
379  DO ibox = 1, nboxes
380  WRITE (com_mol, *) 'Initial Molecules: ', &
381  nchains(:, ibox)
382  END DO
383  DO ibox = 1, nboxes
384  WRITE (cell_unit, *) 'Initial: ', &
385  abc(1:3, ibox)*angstrom
386  WRITE (cbox, '(I4)') ibox
387  CALL open_file(file_name='energy_differences_box'// &
388  trim(adjustl(cbox)), &
389  unit_number=diff(ibox), file_position='APPEND', &
390  file_action='WRITE', file_status='UNKNOWN')
391  IF (sum(nchains(:, ibox)) == 0) THEN
392  WRITE (com_crd, *) ' 0'
393  WRITE (com_crd, *) 'INITIAL BOX '//trim(adjustl(cbox))
394  ELSE
395  CALL write_particle_coordinates(particles_old(ibox)%list%els, &
396  com_crd, dump_xmol, 'POS', 'INITIAL BOX '//trim(adjustl(cbox)), &
397  unit_conv=unit_conv, print_kind=print_kind)
398  END IF
399  CALL open_file(file_name=data_file(ibox), &
400  unit_number=data_unit(ibox), file_position='APPEND', &
401  file_action='WRITE', file_status='UNKNOWN')
402  CALL open_file(file_name=moves_file(ibox), &
403  unit_number=move_unit(ibox), file_position='APPEND', &
404  file_action='WRITE', file_status='UNKNOWN')
405  CALL open_file(file_name=displacement_file(ibox), &
406  unit_number=rm(ibox), file_position='APPEND', &
407  file_action='WRITE', file_status='UNKNOWN')
408  CALL open_file(file_name=cell_file(ibox), &
409  unit_number=cl(ibox), file_position='APPEND', &
410  file_action='WRITE', file_status='UNKNOWN')
411 
412  END DO
413 
414  ! back to parallel mode
415  END IF
416 
417  DO ibox = 1, nboxes
418  CALL group%bcast(cl(ibox), source)
419  CALL group%bcast(rm(ibox), source)
420  CALL group%bcast(diff(ibox), source)
421  ! set all the units numbers that we just opened in the respective mc_par
422  CALL set_mc_par(mc_par(ibox)%mc_par, cl=cl(ibox), rm=rm(ibox), &
423  diff=diff(ibox))
424  END DO
425 
426  ! if we're doing a discrete volume move, we need to set up the array
427  ! that keeps track of which direction we can move in
428  IF (ldiscrete) THEN
429  IF (nboxes .NE. 1) &
430  cpabort('ldiscrete=.true. ONLY for systems with 1 box')
431  CALL create_discrete_array(abc(:, 1), discrete_array(:, :), &
432  discrete_step)
433  END IF
434 
435  ! find out how many steps we're doing...change the updates to be in cycles
436  ! if the total number of steps is measured in cycles
437  IF (.NOT. lstop) THEN
438  nstep = nstep*nchain_total
439  iuptrans = iuptrans*nchain_total
440  iupvolume = iupvolume*nchain_total
441  END IF
442 
443  DO nnstep = nstart + 1, nstart + nstep
444 
445  IF (mod(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
446  WRITE (iw, *)
447  WRITE (iw, *) "------- On Monte Carlo Step ", nnstep
448  END IF
449 
450  IF (ionode) rand = rng_stream%next()
451  ! broadcast the random number, to make sure we're on the same move
452  CALL group%bcast(rand, source)
453 
454  IF (rand .LT. pmvolume) THEN
455 
456  IF (mod(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
457  WRITE (iw, *) "Attempting a volume move"
458  WRITE (iw, *)
459  END IF
460 
461  SELECT CASE (ensemble)
462  CASE ("TRADITIONAL")
463  CALL mc_volume_move(mc_par(1)%mc_par, &
464  force_env(1)%force_env, &
465  moves(1, 1)%moves, move_updates(1, 1)%moves, &
466  old_energy(1), 1, &
467  energy_check(1), r_old(:, :, 1), iw, discrete_array(:, :), &
468  rng_stream)
469  CASE ("GEMC_NVT")
470  CALL mc_ge_volume_move(mc_par, force_env, moves, &
471  move_updates, nnstep, old_energy, energy_check, &
472  r_old, rng_stream)
473  CASE ("GEMC_NPT")
474  ! we need to select a box based on the probability given in the input file
475  IF (ionode) rand = rng_stream%next()
476  CALL group%bcast(rand, source)
477 
478  DO ibox = 1, nboxes
479  IF (rand .LE. pmvol_box(ibox)) THEN
480  box_number = ibox
481  EXIT
482  END IF
483  END DO
484 
485  CALL mc_volume_move(mc_par(box_number)%mc_par, &
486  force_env(box_number)%force_env, &
487  moves(1, box_number)%moves, &
488  move_updates(1, box_number)%moves, &
489  old_energy(box_number), box_number, &
490  energy_check(box_number), r_old(:, :, box_number), iw, &
491  discrete_array(:, :), &
492  rng_stream)
493  END SELECT
494 
495 ! update all the pointers here, because otherwise we may pass wrong information when we're making a bias environment
496  DO ibox = 1, nboxes
497  CALL force_env_get(force_env(ibox)%force_env, &
498  subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
499  CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
500  CALL cp_subsys_get(oldsys(ibox)%subsys, &
501  particles=particles_old(ibox)%list)
502  END DO
503 
504  ! we need a new biasing environment now, if we're into that sort of thing
505  IF (lbias) THEN
506  DO ibox = 1, nboxes
507  CALL force_env_release(bias_env(ibox)%force_env)
508  ! determine the atom names of every particle
509  ALLOCATE (atom_names_box(1:nunits_tot(ibox)))
510  start_mol = 1
511  DO jbox = 1, ibox - 1
512  start_mol = start_mol + sum(nchains(:, jbox))
513  END DO
514  end_mol = start_mol + sum(nchains(:, ibox)) - 1
515  atom_number = 1
516  DO imolecule = 1, sum(nchains(:, ibox))
517  DO iunit = 1, nunits(mol_type(imolecule + start_mol - 1))
518  atom_names_box(atom_number) = &
519  atom_names(iunit, mol_type(imolecule + start_mol - 1))
520  atom_number = atom_number + 1
521  END DO
522  END DO
523 
524 ! need to find out what the cell lengths are
525  CALL force_env_get(force_env(ibox)%force_env, &
526  subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
527  CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
528 
529  CALL get_mc_par(mc_par(ibox)%mc_par, &
530  mc_bias_file=mc_bias_file)
531  nchains_box => nchains(:, ibox)
532 
533  CALL mc_create_bias_force_env(bias_env(ibox)%force_env, &
534  r_old(:, :, ibox), atom_names_box(:), nunits_tot(ibox), &
535  para_env, abc(:, ibox), nchains_box, input_declaration, &
536  mc_bias_file, ionode)
537 
538  IF (sum(nchains(:, ibox)) .NE. 0) THEN
539  CALL force_env_calc_energy_force( &
540  bias_env(ibox)%force_env, &
541  calc_force=.false.)
542  CALL force_env_get(bias_env(ibox)%force_env, &
543  potential_energy=last_bias_energy(ibox))
544  ELSE
545  last_bias_energy(ibox) = 0.0e0_dp
546  END IF
547  bias_energy(ibox) = last_bias_energy(ibox)
548  DEALLOCATE (atom_names_box)
549  END DO
550  END IF
551 
552  ELSEIF (rand .LT. pmswap) THEN
553 
554  ! try a swap move
555  IF (mod(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
556  WRITE (iw, *) "Attempting a swap move"
557  WRITE (iw, *)
558  END IF
559 
560  CALL mc_ge_swap_move(mc_par, force_env, bias_env, moves, &
561  energy_check(:), r_old(:, :, :), old_energy(:), input_declaration, &
562  para_env, bias_energy(:), last_bias_energy(:), rng_stream)
563 
564  ! the number of molecules may have changed, which deallocated the whole
565  ! mc_molecule_info structure
566  CALL get_mc_par(mc_par(1)%mc_par, mc_molecule_info=mc_molecule_info)
567  CALL get_mc_molecule_info(mc_molecule_info, conf_prob=conf_prob, &
568  nchains=nchains, nmol_types=nmol_types, nunits_tot=nunits_tot, &
569  mol_type=mol_type, nchain_total=nchain_total, nunits=nunits, &
570  atom_names=atom_names, mass=mass)
571 
572  ELSEIF (rand .LT. pmhmc) THEN
573 ! try hybrid Monte Carlo
574  IF (mod(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
575  WRITE (iw, *) "Attempting a hybrid Monte Carlo move"
576  WRITE (iw, *)
577  END IF
578 
579 ! pick a box at random
580  IF (ionode) rand = rng_stream%next()
581  CALL group%bcast(rand, source)
582 
583  DO ibox = 1, nboxes
584  IF (rand .LE. pmhmc_box(ibox)) THEN
585  box_number = ibox
586  EXIT
587  END IF
588  END DO
589 
590  CALL mc_hmc_move(mc_par(box_number)%mc_par, &
591  force_env(box_number)%force_env, globenv, &
592  moves(1, box_number)%moves, &
593  move_updates(1, box_number)%moves, &
594  old_energy(box_number), box_number, &
595  energy_check(box_number), r_old(:, :, box_number), &
596  rng_stream)
597 
598  ELSEIF (rand .LT. pmavbmc) THEN
599  ! try an AVBMC move
600  IF (mod(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
601  WRITE (iw, *) "Attempting an AVBMC1 move"
602  WRITE (iw, *)
603  END IF
604 
605  ! first, pick a box to do it for
606  IF (ionode) rand = rng_stream%next()
607  CALL group%bcast(rand, source)
608 
609  IF (nboxes .EQ. 2) THEN
610  IF (rand .LT. 0.1e0_dp) THEN
611  box_number = 1
612  ELSE
613  box_number = 2
614  END IF
615  ELSE
616  box_number = 1
617  END IF
618 
619  ! now pick a molecule type to do it for
620  IF (ionode) rand = rng_stream%next()
621  CALL group%bcast(rand, source)
622  molecule_type_swap = 0
623  DO imol_type = 1, nmol_types
624  IF (rand .LT. pmavbmc_mol(imol_type)) THEN
625  molecule_type_swap = imol_type
626  EXIT
627  END IF
628  END DO
629  IF (molecule_type_swap == 0) &
630  cpabort('Did not choose a molecule type to swap...check AVBMC input')
631 
632  ! now pick a molecule, automatically rejecting the move if the
633  ! box is empty or only has one molecule
634  IF (sum(nchains(:, box_number)) .LE. 1) THEN
635  ! indicate that we tried a move
636  moves(molecule_type_swap, box_number)%moves%empty_avbmc = &
637  moves(molecule_type_swap, box_number)%moves%empty_avbmc + 1
638  ELSE
639 
640  ! pick a molecule to be swapped in the box
641  IF (ionode) THEN
642  CALL find_mc_test_molecule(mc_molecule_info, &
643  start_atom_swap, idum, jdum, rng_stream, &
644  box=box_number, molecule_type_old=molecule_type_swap)
645 
646  ! pick a molecule to act as the target in the box...we don't care what type
647  DO
648  CALL find_mc_test_molecule(mc_molecule_info, &
649  start_atom_target, idum, molecule_type_target, &
650  rng_stream, box=box_number)
651  IF (start_atom_swap .NE. start_atom_target) THEN
652  start_atom_target = start_atom_target + &
653  avbmc_atom(molecule_type_target) - 1
654  EXIT
655  END IF
656  END DO
657 
658  ! choose if we're swapping into the bonded region of mol_target, or
659  ! into the nonbonded region
660  rand = rng_stream%next()
661 
662  END IF
663  CALL group%bcast(start_atom_swap, source)
664  CALL group%bcast(box_number, source)
665  CALL group%bcast(start_atom_target, source)
666  CALL group%bcast(rand, source)
667 
668  IF (rand .LT. pbias(molecule_type_swap)) THEN
669  move_type_avbmc = 'in'
670  ELSE
671  move_type_avbmc = 'out'
672  END IF
673 
674  CALL mc_avbmc_move(mc_par(box_number)%mc_par, &
675  force_env(box_number)%force_env, &
676  bias_env(box_number)%force_env, &
677  moves(molecule_type_swap, box_number)%moves, &
678  energy_check(box_number), &
679  r_old(:, :, box_number), old_energy(box_number), &
680  start_atom_swap, start_atom_target, molecule_type_swap, &
681  box_number, bias_energy(box_number), &
682  last_bias_energy(box_number), &
683  move_type_avbmc, rng_stream)
684 
685  END IF
686 
687  ELSE
688 
689  IF (mod(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
690  WRITE (iw, *) "Attempting an inner move"
691  WRITE (iw, *)
692  END IF
693 
694  DO imove = 1, nmoves
695 
696  IF (ionode) rand = rng_stream%next()
697  CALL group%bcast(rand, source)
698  IF (rand .LT. pmtraion) THEN
699  ! change molecular conformation
700  ! first, pick a box to do it for
701  IF (ionode) rand = rng_stream%next()
702  CALL group%bcast(rand, source)
703  IF (nboxes .EQ. 2) THEN
704  IF (rand .LT. 0.75e0_dp) THEN
705  box_number = 1
706  ELSE
707  box_number = 2
708  END IF
709  ELSE
710  box_number = 1
711  END IF
712 
713  ! figure out which molecule type we're looking for
714  IF (ionode) rand = rng_stream%next()
715  CALL group%bcast(rand, source)
716  molecule_type = 0
717  DO imol_type = 1, nmol_types
718  IF (rand .LT. pmtraion_mol(imol_type)) THEN
719  molecule_type = imol_type
720  EXIT
721  END IF
722  END DO
723  IF (molecule_type == 0) CALL cp_abort( &
724  __location__, &
725  'Did not choose a molecule type to conf change...PMTRAION_MOL should not be all 0.0')
726 
727  ! now pick a molecule, automatically rejecting the move if the
728  ! box is empty
729  IF (nchains(molecule_type, box_number) == 0) THEN
730  ! indicate that we tried a move
731  moves(molecule_type, box_number)%moves%empty_conf = &
732  moves(molecule_type, box_number)%moves%empty_conf + 1
733  ELSE
734  ! pick a molecule in the box
735  IF (ionode) THEN
736  CALL find_mc_test_molecule(mc_molecule_info, &
737  start_atom, idum, &
738  jdum, rng_stream, &
739  box=box_number, molecule_type_old=molecule_type)
740 
741  ! choose if we're changing a bond length or an angle
742  rand = rng_stream%next()
743  END IF
744  CALL group%bcast(rand, source)
745  CALL group%bcast(start_atom, source)
746  CALL group%bcast(box_number, source)
747  CALL group%bcast(molecule_type, source)
748 
749  ! figure out what kind of move we're doing
750  IF (rand .LT. conf_prob(1, molecule_type)) THEN
751  move_type = 'bond'
752  ELSEIF (rand .LT. (conf_prob(1, molecule_type) + &
753  conf_prob(2, molecule_type))) THEN
754  move_type = 'angle'
755  ELSE
756  move_type = 'dihedral'
757  END IF
758  box_flag(box_number) = 1
759  CALL mc_conformation_change(mc_par(box_number)%mc_par, &
760  force_env(box_number)%force_env, &
761  bias_env(box_number)%force_env, &
762  moves(molecule_type, box_number)%moves, &
763  move_updates(molecule_type, box_number)%moves, &
764  start_atom, molecule_type, box_number, &
765  bias_energy(box_number), &
766  move_type, lreject, rng_stream)
767  IF (lreject) EXIT
768  END IF
769  ELSEIF (rand .LT. pmtrans) THEN
770  ! translate a whole molecule in the system
771  ! pick a molecule type
772  IF (ionode) rand = rng_stream%next()
773  CALL group%bcast(rand, source)
774  molecule_type = 0
775  DO imol_type = 1, nmol_types
776  IF (rand .LT. pmtrans_mol(imol_type)) THEN
777  molecule_type = imol_type
778  EXIT
779  END IF
780  END DO
781  IF (molecule_type == 0) CALL cp_abort( &
782  __location__, &
783  'Did not choose a molecule type to translate...PMTRANS_MOL should not be all 0.0')
784 
785  ! now pick a molecule of that type
786  IF (ionode) &
787  CALL find_mc_test_molecule(mc_molecule_info, &
788  start_atom, box_number, idum, rng_stream, &
789  molecule_type_old=molecule_type)
790  CALL group%bcast(start_atom, source)
791  CALL group%bcast(box_number, source)
792  box_flag(box_number) = 1
793  CALL mc_molecule_translation(mc_par(box_number)%mc_par, &
794  force_env(box_number)%force_env, &
795  bias_env(box_number)%force_env, &
796  moves(molecule_type, box_number)%moves, &
797  move_updates(molecule_type, box_number)%moves, &
798  start_atom, box_number, bias_energy(box_number), &
799  molecule_type, lreject, rng_stream)
800  IF (lreject) EXIT
801  ELSEIF (rand .LT. pmcltrans) THEN
802  ! translate a whole cluster in the system
803  ! first, pick a box to do it for
804  IF (ionode) rand = rng_stream%next()
805  CALL group%bcast(rand, source)
806 
807  DO ibox = 1, nboxes
808  IF (rand .LE. pmclus_box(ibox)) THEN
809  box_number = ibox
810  EXIT
811  END IF
812  END DO
813  box_flag(box_number) = 1
814  CALL mc_cluster_translation(mc_par(box_number)%mc_par, &
815  force_env(box_number)%force_env, &
816  bias_env(box_number)%force_env, &
817  moves(1, box_number)%moves, &
818  move_updates(1, box_number)%moves, &
819  box_number, bias_energy(box_number), &
820  lreject, rng_stream)
821  IF (lreject) EXIT
822  ELSE
823  ! rotate a whole molecule in the system
824  ! pick a molecule type
825  IF (ionode) rand = rng_stream%next()
826  CALL group%bcast(rand, source)
827  molecule_type = 0
828  DO imol_type = 1, nmol_types
829  IF (rand .LT. pmrot_mol(imol_type)) THEN
830  molecule_type = imol_type
831  EXIT
832  END IF
833  END DO
834  IF (molecule_type == 0) CALL cp_abort( &
835  __location__, &
836  'Did not choose a molecule type to rotate...PMROT_MOL should not be all 0.0')
837 
838  IF (ionode) &
839  CALL find_mc_test_molecule(mc_molecule_info, &
840  start_atom, box_number, idum, rng_stream, &
841  molecule_type_old=molecule_type)
842  CALL group%bcast(start_atom, source)
843  CALL group%bcast(box_number, source)
844  box_flag(box_number) = 1
845  CALL mc_molecule_rotation(mc_par(box_number)%mc_par, &
846  force_env(box_number)%force_env, &
847  bias_env(box_number)%force_env, &
848  moves(molecule_type, box_number)%moves, &
849  move_updates(molecule_type, box_number)%moves, &
850  box_number, start_atom, &
851  molecule_type, bias_energy(box_number), &
852  lreject, rng_stream)
853  IF (lreject) EXIT
854  END IF
855 
856  END DO
857 
858  ! now do a Quickstep calculation to see if we accept the sequence
859  CALL mc_quickstep_move(mc_par, force_env, bias_env, &
860  moves, lreject, move_updates, energy_check(:), r_old(:, :, :), &
861  nnstep, old_energy(:), bias_energy(:), last_bias_energy(:), &
862  nboxes, box_flag(:), oldsys, particles_old, &
863  rng_stream, unit_conv)
864 
865  END IF
866 
867  ! make sure the pointers are pointing correctly since the subsys may
868  ! have changed
869  DO ibox = 1, nboxes
870  CALL force_env_get(force_env(ibox)%force_env, &
871  subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
872  CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
873  CALL cp_subsys_get(oldsys(ibox)%subsys, &
874  particles=particles_old(ibox)%list)
875  END DO
876 
877  IF (ionode) THEN
878 
879  IF (mod(nnstep, iprint) == 0) THEN
880  WRITE (com_ene, *) nnstep, old_energy(1:nboxes)
881 
882  DO ibox = 1, nboxes
883 
884  ! write the molecule information
885  WRITE (com_mol, *) nnstep, nchains(:, ibox)
886 
887  ! write the move statistics to file
888  DO itype = 1, nmol_types
889  CALL write_move_stats(moves(itype, ibox)%moves, &
890  nnstep, move_unit(ibox))
891  END DO
892 
893  ! write a restart file
894  CALL write_mc_restart(nnstep, mc_par(ibox)%mc_par, &
895  nchains(:, ibox), force_env(ibox)%force_env)
896 
897  ! write cell lengths
898  WRITE (cell_unit, *) nnstep, abc(1:3, ibox)*angstrom
899 
900  ! write particle coordinates
901  WRITE (cbox, '(I4)') ibox
902  WRITE (cstep, '(I8)') nnstep
903  IF (sum(nchains(:, ibox)) == 0) THEN
904  WRITE (com_crd, *) ' 0'
905  WRITE (com_crd, *) 'BOX '//trim(adjustl(cbox))// &
906  ', STEP '//trim(adjustl(cstep))
907  ELSE
908  CALL write_particle_coordinates( &
909  particles_old(ibox)%list%els, &
910  com_crd, dump_xmol, 'POS', &
911  'BOX '//trim(adjustl(cbox))// &
912  ', STEP '//trim(adjustl(cstep)), &
913  unit_conv=unit_conv)
914  END IF
915  END DO
916  END IF ! end the things we only do every iprint moves
917 
918  DO ibox = 1, nboxes
919  ! compute some averages
920  averages(ibox)%averages%ave_energy = &
921  averages(ibox)%averages%ave_energy*real(nnstep - &
922  nstart - 1, dp)/real(nnstep - nstart, dp) + &
923  old_energy(ibox)/real(nnstep - nstart, dp)
924  averages(ibox)%averages%molecules = &
925  averages(ibox)%averages%molecules*real(nnstep - &
926  nstart - 1, dp)/real(nnstep - nstart, dp) + &
927  REAL(sum(nchains(:, ibox)), dp)/REAL(nnstep - nstart, dp)
928  averages(ibox)%averages%ave_volume = &
929  averages(ibox)%averages%ave_volume* &
930  REAL(nnstep - nstart - 1, dp)/REAL(nnstep - nstart, dp) + &
931  abc(1, ibox)*abc(2, ibox)*abc(3, ibox)/ &
932  REAL(nnstep - nstart, dp)
933 
934  ! flush the buffers to the files
935  CALL m_flush(data_unit(ibox))
936  CALL m_flush(diff(ibox))
937  CALL m_flush(move_unit(ibox))
938  CALL m_flush(cl(ibox))
939  CALL m_flush(rm(ibox))
940 
941  END DO
942 
943  ! flush more buffers to the files
944  CALL m_flush(cell_unit)
945  CALL m_flush(com_ene)
946  CALL m_flush(com_crd)
947  CALL m_flush(com_mol)
948 
949  END IF
950 
951  ! reset the box flags
952  box_flag(:) = 0
953 
954  ! check to see if EXIT file exists...if so, end the calculation
955  CALL external_control(should_stop, "MC", globenv=globenv)
956  IF (should_stop) EXIT
957 
958  ! update the move displacements, if necessary
959  DO ibox = 1, nboxes
960  IF (mod(nnstep - nstart, iuptrans) == 0) THEN
961  DO itype = 1, nmol_types
962  CALL mc_move_update(mc_par(ibox)%mc_par, &
963  move_updates(itype, ibox)%moves, itype, &
964  "trans", nnstep, ionode)
965  END DO
966  END IF
967 
968  IF (mod(nnstep - nstart, iupvolume) == 0) THEN
969  CALL mc_move_update(mc_par(ibox)%mc_par, &
970  move_updates(1, ibox)%moves, 1337, &
971  "volume", nnstep, ionode)
972  END IF
973  END DO
974 
975  ! check to see if there are any overlaps in the boxes, and fold coordinates
976 ! don't care about overlaps if we're only doing HMC
977  IF (.NOT. lhmc) THEN
978  DO ibox = 1, nboxes
979  IF (sum(nchains(:, ibox)) .NE. 0) THEN
980  start_mol = 1
981  DO jbox = 1, ibox - 1
982  start_mol = start_mol + sum(nchains(:, jbox))
983  END DO
984  end_mol = start_mol + sum(nchains(:, ibox)) - 1
985  CALL check_for_overlap(force_env(ibox)%force_env, &
986  nchains(:, ibox), nunits, loverlap, &
987  mol_type(start_mol:end_mol))
988  IF (loverlap) THEN
989  IF (iw > 0) WRITE (iw, *) nnstep
990  cpabort('coordinate overlap at the end of the above step')
991  ! now fold the coordinates...don't do this anywhere but here, because
992  ! we can get screwed up with the mc_molecule_info stuff (like in swap move)...
993  ! this is kind of ugly, with allocated and deallocating every time
994  ALLOCATE (r_temp(1:3, 1:nunits_tot(ibox)))
995 
996  DO iunit = 1, nunits_tot(ibox)
997  r_temp(1:3, iunit) = &
998  particles_old(ibox)%list%els(iunit)%r(1:3)
999  END DO
1000 
1001  CALL mc_coordinate_fold(r_temp(:, :), &
1002  sum(nchains(:, ibox)), mol_type(start_mol:end_mol), &
1003  mass, nunits, abc(1:3, ibox))
1004 
1005  ! save the folded coordinates
1006  DO iunit = 1, nunits_tot(ibox)
1007  r_old(1:3, iunit, ibox) = r_temp(1:3, iunit)
1008  particles_old(ibox)%list%els(iunit)%r(1:3) = &
1009  r_temp(1:3, iunit)
1010  END DO
1011 
1012  ! if we're biasing, we need to do the same
1013  IF (lbias) THEN
1014  CALL force_env_get(bias_env(ibox)%force_env, &
1015  subsys=biassys)
1016  CALL cp_subsys_get(biassys, &
1017  particles=particles_bias)
1018 
1019  DO iunit = 1, nunits_tot(ibox)
1020  particles_bias%els(iunit)%r(1:3) = &
1021  r_temp(1:3, iunit)
1022  END DO
1023  END IF
1024 
1025  DEALLOCATE (r_temp)
1026  END IF
1027  END IF
1028  END DO
1029  END IF
1030 
1031  !debug code
1032  IF (debug_this_module) THEN
1033  DO ibox = 1, nboxes
1034  IF (sum(nchains(:, ibox)) .NE. 0) THEN
1035  CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
1036  calc_force=.false.)
1037  CALL force_env_get(force_env(ibox)%force_env, &
1038  potential_energy=test_energy)
1039  ELSE
1040  test_energy = 0.0e0_dp
1041  END IF
1042 
1043  IF (abs(initial_energy(ibox) + energy_check(ibox) - &
1044  test_energy) .GT. 0.0000001e0_dp) THEN
1045  IF (iw > 0) THEN
1046  WRITE (iw, *) '!!!!!!! We have an energy problem. !!!!!!!!'
1047  WRITE (iw, '(A,T64,F16.10)') 'Final Energy = ', test_energy
1048  WRITE (iw, '(A,T64,F16.10)') 'Initial Energy+energy_check=', &
1049  initial_energy(ibox) + energy_check(ibox)
1050  WRITE (iw, *) 'Box ', ibox
1051  WRITE (iw, *) 'nchains ', nchains(:, ibox)
1052  END IF
1053  cpabort('!!!!!!! We have an energy problem. !!!!!!!!')
1054  END IF
1055  END DO
1056  END IF
1057  END DO
1058 
1059  ! write a restart file
1060  IF (ionode) THEN
1061  DO ibox = 1, nboxes
1062  CALL write_mc_restart(nnstep, mc_par(ibox)%mc_par, &
1063  nchains(:, ibox), force_env(ibox)%force_env)
1064  END DO
1065  END IF
1066 
1067  ! calculate the final energy
1068  DO ibox = 1, nboxes
1069  IF (sum(nchains(:, ibox)) .NE. 0) THEN
1070  CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
1071  calc_force=.false.)
1072  CALL force_env_get(force_env(ibox)%force_env, &
1073  potential_energy=final_energy(ibox))
1074  ELSE
1075  final_energy(ibox) = 0.0e0_dp
1076  END IF
1077  IF (lbias) THEN
1078  CALL force_env_release(bias_env(ibox)%force_env)
1079  END IF
1080  END DO
1081 
1082  ! do some stuff in serial
1083  IF (ionode .OR. (iw > 0)) THEN
1084 
1085  WRITE (com_ene, *) 'Final Energies: ', &
1086  final_energy(1:nboxes)
1087 
1088  DO ibox = 1, nboxes
1089  WRITE (cbox, '(I4)') ibox
1090  IF (sum(nchains(:, ibox)) == 0) THEN
1091  WRITE (com_crd, *) ' 0'
1092  WRITE (com_crd, *) 'BOX '//trim(adjustl(cbox))
1093  ELSE
1094  CALL write_particle_coordinates( &
1095  particles_old(ibox)%list%els, &
1096  com_crd, dump_xmol, 'POS', &
1097  'FINAL BOX '//trim(adjustl(cbox)), unit_conv=unit_conv)
1098  END IF
1099 
1100  ! write a bunch of data to the screen
1101  WRITE (iw, '(A)') &
1102  '------------------------------------------------'
1103  WRITE (iw, '(A,I1,A)') &
1104  '| BOX ', ibox, &
1105  ' |'
1106  WRITE (iw, '(A)') &
1107  '------------------------------------------------'
1108  test_moves => moves(:, ibox)
1109  CALL final_mc_write(mc_par(ibox)%mc_par, test_moves, &
1110  iw, energy_check(ibox), &
1111  initial_energy(ibox), final_energy(ibox), &
1112  averages(ibox)%averages)
1113 
1114  ! close any open files
1115  CALL close_file(unit_number=diff(ibox))
1116  CALL close_file(unit_number=data_unit(ibox))
1117  CALL close_file(unit_number=move_unit(ibox))
1118  CALL close_file(unit_number=cl(ibox))
1119  CALL close_file(unit_number=rm(ibox))
1120  END DO
1121 
1122  ! close some more files
1123  CALL close_file(unit_number=cell_unit)
1124  CALL close_file(unit_number=com_ene)
1125  CALL close_file(unit_number=com_crd)
1126  CALL close_file(unit_number=com_mol)
1127  END IF
1128 
1129  DO ibox = 1, nboxes
1130  CALL set_mc_env(mc_env(ibox)%mc_env, &
1131  mc_par=mc_par(ibox)%mc_par, &
1132  force_env=force_env(ibox)%force_env)
1133 
1134  ! deallocate some stuff
1135  DO itype = 1, nmol_types
1136  CALL mc_moves_release(move_updates(itype, ibox)%moves)
1137  CALL mc_moves_release(moves(itype, ibox)%moves)
1138  END DO
1139  CALL mc_averages_release(averages(ibox)%averages)
1140  END DO
1141 
1142  DEALLOCATE (pmhmc_box)
1143  DEALLOCATE (pmvol_box)
1144  DEALLOCATE (pmclus_box)
1145  DEALLOCATE (r_old)
1146  DEALLOCATE (force_env)
1147  DEALLOCATE (bias_env)
1148  DEALLOCATE (cell)
1149  DEALLOCATE (particles_old)
1150  DEALLOCATE (oldsys)
1151  DEALLOCATE (averages)
1152  DEALLOCATE (moves)
1153  DEALLOCATE (move_updates)
1154  DEALLOCATE (mc_par)
1155 
1156  ! end the timing
1157  CALL timestop(handle)
1158 
1159  END SUBROUTINE mc_run_ensemble
1160 
1161 ! **************************************************************************************************
1162 !> \brief Computes the second virial coefficient of a molecule by using the integral form
1163 !> of the second virial coefficient found in McQuarrie "Statistical Thermodynamics",
1164 !> B2(T) = -2Pi Int 0toInf [ Exp[-beta*u(r)] -1] r^2 dr Eq. 15-25
1165 !> I use trapazoidal integration with various step sizes
1166 !> (the integral is broken up into three parts, currently, but that's easily
1167 !> changed by the first variables found below). It generates nvirial configurations,
1168 !> doing the integration for each one, and then averages all the B2(T) to produce
1169 !> the final answer.
1170 !> \param mc_env a pointer that contains all mc_env for all the simulation
1171 !> boxes
1172 !> \param rng_stream the stream we pull random numbers from
1173 !>
1174 !> Suitable for parallel.
1175 !> \author MJM
1176 ! **************************************************************************************************
1177  SUBROUTINE mc_compute_virial(mc_env, rng_stream)
1178 
1179  TYPE(mc_environment_p_type), DIMENSION(:), POINTER :: mc_env
1180  TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
1181 
1182  INTEGER :: current_division, end_atom, ibin, idivision, iparticle, iprint, itemp, iunit, &
1183  ivirial, iw, nbins, nchain_total, nintegral_divisions, nmol_types, nvirial, &
1184  nvirial_temps, source, start_atom
1185  INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot
1186  INTEGER, DIMENSION(:, :), POINTER :: nchains
1187  LOGICAL :: ionode, loverlap
1188  REAL(dp), DIMENSION(:), POINTER :: beta, virial_cutoffs, virial_stepsize, &
1189  virial_temps
1190  REAL(dp), DIMENSION(:, :), POINTER :: mass, mayer, r_old
1191  REAL(kind=dp) :: ave_virial, current_value, distance, exp_max_val, exp_min_val, exponent, &
1192  integral, previous_value, square_value, trial_energy, triangle_value
1193  REAL(kind=dp), DIMENSION(1:3) :: abc, center_of_mass
1194  TYPE(cell_p_type), DIMENSION(:), POINTER :: cell
1195  TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsys
1196  TYPE(force_env_p_type), DIMENSION(:), POINTER :: force_env
1197  TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
1198  TYPE(mc_simulation_parameters_p_type), &
1199  DIMENSION(:), POINTER :: mc_par
1200  TYPE(mp_comm_type) :: group
1201  TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
1202 
1203 ! these are current magic numbers for how we compute the virial...
1204 ! we break it up into three parts to integrate the function so provide
1205 ! better statistics
1206 
1207  nintegral_divisions = 3
1208  ALLOCATE (virial_cutoffs(1:nintegral_divisions))
1209  ALLOCATE (virial_stepsize(1:nintegral_divisions))
1210  virial_cutoffs(1) = 8.0 ! first distance, in bohr
1211  virial_cutoffs(2) = 13.0 ! second distance, in bohr
1212  virial_cutoffs(3) = 22.0 ! maximum distance, in bohr
1213  virial_stepsize(1) = 0.04 ! stepsize from 0 to virial_cutoffs(1)
1214  virial_stepsize(2) = 0.1
1215  virial_stepsize(3) = 0.2
1216 
1217  nbins = ceiling(virial_cutoffs(1)/virial_stepsize(1) + (virial_cutoffs(2) - virial_cutoffs(1))/ &
1218  virial_stepsize(2) + (virial_cutoffs(3) - virial_cutoffs(2))/virial_stepsize(3))
1219 
1220  ! figure out what the default write unit is
1221  iw = cp_logger_get_default_io_unit()
1222 
1223  ! allocate a whole bunch of stuff based on how many boxes we have
1224  ALLOCATE (force_env(1:1))
1225  ALLOCATE (cell(1:1))
1226  ALLOCATE (particles(1:1))
1227  ALLOCATE (subsys(1:1))
1228  ALLOCATE (mc_par(1:1))
1229 
1230  CALL get_mc_env(mc_env(1)%mc_env, &
1231  mc_par=mc_par(1)%mc_par, &
1232  force_env=force_env(1)%force_env)
1233 
1234  ! get some data out of mc_par
1235  CALL get_mc_par(mc_par(1)%mc_par, &
1236  exp_max_val=exp_max_val, &
1237  exp_min_val=exp_min_val, nvirial=nvirial, &
1238  ionode=ionode, source=source, group=group, &
1239  mc_molecule_info=mc_molecule_info, virial_temps=virial_temps)
1240 
1241  IF (iw > 0) THEN
1242  WRITE (iw, *)
1243  WRITE (iw, *)
1244  WRITE (iw, *) 'Beginning the calculation of the second virial coefficient'
1245  WRITE (iw, *)
1246  WRITE (iw, *)
1247  END IF
1248 
1249  ! get some data from the molecule types
1250  CALL get_mc_molecule_info(mc_molecule_info, &
1251  nchains=nchains, nmol_types=nmol_types, nunits_tot=nunits_tot, &
1252  mol_type=mol_type, nchain_total=nchain_total, nunits=nunits, &
1253  mass=mass)
1254 
1255  nvirial_temps = SIZE(virial_temps)
1256  ALLOCATE (beta(1:nvirial_temps))
1257 
1258  DO itemp = 1, nvirial_temps
1259  beta(itemp) = 1/virial_temps(itemp)/boltzmann*joule
1260  END DO
1261 
1262  ! get the subsystems and the cell information
1263  CALL force_env_get(force_env(1)%force_env, &
1264  subsys=subsys(1)%subsys, cell=cell(1)%cell)
1265  CALL get_cell(cell(1)%cell, abc=abc(:))
1266  CALL cp_subsys_get(subsys(1)%subsys, &
1267  particles=particles(1)%list)
1268 
1269  ! check and make sure the box is big enough
1270  IF (abc(1) .NE. abc(2) .OR. abc(2) .NE. abc(3)) THEN
1271  cpabort('The box needs to be cubic for a virial calculation (it is easiest).')
1272  END IF
1273  IF (virial_cutoffs(nintegral_divisions) .GT. abc(1)/2.0e0_dp) THEN
1274  IF (iw > 0) &
1275  WRITE (iw, *) "Box length ", abc(1)*angstrom, " virial cutoff ", &
1276  virial_cutoffs(nintegral_divisions)*angstrom
1277  cpabort('You need a bigger box to deal with this virial cutoff (see above).')
1278  END IF
1279 
1280  ! store the coordinates of the molecules in an array so we can work with it
1281  ALLOCATE (r_old(1:3, 1:nunits_tot(1)))
1282 
1283  DO iparticle = 1, nunits_tot(1)
1284  r_old(1:3, iparticle) = &
1285  particles(1)%list%els(iparticle)%r(1:3)
1286  END DO
1287 
1288  ! move the center of mass of molecule 1 to the origin
1289  start_atom = 1
1290  end_atom = nunits(mol_type(1))
1291  CALL get_center_of_mass(r_old(:, start_atom:end_atom), nunits(mol_type(1)), &
1292  center_of_mass(:), mass(1:nunits(mol_type(1)), mol_type(1)))
1293  DO iunit = start_atom, end_atom
1294  r_old(:, iunit) = r_old(:, iunit) - center_of_mass(:)
1295  END DO
1296  ! set them in the force_env, so the first molecule is ready for the energy calc
1297  DO iparticle = start_atom, end_atom
1298  particles(1)%list%els(iparticle)%r(1:3) = r_old(1:3, iparticle)
1299  END DO
1300 
1301  ! print out a notice every 1%
1302  iprint = floor(real(nvirial, kind=dp)/100.0_dp)
1303  IF (iprint == 0) iprint = 1
1304 
1305  ! we'll compute the average potential, and then integrate that, as opposed to
1306  ! integrating every orientation and then averaging
1307  ALLOCATE (mayer(1:nvirial_temps, 1:nbins))
1308 
1309  mayer(:, :) = 0.0_dp
1310 
1311  ! loop over all nvirial random configurations
1312  DO ivirial = 1, nvirial
1313 
1314  ! move molecule two back to the origin
1315  start_atom = nunits(mol_type(1)) + 1
1316  end_atom = nunits_tot(1)
1317  CALL get_center_of_mass(r_old(:, start_atom:end_atom), nunits(mol_type(2)), &
1318  center_of_mass(:), mass(1:nunits(mol_type(2)), mol_type(2)))
1319  DO iunit = start_atom, end_atom
1320  r_old(:, iunit) = r_old(:, iunit) - center_of_mass(:)
1321  END DO
1322 
1323  ! now we need a random orientation for molecule 2...this routine is
1324  ! only done in serial since it calls a random number
1325  IF (ionode) THEN
1326  CALL rotate_molecule(r_old(:, start_atom:end_atom), &
1327  mass(1:nunits(mol_type(2)), mol_type(2)), &
1328  nunits(mol_type(2)), rng_stream)
1329  END IF
1330  CALL group%bcast(r_old(:, :), source)
1331 
1332  distance = 0.0e0_dp
1333  ibin = 1
1334  DO
1335  ! find out what our stepsize is
1336  current_division = 0
1337  DO idivision = 1, nintegral_divisions
1338  IF (distance .LT. virial_cutoffs(idivision) - virial_stepsize(idivision)/2.0e0_dp) THEN
1339  current_division = idivision
1340  EXIT
1341  END IF
1342  END DO
1343  IF (current_division == 0) EXIT
1344  distance = distance + virial_stepsize(current_division)
1345 
1346  ! move the second molecule only along the x direction
1347  DO iparticle = start_atom, end_atom
1348  particles(1)%list%els(iparticle)%r(1) = r_old(1, iparticle) + distance
1349  particles(1)%list%els(iparticle)%r(2) = r_old(2, iparticle)
1350  particles(1)%list%els(iparticle)%r(3) = r_old(3, iparticle)
1351  END DO
1352 
1353  ! check for overlaps
1354  CALL check_for_overlap(force_env(1)%force_env, nchains(:, 1), nunits, loverlap, mol_type)
1355 
1356  ! compute the energy if there is no overlap
1357  ! exponent is exp(-beta*energy)-1, also called the Mayer term
1358  IF (loverlap) THEN
1359  DO itemp = 1, nvirial_temps
1360  mayer(itemp, ibin) = mayer(itemp, ibin) - 1.0_dp
1361  END DO
1362  ELSE
1363  CALL force_env_calc_energy_force(force_env(1)%force_env, &
1364  calc_force=.false.)
1365  CALL force_env_get(force_env(1)%force_env, &
1366  potential_energy=trial_energy)
1367 
1368  DO itemp = 1, nvirial_temps
1369 
1370  exponent = -beta(itemp)*trial_energy
1371 
1372  IF (exponent .GT. exp_max_val) THEN
1373  exponent = exp_max_val
1374  ELSEIF (exponent .LT. exp_min_val) THEN
1375  exponent = exp_min_val
1376  END IF
1377  mayer(itemp, ibin) = mayer(itemp, ibin) + exp(exponent) - 1.0_dp
1378  END DO
1379  END IF
1380 
1381  ibin = ibin + 1
1382  END DO
1383  ! write out some info that keeps track of where we are
1384  IF (iw > 0) THEN
1385  IF (mod(ivirial, iprint) == 0) &
1386  WRITE (iw, '(A,I6,A,I6)') ' Done with config ', ivirial, ' out of ', nvirial
1387  END IF
1388  END DO
1389 
1390  ! now we integrate this average potential
1391  mayer(:, :) = mayer(:, :)/real(nvirial, dp)
1392 
1393  DO itemp = 1, nvirial_temps
1394  integral = 0.0_dp
1395  previous_value = 0.0_dp
1396  distance = 0.0e0_dp
1397  ibin = 1
1398  DO
1399  current_division = 0
1400  DO idivision = 1, nintegral_divisions
1401  IF (distance .LT. virial_cutoffs(idivision) - virial_stepsize(idivision)/2.0e0_dp) THEN
1402  current_division = idivision
1403  EXIT
1404  END IF
1405  END DO
1406  IF (current_division == 0) EXIT
1407  distance = distance + virial_stepsize(current_division)
1408 
1409  ! now we need to integrate, using the trapazoidal method
1410  ! first, find the value of the square
1411  current_value = mayer(itemp, ibin)*distance**2
1412  square_value = previous_value*virial_stepsize(current_division)
1413  ! now the triangle that sits on top of it, which is half the size of this square...
1414  ! notice this is negative if the current value is less than the previous value
1415  triangle_value = 0.5e0_dp*((current_value - previous_value)*virial_stepsize(current_division))
1416 
1417  integral = integral + square_value + triangle_value
1418  previous_value = current_value
1419  ibin = ibin + 1
1420  END DO
1421 
1422  ! now that the integration is done, compute the second virial that results
1423  ave_virial = -2.0e0_dp*pi*integral
1424 
1425  ! convert from CP2K units to something else
1426  ave_virial = ave_virial*n_avogadro*angstrom**3/1.0e8_dp**3
1427 
1428  IF (iw > 0) THEN
1429  WRITE (iw, *)
1430  WRITE (iw, *) '*********************************************************************'
1431  WRITE (iw, '(A,F12.6,A)') ' *** Temperature = ', virial_temps(itemp), &
1432  ' ***'
1433  WRITE (iw, *) '*** ***'
1434  WRITE (iw, '(A,E12.6,A)') ' *** B2(T) = ', ave_virial, &
1435  ' cm**3/mol ***'
1436  WRITE (iw, *) '*********************************************************************'
1437  WRITE (iw, *)
1438  END IF
1439  END DO
1440 
1441  ! deallocate some stuff
1442  DEALLOCATE (mc_par)
1443  DEALLOCATE (subsys)
1444  DEALLOCATE (force_env)
1445  DEALLOCATE (particles)
1446  DEALLOCATE (cell)
1447  DEALLOCATE (virial_cutoffs)
1448  DEALLOCATE (virial_stepsize)
1449  DEALLOCATE (r_old)
1450  DEALLOCATE (mayer)
1451  DEALLOCATE (beta)
1452 
1453  END SUBROUTINE mc_compute_virial
1454 
1455 END MODULE mc_ensembles
1456 
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cell_types::get_cell
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
cp_external_control
Routines to handle the external control of CP2K.
Definition: cp_external_control.F:15
cp_external_control::external_control
subroutine, public external_control(should_stop, flag, globenv, target_time, start_time, force_check)
External manipulations during a run : when the <PROJECT_NAME>.EXIT_$runtype command is sent the progr...
Definition: cp_external_control.F:90
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition: cp_files.F:16
cp_files::open_file
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
cp_files::close_file
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
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
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...
Definition: cp_log_handling.F:515
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition: cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
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
Definition: cp_subsys_types.F:345
cp_units
unit conversion facility
Definition: cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition: cp_units.F:1179
fft_lib
Definition: fft_lib.F:7
force_env_methods
Interface for the force calculations.
Definition: force_env_methods.F:17
force_env_methods::force_env_calc_energy_force
recursive subroutine, public force_env_calc_energy_force(force_env, calc_force, consistent_energies, skip_external_control, eval_energy_forces, require_consistent_energy_force, linres, calc_stress_tensor)
Interface routine for force and energy calculations.
Definition: force_env_methods.F:193
force_env_types
Interface for the force calculations.
Definition: force_env_types.F:18
force_env_types::force_env_get
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
Definition: force_env_types.F:329
force_env_types::force_env_release
recursive subroutine, public force_env_release(force_env)
releases the given force env
Definition: force_env_types.F:200
global_types
Define type storing the global information of a run. Keep the amount of stored data small....
Definition: global_types.F:21
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::dump_xmol
integer, parameter, public dump_xmol
Definition: input_constants.F:375
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_val_get
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
Definition: input_section_types.F:1040
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
machine
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition: machine.F:106
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
mc_control
contains some general routines for dealing with the restart files and creating force_env for MC use
Definition: mc_control.F:15
mc_control::write_mc_restart
subroutine, public write_mc_restart(nnstep, mc_par, nchains, force_env)
writes the coordinates of the current step to a file that can be read in at the start of the next sim...
Definition: mc_control.F:79
mc_control::mc_create_bias_force_env
subroutine, public mc_create_bias_force_env(bias_env, r, atom_symbols, nunits_tot, para_env, box_length, nchains, input_declaration, mc_input_file, ionode)
essentially copies the cell size and coordinates of one force env to another that we will use to bias...
Definition: mc_control.F:395
mc_coordinates
contains miscellaneous subroutines used in the Monte Carlo runs,mostly geared towards changes in coor...
Definition: mc_coordinates.F:13
mc_coordinates::get_center_of_mass
subroutine, public get_center_of_mass(coordinates, natom, center_of_mass, mass)
calculates the center of mass of a given molecule
Definition: mc_coordinates.F:190
mc_coordinates::check_for_overlap
subroutine, public check_for_overlap(force_env, nchains, nunits, loverlap, mol_type, cell_length, molecule_number)
looks for overlaps (intermolecular distances less than rmin)
Definition: mc_coordinates.F:69
mc_coordinates::mc_coordinate_fold
subroutine, public mc_coordinate_fold(coordinates, nchains_tot, mol_type, mass, nunits, box_length)
folds all the coordinates into the center simulation box using a center of mass cutoff
Definition: mc_coordinates.F:236
mc_coordinates::find_mc_test_molecule
subroutine, public find_mc_test_molecule(mc_molecule_info, start_atom, box_number, molecule_type, rng_stream, box, molecule_type_old)
selects a molecule at random to perform a MC move on...you can specify the box the molecule should be...
Definition: mc_coordinates.F:716
mc_coordinates::rotate_molecule
subroutine, public rotate_molecule(r, mass, natoms, rng_stream)
rotates a molecule randomly around the center of mass, sequentially in x, y, and z directions
Definition: mc_coordinates.F:627
mc_coordinates::create_discrete_array
subroutine, public create_discrete_array(cell, discrete_array, step_size)
generates an array that tells us which sides of the simulation cell we can increase or decrease using...
Definition: mc_coordinates.F:867
mc_ensembles
Used to run the bulk of the MC simulation, doing things like choosing move types and writing data to ...
Definition: mc_ensembles.F:16
mc_ensembles::mc_compute_virial
subroutine, public mc_compute_virial(mc_env, rng_stream)
Computes the second virial coefficient of a molecule by using the integral form of the second virial ...
Definition: mc_ensembles.F:1178
mc_ensembles::mc_run_ensemble
subroutine, public mc_run_ensemble(mc_env, para_env, globenv, input_declaration, nboxes, rng_stream)
directs the program in running one or two box MC simulations
Definition: mc_ensembles.F:115
mc_environment_types
contains the subroutines for dealing with the mc_env
Definition: mc_environment_types.F:12
mc_environment_types::get_mc_env
subroutine, public get_mc_env(mc_env, mc_par, force_env)
provides a method for getting the various structures attached to an mc_env
Definition: mc_environment_types.F:95
mc_environment_types::set_mc_env
subroutine, public set_mc_env(mc_env, mc_par, force_env)
provides a method for attaching various structures to an mc_env
Definition: mc_environment_types.F:70
mc_ge_moves
contains the Monte Carlo moves that can handle more than one box, including the Quickstep move,...
Definition: mc_ge_moves.F:18
mc_ge_moves::mc_ge_volume_move
subroutine, public mc_ge_volume_move(mc_par, force_env, moves, move_updates, nnstep, old_energy, energy_check, r_old, rng_stream)
performs a Monte Carlo move that alters the volume of the simulation boxes, keeping the total volume ...
Definition: mc_ge_moves.F:1134
mc_ge_moves::mc_quickstep_move
subroutine, public mc_quickstep_move(mc_par, force_env, bias_env, moves, lreject, move_updates, energy_check, r_old, nnstep, old_energy, bias_energy_new, last_bias_energy, nboxes, box_flag, subsys, particles, rng_stream, unit_conv)
computes the acceptance of a series of biased or unbiased moves (translation, rotation,...
Definition: mc_ge_moves.F:105
mc_ge_moves::mc_ge_swap_move
subroutine, public mc_ge_swap_move(mc_par, force_env, bias_env, moves, energy_check, r_old, old_energy, input_declaration, para_env, bias_energy_old, last_bias_energy, rng_stream)
attempts a swap move between two simulation boxes
Definition: mc_ge_moves.F:428
mc_misc
contains miscellaneous subroutines used in the Monte Carlo runs, mostly I/O stuff
Definition: mc_misc.F:13
mc_misc::mc_averages_release
subroutine, public mc_averages_release(averages)
deallocates the structure that holds running averages of MC variables
Definition: mc_misc.F:77
mc_misc::final_mc_write
subroutine, public final_mc_write(mc_par, all_moves, iw, energy_check, initial_energy, final_energy, averages)
writes a bunch of simulation data to the specified unit
Definition: mc_misc.F:115
mc_misc::mc_averages_create
subroutine, public mc_averages_create(averages)
initializes the structure that holds running averages of MC variables
Definition: mc_misc.F:46
mc_move_control
control the handling of the move data in Monte Carlo (MC) simulations
Definition: mc_move_control.F:14
mc_move_control::mc_moves_release
subroutine, public mc_moves_release(moves)
deallocates all the structures and nullifies the pointer
Definition: mc_move_control.F:97
mc_move_control::write_move_stats
subroutine, public write_move_stats(moves, nnstep, unit)
writes the number of accepted and attempted moves to a file for the various move types
Definition: mc_move_control.F:273
mc_move_control::mc_move_update
subroutine, public mc_move_update(mc_par, move_updates, molecule_type, flag, nnstep, ionode)
updates the maximum displacements of a Monte Carlo simulation, based on the ratio of successful moves...
Definition: mc_move_control.F:360
mc_move_control::init_mc_moves
subroutine, public init_mc_moves(moves)
allocates and initializes the structure to record all move attempts/successes
Definition: mc_move_control.F:50
mc_moves
the various moves in Monte Carlo (MC) simulations, including change of internal conformation,...
Definition: mc_moves.F:16
mc_moves::mc_molecule_rotation
subroutine, public mc_molecule_rotation(mc_par, force_env, bias_env, moves, move_updates, box_number, start_atom, molecule_type, bias_energy, lreject, rng_stream)
rotates the given molecule randomly around the x,y, or z axis... only works for water at the moment
Definition: mc_moves.F:660
mc_moves::mc_avbmc_move
subroutine, public mc_avbmc_move(mc_par, force_env, bias_env, moves, energy_check, r_old, old_energy, start_atom_swap, target_atom, molecule_type, box_number, bias_energy_old, last_bias_energy, move_type, rng_stream)
performs either a bond or angle change move for a given molecule
Definition: mc_moves.F:1950
mc_moves::mc_volume_move
subroutine, public mc_volume_move(mc_par, force_env, moves, move_updates, old_energy, box_number, energy_check, r_old, iw, discrete_array, rng_stream)
performs a Monte Carlo move that alters the volume of the simulation box
Definition: mc_moves.F:962
mc_moves::mc_molecule_translation
subroutine, public mc_molecule_translation(mc_par, force_env, bias_env, moves, move_updates, start_atom, box_number, bias_energy, molecule_type, lreject, rng_stream)
translates the given molecule randomly in either the x,y, or z direction
Definition: mc_moves.F:436
mc_moves::mc_cluster_translation
subroutine, public mc_cluster_translation(mc_par, force_env, bias_env, moves, move_updates, box_number, bias_energy, lreject, rng_stream)
translates the cluster randomly in either the x,y, or z direction
Definition: mc_moves.F:2487
mc_moves::mc_hmc_move
subroutine, public mc_hmc_move(mc_par, force_env, globenv, moves, move_updates, old_energy, box_number, energy_check, r_old, rng_stream)
performs a hybrid Monte Carlo move that runs a short MD sequence
Definition: mc_moves.F:2349
mc_moves::mc_conformation_change
subroutine, public mc_conformation_change(mc_par, force_env, bias_env, moves, move_updates, start_atom, molecule_type, box_number, bias_energy, move_type, lreject, rng_stream)
performs either a bond or angle change move for a given molecule
Definition: mc_moves.F:140
mc_types
holds all the structure types needed for Monte Carlo, except the mc_environment_type
Definition: mc_types.F:15
mc_types::get_mc_molecule_info
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
mc_types::set_mc_par
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
mc_types::get_mc_par
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
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
parallel_rng_types
Parallel (pseudo)random number generator (RNG) for multiple streams and substreams of random numbers.
Definition: parallel_rng_types.F:51
particle_list_types
represent a simple array based list of the given type
Definition: particle_list_types.F:15
particle_methods
Define methods related to particle_type.
Definition: particle_methods.F:14
particle_methods::write_particle_coordinates
subroutine, public write_particle_coordinates(particle_set, iunit, output_format, content, title, cell, array, unit_conv, charge_occup, charge_beta, charge_extended, print_kind)
Should be able to write a few formats e.g. xmol, and some binary format (dcd) some format can be used...
Definition: particle_methods.F:183
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::boltzmann
real(kind=dp), parameter, public boltzmann
Definition: physcon.F:129
physcon::n_avogadro
real(kind=dp), parameter, public n_avogadro
Definition: physcon.F:126
physcon::joule
real(kind=dp), parameter, public joule
Definition: physcon.F:159
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition: physcon.F:144