20 mc_molecule_info_type,&
25 #include "../../base/base_uses.f90"
33 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'mc_move_control'
51 TYPE(mc_moves_type),
POINTER :: moves
53 CHARACTER(len=*),
PARAMETER :: routinen =
'init_mc_moves'
59 CALL timeset(routinen, handle)
64 ALLOCATE (moves%angle)
65 ALLOCATE (moves%dihedral)
66 ALLOCATE (moves%trans)
67 ALLOCATE (moves%cltrans)
69 ALLOCATE (moves%bias_bond)
70 ALLOCATE (moves%bias_angle)
71 ALLOCATE (moves%bias_dihedral)
72 ALLOCATE (moves%bias_trans)
73 ALLOCATE (moves%bias_cltrans)
74 ALLOCATE (moves%bias_rot)
75 ALLOCATE (moves%volume)
77 ALLOCATE (moves%avbmc_inin)
78 ALLOCATE (moves%avbmc_inout)
79 ALLOCATE (moves%avbmc_outin)
80 ALLOCATE (moves%avbmc_outout)
82 ALLOCATE (moves%Quickstep)
98 TYPE(mc_moves_type),
POINTER :: moves
100 CHARACTER(len=*),
PARAMETER :: routinen =
'mc_moves_release'
106 CALL timeset(routinen, handle)
109 DEALLOCATE (moves%bond)
110 DEALLOCATE (moves%angle)
111 DEALLOCATE (moves%dihedral)
112 DEALLOCATE (moves%trans)
113 DEALLOCATE (moves%cltrans)
114 DEALLOCATE (moves%rot)
115 DEALLOCATE (moves%bias_bond)
116 DEALLOCATE (moves%bias_angle)
117 DEALLOCATE (moves%bias_dihedral)
118 DEALLOCATE (moves%bias_trans)
119 DEALLOCATE (moves%bias_cltrans)
120 DEALLOCATE (moves%bias_rot)
121 DEALLOCATE (moves%volume)
122 DEALLOCATE (moves%hmc)
123 DEALLOCATE (moves%avbmc_inin)
124 DEALLOCATE (moves%avbmc_inout)
125 DEALLOCATE (moves%avbmc_outin)
126 DEALLOCATE (moves%avbmc_outout)
127 DEALLOCATE (moves%swap)
128 DEALLOCATE (moves%Quickstep)
136 CALL timestop(handle)
151 TYPE(mc_moves_type),
POINTER :: moves
152 LOGICAL,
INTENT(IN) :: lbias
154 CHARACTER(len=*),
PARAMETER :: routinen =
'move_q_reinit'
160 CALL timeset(routinen, handle)
164 moves%bias_bond%qsuccesses = 0
165 moves%bias_angle%qsuccesses = 0
166 moves%bias_dihedral%qsuccesses = 0
167 moves%bias_trans%qsuccesses = 0
168 moves%bias_cltrans%qsuccesses = 0
169 moves%bias_rot%qsuccesses = 0
171 moves%bond%qsuccesses = 0
172 moves%angle%qsuccesses = 0
173 moves%dihedral%qsuccesses = 0
174 moves%trans%qsuccesses = 0
175 moves%cltrans%qsuccesses = 0
176 moves%rot%qsuccesses = 0
177 moves%volume%qsuccesses = 0
178 moves%hmc%qsuccesses = 0
179 moves%qtrans_dis = 0.0e0_dp
180 moves%qcltrans_dis = 0.0e0_dp
184 CALL timestop(handle)
201 TYPE(mc_moves_type),
POINTER :: moves
202 LOGICAL,
INTENT(IN) :: lbias
204 CHARACTER(len=*),
PARAMETER :: routinen =
'q_move_accept'
210 CALL timeset(routinen, handle)
214 moves%bias_bond%successes = moves%bias_bond%successes &
215 + moves%bias_bond%qsuccesses
216 moves%bias_angle%successes = moves%bias_angle%successes &
217 + moves%bias_angle%qsuccesses
218 moves%bias_dihedral%successes = moves%bias_dihedral%successes &
219 + moves%bias_dihedral%qsuccesses
220 moves%bias_trans%successes = moves%bias_trans%successes &
221 + moves%bias_trans%qsuccesses
222 moves%bias_cltrans%successes = moves%bias_cltrans%successes &
223 + moves%bias_cltrans%qsuccesses
224 moves%bias_rot%successes = moves%bias_rot%successes &
225 + moves%bias_rot%qsuccesses
228 moves%bond%successes = moves%bond%successes &
229 + moves%bond%qsuccesses
230 moves%angle%successes = moves%angle%successes &
231 + moves%angle%qsuccesses
232 moves%dihedral%successes = moves%dihedral%successes &
233 + moves%dihedral%qsuccesses
234 moves%trans%successes = moves%trans%successes &
235 + moves%trans%qsuccesses
236 moves%cltrans%successes = moves%cltrans%successes &
237 + moves%cltrans%qsuccesses
238 moves%rot%successes = moves%rot%successes &
239 + moves%rot%qsuccesses
240 moves%hmc%successes = moves%hmc%successes &
241 + moves%hmc%qsuccesses
242 moves%volume%successes = moves%volume%successes &
243 + moves%volume%qsuccesses
244 moves%avbmc_inin%successes = moves%avbmc_inin%successes &
245 + moves%avbmc_inin%qsuccesses
246 moves%avbmc_inout%successes = moves%avbmc_inout%successes &
247 + moves%avbmc_inout%qsuccesses
248 moves%avbmc_outin%successes = moves%avbmc_outin%successes &
249 + moves%avbmc_outin%qsuccesses
250 moves%avbmc_outout%successes = moves%avbmc_outout%successes &
251 + moves%avbmc_outout%qsuccesses
253 moves%trans_dis = moves%trans_dis + moves%qtrans_dis
254 moves%cltrans_dis = moves%cltrans_dis + moves%qcltrans_dis
258 CALL timestop(handle)
274 TYPE(mc_moves_type),
POINTER :: moves
275 INTEGER,
INTENT(IN) :: nnstep, unit
277 CHARACTER(len=*),
PARAMETER :: routinen =
'write_move_stats'
283 CALL timeset(routinen, handle)
285 WRITE (unit, 1000) nnstep,
' bias_bond ', &
286 moves%bias_bond%successes, moves%bias_bond%attempts
287 WRITE (unit, 1000) nnstep,
' bias_angle ', &
288 moves%bias_angle%successes, moves%bias_angle%attempts
289 WRITE (unit, 1000) nnstep,
' bias_dihedral ', &
290 moves%bias_dihedral%successes, moves%bias_dihedral%attempts
291 WRITE (unit, 1000) nnstep,
' bias_trans ', &
292 moves%bias_trans%successes, moves%bias_trans%attempts
293 WRITE (unit, 1000) nnstep,
' bias_cltrans ', &
294 moves%bias_cltrans%successes, moves%bias_cltrans%attempts
295 WRITE (unit, 1000) nnstep,
' bias_rot ', &
296 moves%bias_rot%successes, moves%bias_rot%attempts
298 WRITE (unit, 1000) nnstep,
' bond ', &
299 moves%bond%successes, moves%bond%attempts
300 WRITE (unit, 1000) nnstep,
' angle ', &
301 moves%angle%successes, moves%angle%attempts
302 WRITE (unit, 1000) nnstep,
' dihedral ', &
303 moves%dihedral%successes, moves%dihedral%attempts
304 WRITE (unit, 1000) nnstep,
' trans ', &
305 moves%trans%successes, moves%trans%attempts
306 WRITE (unit, 1000) nnstep,
' cltrans ', &
307 moves%cltrans%successes, moves%cltrans%attempts
308 WRITE (unit, 1000) nnstep,
' rot ', &
309 moves%rot%successes, moves%rot%attempts
310 WRITE (unit, 1000) nnstep,
' swap ', &
311 moves%swap%successes, moves%swap%attempts
312 WRITE (unit, 1001) nnstep,
' grown ', &
314 WRITE (unit, 1001) nnstep,
' empty_swap ', &
316 WRITE (unit, 1001) nnstep,
' empty_conf ', &
318 WRITE (unit, 1000) nnstep,
' volume ', &
319 moves%volume%successes, moves%volume%attempts
320 WRITE (unit, 1000) nnstep,
' HMC ', &
321 moves%hmc%successes, moves%hmc%attempts
322 WRITE (unit, 1000) nnstep,
' avbmc_inin ', &
323 moves%avbmc_inin%successes, moves%avbmc_inin%attempts
324 WRITE (unit, 1000) nnstep,
' avbmc_inout ', &
325 moves%avbmc_inout%successes, moves%avbmc_inout%attempts
326 WRITE (unit, 1000) nnstep,
' avbmc_outin ', &
327 moves%avbmc_outin%successes, moves%avbmc_outin%attempts
328 WRITE (unit, 1000) nnstep,
' avbmc_outout ', &
329 moves%avbmc_outout%successes, moves%avbmc_outout%attempts
330 WRITE (unit, 1001) nnstep,
' empty_avbmc ', &
332 WRITE (unit, 1000) nnstep,
' Quickstep ', &
333 moves%quickstep%successes, moves%quickstep%attempts
335 1000
FORMAT(i10, 2x, a, 2x, i10, 2x, i10)
336 1001
FORMAT(i10, 2x, a, 2x, i10)
338 CALL timestop(handle)
361 TYPE(mc_simpar_type),
POINTER :: mc_par
362 TYPE(mc_moves_type),
POINTER :: move_updates
363 INTEGER,
INTENT(IN) :: molecule_type
364 CHARACTER(LEN=*),
INTENT(IN) :: flag
365 INTEGER,
INTENT(IN) :: nnstep
366 LOGICAL,
INTENT(IN) :: ionode
368 CHARACTER(len=*),
PARAMETER :: routinen =
'mc_move_update'
370 INTEGER :: handle, nmol_types, rm
371 REAL(
dp),
DIMENSION(:),
POINTER :: rmangle, rmbond, rmdihedral, rmrot, &
373 REAL(kind=
dp) :: rmcltrans, rmvolume, test_ratio
374 TYPE(mc_molecule_info_type),
POINTER :: mc_molecule_info
378 CALL timeset(routinen, handle)
380 NULLIFY (rmangle, rmbond, rmdihedral, rmrot, rmtrans)
383 CALL get_mc_par(mc_par, rmbond=rmbond, rmangle=rmangle, rmrot=rmrot, &
384 rmtrans=rmtrans, rmcltrans=rmcltrans, rmvolume=rmvolume, rm=rm, rmdihedral=rmdihedral, &
385 mc_molecule_info=mc_molecule_info)
390 WRITE (*, *)
'flag =', flag
391 cpabort(
"Wrong option passed")
395 IF (ionode)
WRITE (rm, *) nnstep,
' Data for molecule type ', &
399 IF (move_updates%bias_bond%attempts .GT. 0)
THEN
402 IF (move_updates%bias_bond%successes == 0)
THEN
403 rmbond(molecule_type) = rmbond(molecule_type)/2.0e0_dp
404 ELSEIF (move_updates%bias_bond%successes == &
405 move_updates%bias_bond%attempts)
THEN
406 rmbond(molecule_type) = rmbond(molecule_type)*2.0e0_dp
409 test_ratio = real(move_updates%bias_bond%successes,
dp) &
410 /real(move_updates%bias_bond%attempts,
dp)/0.5e0_dp
411 IF (test_ratio .GT. 2.0e0_dp) test_ratio = 2.0e0_dp
412 IF (test_ratio .LT. 0.5e0_dp) test_ratio = 0.5e0_dp
413 rmbond(molecule_type) = rmbond(molecule_type)*test_ratio
417 move_updates%bias_bond%attempts = 0
418 move_updates%bias_bond%successes = 0
421 IF (ionode)
WRITE (rm, *) nnstep,
' rmbond = ', &
422 rmbond(molecule_type)*
angstrom,
' angstroms'
427 IF (move_updates%bias_angle%attempts .GT. 0)
THEN
430 IF (move_updates%bias_angle%successes == 0)
THEN
431 rmangle(molecule_type) = rmangle(molecule_type)/2.0e0_dp
432 ELSEIF (move_updates%bias_angle%successes == &
433 move_updates%bias_angle%attempts)
THEN
434 rmangle(molecule_type) = rmangle(molecule_type)*2.0e0_dp
437 test_ratio = real(move_updates%bias_angle%successes,
dp) &
438 /real(move_updates%bias_angle%attempts,
dp)/0.5e0_dp
439 IF (test_ratio .GT. 2.0e0_dp) test_ratio = 2.0e0_dp
440 IF (test_ratio .LT. 0.5e0_dp) test_ratio = 0.5e0_dp
441 rmangle(molecule_type) = rmangle(molecule_type)*test_ratio
445 IF (rmangle(molecule_type) .GT.
pi) rmangle(molecule_type) =
pi
448 move_updates%bias_angle%attempts = 0
449 move_updates%bias_angle%successes = 0
452 IF (ionode)
WRITE (rm, *) nnstep,
' rmangle = ', &
453 rmangle(molecule_type)/
pi*180.0e0_dp,
' degrees'
457 IF (move_updates%bias_dihedral%attempts .GT. 0)
THEN
460 IF (move_updates%bias_dihedral%successes == 0)
THEN
461 rmdihedral(molecule_type) = rmdihedral(molecule_type)/2.0e0_dp
462 ELSEIF (move_updates%bias_dihedral%successes == &
463 move_updates%bias_dihedral%attempts)
THEN
464 rmdihedral(molecule_type) = rmdihedral(molecule_type)*2.0e0_dp
467 test_ratio = real(move_updates%bias_dihedral%successes,
dp) &
468 /real(move_updates%bias_dihedral%attempts,
dp)/0.5e0_dp
469 IF (test_ratio .GT. 2.0e0_dp) test_ratio = 2.0e0_dp
470 IF (test_ratio .LT. 0.5e0_dp) test_ratio = 0.5e0_dp
471 rmdihedral(molecule_type) = rmdihedral(molecule_type)*test_ratio
475 IF (rmdihedral(molecule_type) .GT.
pi) rmdihedral(molecule_type) =
pi
478 move_updates%bias_dihedral%attempts = 0
479 move_updates%bias_dihedral%successes = 0
482 IF (ionode)
WRITE (rm, *) nnstep,
' rmdihedral = ', &
483 rmdihedral(molecule_type)/
pi*180.0e0_dp,
' degrees'
487 IF (move_updates%bias_trans%attempts .GT. 0)
THEN
490 IF (move_updates%bias_trans%successes == 0)
THEN
491 rmtrans(molecule_type) = rmtrans(molecule_type)/2.0e0_dp
492 ELSEIF (move_updates%bias_trans%successes == &
493 move_updates%bias_trans%attempts)
THEN
494 rmtrans(molecule_type) = rmtrans(molecule_type)*2.0e0_dp
497 test_ratio = real(move_updates%bias_trans%successes,
dp) &
498 /real(move_updates%bias_trans%attempts,
dp)/0.5e0_dp
499 IF (test_ratio .GT. 2.0e0_dp) test_ratio = 2.0e0_dp
500 IF (test_ratio .LT. 0.5e0_dp) test_ratio = 0.5e0_dp
501 rmtrans(molecule_type) = rmtrans(molecule_type)*test_ratio
505 IF (rmtrans(molecule_type) .GT. 10.0e0_dp) &
506 rmtrans(molecule_type) = 10.0e0_dp
509 move_updates%bias_trans%attempts = 0
510 move_updates%bias_trans%successes = 0
513 IF (ionode)
WRITE (rm, *) nnstep,
' rmtrans = ', &
514 rmtrans(molecule_type)*
angstrom,
' angstroms'
518 IF (move_updates%bias_cltrans%attempts .GT. 0)
THEN
521 IF (move_updates%bias_cltrans%successes == 0)
THEN
522 rmcltrans = rmcltrans/2.0e0_dp
523 ELSEIF (move_updates%bias_cltrans%successes == &
524 move_updates%bias_cltrans%attempts)
THEN
525 rmcltrans = rmcltrans*2.0e0_dp
528 test_ratio = real(move_updates%bias_cltrans%successes,
dp) &
529 /real(move_updates%bias_cltrans%attempts,
dp)/0.5e0_dp
530 IF (test_ratio .GT. 2.0e0_dp) test_ratio = 2.0e0_dp
531 IF (test_ratio .LT. 0.5e0_dp) test_ratio = 0.5e0_dp
532 rmcltrans = rmcltrans*test_ratio
536 IF (rmcltrans .GT. 10.0e0_dp) &
537 rmcltrans = 10.0e0_dp
540 move_updates%bias_cltrans%attempts = 0
541 move_updates%bias_cltrans%successes = 0
544 IF (ionode)
WRITE (rm, *) nnstep,
' rmcltrans = ', &
549 IF (move_updates%bias_rot%attempts .GT. 0)
THEN
552 IF (move_updates%bias_rot%successes == 0)
THEN
553 rmrot = rmrot/2.0e0_dp
555 IF (rmrot(molecule_type) .GT.
pi) rmrot(molecule_type) =
pi
557 ELSEIF (move_updates%bias_rot%successes == &
558 move_updates%bias_rot%attempts)
THEN
559 rmrot(molecule_type) = rmrot(molecule_type)*2.0e0_dp
562 IF (rmrot(molecule_type) .GT.
pi) rmrot(molecule_type) =
pi
566 test_ratio = real(move_updates%bias_rot%successes,
dp) &
567 /real(move_updates%bias_rot%attempts,
dp)/0.5e0_dp
568 IF (test_ratio .GT. 2.0e0_dp) test_ratio = 2.0e0_dp
569 IF (test_ratio .LT. 0.5e0_dp) test_ratio = 0.5e0_dp
570 rmrot(molecule_type) = rmrot(molecule_type)*test_ratio
573 IF (rmrot(molecule_type) .GT.
pi) rmrot(molecule_type) =
pi
578 move_updates%bias_rot%attempts = 0
579 move_updates%bias_rot%successes = 0
582 IF (ionode)
WRITE (rm, *) nnstep,
' rmrot = ', &
583 rmrot(molecule_type)/
pi*180.0e0_dp,
' degrees'
589 IF (move_updates%volume%attempts .NE. 0)
THEN
592 IF (move_updates%volume%successes == 0)
THEN
593 rmvolume = rmvolume/2.0e0_dp
595 ELSEIF (move_updates%volume%successes == &
596 move_updates%volume%attempts)
THEN
597 rmvolume = rmvolume*2.0e0_dp
600 test_ratio = real(move_updates%volume%successes,
dp)/ &
601 REAL(move_updates%volume%attempts,
dp)/0.5e0_dp
602 IF (test_ratio .GT. 2.0e0_dp) test_ratio = 2.0e0_dp
603 IF (test_ratio .LT. 0.5e0_dp) test_ratio = 0.5e0_dp
604 rmvolume = rmvolume*test_ratio
609 move_updates%volume%attempts = 0
610 move_updates%volume%successes = 0
613 IF (ionode)
WRITE (rm, *) nnstep,
' rmvolume = ', &
614 rmvolume*
angstrom**3,
' angstroms^3'
621 CALL set_mc_par(mc_par, rmbond=rmbond, rmangle=rmangle, rmrot=rmrot, &
622 rmtrans=rmtrans, rmcltrans=rmcltrans, rmvolume=rmvolume, rmdihedral=rmdihedral)
625 CALL timestop(handle)
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
control the handling of the move data in Monte Carlo (MC) simulations
subroutine, public mc_moves_release(moves)
deallocates all the structures and nullifies the pointer
subroutine, public q_move_accept(moves, lbias)
updates accepted moves in the given structure...assumes you've been recording all successful moves in...
subroutine, public write_move_stats(moves, nnstep, unit)
writes the number of accepted and attempted moves to a file for the various move types
subroutine, public move_q_reinit(moves, lbias)
sets all qsuccess counters back to zero
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...
subroutine, public init_mc_moves(moves)
allocates and initializes the structure to record all move attempts/successes
holds all the structure types needed for Monte Carlo, except the mc_environment_type
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)
...
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
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 of physical constants:
real(kind=dp), parameter, public angstrom