(git:6344619)
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 CALL cp_abort(__location__, &
391 "<trans>, <volume> are supported as the <flag> for mc_move_update, "// &
392 "found unknown option <"//trim(flag)//">")
393 CASE ("trans")
394
395! we need to update all the displacements for every molecule type
396 IF (ionode) WRITE (rm, *) nnstep, ' Data for molecule type ', &
397 molecule_type
398
399! update the maximum displacement for bond length change
400 IF (move_updates%bias_bond%attempts > 0) THEN
401
402! first account for the extreme cases
403 IF (move_updates%bias_bond%successes == 0) THEN
404 rmbond(molecule_type) = rmbond(molecule_type)/2.0e0_dp
405 ELSE IF (move_updates%bias_bond%successes == &
406 move_updates%bias_bond%attempts) THEN
407 rmbond(molecule_type) = rmbond(molecule_type)*2.0e0_dp
408 ELSE
409! now for the middle case
410 test_ratio = real(move_updates%bias_bond%successes, dp) &
411 /real(move_updates%bias_bond%attempts, dp)/0.5e0_dp
412 IF (test_ratio > 2.0e0_dp) test_ratio = 2.0e0_dp
413 IF (test_ratio < 0.5e0_dp) test_ratio = 0.5e0_dp
414 rmbond(molecule_type) = rmbond(molecule_type)*test_ratio
415 END IF
416
417! update and clear the counters
418 move_updates%bias_bond%attempts = 0
419 move_updates%bias_bond%successes = 0
420
421! write the new displacement to a file
422 IF (ionode) WRITE (rm, *) nnstep, ' rmbond = ', &
423 rmbond(molecule_type)*angstrom, ' angstroms'
424
425 END IF
426
427! update the maximum displacement for bond angle change
428 IF (move_updates%bias_angle%attempts > 0) THEN
429
430! first account for the extreme cases
431 IF (move_updates%bias_angle%successes == 0) THEN
432 rmangle(molecule_type) = rmangle(molecule_type)/2.0e0_dp
433 ELSE IF (move_updates%bias_angle%successes == &
434 move_updates%bias_angle%attempts) THEN
435 rmangle(molecule_type) = rmangle(molecule_type)*2.0e0_dp
436 ELSE
437! now for the middle case
438 test_ratio = real(move_updates%bias_angle%successes, dp) &
439 /real(move_updates%bias_angle%attempts, dp)/0.5e0_dp
440 IF (test_ratio > 2.0e0_dp) test_ratio = 2.0e0_dp
441 IF (test_ratio < 0.5e0_dp) test_ratio = 0.5e0_dp
442 rmangle(molecule_type) = rmangle(molecule_type)*test_ratio
443 END IF
444
445! more than pi changes meaningless
446 IF (rmangle(molecule_type) > pi) rmangle(molecule_type) = pi
447
448! clear the counters
449 move_updates%bias_angle%attempts = 0
450 move_updates%bias_angle%successes = 0
451
452! write the new displacement to a file
453 IF (ionode) WRITE (rm, *) nnstep, ' rmangle = ', &
454 rmangle(molecule_type)/pi*180.0e0_dp, ' degrees'
455 END IF
456
457! update the maximum displacement for a dihedral change
458 IF (move_updates%bias_dihedral%attempts > 0) THEN
459
460! first account for the extreme cases
461 IF (move_updates%bias_dihedral%successes == 0) THEN
462 rmdihedral(molecule_type) = rmdihedral(molecule_type)/2.0e0_dp
463 ELSE IF (move_updates%bias_dihedral%successes == &
464 move_updates%bias_dihedral%attempts) THEN
465 rmdihedral(molecule_type) = rmdihedral(molecule_type)*2.0e0_dp
466 ELSE
467! now for the middle case
468 test_ratio = real(move_updates%bias_dihedral%successes, dp) &
469 /real(move_updates%bias_dihedral%attempts, dp)/0.5e0_dp
470 IF (test_ratio > 2.0e0_dp) test_ratio = 2.0e0_dp
471 IF (test_ratio < 0.5e0_dp) test_ratio = 0.5e0_dp
472 rmdihedral(molecule_type) = rmdihedral(molecule_type)*test_ratio
473 END IF
474
475! more than pi changes meaningless
476 IF (rmdihedral(molecule_type) > pi) rmdihedral(molecule_type) = pi
477
478! clear the counters
479 move_updates%bias_dihedral%attempts = 0
480 move_updates%bias_dihedral%successes = 0
481
482! write the new displacement to a file
483 IF (ionode) WRITE (rm, *) nnstep, ' rmdihedral = ', &
484 rmdihedral(molecule_type)/pi*180.0e0_dp, ' degrees'
485 END IF
486
487! update the maximum displacement for molecule translation
488 IF (move_updates%bias_trans%attempts > 0) THEN
489
490! first account for the extreme cases
491 IF (move_updates%bias_trans%successes == 0) THEN
492 rmtrans(molecule_type) = rmtrans(molecule_type)/2.0e0_dp
493 ELSE IF (move_updates%bias_trans%successes == &
494 move_updates%bias_trans%attempts) THEN
495 rmtrans(molecule_type) = rmtrans(molecule_type)*2.0e0_dp
496 ELSE
497! now for the middle case
498 test_ratio = real(move_updates%bias_trans%successes, dp) &
499 /real(move_updates%bias_trans%attempts, dp)/0.5e0_dp
500 IF (test_ratio > 2.0e0_dp) test_ratio = 2.0e0_dp
501 IF (test_ratio < 0.5e0_dp) test_ratio = 0.5e0_dp
502 rmtrans(molecule_type) = rmtrans(molecule_type)*test_ratio
503 END IF
504
505 ! make an upper bound...10 a.u.
506 IF (rmtrans(molecule_type) > 10.0e0_dp) THEN
507 rmtrans(molecule_type) = 10.0e0_dp
508 END IF
509
510 ! clear the counters
511 move_updates%bias_trans%attempts = 0
512 move_updates%bias_trans%successes = 0
513
514! write the new displacement to a file
515 IF (ionode) WRITE (rm, *) nnstep, ' rmtrans = ', &
516 rmtrans(molecule_type)*angstrom, ' angstroms'
517 END IF
518
519! update the maximum displacement for cluster translation
520 IF (move_updates%bias_cltrans%attempts > 0) THEN
521
522! first account for the extreme cases
523 IF (move_updates%bias_cltrans%successes == 0) THEN
524 rmcltrans = rmcltrans/2.0e0_dp
525 ELSE IF (move_updates%bias_cltrans%successes == &
526 move_updates%bias_cltrans%attempts) THEN
527 rmcltrans = rmcltrans*2.0e0_dp
528 ELSE
529! now for the middle case
530 test_ratio = real(move_updates%bias_cltrans%successes, dp) &
531 /real(move_updates%bias_cltrans%attempts, dp)/0.5e0_dp
532 IF (test_ratio > 2.0e0_dp) test_ratio = 2.0e0_dp
533 IF (test_ratio < 0.5e0_dp) test_ratio = 0.5e0_dp
534 rmcltrans = rmcltrans*test_ratio
535 END IF
536
537 ! make an upper bound...10 a.u.
538 IF (rmcltrans > 10.0e0_dp) THEN
539 rmcltrans = 10.0e0_dp
540 END IF
541
542 ! clear the counters
543 move_updates%bias_cltrans%attempts = 0
544 move_updates%bias_cltrans%successes = 0
545
546! write the new displacement to a file
547 IF (ionode) WRITE (rm, *) nnstep, ' rmcltrans = ', &
548 rmcltrans*angstrom, ' angstroms'
549 END IF
550
551! update the maximum displacement for molecule rotation
552 IF (move_updates%bias_rot%attempts > 0) THEN
553
554! first account for the extreme cases
555 IF (move_updates%bias_rot%successes == 0) THEN
556 rmrot = rmrot/2.0e0_dp
557
558 IF (rmrot(molecule_type) > pi) rmrot(molecule_type) = pi
559
560 ELSE IF (move_updates%bias_rot%successes == &
561 move_updates%bias_rot%attempts) THEN
562 rmrot(molecule_type) = rmrot(molecule_type)*2.0e0_dp
563
564! more than pi rotation is meaningless
565 IF (rmrot(molecule_type) > pi) rmrot(molecule_type) = pi
566
567 ELSE
568! now for the middle case
569 test_ratio = real(move_updates%bias_rot%successes, dp) &
570 /real(move_updates%bias_rot%attempts, dp)/0.5e0_dp
571 IF (test_ratio > 2.0e0_dp) test_ratio = 2.0e0_dp
572 IF (test_ratio < 0.5e0_dp) test_ratio = 0.5e0_dp
573 rmrot(molecule_type) = rmrot(molecule_type)*test_ratio
574
575! more than pi rotation is meaningless
576 IF (rmrot(molecule_type) > pi) rmrot(molecule_type) = pi
577
578 END IF
579
580! clear the counters
581 move_updates%bias_rot%attempts = 0
582 move_updates%bias_rot%successes = 0
583
584! write the new displacement to a file
585 IF (ionode) WRITE (rm, *) nnstep, ' rmrot = ', &
586 rmrot(molecule_type)/pi*180.0e0_dp, ' degrees'
587 END IF
588
589 CASE ("volume")
590
591! update the maximum displacement for volume displacement
592 IF (move_updates%volume%attempts /= 0) THEN
593
594! first account for the extreme cases
595 IF (move_updates%volume%successes == 0) THEN
596 rmvolume = rmvolume/2.0e0_dp
597
598 ELSE IF (move_updates%volume%successes == &
599 move_updates%volume%attempts) THEN
600 rmvolume = rmvolume*2.0e0_dp
601 ELSE
602! now for the middle case
603 test_ratio = real(move_updates%volume%successes, dp)/ &
604 REAL(move_updates%volume%attempts, dp)/0.5e0_dp
605 IF (test_ratio > 2.0e0_dp) test_ratio = 2.0e0_dp
606 IF (test_ratio < 0.5e0_dp) test_ratio = 0.5e0_dp
607 rmvolume = rmvolume*test_ratio
608
609 END IF
610
611! clear the counters
612 move_updates%volume%attempts = 0
613 move_updates%volume%successes = 0
614
615! write the new displacement to a file
616 IF (ionode) WRITE (rm, *) nnstep, ' rmvolume = ', &
617 rmvolume*angstrom**3, ' angstroms^3'
618
619 END IF
620
621 END SELECT
622
623! set some of the MC parameters
624 CALL set_mc_par(mc_par, rmbond=rmbond, rmangle=rmangle, rmrot=rmrot, &
625 rmtrans=rmtrans, rmcltrans=rmcltrans, rmvolume=rmvolume, rmdihedral=rmdihedral)
626
627! end the timing
628 CALL timestop(handle)
629
630 END SUBROUTINE mc_move_update
631
632END MODULE mc_move_control
633
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