(git:71c3ab0)
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-2026 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 > 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) < 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) < 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) < 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)) /= 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)) /= 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 /= 1) THEN
430 cpabort('ldiscrete=.true. ONLY for systems with 1 box')
431 END IF
432 CALL create_discrete_array(abc(:, 1), discrete_array(:, :), &
433 discrete_step)
434 END IF
435
436 ! find out how many steps we're doing...change the updates to be in cycles
437 ! if the total number of steps is measured in cycles
438 IF (.NOT. lstop) THEN
439 nstep = nstep*nchain_total
440 iuptrans = iuptrans*nchain_total
441 iupvolume = iupvolume*nchain_total
442 END IF
443
444 DO nnstep = nstart + 1, nstart + nstep
445
446 IF (mod(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
447 WRITE (iw, *)
448 WRITE (iw, *) "------- On Monte Carlo Step ", nnstep
449 END IF
450
451 IF (ionode) rand = rng_stream%next()
452 ! broadcast the random number, to make sure we're on the same move
453 CALL group%bcast(rand, source)
454
455 IF (rand < pmvolume) THEN
456
457 IF (mod(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
458 WRITE (iw, *) "Attempting a volume move"
459 WRITE (iw, *)
460 END IF
461
462 SELECT CASE (ensemble)
463 CASE ("TRADITIONAL")
464 CALL mc_volume_move(mc_par(1)%mc_par, &
465 force_env(1)%force_env, &
466 moves(1, 1)%moves, move_updates(1, 1)%moves, &
467 old_energy(1), 1, &
468 energy_check(1), r_old(:, :, 1), iw, discrete_array(:, :), &
469 rng_stream)
470 CASE ("GEMC_NVT")
471 CALL mc_ge_volume_move(mc_par, force_env, moves, &
472 move_updates, nnstep, old_energy, energy_check, &
473 r_old, rng_stream)
474 CASE ("GEMC_NPT")
475 ! we need to select a box based on the probability given in the input file
476 IF (ionode) rand = rng_stream%next()
477 CALL group%bcast(rand, source)
478
479 DO ibox = 1, nboxes
480 IF (rand <= pmvol_box(ibox)) THEN
481 box_number = ibox
482 EXIT
483 END IF
484 END DO
485
486 CALL mc_volume_move(mc_par(box_number)%mc_par, &
487 force_env(box_number)%force_env, &
488 moves(1, box_number)%moves, &
489 move_updates(1, box_number)%moves, &
490 old_energy(box_number), box_number, &
491 energy_check(box_number), r_old(:, :, box_number), iw, &
492 discrete_array(:, :), &
493 rng_stream)
494 END SELECT
495
496! update all the pointers here, because otherwise we may pass wrong information when we're making a bias environment
497 DO ibox = 1, nboxes
498 CALL force_env_get(force_env(ibox)%force_env, &
499 subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
500 CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
501 CALL cp_subsys_get(oldsys(ibox)%subsys, &
502 particles=particles_old(ibox)%list)
503 END DO
504
505 ! we need a new biasing environment now, if we're into that sort of thing
506 IF (lbias) THEN
507 DO ibox = 1, nboxes
508 CALL force_env_release(bias_env(ibox)%force_env)
509 ! determine the atom names of every particle
510 ALLOCATE (atom_names_box(1:nunits_tot(ibox)))
511 start_mol = 1
512 DO jbox = 1, ibox - 1
513 start_mol = start_mol + sum(nchains(:, jbox))
514 END DO
515 end_mol = start_mol + sum(nchains(:, ibox)) - 1
516 atom_number = 1
517 DO imolecule = 1, sum(nchains(:, ibox))
518 DO iunit = 1, nunits(mol_type(imolecule + start_mol - 1))
519 atom_names_box(atom_number) = &
520 atom_names(iunit, mol_type(imolecule + start_mol - 1))
521 atom_number = atom_number + 1
522 END DO
523 END DO
524
525! need to find out what the cell lengths are
526 CALL force_env_get(force_env(ibox)%force_env, &
527 subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
528 CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
529
530 CALL get_mc_par(mc_par(ibox)%mc_par, &
531 mc_bias_file=mc_bias_file)
532 nchains_box => nchains(:, ibox)
533
534 CALL mc_create_bias_force_env(bias_env(ibox)%force_env, &
535 r_old(:, :, ibox), atom_names_box(:), nunits_tot(ibox), &
536 para_env, abc(:, ibox), nchains_box, input_declaration, &
537 mc_bias_file, ionode)
538
539 IF (sum(nchains(:, ibox)) /= 0) THEN
540 CALL force_env_calc_energy_force( &
541 bias_env(ibox)%force_env, &
542 calc_force=.false.)
543 CALL force_env_get(bias_env(ibox)%force_env, &
544 potential_energy=last_bias_energy(ibox))
545 ELSE
546 last_bias_energy(ibox) = 0.0e0_dp
547 END IF
548 bias_energy(ibox) = last_bias_energy(ibox)
549 DEALLOCATE (atom_names_box)
550 END DO
551 END IF
552
553 ELSE IF (rand < pmswap) THEN
554
555 ! try a swap move
556 IF (mod(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
557 WRITE (iw, *) "Attempting a swap move"
558 WRITE (iw, *)
559 END IF
560
561 CALL mc_ge_swap_move(mc_par, force_env, bias_env, moves, &
562 energy_check(:), r_old(:, :, :), old_energy(:), input_declaration, &
563 para_env, bias_energy(:), last_bias_energy(:), rng_stream)
564
565 ! the number of molecules may have changed, which deallocated the whole
566 ! mc_molecule_info structure
567 CALL get_mc_par(mc_par(1)%mc_par, mc_molecule_info=mc_molecule_info)
568 CALL get_mc_molecule_info(mc_molecule_info, conf_prob=conf_prob, &
569 nchains=nchains, nmol_types=nmol_types, nunits_tot=nunits_tot, &
570 mol_type=mol_type, nchain_total=nchain_total, nunits=nunits, &
571 atom_names=atom_names, mass=mass)
572
573 ELSE IF (rand < pmhmc) THEN
574! try hybrid Monte Carlo
575 IF (mod(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
576 WRITE (iw, *) "Attempting a hybrid Monte Carlo move"
577 WRITE (iw, *)
578 END IF
579
580! pick a box at random
581 IF (ionode) rand = rng_stream%next()
582 CALL group%bcast(rand, source)
583
584 DO ibox = 1, nboxes
585 IF (rand <= pmhmc_box(ibox)) THEN
586 box_number = ibox
587 EXIT
588 END IF
589 END DO
590
591 CALL mc_hmc_move(mc_par(box_number)%mc_par, &
592 force_env(box_number)%force_env, globenv, &
593 moves(1, box_number)%moves, &
594 move_updates(1, box_number)%moves, &
595 old_energy(box_number), box_number, &
596 energy_check(box_number), r_old(:, :, box_number), &
597 rng_stream)
598
599 ELSE IF (rand < pmavbmc) THEN
600 ! try an AVBMC move
601 IF (mod(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
602 WRITE (iw, *) "Attempting an AVBMC1 move"
603 WRITE (iw, *)
604 END IF
605
606 ! first, pick a box to do it for
607 IF (ionode) rand = rng_stream%next()
608 CALL group%bcast(rand, source)
609
610 IF (nboxes == 2) THEN
611 IF (rand < 0.1e0_dp) THEN
612 box_number = 1
613 ELSE
614 box_number = 2
615 END IF
616 ELSE
617 box_number = 1
618 END IF
619
620 ! now pick a molecule type to do it for
621 IF (ionode) rand = rng_stream%next()
622 CALL group%bcast(rand, source)
623 molecule_type_swap = 0
624 DO imol_type = 1, nmol_types
625 IF (rand < pmavbmc_mol(imol_type)) THEN
626 molecule_type_swap = imol_type
627 EXIT
628 END IF
629 END DO
630 IF (molecule_type_swap == 0) THEN
631 cpabort('Did not choose a molecule type to swap...check AVBMC input')
632 END IF
633
634 ! now pick a molecule, automatically rejecting the move if the
635 ! box is empty or only has one molecule
636 IF (sum(nchains(:, box_number)) <= 1) THEN
637 ! indicate that we tried a move
638 moves(molecule_type_swap, box_number)%moves%empty_avbmc = &
639 moves(molecule_type_swap, box_number)%moves%empty_avbmc + 1
640 ELSE
641
642 ! pick a molecule to be swapped in the box
643 IF (ionode) THEN
644 CALL find_mc_test_molecule(mc_molecule_info, &
645 start_atom_swap, idum, jdum, rng_stream, &
646 box=box_number, molecule_type_old=molecule_type_swap)
647
648 ! pick a molecule to act as the target in the box...we don't care what type
649 DO
650 CALL find_mc_test_molecule(mc_molecule_info, &
651 start_atom_target, idum, molecule_type_target, &
652 rng_stream, box=box_number)
653 IF (start_atom_swap /= start_atom_target) THEN
654 start_atom_target = start_atom_target + &
655 avbmc_atom(molecule_type_target) - 1
656 EXIT
657 END IF
658 END DO
659
660 ! choose if we're swapping into the bonded region of mol_target, or
661 ! into the nonbonded region
662 rand = rng_stream%next()
663
664 END IF
665 CALL group%bcast(start_atom_swap, source)
666 CALL group%bcast(box_number, source)
667 CALL group%bcast(start_atom_target, source)
668 CALL group%bcast(rand, source)
669
670 IF (rand < pbias(molecule_type_swap)) THEN
671 move_type_avbmc = 'in'
672 ELSE
673 move_type_avbmc = 'out'
674 END IF
675
676 CALL mc_avbmc_move(mc_par(box_number)%mc_par, &
677 force_env(box_number)%force_env, &
678 bias_env(box_number)%force_env, &
679 moves(molecule_type_swap, box_number)%moves, &
680 energy_check(box_number), &
681 r_old(:, :, box_number), old_energy(box_number), &
682 start_atom_swap, start_atom_target, molecule_type_swap, &
683 box_number, bias_energy(box_number), &
684 last_bias_energy(box_number), &
685 move_type_avbmc, rng_stream)
686
687 END IF
688
689 ELSE
690
691 IF (mod(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
692 WRITE (iw, *) "Attempting an inner move"
693 WRITE (iw, *)
694 END IF
695
696 DO imove = 1, nmoves
697
698 IF (ionode) rand = rng_stream%next()
699 CALL group%bcast(rand, source)
700 IF (rand < pmtraion) THEN
701 ! change molecular conformation
702 ! first, pick a box to do it for
703 IF (ionode) rand = rng_stream%next()
704 CALL group%bcast(rand, source)
705 IF (nboxes == 2) THEN
706 IF (rand < 0.75e0_dp) THEN
707 box_number = 1
708 ELSE
709 box_number = 2
710 END IF
711 ELSE
712 box_number = 1
713 END IF
714
715 ! figure out which molecule type we're looking for
716 IF (ionode) rand = rng_stream%next()
717 CALL group%bcast(rand, source)
718 molecule_type = 0
719 DO imol_type = 1, nmol_types
720 IF (rand < pmtraion_mol(imol_type)) THEN
721 molecule_type = imol_type
722 EXIT
723 END IF
724 END DO
725 IF (molecule_type == 0) CALL cp_abort( &
726 __location__, &
727 'Did not choose a molecule type to conf change...PMTRAION_MOL should not be all 0.0')
728
729 ! now pick a molecule, automatically rejecting the move if the
730 ! box is empty
731 IF (nchains(molecule_type, box_number) == 0) THEN
732 ! indicate that we tried a move
733 moves(molecule_type, box_number)%moves%empty_conf = &
734 moves(molecule_type, box_number)%moves%empty_conf + 1
735 ELSE
736 ! pick a molecule in the box
737 IF (ionode) THEN
738 CALL find_mc_test_molecule(mc_molecule_info, &
739 start_atom, idum, &
740 jdum, rng_stream, &
741 box=box_number, molecule_type_old=molecule_type)
742
743 ! choose if we're changing a bond length or an angle
744 rand = rng_stream%next()
745 END IF
746 CALL group%bcast(rand, source)
747 CALL group%bcast(start_atom, source)
748 CALL group%bcast(box_number, source)
749 CALL group%bcast(molecule_type, source)
750
751 ! figure out what kind of move we're doing
752 IF (rand < conf_prob(1, molecule_type)) THEN
753 move_type = 'bond'
754 ELSE IF (rand < (conf_prob(1, molecule_type) + &
755 conf_prob(2, molecule_type))) THEN
756 move_type = 'angle'
757 ELSE
758 move_type = 'dihedral'
759 END IF
760 box_flag(box_number) = 1
761 CALL mc_conformation_change(mc_par(box_number)%mc_par, &
762 force_env(box_number)%force_env, &
763 bias_env(box_number)%force_env, &
764 moves(molecule_type, box_number)%moves, &
765 move_updates(molecule_type, box_number)%moves, &
766 start_atom, molecule_type, box_number, &
767 bias_energy(box_number), &
768 move_type, lreject, rng_stream)
769 IF (lreject) EXIT
770 END IF
771 ELSE IF (rand < pmtrans) THEN
772 ! translate a whole molecule in the system
773 ! pick a molecule type
774 IF (ionode) rand = rng_stream%next()
775 CALL group%bcast(rand, source)
776 molecule_type = 0
777 DO imol_type = 1, nmol_types
778 IF (rand < pmtrans_mol(imol_type)) THEN
779 molecule_type = imol_type
780 EXIT
781 END IF
782 END DO
783 IF (molecule_type == 0) CALL cp_abort( &
784 __location__, &
785 'Did not choose a molecule type to translate...PMTRANS_MOL should not be all 0.0')
786
787 ! now pick a molecule of that type
788 IF (ionode) THEN
789 CALL find_mc_test_molecule(mc_molecule_info, &
790 start_atom, box_number, idum, rng_stream, &
791 molecule_type_old=molecule_type)
792 END IF
793 CALL group%bcast(start_atom, source)
794 CALL group%bcast(box_number, source)
795 box_flag(box_number) = 1
796 CALL mc_molecule_translation(mc_par(box_number)%mc_par, &
797 force_env(box_number)%force_env, &
798 bias_env(box_number)%force_env, &
799 moves(molecule_type, box_number)%moves, &
800 move_updates(molecule_type, box_number)%moves, &
801 start_atom, box_number, bias_energy(box_number), &
802 molecule_type, lreject, rng_stream)
803 IF (lreject) EXIT
804 ELSE IF (rand < pmcltrans) THEN
805 ! translate a whole cluster in the system
806 ! first, pick a box to do it for
807 IF (ionode) rand = rng_stream%next()
808 CALL group%bcast(rand, source)
809
810 DO ibox = 1, nboxes
811 IF (rand <= pmclus_box(ibox)) THEN
812 box_number = ibox
813 EXIT
814 END IF
815 END DO
816 box_flag(box_number) = 1
817 CALL mc_cluster_translation(mc_par(box_number)%mc_par, &
818 force_env(box_number)%force_env, &
819 bias_env(box_number)%force_env, &
820 moves(1, box_number)%moves, &
821 move_updates(1, box_number)%moves, &
822 box_number, bias_energy(box_number), &
823 lreject, rng_stream)
824 IF (lreject) EXIT
825 ELSE
826 ! rotate a whole molecule in the system
827 ! pick a molecule type
828 IF (ionode) rand = rng_stream%next()
829 CALL group%bcast(rand, source)
830 molecule_type = 0
831 DO imol_type = 1, nmol_types
832 IF (rand < pmrot_mol(imol_type)) THEN
833 molecule_type = imol_type
834 EXIT
835 END IF
836 END DO
837 IF (molecule_type == 0) CALL cp_abort( &
838 __location__, &
839 'Did not choose a molecule type to rotate...PMROT_MOL should not be all 0.0')
840
841 IF (ionode) THEN
842 CALL find_mc_test_molecule(mc_molecule_info, &
843 start_atom, box_number, idum, rng_stream, &
844 molecule_type_old=molecule_type)
845 END IF
846 CALL group%bcast(start_atom, source)
847 CALL group%bcast(box_number, source)
848 box_flag(box_number) = 1
849 CALL mc_molecule_rotation(mc_par(box_number)%mc_par, &
850 force_env(box_number)%force_env, &
851 bias_env(box_number)%force_env, &
852 moves(molecule_type, box_number)%moves, &
853 move_updates(molecule_type, box_number)%moves, &
854 box_number, start_atom, &
855 molecule_type, bias_energy(box_number), &
856 lreject, rng_stream)
857 IF (lreject) EXIT
858 END IF
859
860 END DO
861
862 ! now do a Quickstep calculation to see if we accept the sequence
863 CALL mc_quickstep_move(mc_par, force_env, bias_env, &
864 moves, lreject, move_updates, energy_check(:), r_old(:, :, :), &
865 nnstep, old_energy(:), bias_energy(:), last_bias_energy(:), &
866 nboxes, box_flag(:), oldsys, particles_old, &
867 rng_stream, unit_conv)
868
869 END IF
870
871 ! make sure the pointers are pointing correctly since the subsys may
872 ! have changed
873 DO ibox = 1, nboxes
874 CALL force_env_get(force_env(ibox)%force_env, &
875 subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
876 CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
877 CALL cp_subsys_get(oldsys(ibox)%subsys, &
878 particles=particles_old(ibox)%list)
879 END DO
880
881 IF (ionode) THEN
882
883 IF (mod(nnstep, iprint) == 0) THEN
884 WRITE (com_ene, *) nnstep, old_energy(1:nboxes)
885
886 DO ibox = 1, nboxes
887
888 ! write the molecule information
889 WRITE (com_mol, *) nnstep, nchains(:, ibox)
890
891 ! write the move statistics to file
892 DO itype = 1, nmol_types
893 CALL write_move_stats(moves(itype, ibox)%moves, &
894 nnstep, move_unit(ibox))
895 END DO
896
897 ! write a restart file
898 CALL write_mc_restart(nnstep, mc_par(ibox)%mc_par, &
899 nchains(:, ibox), force_env(ibox)%force_env)
900
901 ! write cell lengths
902 WRITE (cell_unit, *) nnstep, abc(1:3, ibox)*angstrom
903
904 ! write particle coordinates
905 WRITE (cbox, '(I4)') ibox
906 WRITE (cstep, '(I8)') nnstep
907 IF (sum(nchains(:, ibox)) == 0) THEN
908 WRITE (com_crd, *) ' 0'
909 WRITE (com_crd, *) 'BOX '//trim(adjustl(cbox))// &
910 ', STEP '//trim(adjustl(cstep))
911 ELSE
912 CALL write_particle_coordinates( &
913 particles_old(ibox)%list%els, &
914 com_crd, dump_xmol, 'POS', &
915 'BOX '//trim(adjustl(cbox))// &
916 ', STEP '//trim(adjustl(cstep)), &
917 unit_conv=unit_conv)
918 END IF
919 END DO
920 END IF ! end the things we only do every iprint moves
921
922 DO ibox = 1, nboxes
923 ! compute some averages
924 averages(ibox)%averages%ave_energy = &
925 averages(ibox)%averages%ave_energy*real(nnstep - &
926 nstart - 1, dp)/real(nnstep - nstart, dp) + &
927 old_energy(ibox)/real(nnstep - nstart, dp)
928 averages(ibox)%averages%molecules = &
929 averages(ibox)%averages%molecules*real(nnstep - &
930 nstart - 1, dp)/real(nnstep - nstart, dp) + &
931 REAL(sum(nchains(:, ibox)), dp)/REAL(nnstep - nstart, dp)
932 averages(ibox)%averages%ave_volume = &
933 averages(ibox)%averages%ave_volume* &
934 REAL(nnstep - nstart - 1, dp)/REAL(nnstep - nstart, dp) + &
935 abc(1, ibox)*abc(2, ibox)*abc(3, ibox)/ &
936 REAL(nnstep - nstart, dp)
937
938 ! flush the buffers to the files
939 CALL m_flush(data_unit(ibox))
940 CALL m_flush(diff(ibox))
941 CALL m_flush(move_unit(ibox))
942 CALL m_flush(cl(ibox))
943 CALL m_flush(rm(ibox))
944
945 END DO
946
947 ! flush more buffers to the files
948 CALL m_flush(cell_unit)
949 CALL m_flush(com_ene)
950 CALL m_flush(com_crd)
951 CALL m_flush(com_mol)
952
953 END IF
954
955 ! reset the box flags
956 box_flag(:) = 0
957
958 ! check to see if EXIT file exists...if so, end the calculation
959 CALL external_control(should_stop, "MC", globenv=globenv)
960 IF (should_stop) EXIT
961
962 ! update the move displacements, if necessary
963 DO ibox = 1, nboxes
964 IF (mod(nnstep - nstart, iuptrans) == 0) THEN
965 DO itype = 1, nmol_types
966 CALL mc_move_update(mc_par(ibox)%mc_par, &
967 move_updates(itype, ibox)%moves, itype, &
968 "trans", nnstep, ionode)
969 END DO
970 END IF
971
972 IF (mod(nnstep - nstart, iupvolume) == 0) THEN
973 CALL mc_move_update(mc_par(ibox)%mc_par, &
974 move_updates(1, ibox)%moves, 1337, &
975 "volume", nnstep, ionode)
976 END IF
977 END DO
978
979 ! check to see if there are any overlaps in the boxes, and fold coordinates
980! don't care about overlaps if we're only doing HMC
981 IF (.NOT. lhmc) THEN
982 DO ibox = 1, nboxes
983 IF (sum(nchains(:, ibox)) /= 0) THEN
984 start_mol = 1
985 DO jbox = 1, ibox - 1
986 start_mol = start_mol + sum(nchains(:, jbox))
987 END DO
988 end_mol = start_mol + sum(nchains(:, ibox)) - 1
989 CALL check_for_overlap(force_env(ibox)%force_env, &
990 nchains(:, ibox), nunits, loverlap, &
991 mol_type(start_mol:end_mol))
992 IF (loverlap) THEN
993 IF (iw > 0) WRITE (iw, *) nnstep
994 cpabort('coordinate overlap at the end of the above step')
995 ! now fold the coordinates...don't do this anywhere but here, because
996 ! we can get screwed up with the mc_molecule_info stuff (like in swap move)...
997 ! this is kind of ugly, with allocated and deallocating every time
998 ALLOCATE (r_temp(1:3, 1:nunits_tot(ibox)))
999
1000 DO iunit = 1, nunits_tot(ibox)
1001 r_temp(1:3, iunit) = &
1002 particles_old(ibox)%list%els(iunit)%r(1:3)
1003 END DO
1004
1005 CALL mc_coordinate_fold(r_temp(:, :), &
1006 sum(nchains(:, ibox)), mol_type(start_mol:end_mol), &
1007 mass, nunits, abc(1:3, ibox))
1008
1009 ! save the folded coordinates
1010 DO iunit = 1, nunits_tot(ibox)
1011 r_old(1:3, iunit, ibox) = r_temp(1:3, iunit)
1012 particles_old(ibox)%list%els(iunit)%r(1:3) = &
1013 r_temp(1:3, iunit)
1014 END DO
1015
1016 ! if we're biasing, we need to do the same
1017 IF (lbias) THEN
1018 CALL force_env_get(bias_env(ibox)%force_env, &
1019 subsys=biassys)
1020 CALL cp_subsys_get(biassys, &
1021 particles=particles_bias)
1022
1023 DO iunit = 1, nunits_tot(ibox)
1024 particles_bias%els(iunit)%r(1:3) = &
1025 r_temp(1:3, iunit)
1026 END DO
1027 END IF
1028
1029 DEALLOCATE (r_temp)
1030 END IF
1031 END IF
1032 END DO
1033 END IF
1034
1035 !debug code
1036 IF (debug_this_module) THEN
1037 DO ibox = 1, nboxes
1038 IF (sum(nchains(:, ibox)) /= 0) THEN
1039 CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
1040 calc_force=.false.)
1041 CALL force_env_get(force_env(ibox)%force_env, &
1042 potential_energy=test_energy)
1043 ELSE
1044 test_energy = 0.0e0_dp
1045 END IF
1046
1047 IF (abs(initial_energy(ibox) + energy_check(ibox) - &
1048 test_energy) > 0.0000001e0_dp) THEN
1049 IF (iw > 0) THEN
1050 WRITE (iw, *) '!!!!!!! We have an energy problem. !!!!!!!!'
1051 WRITE (iw, '(A,T64,F16.10)') 'Final Energy = ', test_energy
1052 WRITE (iw, '(A,T64,F16.10)') 'Initial Energy+energy_check=', &
1053 initial_energy(ibox) + energy_check(ibox)
1054 WRITE (iw, *) 'Box ', ibox
1055 WRITE (iw, *) 'nchains ', nchains(:, ibox)
1056 END IF
1057 cpabort('!!!!!!! We have an energy problem. !!!!!!!!')
1058 END IF
1059 END DO
1060 END IF
1061 END DO
1062
1063 ! write a restart file
1064 IF (ionode) THEN
1065 DO ibox = 1, nboxes
1066 CALL write_mc_restart(nnstep, mc_par(ibox)%mc_par, &
1067 nchains(:, ibox), force_env(ibox)%force_env)
1068 END DO
1069 END IF
1070
1071 ! calculate the final energy
1072 DO ibox = 1, nboxes
1073 IF (sum(nchains(:, ibox)) /= 0) THEN
1074 CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
1075 calc_force=.false.)
1076 CALL force_env_get(force_env(ibox)%force_env, &
1077 potential_energy=final_energy(ibox))
1078 ELSE
1079 final_energy(ibox) = 0.0e0_dp
1080 END IF
1081 IF (lbias) THEN
1082 CALL force_env_release(bias_env(ibox)%force_env)
1083 END IF
1084 END DO
1085
1086 ! do some stuff in serial
1087 IF (ionode .OR. (iw > 0)) THEN
1088
1089 WRITE (com_ene, *) 'Final Energies: ', &
1090 final_energy(1:nboxes)
1091
1092 DO ibox = 1, nboxes
1093 WRITE (cbox, '(I4)') ibox
1094 IF (sum(nchains(:, ibox)) == 0) THEN
1095 WRITE (com_crd, *) ' 0'
1096 WRITE (com_crd, *) 'BOX '//trim(adjustl(cbox))
1097 ELSE
1098 CALL write_particle_coordinates( &
1099 particles_old(ibox)%list%els, &
1100 com_crd, dump_xmol, 'POS', &
1101 'FINAL BOX '//trim(adjustl(cbox)), unit_conv=unit_conv)
1102 END IF
1103
1104 ! write a bunch of data to the screen
1105 WRITE (iw, '(A)') &
1106 '------------------------------------------------'
1107 WRITE (iw, '(A,I1,A)') &
1108 '| BOX ', ibox, &
1109 ' |'
1110 WRITE (iw, '(A)') &
1111 '------------------------------------------------'
1112 test_moves => moves(:, ibox)
1113 CALL final_mc_write(mc_par(ibox)%mc_par, test_moves, &
1114 iw, energy_check(ibox), &
1115 initial_energy(ibox), final_energy(ibox), &
1116 averages(ibox)%averages)
1117
1118 ! close any open files
1119 CALL close_file(unit_number=diff(ibox))
1120 CALL close_file(unit_number=data_unit(ibox))
1121 CALL close_file(unit_number=move_unit(ibox))
1122 CALL close_file(unit_number=cl(ibox))
1123 CALL close_file(unit_number=rm(ibox))
1124 END DO
1125
1126 ! close some more files
1127 CALL close_file(unit_number=cell_unit)
1128 CALL close_file(unit_number=com_ene)
1129 CALL close_file(unit_number=com_crd)
1130 CALL close_file(unit_number=com_mol)
1131 END IF
1132
1133 DO ibox = 1, nboxes
1134 CALL set_mc_env(mc_env(ibox)%mc_env, &
1135 mc_par=mc_par(ibox)%mc_par, &
1136 force_env=force_env(ibox)%force_env)
1137
1138 ! deallocate some stuff
1139 DO itype = 1, nmol_types
1140 CALL mc_moves_release(move_updates(itype, ibox)%moves)
1141 CALL mc_moves_release(moves(itype, ibox)%moves)
1142 END DO
1143 CALL mc_averages_release(averages(ibox)%averages)
1144 END DO
1145
1146 DEALLOCATE (pmhmc_box)
1147 DEALLOCATE (pmvol_box)
1148 DEALLOCATE (pmclus_box)
1149 DEALLOCATE (r_old)
1150 DEALLOCATE (force_env)
1151 DEALLOCATE (bias_env)
1152 DEALLOCATE (cell)
1153 DEALLOCATE (particles_old)
1154 DEALLOCATE (oldsys)
1155 DEALLOCATE (averages)
1156 DEALLOCATE (moves)
1157 DEALLOCATE (move_updates)
1158 DEALLOCATE (mc_par)
1159
1160 ! end the timing
1161 CALL timestop(handle)
1162
1163 END SUBROUTINE mc_run_ensemble
1164
1165! **************************************************************************************************
1166!> \brief Computes the second virial coefficient of a molecule by using the integral form
1167!> of the second virial coefficient found in McQuarrie "Statistical Thermodynamics",
1168!> B2(T) = -2Pi Int 0toInf [ Exp[-beta*u(r)] -1] r^2 dr Eq. 15-25
1169!> I use trapazoidal integration with various step sizes
1170!> (the integral is broken up into three parts, currently, but that's easily
1171!> changed by the first variables found below). It generates nvirial configurations,
1172!> doing the integration for each one, and then averages all the B2(T) to produce
1173!> the final answer.
1174!> \param mc_env a pointer that contains all mc_env for all the simulation
1175!> boxes
1176!> \param rng_stream the stream we pull random numbers from
1177!>
1178!> Suitable for parallel.
1179!> \author MJM
1180! **************************************************************************************************
1181 SUBROUTINE mc_compute_virial(mc_env, rng_stream)
1182
1183 TYPE(mc_environment_p_type), DIMENSION(:), POINTER :: mc_env
1184 TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
1185
1186 INTEGER :: current_division, end_atom, ibin, idivision, iparticle, iprint, itemp, iunit, &
1187 ivirial, iw, nbins, nchain_total, nintegral_divisions, nmol_types, nvirial, &
1188 nvirial_temps, source, start_atom
1189 INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot
1190 INTEGER, DIMENSION(:, :), POINTER :: nchains
1191 LOGICAL :: ionode, loverlap
1192 REAL(dp), DIMENSION(:), POINTER :: beta, virial_cutoffs, virial_stepsize, &
1193 virial_temps
1194 REAL(dp), DIMENSION(:, :), POINTER :: mass, mayer, r_old
1195 REAL(kind=dp) :: ave_virial, current_value, distance, exp_max_val, exp_min_val, exponent, &
1196 integral, previous_value, square_value, trial_energy, triangle_value
1197 REAL(kind=dp), DIMENSION(1:3) :: abc, center_of_mass
1198 TYPE(cell_p_type), DIMENSION(:), POINTER :: cell
1199 TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsys
1200 TYPE(force_env_p_type), DIMENSION(:), POINTER :: force_env
1201 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
1202 TYPE(mc_simulation_parameters_p_type), &
1203 DIMENSION(:), POINTER :: mc_par
1204 TYPE(mp_comm_type) :: group
1205 TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
1206
1207! these are current magic numbers for how we compute the virial...
1208! we break it up into three parts to integrate the function so provide
1209! better statistics
1210
1211 nintegral_divisions = 3
1212 ALLOCATE (virial_cutoffs(1:nintegral_divisions))
1213 ALLOCATE (virial_stepsize(1:nintegral_divisions))
1214 virial_cutoffs(1) = 8.0 ! first distance, in bohr
1215 virial_cutoffs(2) = 13.0 ! second distance, in bohr
1216 virial_cutoffs(3) = 22.0 ! maximum distance, in bohr
1217 virial_stepsize(1) = 0.04 ! stepsize from 0 to virial_cutoffs(1)
1218 virial_stepsize(2) = 0.1
1219 virial_stepsize(3) = 0.2
1220
1221 nbins = ceiling(virial_cutoffs(1)/virial_stepsize(1) + (virial_cutoffs(2) - virial_cutoffs(1))/ &
1222 virial_stepsize(2) + (virial_cutoffs(3) - virial_cutoffs(2))/virial_stepsize(3))
1223
1224 ! figure out what the default write unit is
1225 iw = cp_logger_get_default_io_unit()
1226
1227 ! allocate a whole bunch of stuff based on how many boxes we have
1228 ALLOCATE (force_env(1:1))
1229 ALLOCATE (cell(1:1))
1230 ALLOCATE (particles(1:1))
1231 ALLOCATE (subsys(1:1))
1232 ALLOCATE (mc_par(1:1))
1233
1234 CALL get_mc_env(mc_env(1)%mc_env, &
1235 mc_par=mc_par(1)%mc_par, &
1236 force_env=force_env(1)%force_env)
1237
1238 ! get some data out of mc_par
1239 CALL get_mc_par(mc_par(1)%mc_par, &
1240 exp_max_val=exp_max_val, &
1241 exp_min_val=exp_min_val, nvirial=nvirial, &
1242 ionode=ionode, source=source, group=group, &
1243 mc_molecule_info=mc_molecule_info, virial_temps=virial_temps)
1244
1245 IF (iw > 0) THEN
1246 WRITE (iw, *)
1247 WRITE (iw, *)
1248 WRITE (iw, *) 'Beginning the calculation of the second virial coefficient'
1249 WRITE (iw, *)
1250 WRITE (iw, *)
1251 END IF
1252
1253 ! get some data from the molecule types
1254 CALL get_mc_molecule_info(mc_molecule_info, &
1255 nchains=nchains, nmol_types=nmol_types, nunits_tot=nunits_tot, &
1256 mol_type=mol_type, nchain_total=nchain_total, nunits=nunits, &
1257 mass=mass)
1258
1259 nvirial_temps = SIZE(virial_temps)
1260 ALLOCATE (beta(1:nvirial_temps))
1261
1262 DO itemp = 1, nvirial_temps
1263 beta(itemp) = 1/virial_temps(itemp)/boltzmann*joule
1264 END DO
1265
1266 ! get the subsystems and the cell information
1267 CALL force_env_get(force_env(1)%force_env, &
1268 subsys=subsys(1)%subsys, cell=cell(1)%cell)
1269 CALL get_cell(cell(1)%cell, abc=abc(:))
1270 CALL cp_subsys_get(subsys(1)%subsys, &
1271 particles=particles(1)%list)
1272
1273 ! check and make sure the box is big enough
1274 IF (abc(1) /= abc(2) .OR. abc(2) /= abc(3)) THEN
1275 cpabort('The box needs to be cubic for a virial calculation (it is easiest).')
1276 END IF
1277 IF (virial_cutoffs(nintegral_divisions) > abc(1)/2.0e0_dp) THEN
1278 IF (iw > 0) THEN
1279 WRITE (iw, *) "Box length ", abc(1)*angstrom, " virial cutoff ", &
1280 virial_cutoffs(nintegral_divisions)*angstrom
1281 END IF
1282 cpabort('You need a bigger box to deal with this virial cutoff (see above).')
1283 END IF
1284
1285 ! store the coordinates of the molecules in an array so we can work with it
1286 ALLOCATE (r_old(1:3, 1:nunits_tot(1)))
1287
1288 DO iparticle = 1, nunits_tot(1)
1289 r_old(1:3, iparticle) = &
1290 particles(1)%list%els(iparticle)%r(1:3)
1291 END DO
1292
1293 ! move the center of mass of molecule 1 to the origin
1294 start_atom = 1
1295 end_atom = nunits(mol_type(1))
1296 CALL get_center_of_mass(r_old(:, start_atom:end_atom), nunits(mol_type(1)), &
1297 center_of_mass(:), mass(1:nunits(mol_type(1)), mol_type(1)))
1298 DO iunit = start_atom, end_atom
1299 r_old(:, iunit) = r_old(:, iunit) - center_of_mass(:)
1300 END DO
1301 ! set them in the force_env, so the first molecule is ready for the energy calc
1302 DO iparticle = start_atom, end_atom
1303 particles(1)%list%els(iparticle)%r(1:3) = r_old(1:3, iparticle)
1304 END DO
1305
1306 ! print out a notice every 1%
1307 iprint = floor(real(nvirial, kind=dp)/100.0_dp)
1308 IF (iprint == 0) iprint = 1
1309
1310 ! we'll compute the average potential, and then integrate that, as opposed to
1311 ! integrating every orientation and then averaging
1312 ALLOCATE (mayer(1:nvirial_temps, 1:nbins))
1313
1314 mayer(:, :) = 0.0_dp
1315
1316 ! loop over all nvirial random configurations
1317 DO ivirial = 1, nvirial
1318
1319 ! move molecule two back to the origin
1320 start_atom = nunits(mol_type(1)) + 1
1321 end_atom = nunits_tot(1)
1322 CALL get_center_of_mass(r_old(:, start_atom:end_atom), nunits(mol_type(2)), &
1323 center_of_mass(:), mass(1:nunits(mol_type(2)), mol_type(2)))
1324 DO iunit = start_atom, end_atom
1325 r_old(:, iunit) = r_old(:, iunit) - center_of_mass(:)
1326 END DO
1327
1328 ! now we need a random orientation for molecule 2...this routine is
1329 ! only done in serial since it calls a random number
1330 IF (ionode) THEN
1331 CALL rotate_molecule(r_old(:, start_atom:end_atom), &
1332 mass(1:nunits(mol_type(2)), mol_type(2)), &
1333 nunits(mol_type(2)), rng_stream)
1334 END IF
1335 CALL group%bcast(r_old(:, :), source)
1336
1337 distance = 0.0e0_dp
1338 ibin = 1
1339 DO
1340 ! find out what our stepsize is
1341 current_division = 0
1342 DO idivision = 1, nintegral_divisions
1343 IF (distance < virial_cutoffs(idivision) - virial_stepsize(idivision)/2.0e0_dp) THEN
1344 current_division = idivision
1345 EXIT
1346 END IF
1347 END DO
1348 IF (current_division == 0) EXIT
1349 distance = distance + virial_stepsize(current_division)
1350
1351 ! move the second molecule only along the x direction
1352 DO iparticle = start_atom, end_atom
1353 particles(1)%list%els(iparticle)%r(1) = r_old(1, iparticle) + distance
1354 particles(1)%list%els(iparticle)%r(2) = r_old(2, iparticle)
1355 particles(1)%list%els(iparticle)%r(3) = r_old(3, iparticle)
1356 END DO
1357
1358 ! check for overlaps
1359 CALL check_for_overlap(force_env(1)%force_env, nchains(:, 1), nunits, loverlap, mol_type)
1360
1361 ! compute the energy if there is no overlap
1362 ! exponent is exp(-beta*energy)-1, also called the Mayer term
1363 IF (loverlap) THEN
1364 DO itemp = 1, nvirial_temps
1365 mayer(itemp, ibin) = mayer(itemp, ibin) - 1.0_dp
1366 END DO
1367 ELSE
1368 CALL force_env_calc_energy_force(force_env(1)%force_env, &
1369 calc_force=.false.)
1370 CALL force_env_get(force_env(1)%force_env, &
1371 potential_energy=trial_energy)
1372
1373 DO itemp = 1, nvirial_temps
1374
1375 exponent = -beta(itemp)*trial_energy
1376
1377 IF (exponent > exp_max_val) THEN
1378 exponent = exp_max_val
1379 ELSE IF (exponent < exp_min_val) THEN
1380 exponent = exp_min_val
1381 END IF
1382 mayer(itemp, ibin) = mayer(itemp, ibin) + exp(exponent) - 1.0_dp
1383 END DO
1384 END IF
1385
1386 ibin = ibin + 1
1387 END DO
1388 ! write out some info that keeps track of where we are
1389 IF (iw > 0) THEN
1390 IF (mod(ivirial, iprint) == 0) THEN
1391 WRITE (iw, '(A,I6,A,I6)') ' Done with config ', ivirial, ' out of ', nvirial
1392 END IF
1393 END IF
1394 END DO
1395
1396 ! now we integrate this average potential
1397 mayer(:, :) = mayer(:, :)/real(nvirial, dp)
1398
1399 DO itemp = 1, nvirial_temps
1400 integral = 0.0_dp
1401 previous_value = 0.0_dp
1402 distance = 0.0e0_dp
1403 ibin = 1
1404 DO
1405 current_division = 0
1406 DO idivision = 1, nintegral_divisions
1407 IF (distance < virial_cutoffs(idivision) - virial_stepsize(idivision)/2.0e0_dp) THEN
1408 current_division = idivision
1409 EXIT
1410 END IF
1411 END DO
1412 IF (current_division == 0) EXIT
1413 distance = distance + virial_stepsize(current_division)
1414
1415 ! now we need to integrate, using the trapazoidal method
1416 ! first, find the value of the square
1417 current_value = mayer(itemp, ibin)*distance**2
1418 square_value = previous_value*virial_stepsize(current_division)
1419 ! now the triangle that sits on top of it, which is half the size of this square...
1420 ! notice this is negative if the current value is less than the previous value
1421 triangle_value = 0.5e0_dp*((current_value - previous_value)*virial_stepsize(current_division))
1422
1423 integral = integral + square_value + triangle_value
1424 previous_value = current_value
1425 ibin = ibin + 1
1426 END DO
1427
1428 ! now that the integration is done, compute the second virial that results
1429 ave_virial = -2.0e0_dp*pi*integral
1430
1431 ! convert from CP2K units to something else
1432 ave_virial = ave_virial*n_avogadro*angstrom**3/1.0e8_dp**3
1433
1434 IF (iw > 0) THEN
1435 WRITE (iw, *)
1436 WRITE (iw, *) '*********************************************************************'
1437 WRITE (iw, '(A,F12.6,A)') ' *** Temperature = ', virial_temps(itemp), &
1438 ' ***'
1439 WRITE (iw, *) '*** ***'
1440 WRITE (iw, '(A,E12.6,A)') ' *** B2(T) = ', ave_virial, &
1441 ' cm**3/mol ***'
1442 WRITE (iw, *) '*********************************************************************'
1443 WRITE (iw, *)
1444 END IF
1445 END DO
1446
1447 ! deallocate some stuff
1448 DEALLOCATE (mc_par)
1449 DEALLOCATE (subsys)
1450 DEALLOCATE (force_env)
1451 DEALLOCATE (particles)
1452 DEALLOCATE (cell)
1453 DEALLOCATE (virial_cutoffs)
1454 DEALLOCATE (virial_stepsize)
1455 DEALLOCATE (r_old)
1456 DEALLOCATE (mayer)
1457 DEALLOCATE (beta)
1458
1459 END SUBROUTINE mc_compute_virial
1460
1461END MODULE mc_ensembles
1462
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:210
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:91
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:311
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:122
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:523
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, cell_ref, use_ref_cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:364
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:1251
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:218
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, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:344
force_env_types::force_env_release
recursive subroutine, public force_env_release(force_env)
releases the given force env
Definition force_env_types.F:207
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:18
input_constants::dump_xmol
integer, parameter, public dump_xmol
Definition input_constants.F:386
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:1052
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:124
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:399
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:717
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:628
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:868
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:1182
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:1137
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:431
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:1951
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:2488
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:2350
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:226
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:75
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:123
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:174
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:127
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:158
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:746
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