Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 : ! **************************************************************************************************
14 : MODULE 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 :
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 84 : 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 84 : CALL timeset(routineN, handle)
60 :
61 : ! allocate all the structures
62 84 : ALLOCATE (moves)
63 84 : ALLOCATE (moves%bond)
64 84 : ALLOCATE (moves%angle)
65 84 : ALLOCATE (moves%dihedral)
66 84 : ALLOCATE (moves%trans)
67 84 : ALLOCATE (moves%cltrans)
68 84 : ALLOCATE (moves%rot)
69 84 : ALLOCATE (moves%bias_bond)
70 84 : ALLOCATE (moves%bias_angle)
71 84 : ALLOCATE (moves%bias_dihedral)
72 84 : ALLOCATE (moves%bias_trans)
73 84 : ALLOCATE (moves%bias_cltrans)
74 84 : ALLOCATE (moves%bias_rot)
75 84 : ALLOCATE (moves%volume)
76 84 : ALLOCATE (moves%hmc)
77 84 : ALLOCATE (moves%avbmc_inin)
78 84 : ALLOCATE (moves%avbmc_inout)
79 84 : ALLOCATE (moves%avbmc_outin)
80 84 : ALLOCATE (moves%avbmc_outout)
81 84 : ALLOCATE (moves%swap)
82 84 : ALLOCATE (moves%Quickstep)
83 :
84 : ! end the timing
85 84 : CALL timestop(handle)
86 :
87 84 : 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 84 : 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 84 : CALL timeset(routineN, handle)
107 :
108 : ! allocate all the structures
109 84 : DEALLOCATE (moves%bond)
110 84 : DEALLOCATE (moves%angle)
111 84 : DEALLOCATE (moves%dihedral)
112 84 : DEALLOCATE (moves%trans)
113 84 : DEALLOCATE (moves%cltrans)
114 84 : DEALLOCATE (moves%rot)
115 84 : DEALLOCATE (moves%bias_bond)
116 84 : DEALLOCATE (moves%bias_angle)
117 84 : DEALLOCATE (moves%bias_dihedral)
118 84 : DEALLOCATE (moves%bias_trans)
119 84 : DEALLOCATE (moves%bias_cltrans)
120 84 : DEALLOCATE (moves%bias_rot)
121 84 : DEALLOCATE (moves%volume)
122 84 : DEALLOCATE (moves%hmc)
123 84 : DEALLOCATE (moves%avbmc_inin)
124 84 : DEALLOCATE (moves%avbmc_inout)
125 84 : DEALLOCATE (moves%avbmc_outin)
126 84 : DEALLOCATE (moves%avbmc_outout)
127 84 : DEALLOCATE (moves%swap)
128 84 : DEALLOCATE (moves%Quickstep)
129 :
130 84 : DEALLOCATE (moves)
131 :
132 : ! now nullify the moves
133 : NULLIFY (moves)
134 :
135 : ! end the timing
136 84 : CALL timestop(handle)
137 :
138 84 : 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 2800 : 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 2800 : CALL timeset(routineN, handle)
161 :
162 : ! set all the counters equal to zero
163 2800 : IF (lbias) THEN
164 1060 : moves%bias_bond%qsuccesses = 0
165 1060 : moves%bias_angle%qsuccesses = 0
166 1060 : moves%bias_dihedral%qsuccesses = 0
167 1060 : moves%bias_trans%qsuccesses = 0
168 1060 : moves%bias_cltrans%qsuccesses = 0
169 1060 : moves%bias_rot%qsuccesses = 0
170 : ELSE
171 1740 : moves%bond%qsuccesses = 0
172 1740 : moves%angle%qsuccesses = 0
173 1740 : moves%dihedral%qsuccesses = 0
174 1740 : moves%trans%qsuccesses = 0
175 1740 : moves%cltrans%qsuccesses = 0
176 1740 : moves%rot%qsuccesses = 0
177 1740 : moves%volume%qsuccesses = 0
178 1740 : moves%hmc%qsuccesses = 0
179 1740 : moves%qtrans_dis = 0.0E0_dp
180 1740 : moves%qcltrans_dis = 0.0E0_dp
181 : END IF
182 :
183 : ! end the timing
184 2800 : CALL timestop(handle)
185 :
186 2800 : 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 2144 : 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 2144 : CALL timeset(routineN, handle)
211 :
212 2144 : 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 732 : + moves%bias_bond%qsuccesses
216 : moves%bias_angle%successes = moves%bias_angle%successes &
217 732 : + moves%bias_angle%qsuccesses
218 : moves%bias_dihedral%successes = moves%bias_dihedral%successes &
219 732 : + moves%bias_dihedral%qsuccesses
220 : moves%bias_trans%successes = moves%bias_trans%successes &
221 732 : + moves%bias_trans%qsuccesses
222 : moves%bias_cltrans%successes = moves%bias_cltrans%successes &
223 732 : + moves%bias_cltrans%qsuccesses
224 : moves%bias_rot%successes = moves%bias_rot%successes &
225 732 : + 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 1412 : + moves%bond%qsuccesses
230 : moves%angle%successes = moves%angle%successes &
231 1412 : + moves%angle%qsuccesses
232 : moves%dihedral%successes = moves%dihedral%successes &
233 1412 : + moves%dihedral%qsuccesses
234 : moves%trans%successes = moves%trans%successes &
235 1412 : + moves%trans%qsuccesses
236 : moves%cltrans%successes = moves%cltrans%successes &
237 1412 : + moves%cltrans%qsuccesses
238 : moves%rot%successes = moves%rot%successes &
239 1412 : + moves%rot%qsuccesses
240 : moves%hmc%successes = moves%hmc%successes &
241 1412 : + moves%hmc%qsuccesses
242 : moves%volume%successes = moves%volume%successes &
243 1412 : + moves%volume%qsuccesses
244 : moves%avbmc_inin%successes = moves%avbmc_inin%successes &
245 1412 : + moves%avbmc_inin%qsuccesses
246 : moves%avbmc_inout%successes = moves%avbmc_inout%successes &
247 1412 : + moves%avbmc_inout%qsuccesses
248 : moves%avbmc_outin%successes = moves%avbmc_outin%successes &
249 1412 : + moves%avbmc_outin%qsuccesses
250 : moves%avbmc_outout%successes = moves%avbmc_outout%successes &
251 1412 : + moves%avbmc_outout%qsuccesses
252 :
253 1412 : moves%trans_dis = moves%trans_dis + moves%qtrans_dis
254 1412 : moves%cltrans_dis = moves%cltrans_dis + moves%qcltrans_dis
255 : END IF
256 :
257 : ! end the timing
258 2144 : CALL timestop(handle)
259 :
260 2144 : 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 29 : 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 29 : CALL timeset(routineN, handle)
284 :
285 29 : WRITE (unit, 1000) nnstep, ' bias_bond ', &
286 58 : moves%bias_bond%successes, moves%bias_bond%attempts
287 29 : WRITE (unit, 1000) nnstep, ' bias_angle ', &
288 58 : moves%bias_angle%successes, moves%bias_angle%attempts
289 29 : WRITE (unit, 1000) nnstep, ' bias_dihedral ', &
290 58 : moves%bias_dihedral%successes, moves%bias_dihedral%attempts
291 29 : WRITE (unit, 1000) nnstep, ' bias_trans ', &
292 58 : moves%bias_trans%successes, moves%bias_trans%attempts
293 29 : WRITE (unit, 1000) nnstep, ' bias_cltrans ', &
294 58 : moves%bias_cltrans%successes, moves%bias_cltrans%attempts
295 29 : WRITE (unit, 1000) nnstep, ' bias_rot ', &
296 58 : moves%bias_rot%successes, moves%bias_rot%attempts
297 :
298 29 : WRITE (unit, 1000) nnstep, ' bond ', &
299 58 : moves%bond%successes, moves%bond%attempts
300 29 : WRITE (unit, 1000) nnstep, ' angle ', &
301 58 : moves%angle%successes, moves%angle%attempts
302 29 : WRITE (unit, 1000) nnstep, ' dihedral ', &
303 58 : moves%dihedral%successes, moves%dihedral%attempts
304 29 : WRITE (unit, 1000) nnstep, ' trans ', &
305 58 : moves%trans%successes, moves%trans%attempts
306 29 : WRITE (unit, 1000) nnstep, ' cltrans ', &
307 58 : moves%cltrans%successes, moves%cltrans%attempts
308 29 : WRITE (unit, 1000) nnstep, ' rot ', &
309 58 : moves%rot%successes, moves%rot%attempts
310 29 : WRITE (unit, 1000) nnstep, ' swap ', &
311 58 : moves%swap%successes, moves%swap%attempts
312 29 : WRITE (unit, 1001) nnstep, ' grown ', &
313 58 : moves%grown
314 29 : WRITE (unit, 1001) nnstep, ' empty_swap ', &
315 58 : moves%empty
316 29 : WRITE (unit, 1001) nnstep, ' empty_conf ', &
317 58 : moves%empty_conf
318 29 : WRITE (unit, 1000) nnstep, ' volume ', &
319 58 : moves%volume%successes, moves%volume%attempts
320 29 : WRITE (unit, 1000) nnstep, ' HMC ', &
321 58 : moves%hmc%successes, moves%hmc%attempts
322 29 : WRITE (unit, 1000) nnstep, ' avbmc_inin ', &
323 58 : moves%avbmc_inin%successes, moves%avbmc_inin%attempts
324 29 : WRITE (unit, 1000) nnstep, ' avbmc_inout ', &
325 58 : moves%avbmc_inout%successes, moves%avbmc_inout%attempts
326 29 : WRITE (unit, 1000) nnstep, ' avbmc_outin ', &
327 58 : moves%avbmc_outin%successes, moves%avbmc_outin%attempts
328 29 : WRITE (unit, 1000) nnstep, ' avbmc_outout ', &
329 58 : moves%avbmc_outout%successes, moves%avbmc_outout%attempts
330 29 : WRITE (unit, 1001) nnstep, ' empty_avbmc ', &
331 58 : moves%empty_avbmc
332 29 : WRITE (unit, 1000) nnstep, ' Quickstep ', &
333 58 : 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 29 : CALL timestop(handle)
339 :
340 29 : 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 0 : 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 0 : REAL(dp), DIMENSION(:), POINTER :: rmangle, rmbond, rmdihedral, rmrot, &
372 0 : 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 0 : CALL timeset(routineN, handle)
379 :
380 0 : 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 0 : mc_molecule_info=mc_molecule_info)
386 0 : CALL get_mc_molecule_info(mc_molecule_info, nmol_types=nmol_types)
387 :
388 0 : SELECT CASE (flag)
389 : CASE DEFAULT
390 0 : WRITE (*, *) 'flag =', flag
391 0 : CPABORT("Wrong option passed")
392 : CASE ("trans")
393 :
394 : ! we need to update all the displacements for every molecule type
395 0 : IF (ionode) WRITE (rm, *) nnstep, ' Data for molecule type ', &
396 0 : molecule_type
397 :
398 : ! update the maximum displacement for bond length change
399 0 : IF (move_updates%bias_bond%attempts .GT. 0) THEN
400 :
401 : ! first account for the extreme cases
402 0 : IF (move_updates%bias_bond%successes == 0) THEN
403 0 : rmbond(molecule_type) = rmbond(molecule_type)/2.0E0_dp
404 0 : ELSEIF (move_updates%bias_bond%successes == &
405 : move_updates%bias_bond%attempts) THEN
406 0 : 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 0 : /REAL(move_updates%bias_bond%attempts, dp)/0.5E0_dp
411 0 : IF (test_ratio .GT. 2.0E0_dp) test_ratio = 2.0E0_dp
412 0 : IF (test_ratio .LT. 0.5E0_dp) test_ratio = 0.5E0_dp
413 0 : rmbond(molecule_type) = rmbond(molecule_type)*test_ratio
414 : END IF
415 :
416 : ! update and clear the counters
417 0 : move_updates%bias_bond%attempts = 0
418 0 : move_updates%bias_bond%successes = 0
419 :
420 : ! write the new displacement to a file
421 0 : IF (ionode) WRITE (rm, *) nnstep, ' rmbond = ', &
422 0 : rmbond(molecule_type)*angstrom, ' angstroms'
423 :
424 : END IF
425 :
426 : ! update the maximum displacement for bond angle change
427 0 : IF (move_updates%bias_angle%attempts .GT. 0) THEN
428 :
429 : ! first account for the extreme cases
430 0 : IF (move_updates%bias_angle%successes == 0) THEN
431 0 : rmangle(molecule_type) = rmangle(molecule_type)/2.0E0_dp
432 0 : ELSEIF (move_updates%bias_angle%successes == &
433 : move_updates%bias_angle%attempts) THEN
434 0 : 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 0 : /REAL(move_updates%bias_angle%attempts, dp)/0.5E0_dp
439 0 : IF (test_ratio .GT. 2.0E0_dp) test_ratio = 2.0E0_dp
440 0 : IF (test_ratio .LT. 0.5E0_dp) test_ratio = 0.5E0_dp
441 0 : rmangle(molecule_type) = rmangle(molecule_type)*test_ratio
442 : END IF
443 :
444 : ! more than pi changes meaningless
445 0 : IF (rmangle(molecule_type) .GT. pi) rmangle(molecule_type) = pi
446 :
447 : ! clear the counters
448 0 : move_updates%bias_angle%attempts = 0
449 0 : move_updates%bias_angle%successes = 0
450 :
451 : ! write the new displacement to a file
452 0 : IF (ionode) WRITE (rm, *) nnstep, ' rmangle = ', &
453 0 : rmangle(molecule_type)/pi*180.0E0_dp, ' degrees'
454 : END IF
455 :
456 : ! update the maximum displacement for a dihedral change
457 0 : IF (move_updates%bias_dihedral%attempts .GT. 0) THEN
458 :
459 : ! first account for the extreme cases
460 0 : IF (move_updates%bias_dihedral%successes == 0) THEN
461 0 : rmdihedral(molecule_type) = rmdihedral(molecule_type)/2.0E0_dp
462 0 : ELSEIF (move_updates%bias_dihedral%successes == &
463 : move_updates%bias_dihedral%attempts) THEN
464 0 : 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 0 : /REAL(move_updates%bias_dihedral%attempts, dp)/0.5E0_dp
469 0 : IF (test_ratio .GT. 2.0E0_dp) test_ratio = 2.0E0_dp
470 0 : IF (test_ratio .LT. 0.5E0_dp) test_ratio = 0.5E0_dp
471 0 : rmdihedral(molecule_type) = rmdihedral(molecule_type)*test_ratio
472 : END IF
473 :
474 : ! more than pi changes meaningless
475 0 : IF (rmdihedral(molecule_type) .GT. pi) rmdihedral(molecule_type) = pi
476 :
477 : ! clear the counters
478 0 : move_updates%bias_dihedral%attempts = 0
479 0 : move_updates%bias_dihedral%successes = 0
480 :
481 : ! write the new displacement to a file
482 0 : IF (ionode) WRITE (rm, *) nnstep, ' rmdihedral = ', &
483 0 : rmdihedral(molecule_type)/pi*180.0E0_dp, ' degrees'
484 : END IF
485 :
486 : ! update the maximum displacement for molecule translation
487 0 : IF (move_updates%bias_trans%attempts .GT. 0) THEN
488 :
489 : ! first account for the extreme cases
490 0 : IF (move_updates%bias_trans%successes == 0) THEN
491 0 : rmtrans(molecule_type) = rmtrans(molecule_type)/2.0E0_dp
492 0 : ELSEIF (move_updates%bias_trans%successes == &
493 : move_updates%bias_trans%attempts) THEN
494 0 : 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 0 : /REAL(move_updates%bias_trans%attempts, dp)/0.5E0_dp
499 0 : IF (test_ratio .GT. 2.0E0_dp) test_ratio = 2.0E0_dp
500 0 : IF (test_ratio .LT. 0.5E0_dp) test_ratio = 0.5E0_dp
501 0 : rmtrans(molecule_type) = rmtrans(molecule_type)*test_ratio
502 : END IF
503 :
504 : ! make an upper bound...10 a.u.
505 0 : IF (rmtrans(molecule_type) .GT. 10.0E0_dp) &
506 0 : rmtrans(molecule_type) = 10.0E0_dp
507 :
508 : ! clear the counters
509 0 : move_updates%bias_trans%attempts = 0
510 0 : move_updates%bias_trans%successes = 0
511 :
512 : ! write the new displacement to a file
513 0 : IF (ionode) WRITE (rm, *) nnstep, ' rmtrans = ', &
514 0 : rmtrans(molecule_type)*angstrom, ' angstroms'
515 : END IF
516 :
517 : ! update the maximum displacement for cluster translation
518 0 : IF (move_updates%bias_cltrans%attempts .GT. 0) THEN
519 :
520 : ! first account for the extreme cases
521 0 : IF (move_updates%bias_cltrans%successes == 0) THEN
522 0 : rmcltrans = rmcltrans/2.0E0_dp
523 0 : ELSEIF (move_updates%bias_cltrans%successes == &
524 : move_updates%bias_cltrans%attempts) THEN
525 0 : rmcltrans = rmcltrans*2.0E0_dp
526 : ELSE
527 : ! now for the middle case
528 : test_ratio = REAL(move_updates%bias_cltrans%successes, dp) &
529 0 : /REAL(move_updates%bias_cltrans%attempts, dp)/0.5E0_dp
530 0 : IF (test_ratio .GT. 2.0E0_dp) test_ratio = 2.0E0_dp
531 0 : IF (test_ratio .LT. 0.5E0_dp) test_ratio = 0.5E0_dp
532 0 : rmcltrans = rmcltrans*test_ratio
533 : END IF
534 :
535 : ! make an upper bound...10 a.u.
536 0 : IF (rmcltrans .GT. 10.0E0_dp) &
537 0 : rmcltrans = 10.0E0_dp
538 :
539 : ! clear the counters
540 0 : move_updates%bias_cltrans%attempts = 0
541 0 : move_updates%bias_cltrans%successes = 0
542 :
543 : ! write the new displacement to a file
544 0 : IF (ionode) WRITE (rm, *) nnstep, ' rmcltrans = ', &
545 0 : rmcltrans*angstrom, ' angstroms'
546 : END IF
547 :
548 : ! update the maximum displacement for molecule rotation
549 0 : IF (move_updates%bias_rot%attempts .GT. 0) THEN
550 :
551 : ! first account for the extreme cases
552 0 : IF (move_updates%bias_rot%successes == 0) THEN
553 0 : rmrot = rmrot/2.0E0_dp
554 :
555 0 : IF (rmrot(molecule_type) .GT. pi) rmrot(molecule_type) = pi
556 :
557 0 : ELSEIF (move_updates%bias_rot%successes == &
558 : move_updates%bias_rot%attempts) THEN
559 0 : rmrot(molecule_type) = rmrot(molecule_type)*2.0E0_dp
560 :
561 : ! more than pi rotation is meaningless
562 0 : 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 0 : /REAL(move_updates%bias_rot%attempts, dp)/0.5E0_dp
568 0 : IF (test_ratio .GT. 2.0E0_dp) test_ratio = 2.0E0_dp
569 0 : IF (test_ratio .LT. 0.5E0_dp) test_ratio = 0.5E0_dp
570 0 : rmrot(molecule_type) = rmrot(molecule_type)*test_ratio
571 :
572 : ! more than pi rotation is meaningless
573 0 : IF (rmrot(molecule_type) .GT. pi) rmrot(molecule_type) = pi
574 :
575 : END IF
576 :
577 : ! clear the counters
578 0 : move_updates%bias_rot%attempts = 0
579 0 : move_updates%bias_rot%successes = 0
580 :
581 : ! write the new displacement to a file
582 0 : IF (ionode) WRITE (rm, *) nnstep, ' rmrot = ', &
583 0 : 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 0 : IF (move_updates%volume%attempts .NE. 0) THEN
590 :
591 : ! first account for the extreme cases
592 0 : IF (move_updates%volume%successes == 0) THEN
593 0 : rmvolume = rmvolume/2.0E0_dp
594 :
595 0 : ELSEIF (move_updates%volume%successes == &
596 : move_updates%volume%attempts) THEN
597 0 : rmvolume = rmvolume*2.0E0_dp
598 : ELSE
599 : ! now for the middle case
600 : test_ratio = REAL(move_updates%volume%successes, dp)/ &
601 0 : REAL(move_updates%volume%attempts, dp)/0.5E0_dp
602 0 : IF (test_ratio .GT. 2.0E0_dp) test_ratio = 2.0E0_dp
603 0 : IF (test_ratio .LT. 0.5E0_dp) test_ratio = 0.5E0_dp
604 0 : rmvolume = rmvolume*test_ratio
605 :
606 : END IF
607 :
608 : ! clear the counters
609 0 : move_updates%volume%attempts = 0
610 0 : move_updates%volume%successes = 0
611 :
612 : ! write the new displacement to a file
613 0 : IF (ionode) WRITE (rm, *) nnstep, ' rmvolume = ', &
614 0 : 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 0 : rmtrans=rmtrans, rmcltrans=rmcltrans, rmvolume=rmvolume, rmdihedral=rmdihedral)
623 :
624 : ! end the timing
625 0 : CALL timestop(handle)
626 :
627 0 : END SUBROUTINE mc_move_update
628 :
629 : END MODULE mc_move_control
630 :
|
