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