(git:71c3ab0)
Loading...
Searching...
No Matches
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-2026 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! **************************************************************************************************
14MODULE mc_move_control
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,&
23 set_mc_par
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, &
36 mc_move_update, move_q_reinit, q_move_accept, mc_moves_release, &
37 write_move_stats
38
39CONTAINS
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
3351000 FORMAT(i10, 2x, a, 2x, i10, 2x, i10)
3361001 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 > 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 ELSE IF (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 > 2.0e0_dp) test_ratio = 2.0e0_dp
412 IF (test_ratio < 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 > 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 ELSE IF (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 > 2.0e0_dp) test_ratio = 2.0e0_dp
440 IF (test_ratio < 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) > 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 > 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 ELSE IF (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 > 2.0e0_dp) test_ratio = 2.0e0_dp
470 IF (test_ratio < 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) > 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 > 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 ELSE IF (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 > 2.0e0_dp) test_ratio = 2.0e0_dp
500 IF (test_ratio < 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) > 10.0e0_dp) THEN
506 rmtrans(molecule_type) = 10.0e0_dp
507 END IF
508
509 ! clear the counters
510 move_updates%bias_trans%attempts = 0
511 move_updates%bias_trans%successes = 0
512
513! write the new displacement to a file
514 IF (ionode) WRITE (rm, *) nnstep, ' rmtrans = ', &
515 rmtrans(molecule_type)*angstrom, ' angstroms'
516 END IF
517
518! update the maximum displacement for cluster translation
519 IF (move_updates%bias_cltrans%attempts > 0) THEN
520
521! first account for the extreme cases
522 IF (move_updates%bias_cltrans%successes == 0) THEN
523 rmcltrans = rmcltrans/2.0e0_dp
524 ELSE IF (move_updates%bias_cltrans%successes == &
525 move_updates%bias_cltrans%attempts) THEN
526 rmcltrans = rmcltrans*2.0e0_dp
527 ELSE
528! now for the middle case
529 test_ratio = real(move_updates%bias_cltrans%successes, dp) &
530 /real(move_updates%bias_cltrans%attempts, dp)/0.5e0_dp
531 IF (test_ratio > 2.0e0_dp) test_ratio = 2.0e0_dp
532 IF (test_ratio < 0.5e0_dp) test_ratio = 0.5e0_dp
533 rmcltrans = rmcltrans*test_ratio
534 END IF
535
536 ! make an upper bound...10 a.u.
537 IF (rmcltrans > 10.0e0_dp) THEN
538 rmcltrans = 10.0e0_dp
539 END IF
540
541 ! clear the counters
542 move_updates%bias_cltrans%attempts = 0
543 move_updates%bias_cltrans%successes = 0
544
545! write the new displacement to a file
546 IF (ionode) WRITE (rm, *) nnstep, ' rmcltrans = ', &
547 rmcltrans*angstrom, ' angstroms'
548 END IF
549
550! update the maximum displacement for molecule rotation
551 IF (move_updates%bias_rot%attempts > 0) THEN
552
553! first account for the extreme cases
554 IF (move_updates%bias_rot%successes == 0) THEN
555 rmrot = rmrot/2.0e0_dp
556
557 IF (rmrot(molecule_type) > pi) rmrot(molecule_type) = pi
558
559 ELSE IF (move_updates%bias_rot%successes == &
560 move_updates%bias_rot%attempts) THEN
561 rmrot(molecule_type) = rmrot(molecule_type)*2.0e0_dp
562
563! more than pi rotation is meaningless
564 IF (rmrot(molecule_type) > pi) rmrot(molecule_type) = pi
565
566 ELSE
567! now for the middle case
568 test_ratio = real(move_updates%bias_rot%successes, dp) &
569 /real(move_updates%bias_rot%attempts, dp)/0.5e0_dp
570 IF (test_ratio > 2.0e0_dp) test_ratio = 2.0e0_dp
571 IF (test_ratio < 0.5e0_dp) test_ratio = 0.5e0_dp
572 rmrot(molecule_type) = rmrot(molecule_type)*test_ratio
573
574! more than pi rotation is meaningless
575 IF (rmrot(molecule_type) > pi) rmrot(molecule_type) = pi
576
577 END IF
578
579! clear the counters
580 move_updates%bias_rot%attempts = 0
581 move_updates%bias_rot%successes = 0
582
583! write the new displacement to a file
584 IF (ionode) WRITE (rm, *) nnstep, ' rmrot = ', &
585 rmrot(molecule_type)/pi*180.0e0_dp, ' degrees'
586 END IF
587
588 CASE ("volume")
589
590! update the maximum displacement for volume displacement
591 IF (move_updates%volume%attempts /= 0) THEN
592
593! first account for the extreme cases
594 IF (move_updates%volume%successes == 0) THEN
595 rmvolume = rmvolume/2.0e0_dp
596
597 ELSE IF (move_updates%volume%successes == &
598 move_updates%volume%attempts) THEN
599 rmvolume = rmvolume*2.0e0_dp
600 ELSE
601! now for the middle case
602 test_ratio = real(move_updates%volume%successes, dp)/ &
603 REAL(move_updates%volume%attempts, dp)/0.5e0_dp
604 IF (test_ratio > 2.0e0_dp) test_ratio = 2.0e0_dp
605 IF (test_ratio < 0.5e0_dp) test_ratio = 0.5e0_dp
606 rmvolume = rmvolume*test_ratio
607
608 END IF
609
610! clear the counters
611 move_updates%volume%attempts = 0
612 move_updates%volume%successes = 0
613
614! write the new displacement to a file
615 IF (ionode) WRITE (rm, *) nnstep, ' rmvolume = ', &
616 rmvolume*angstrom**3, ' angstroms^3'
617
618 END IF
619
620 END SELECT
621
622! set some of the MC parameters
623 CALL set_mc_par(mc_par, rmbond=rmbond, rmangle=rmangle, rmrot=rmrot, &
624 rmtrans=rmtrans, rmcltrans=rmcltrans, rmvolume=rmvolume, rmdihedral=rmdihedral)
625
626! end the timing
627 CALL timestop(handle)
628
629 END SUBROUTINE mc_move_update
630
631END MODULE mc_move_control
632
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
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::q_move_accept
subroutine, public q_move_accept(moves, lbias)
updates accepted moves in the given structure...assumes you've been recording all successful moves in...
Definition mc_move_control.F:200
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::move_q_reinit
subroutine, public move_q_reinit(moves, lbias)
sets all qsuccess counters back to zero
Definition mc_move_control.F:150
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_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
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
mc_types::mc_molecule_info_type
Definition mc_types.F:169
mc_types::mc_moves_type
Definition mc_types.F:204
mc_types::mc_simpar_type
Definition mc_types.F:81