(git:58e3e09)
mc_move_control.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 control the handling of the move data in Monte Carlo (MC) simulations
10 !> \par History
11 !> none
12 !> \author Matthew J. McGrath (10.16.2003)
13 ! **************************************************************************************************
15 
16  USE kinds, ONLY: dp
17  USE mathconstants, ONLY: pi
18  USE mc_types, ONLY: get_mc_molecule_info,&
19  get_mc_par,&
20  mc_molecule_info_type,&
21  mc_moves_type,&
22  mc_simpar_type,&
24  USE physcon, ONLY: angstrom
25 #include "../../base/base_uses.f90"
26 
27  IMPLICIT NONE
28 
29  PRIVATE
30 
31 ! *** Global parameters ***
32 
33  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_move_control'
34 
35  PUBLIC :: init_mc_moves, &
38 
39 CONTAINS
40 
41 ! **************************************************************************************************
42 !> \brief allocates and initializes the structure to record all move
43 !> attempts/successes
44 !> \param moves the move structure to update
45 !>
46 !> Suitable for parallel.
47 !> \author MJM
48 ! **************************************************************************************************
49  SUBROUTINE init_mc_moves(moves)
50 
51  TYPE(mc_moves_type), POINTER :: moves
52 
53  CHARACTER(len=*), PARAMETER :: routinen = 'init_mc_moves'
54 
55  INTEGER :: handle
56 
57 ! begin the timing of the subroutine
58 
59  CALL timeset(routinen, handle)
60 
61 ! allocate all the structures
62  ALLOCATE (moves)
63  ALLOCATE (moves%bond)
64  ALLOCATE (moves%angle)
65  ALLOCATE (moves%dihedral)
66  ALLOCATE (moves%trans)
67  ALLOCATE (moves%cltrans)
68  ALLOCATE (moves%rot)
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)
76  ALLOCATE (moves%hmc)
77  ALLOCATE (moves%avbmc_inin)
78  ALLOCATE (moves%avbmc_inout)
79  ALLOCATE (moves%avbmc_outin)
80  ALLOCATE (moves%avbmc_outout)
81  ALLOCATE (moves%swap)
82  ALLOCATE (moves%Quickstep)
83 
84  ! end the timing
85  CALL timestop(handle)
86 
87  END SUBROUTINE init_mc_moves
88 
89 ! **************************************************************************************************
90 !> \brief deallocates all the structures and nullifies the pointer
91 !> \param moves the move structure to release
92 !>
93 !> Suitable for parallel.
94 !> \author MJM
95 ! **************************************************************************************************
96  SUBROUTINE mc_moves_release(moves)
97 
98  TYPE(mc_moves_type), POINTER :: moves
99 
100  CHARACTER(len=*), PARAMETER :: routinen = 'mc_moves_release'
101 
102  INTEGER :: handle
103 
104 ! begin the timing of the subroutine
105 
106  CALL timeset(routinen, handle)
107 
108 ! allocate all the structures
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)
129 
130  DEALLOCATE (moves)
131 
132 ! now nullify the moves
133  NULLIFY (moves)
134 
135  ! end the timing
136  CALL timestop(handle)
137 
138  END SUBROUTINE mc_moves_release
139 
140 ! **************************************************************************************************
141 !> \brief sets all qsuccess counters back to zero
142 !> \param moves the move structure to update
143 !> \param lbias are we biasing translations/rotations/conformational changes
144 !> with a different potential?
145 !>
146 !> Suitable for parallel.
147 !> \author MJM
148 ! **************************************************************************************************
149  SUBROUTINE move_q_reinit(moves, lbias)
150 
151  TYPE(mc_moves_type), POINTER :: moves
152  LOGICAL, INTENT(IN) :: lbias
153 
154  CHARACTER(len=*), PARAMETER :: routinen = 'move_q_reinit'
155 
156  INTEGER :: handle
157 
158 ! begin the timing of the subroutine
159 
160  CALL timeset(routinen, handle)
161 
162 ! set all the counters equal to zero
163  IF (lbias) THEN
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
170  ELSE
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
181  END IF
182 
183  ! end the timing
184  CALL timestop(handle)
185 
186  END SUBROUTINE move_q_reinit
187 
188 ! **************************************************************************************************
189 !> \brief updates accepted moves in the given structure...assumes you've been
190 !> recording all successful moves in "qsuccesses"...this was done to
191 !> compensate for doing multiple inner moves between Quickstep moves
192 !> (which determine ultimate acceptance of moves)
193 !> \param moves the move structure to update
194 !> \param lbias are we biasing non-swap particle moves with a cheaper potential
195 !>
196 !> Suitable for parallel.
197 !> \author MJM
198 ! **************************************************************************************************
199  SUBROUTINE q_move_accept(moves, lbias)
200 
201  TYPE(mc_moves_type), POINTER :: moves
202  LOGICAL, INTENT(IN) :: lbias
203 
204  CHARACTER(len=*), PARAMETER :: routinen = 'q_move_accept'
205 
206  INTEGER :: handle
207 
208 ! begin the timing of the subroutine
209 
210  CALL timeset(routinen, handle)
211 
212  IF (lbias) THEN
213 ! change the number of successful moves for the total move counter
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
226  ELSE
227 ! change the number of successful moves for the total move counter
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
252 
253  moves%trans_dis = moves%trans_dis + moves%qtrans_dis
254  moves%cltrans_dis = moves%cltrans_dis + moves%qcltrans_dis
255  END IF
256 
257 ! end the timing
258  CALL timestop(handle)
259 
260  END SUBROUTINE q_move_accept
261 
262 ! **************************************************************************************************
263 !> \brief writes the number of accepted and attempted moves to a file for
264 !> the various move types
265 !> \param moves the structure containing the move data
266 !> \param nnstep what step we're on
267 !> \param unit the unit of the file we're writing to
268 !>
269 !> Use only in serial.
270 !> \author MJM
271 ! **************************************************************************************************
272  SUBROUTINE write_move_stats(moves, nnstep, unit)
273 
274  TYPE(mc_moves_type), POINTER :: moves
275  INTEGER, INTENT(IN) :: nnstep, unit
276 
277  CHARACTER(len=*), PARAMETER :: routinen = 'write_move_stats'
278 
279  INTEGER :: handle
280 
281 ! begin the timing of the subroutine
282 
283  CALL timeset(routinen, handle)
284 
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
297 
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 ', &
313  moves%grown
314  WRITE (unit, 1001) nnstep, ' empty_swap ', &
315  moves%empty
316  WRITE (unit, 1001) nnstep, ' empty_conf ', &
317  moves%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 ', &
331  moves%empty_avbmc
332  WRITE (unit, 1000) nnstep, ' Quickstep ', &
333  moves%quickstep%successes, moves%quickstep%attempts
334 
335 1000 FORMAT(i10, 2x, a, 2x, i10, 2x, i10)
336 1001 FORMAT(i10, 2x, a, 2x, i10)
337 ! end the timing
338  CALL timestop(handle)
339 
340  END SUBROUTINE write_move_stats
341 
342 ! **************************************************************************************************
343 !> \brief updates the maximum displacements of a Monte Carlo simulation,
344 !> based on the ratio of successful moves to attempts...tries to hit a
345 !> target of 0.5 acceptance ratio
346 !> \param mc_par the mc parameters for the force env
347 !> \param move_updates holds the accepted/attempted moves since the last
348 !> update (or start of simulation)
349 !> \param molecule_type ...
350 !> \param flag indicates which displacements to update..."volume" is for
351 !> volume moves and "trans" is for everything else
352 !> \param nnstep how many steps the simulation has run
353 !> \param ionode is this the main CPU running this job?
354 !>
355 !> Suitable for parallel.
356 !> \author MJM
357 ! **************************************************************************************************
358  SUBROUTINE mc_move_update(mc_par, move_updates, molecule_type, flag, &
359  nnstep, ionode)
360 
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
367 
368  CHARACTER(len=*), PARAMETER :: routinen = 'mc_move_update'
369 
370  INTEGER :: handle, nmol_types, rm
371  REAL(dp), DIMENSION(:), POINTER :: rmangle, rmbond, rmdihedral, rmrot, &
372  rmtrans
373  REAL(kind=dp) :: rmcltrans, rmvolume, test_ratio
374  TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
375 
376 ! begin the timing of the subroutine
377 
378  CALL timeset(routinen, handle)
379 
380  NULLIFY (rmangle, rmbond, rmdihedral, rmrot, rmtrans)
381 
382 ! grab some stuff from mc_par
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)
386  CALL get_mc_molecule_info(mc_molecule_info, nmol_types=nmol_types)
387 
388  SELECT CASE (flag)
389  CASE DEFAULT
390  WRITE (*, *) 'flag =', flag
391  cpabort("Wrong option passed")
392  CASE ("trans")
393 
394 ! we need to update all the displacements for every molecule type
395  IF (ionode) WRITE (rm, *) nnstep, ' Data for molecule type ', &
396  molecule_type
397 
398 ! update the maximum displacement for bond length change
399  IF (move_updates%bias_bond%attempts .GT. 0) THEN
400 
401 ! first account for the extreme cases
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
407  ELSE
408 ! now for the middle case
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
414  END IF
415 
416 ! update and clear the counters
417  move_updates%bias_bond%attempts = 0
418  move_updates%bias_bond%successes = 0
419 
420 ! write the new displacement to a file
421  IF (ionode) WRITE (rm, *) nnstep, ' rmbond = ', &
422  rmbond(molecule_type)*angstrom, ' angstroms'
423 
424  END IF
425 
426 ! update the maximum displacement for bond angle change
427  IF (move_updates%bias_angle%attempts .GT. 0) THEN
428 
429 ! first account for the extreme cases
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
435  ELSE
436 ! now for the middle case
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
442  END IF
443 
444 ! more than pi changes meaningless
445  IF (rmangle(molecule_type) .GT. pi) rmangle(molecule_type) = pi
446 
447 ! clear the counters
448  move_updates%bias_angle%attempts = 0
449  move_updates%bias_angle%successes = 0
450 
451 ! write the new displacement to a file
452  IF (ionode) WRITE (rm, *) nnstep, ' rmangle = ', &
453  rmangle(molecule_type)/pi*180.0e0_dp, ' degrees'
454  END IF
455 
456 ! update the maximum displacement for a dihedral change
457  IF (move_updates%bias_dihedral%attempts .GT. 0) THEN
458 
459 ! first account for the extreme cases
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
465  ELSE
466 ! now for the middle case
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
472  END IF
473 
474 ! more than pi changes meaningless
475  IF (rmdihedral(molecule_type) .GT. pi) rmdihedral(molecule_type) = pi
476 
477 ! clear the counters
478  move_updates%bias_dihedral%attempts = 0
479  move_updates%bias_dihedral%successes = 0
480 
481 ! write the new displacement to a file
482  IF (ionode) WRITE (rm, *) nnstep, ' rmdihedral = ', &
483  rmdihedral(molecule_type)/pi*180.0e0_dp, ' degrees'
484  END IF
485 
486 ! update the maximum displacement for molecule translation
487  IF (move_updates%bias_trans%attempts .GT. 0) THEN
488 
489 ! first account for the extreme cases
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
495  ELSE
496 ! now for the middle case
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
502  END IF
503 
504  ! make an upper bound...10 a.u.
505  IF (rmtrans(molecule_type) .GT. 10.0e0_dp) &
506  rmtrans(molecule_type) = 10.0e0_dp
507 
508  ! clear the counters
509  move_updates%bias_trans%attempts = 0
510  move_updates%bias_trans%successes = 0
511 
512 ! write the new displacement to a file
513  IF (ionode) WRITE (rm, *) nnstep, ' rmtrans = ', &
514  rmtrans(molecule_type)*angstrom, ' angstroms'
515  END IF
516 
517 ! update the maximum displacement for cluster translation
518  IF (move_updates%bias_cltrans%attempts .GT. 0) THEN
519 
520 ! first account for the extreme cases
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
526  ELSE
527 ! now for the middle case
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
533  END IF
534 
535  ! make an upper bound...10 a.u.
536  IF (rmcltrans .GT. 10.0e0_dp) &
537  rmcltrans = 10.0e0_dp
538 
539  ! clear the counters
540  move_updates%bias_cltrans%attempts = 0
541  move_updates%bias_cltrans%successes = 0
542 
543 ! write the new displacement to a file
544  IF (ionode) WRITE (rm, *) nnstep, ' rmcltrans = ', &
545  rmcltrans*angstrom, ' angstroms'
546  END IF
547 
548 ! update the maximum displacement for molecule rotation
549  IF (move_updates%bias_rot%attempts .GT. 0) THEN
550 
551 ! first account for the extreme cases
552  IF (move_updates%bias_rot%successes == 0) THEN
553  rmrot = rmrot/2.0e0_dp
554 
555  IF (rmrot(molecule_type) .GT. pi) rmrot(molecule_type) = pi
556 
557  ELSEIF (move_updates%bias_rot%successes == &
558  move_updates%bias_rot%attempts) THEN
559  rmrot(molecule_type) = rmrot(molecule_type)*2.0e0_dp
560 
561 ! more than pi rotation is meaningless
562  IF (rmrot(molecule_type) .GT. pi) rmrot(molecule_type) = pi
563 
564  ELSE
565 ! now for the middle case
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
571 
572 ! more than pi rotation is meaningless
573  IF (rmrot(molecule_type) .GT. pi) rmrot(molecule_type) = pi
574 
575  END IF
576 
577 ! clear the counters
578  move_updates%bias_rot%attempts = 0
579  move_updates%bias_rot%successes = 0
580 
581 ! write the new displacement to a file
582  IF (ionode) WRITE (rm, *) nnstep, ' rmrot = ', &
583  rmrot(molecule_type)/pi*180.0e0_dp, ' degrees'
584  END IF
585 
586  CASE ("volume")
587 
588 ! update the maximum displacement for volume displacement
589  IF (move_updates%volume%attempts .NE. 0) THEN
590 
591 ! first account for the extreme cases
592  IF (move_updates%volume%successes == 0) THEN
593  rmvolume = rmvolume/2.0e0_dp
594 
595  ELSEIF (move_updates%volume%successes == &
596  move_updates%volume%attempts) THEN
597  rmvolume = rmvolume*2.0e0_dp
598  ELSE
599 ! now for the middle case
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
605 
606  END IF
607 
608 ! clear the counters
609  move_updates%volume%attempts = 0
610  move_updates%volume%successes = 0
611 
612 ! write the new displacement to a file
613  IF (ionode) WRITE (rm, *) nnstep, ' rmvolume = ', &
614  rmvolume*angstrom**3, ' angstroms^3'
615 
616  END IF
617 
618  END SELECT
619 
620 ! set some of the MC parameters
621  CALL set_mc_par(mc_par, rmbond=rmbond, rmangle=rmangle, rmrot=rmrot, &
622  rmtrans=rmtrans, rmcltrans=rmcltrans, rmvolume=rmvolume, rmdihedral=rmdihedral)
623 
624 ! end the timing
625  CALL timestop(handle)
626 
627  END SUBROUTINE mc_move_update
628 
629 END MODULE mc_move_control
630 
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
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
Definition: mc_types.F:15
subroutine, public get_mc_molecule_info(mc_molecule_info, nmol_types, nchain_total, nboxes, names, conf_prob, nchains, nunits, mol_type, nunits_tot, in_box, atom_names, mass)
...
Definition: mc_types.F:554
subroutine, public 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
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
Definition of physical constants:
Definition: physcon.F:68
real(kind=dp), parameter, public angstrom
Definition: physcon.F:144