(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
16MODULE 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
99CONTAINS
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
1455END 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_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
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
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
cell_types::cell_p_type
Definition cell_types.F:66
cp_subsys_types::cp_subsys_p_type
represent a pointer to a subsys, to be able to create arrays of pointers
Definition cp_subsys_types.F:122
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
force_env_types::force_env_p_type
allows for the creation of an array of force_env
Definition force_env_types.F:167
global_types::global_environment_type
contains the initially parsed file and the initial parallel environment
Definition global_types.F:66
input_section_types::section_type
represent a section of the input file
Definition input_section_types.F:102
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
mc_environment_types::mc_environment_p_type
Definition mc_environment_types.F:29
mc_types::mc_averages_p_type
Definition mc_types.F:199
mc_types::mc_input_file_type
Definition mc_types.F:155
mc_types::mc_molecule_info_type
Definition mc_types.F:169
mc_types::mc_moves_p_type
Definition mc_types.F:239
mc_types::mc_simulation_parameters_p_type
Definition mc_types.F:184
message_passing::mp_comm_type
Definition message_passing.F:189
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
parallel_rng_types::rng_stream_type
Definition parallel_rng_types.F:157
particle_list_types::particle_list_p_type
represents a pointer to a list
Definition particle_list_types.F:59
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46