(git:374b731)
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-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! **************************************************************************************************
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 .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
629END MODULE mc_move_control
630
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