|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00012 MODULE mc_move_control 00013 USE f77_blas 00014 USE kinds, ONLY: dp 00015 USE mathconstants, ONLY: pi 00016 USE mc_types, ONLY: get_mc_molecule_info,& 00017 get_mc_par,& 00018 mc_molecule_info_type,& 00019 mc_moves_type,& 00020 mc_simpar_type,& 00021 set_mc_par 00022 USE physcon, ONLY: angstrom 00023 USE termination, ONLY: stop_memory,& 00024 stop_program 00025 USE timings, ONLY: timeset,& 00026 timestop 00027 #include "cp_common_uses.h" 00028 00029 IMPLICIT NONE 00030 00031 PRIVATE 00032 00033 ! *** Global parameters *** 00034 00035 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_move_control' 00036 00037 PUBLIC :: init_mc_moves, & 00038 mc_move_update,move_q_reinit,q_move_accept,mc_moves_release,& 00039 write_move_stats 00040 00041 CONTAINS 00042 00043 ! ***************************************************************************** 00051 SUBROUTINE init_mc_moves ( moves ) 00052 00053 TYPE(mc_moves_type), POINTER :: moves 00054 00055 CHARACTER(len=*), PARAMETER :: routineN = 'init_mc_moves', 00056 routineP = moduleN//':'//routineN 00057 00058 INTEGER :: handle, stat 00059 00060 ! begin the timing of the subroutine 00061 00062 CALL timeset(routineN,handle) 00063 00064 ! allocate all the structures 00065 ALLOCATE (moves,stat=stat) 00066 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00067 "moves",0) 00068 ALLOCATE (moves%bond,stat=stat) 00069 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00070 "moves%bond",0) 00071 ALLOCATE (moves%angle,stat=stat) 00072 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00073 "moves%angle",0) 00074 ALLOCATE (moves%dihedral,stat=stat) 00075 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00076 "moves%dihedral",0) 00077 ALLOCATE (moves%trans,stat=stat) 00078 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00079 "moves%trans",0) 00080 ALLOCATE (moves%rot,stat=stat) 00081 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00082 "moves%rot",0) 00083 ALLOCATE (moves%bias_bond,stat=stat) 00084 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00085 "moves%bias_bond",0) 00086 ALLOCATE (moves%bias_angle,stat=stat) 00087 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00088 "moves%bias_angle",0) 00089 ALLOCATE (moves%bias_dihedral,stat=stat) 00090 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00091 "moves%bias_dihedral",0) 00092 ALLOCATE (moves%bias_trans,stat=stat) 00093 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00094 "moves%bias_trans",0) 00095 ALLOCATE (moves%bias_rot,stat=stat) 00096 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00097 "moves%bias_rot",0) 00098 ALLOCATE (moves%volume,stat=stat) 00099 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00100 "moves%volume",0) 00101 ALLOCATE (moves%hmc,stat=stat) 00102 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00103 "moves%hmc",0) 00104 ALLOCATE (moves%avbmc_inin,stat=stat) 00105 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00106 "moves%avbmc_inin",0) 00107 ALLOCATE (moves%avbmc_inout,stat=stat) 00108 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00109 "moves%avbmc_inout",0) 00110 ALLOCATE (moves%avbmc_outin,stat=stat) 00111 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00112 "moves%avbmc_outin",0) 00113 ALLOCATE (moves%avbmc_outout,stat=stat) 00114 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00115 "moves%avbmc_outout",0) 00116 ALLOCATE (moves%swap,stat=stat) 00117 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00118 "moves%swap",0) 00119 ALLOCATE (moves%Quickstep,stat=stat) 00120 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00121 "moves%Quickstep",0) 00122 ! set all the counters equal to zero 00123 moves%bias_bond%attempts=0 00124 moves%bias_bond%successes=0 00125 moves%bias_bond%qsuccesses=0 00126 moves%bias_angle%attempts=0 00127 moves%bias_angle%successes=0 00128 moves%bias_angle%qsuccesses=0 00129 moves%bias_dihedral%attempts=0 00130 moves%bias_dihedral%successes=0 00131 moves%bias_dihedral%qsuccesses=0 00132 moves%bias_trans%attempts=0 00133 moves%bias_trans%successes=0 00134 moves%bias_trans%qsuccesses=0 00135 moves%bias_rot%attempts=0 00136 moves%bias_rot%successes=0 00137 moves%bias_rot%qsuccesses=0 00138 00139 moves%bond%attempts=0 00140 moves%bond%successes=0 00141 moves%bond%qsuccesses=0 00142 moves%angle%attempts=0 00143 moves%angle%successes=0 00144 moves%angle%qsuccesses=0 00145 moves%dihedral%attempts=0 00146 moves%dihedral%successes=0 00147 moves%dihedral%qsuccesses=0 00148 moves%trans%attempts=0 00149 moves%trans%successes=0 00150 moves%trans%qsuccesses=0 00151 moves%rot%attempts=0 00152 moves%rot%successes=0 00153 moves%rot%qsuccesses=0 00154 moves%avbmc_inin%attempts=0 00155 moves%avbmc_inin%successes=0 00156 moves%avbmc_inin%qsuccesses=0 00157 moves%avbmc_inout%attempts=0 00158 moves%avbmc_inout%successes=0 00159 moves%avbmc_inout%qsuccesses=0 00160 moves%avbmc_outin%attempts=0 00161 moves%avbmc_outin%successes=0 00162 moves%avbmc_outin%qsuccesses=0 00163 moves%avbmc_outout%attempts=0 00164 moves%avbmc_outout%successes=0 00165 moves%avbmc_outout%qsuccesses=0 00166 moves%volume%attempts=0 00167 moves%volume%successes=0 00168 moves%volume%qsuccesses=0 00169 moves%hmc%attempts=0 00170 moves%hmc%successes=0 00171 moves%hmc%qsuccesses=0 00172 moves%swap%attempts=0 00173 moves%swap%successes=0 00174 moves%swap%qsuccesses=0 00175 moves%Quickstep%attempts=0 00176 moves%Quickstep%successes=0 00177 moves%Quickstep%qsuccesses=0 00178 moves%trans_dis=0.0E0_dp 00179 moves%qtrans_dis=0.0E0_dp 00180 moves%empty=0 00181 moves%empty_conf=0 00182 moves%empty_avbmc=0 00183 moves%grown=0 00184 ! moves%force_create=1 00185 00186 ! end the timing 00187 CALL timestop(handle) 00188 00189 END SUBROUTINE init_mc_moves 00190 00191 ! ***************************************************************************** 00198 SUBROUTINE mc_moves_release ( moves ) 00199 00200 TYPE(mc_moves_type), POINTER :: moves 00201 00202 CHARACTER(len=*), PARAMETER :: routineN = 'mc_moves_release', 00203 routineP = moduleN//':'//routineN 00204 00205 INTEGER :: handle, stat 00206 00207 ! begin the timing of the subroutine 00208 00209 CALL timeset(routineN,handle) 00210 00211 ! allocate all the structures 00212 DEALLOCATE (moves%bond,stat=stat) 00213 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00214 "moves%bond",0) 00215 DEALLOCATE (moves%angle,stat=stat) 00216 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00217 "moves%angle",0) 00218 DEALLOCATE (moves%dihedral,stat=stat) 00219 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00220 "moves%dihedral",0) 00221 DEALLOCATE (moves%trans,stat=stat) 00222 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00223 "moves%trans",0) 00224 DEALLOCATE (moves%rot,stat=stat) 00225 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00226 "moves%rot",0) 00227 DEALLOCATE (moves%bias_bond,stat=stat) 00228 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00229 "moves%bias_bond",0) 00230 DEALLOCATE (moves%bias_angle,stat=stat) 00231 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00232 "moves%bias_angle",0) 00233 DEALLOCATE (moves%bias_dihedral,stat=stat) 00234 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00235 "moves%bias_dihedral",0) 00236 DEALLOCATE (moves%bias_trans,stat=stat) 00237 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00238 "moves%bias_trans",0) 00239 DEALLOCATE (moves%bias_rot,stat=stat) 00240 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00241 "moves%bias_rot",0) 00242 DEALLOCATE (moves%volume,stat=stat) 00243 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00244 "moves%volume",0) 00245 DEALLOCATE (moves%hmc,stat=stat) 00246 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00247 "moves%hmc",0) 00248 DEALLOCATE (moves%avbmc_inin,stat=stat) 00249 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00250 "moves%avbmc_inin",0) 00251 DEALLOCATE (moves%avbmc_inout,stat=stat) 00252 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00253 "moves%avbmc_inout",0) 00254 DEALLOCATE (moves%avbmc_outin,stat=stat) 00255 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00256 "moves%avbmc_outin",0) 00257 DEALLOCATE (moves%avbmc_outout,stat=stat) 00258 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00259 "moves%avbmc_outout",0) 00260 DEALLOCATE (moves%swap,stat=stat) 00261 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00262 "moves%swap",0) 00263 DEALLOCATE (moves%Quickstep,stat=stat) 00264 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00265 "moves%Quickstep",0) 00266 00267 DEALLOCATE (moves,stat=stat) 00268 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00269 "moves",0) 00270 00271 ! now nullify the moves 00272 NULLIFY(moves) 00273 00274 ! end the timing 00275 CALL timestop(handle) 00276 00277 END SUBROUTINE mc_moves_release 00278 00279 ! ***************************************************************************** 00288 SUBROUTINE move_q_reinit ( moves , lbias ) 00289 00290 TYPE(mc_moves_type), POINTER :: moves 00291 LOGICAL, INTENT(IN) :: lbias 00292 00293 CHARACTER(len=*), PARAMETER :: routineN = 'move_q_reinit', 00294 routineP = moduleN//':'//routineN 00295 00296 INTEGER :: handle 00297 00298 ! begin the timing of the subroutine 00299 00300 CALL timeset(routineN,handle) 00301 00302 ! set all the counters equal to zero 00303 IF(lbias) THEN 00304 moves%bias_bond%qsuccesses=0 00305 moves%bias_angle%qsuccesses=0 00306 moves%bias_dihedral%qsuccesses=0 00307 moves%bias_trans%qsuccesses=0 00308 moves%bias_rot%qsuccesses=0 00309 ELSE 00310 moves%bond%qsuccesses=0 00311 moves%angle%qsuccesses=0 00312 moves%dihedral%qsuccesses=0 00313 moves%trans%qsuccesses=0 00314 moves%rot%qsuccesses=0 00315 moves%volume%qsuccesses=0 00316 moves%hmc%qsuccesses=0 00317 moves%qtrans_dis=0.0E0_dp 00318 ENDIF 00319 00320 ! end the timing 00321 CALL timestop(handle) 00322 00323 END SUBROUTINE move_q_reinit 00324 00325 ! ***************************************************************************** 00336 SUBROUTINE q_move_accept(moves,lbias) 00337 00338 TYPE(mc_moves_type), POINTER :: moves 00339 LOGICAL, INTENT(IN) :: lbias 00340 00341 CHARACTER(len=*), PARAMETER :: routineN = 'q_move_accept', 00342 routineP = moduleN//':'//routineN 00343 00344 INTEGER :: handle 00345 00346 ! begin the timing of the subroutine 00347 00348 CALL timeset(routineN,handle) 00349 00350 IF(lbias) THEN 00351 ! change the number of successful moves for the total move counter 00352 moves%bias_bond%successes=moves%bias_bond%successes& 00353 +moves%bias_bond%qsuccesses 00354 moves%bias_angle%successes=moves%bias_angle%successes& 00355 +moves%bias_angle%qsuccesses 00356 moves%bias_dihedral%successes=moves%bias_dihedral%successes& 00357 +moves%bias_dihedral%qsuccesses 00358 moves%bias_trans%successes=moves%bias_trans%successes& 00359 +moves%bias_trans%qsuccesses 00360 moves%bias_rot%successes=moves%bias_rot%successes& 00361 +moves%bias_rot%qsuccesses 00362 ELSE 00363 ! change the number of successful moves for the total move counter 00364 moves%bond%successes=moves%bond%successes& 00365 +moves%bond%qsuccesses 00366 moves%angle%successes=moves%angle%successes& 00367 +moves%angle%qsuccesses 00368 moves%dihedral%successes=moves%dihedral%successes& 00369 +moves%dihedral%qsuccesses 00370 moves%trans%successes=moves%trans%successes& 00371 +moves%trans%qsuccesses 00372 moves%rot%successes=moves%rot%successes& 00373 +moves%rot%qsuccesses 00374 moves%hmc%successes=moves%hmc%successes& 00375 +moves%hmc%qsuccesses 00376 moves%volume%successes=moves%volume%successes& 00377 +moves%volume%qsuccesses 00378 moves%avbmc_inin%successes=moves%avbmc_inin%successes& 00379 +moves%avbmc_inin%qsuccesses 00380 moves%avbmc_inout%successes=moves%avbmc_inout%successes& 00381 +moves%avbmc_inout%qsuccesses 00382 moves%avbmc_outin%successes=moves%avbmc_outin%successes& 00383 +moves%avbmc_outin%qsuccesses 00384 moves%avbmc_outout%successes=moves%avbmc_outout%successes& 00385 +moves%avbmc_outout%qsuccesses 00386 00387 moves%trans_dis=moves%trans_dis+moves%qtrans_dis 00388 ENDIF 00389 00390 ! end the timing 00391 CALL timestop(handle) 00392 00393 END SUBROUTINE q_move_accept 00394 00395 ! ***************************************************************************** 00405 SUBROUTINE write_move_stats(moves,nnstep,unit) 00406 00407 TYPE(mc_moves_type), POINTER :: moves 00408 INTEGER, INTENT(IN) :: nnstep, unit 00409 00410 CHARACTER(len=*), PARAMETER :: routineN = 'write_move_stats', 00411 routineP = moduleN//':'//routineN 00412 00413 INTEGER :: handle 00414 00415 ! begin the timing of the subroutine 00416 00417 CALL timeset(routineN,handle) 00418 00419 WRITE(unit,1000) nnstep,' bias_bond ',& 00420 moves%bias_bond%successes,moves%bias_bond%attempts 00421 WRITE(unit,1000) nnstep,' bias_angle ',& 00422 moves%bias_angle%successes,moves%bias_angle%attempts 00423 WRITE(unit,1000) nnstep,' bias_dihedral ',& 00424 moves%bias_dihedral%successes,moves%bias_dihedral%attempts 00425 WRITE(unit,1000) nnstep,' bias_trans ',& 00426 moves%bias_trans%successes,moves%bias_trans%attempts 00427 WRITE(unit,1000) nnstep,' bias_rot ',& 00428 moves%bias_rot%successes,moves%bias_rot%attempts 00429 00430 WRITE(unit,1000) nnstep,' bond ',& 00431 moves%bond%successes,moves%bond%attempts 00432 WRITE(unit,1000) nnstep,' angle ',& 00433 moves%angle%successes,moves%angle%attempts 00434 WRITE(unit,1000) nnstep,' dihedral ',& 00435 moves%dihedral%successes,moves%dihedral%attempts 00436 WRITE(unit,1000) nnstep,' trans ',& 00437 moves%trans%successes,moves%trans%attempts 00438 WRITE(unit,1000) nnstep,' rot ',& 00439 moves%rot%successes,moves%rot%attempts 00440 WRITE(unit,1000) nnstep,' swap ',& 00441 moves%swap%successes,moves%swap%attempts 00442 WRITE(unit,1001) nnstep,' grown ',& 00443 moves%grown 00444 WRITE(unit,1001) nnstep,' empty_swap ',& 00445 moves%empty 00446 WRITE(unit,1001) nnstep,' empty_conf ',& 00447 moves%empty_conf 00448 WRITE(unit,1000) nnstep,' volume ',& 00449 moves%volume%successes,moves%volume%attempts 00450 WRITE(unit,1000) nnstep,' HMC ',& 00451 moves%hmc%successes,moves%hmc%attempts 00452 WRITE(unit,1000) nnstep,' avbmc_inin ',& 00453 moves%avbmc_inin%successes,moves%avbmc_inin%attempts 00454 WRITE(unit,1000) nnstep,' avbmc_inout ',& 00455 moves%avbmc_inout%successes,moves%avbmc_inout%attempts 00456 WRITE(unit,1000) nnstep,' avbmc_outin ',& 00457 moves%avbmc_outin%successes,moves%avbmc_outin%attempts 00458 WRITE(unit,1000) nnstep,' avbmc_outout ',& 00459 moves%avbmc_outout%successes,moves%avbmc_outout%attempts 00460 WRITE(unit,1001) nnstep,' empty_avbmc ',& 00461 moves%empty_avbmc 00462 WRITE(unit,1000) nnstep,' Quickstep ',& 00463 moves%quickstep%successes,moves%quickstep%attempts 00464 00465 1000 FORMAT(I10,2X,A,2X,I10,2X,I10) 00466 1001 FORMAT(I10,2X,A,2X,I10) 00467 ! end the timing 00468 CALL timestop(handle) 00469 00470 END SUBROUTINE write_move_stats 00471 00472 ! ***************************************************************************** 00487 SUBROUTINE mc_move_update ( mc_par,move_updates,molecule_type,flag,& 00488 nnstep,ionode ) 00489 00490 TYPE(mc_simpar_type), POINTER :: mc_par 00491 TYPE(mc_moves_type), POINTER :: move_updates 00492 INTEGER, INTENT(IN) :: molecule_type 00493 CHARACTER(LEN=*), INTENT(IN) :: flag 00494 INTEGER, INTENT(IN) :: nnstep 00495 LOGICAL, INTENT(IN) :: ionode 00496 00497 CHARACTER(len=*), PARAMETER :: routineN = 'mc_move_update', 00498 routineP = moduleN//':'//routineN 00499 00500 INTEGER :: handle, nmol_types, rm 00501 REAL(dp), DIMENSION(:), POINTER :: rmangle, rmbond, rmdihedral, 00502 rmrot, rmtrans 00503 REAL(KIND=dp) :: rmvolume, test_ratio 00504 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 00505 00506 ! begin the timing of the subroutine 00507 00508 CALL timeset(routineN,handle) 00509 00510 NULLIFY(rmangle,rmbond,rmdihedral,rmrot,rmtrans) 00511 00512 ! grab some stuff from mc_par 00513 CALL get_mc_par(mc_par,rmbond=rmbond,rmangle=rmangle,rmrot=rmrot,& 00514 rmtrans=rmtrans,rmvolume=rmvolume,rm=rm,rmdihedral=rmdihedral,& 00515 mc_molecule_info=mc_molecule_info) 00516 CALL get_mc_molecule_info(mc_molecule_info,nmol_types=nmol_types) 00517 00518 SELECT CASE (flag) 00519 CASE DEFAULT 00520 WRITE(6,*) 'flag =',flag 00521 CALL stop_program(routineN,moduleN,__LINE__,"Wrong option passed") 00522 CASE ("trans") 00523 00524 ! we need to update all the displacements for every molecule type 00525 IF(ionode) WRITE(rm,*) nnstep,' Data for molecule type ',& 00526 molecule_type 00527 00528 ! update the maximum displacement for bond length change 00529 IF( move_updates%bias_bond%attempts .GT. 0 ) THEN 00530 00531 ! first account for the extreme cases 00532 IF ( move_updates%bias_bond%successes == 0 ) THEN 00533 rmbond(molecule_type)=rmbond(molecule_type)/2.0E0_dp 00534 ELSEIF ( move_updates%bias_bond%successes == & 00535 move_updates%bias_bond%attempts ) THEN 00536 rmbond(molecule_type)=rmbond(molecule_type)*2.0E0_dp 00537 ELSE 00538 ! now for the middle case 00539 test_ratio=REAL(move_updates%bias_bond%successes,dp) 00540 /REAL(move_updates%bias_bond%attempts,dp)/0.5E0_dp 00541 IF (test_ratio .GT. 2.0E0_dp) test_ratio=2.0E0_dp 00542 IF (test_ratio .LT. 0.5E0_dp) test_ratio=0.5E0_dp 00543 rmbond(molecule_type)=rmbond(molecule_type)*test_ratio 00544 ENDIF 00545 00546 ! update and clear the counters 00547 move_updates%bias_bond%attempts=0 00548 move_updates%bias_bond%successes=0 00549 00550 ! write the new displacement to a file 00551 IF(ionode) WRITE(rm,*) nnstep,' rmbond = ',& 00552 rmbond(molecule_type)*angstrom,' angstroms' 00553 00554 ENDIF 00555 00556 ! update the maximum displacement for bond angle change 00557 IF( move_updates%bias_angle%attempts .GT. 0 ) THEN 00558 00559 ! first account for the extreme cases 00560 IF ( move_updates%bias_angle%successes == 0 ) THEN 00561 rmangle(molecule_type)=rmangle(molecule_type)/2.0E0_dp 00562 ELSEIF ( move_updates%bias_angle%successes == & 00563 move_updates%bias_angle%attempts ) THEN 00564 rmangle(molecule_type)=rmangle(molecule_type)*2.0E0_dp 00565 ELSE 00566 ! now for the middle case 00567 test_ratio=REAL(move_updates%bias_angle%successes,dp) 00568 /REAL(move_updates%bias_angle%attempts,dp)/0.5E0_dp 00569 IF (test_ratio .GT. 2.0E0_dp) test_ratio=2.0E0_dp 00570 IF (test_ratio .LT. 0.5E0_dp) test_ratio=0.5E0_dp 00571 rmangle(molecule_type)=rmangle(molecule_type)*test_ratio 00572 ENDIF 00573 00574 ! more than pi changes meaningless 00575 IF (rmangle(molecule_type) .GT. pi) rmangle(molecule_type)=pi 00576 00577 ! clear the counters 00578 move_updates%bias_angle%attempts=0 00579 move_updates%bias_angle%successes=0 00580 00581 ! write the new displacement to a file 00582 IF(ionode) WRITE(rm,*) nnstep,' rmangle = ',& 00583 rmangle(molecule_type)/pi*180.0E0_dp,' degrees' 00584 ENDIF 00585 00586 ! update the maximum displacement for a dihedral change 00587 IF( move_updates%bias_dihedral%attempts .GT. 0 ) THEN 00588 00589 ! first account for the extreme cases 00590 IF ( move_updates%bias_dihedral%successes == 0 ) THEN 00591 rmdihedral(molecule_type)=rmdihedral(molecule_type)/2.0E0_dp 00592 ELSEIF ( move_updates%bias_dihedral%successes == & 00593 move_updates%bias_dihedral%attempts ) THEN 00594 rmdihedral(molecule_type)=rmdihedral(molecule_type)*2.0E0_dp 00595 ELSE 00596 ! now for the middle case 00597 test_ratio=REAL(move_updates%bias_dihedral%successes,dp) 00598 /REAL(move_updates%bias_dihedral%attempts,dp)/0.5E0_dp 00599 IF (test_ratio .GT. 2.0E0_dp) test_ratio=2.0E0_dp 00600 IF (test_ratio .LT. 0.5E0_dp) test_ratio=0.5E0_dp 00601 rmdihedral(molecule_type)=rmdihedral(molecule_type)*test_ratio 00602 ENDIF 00603 00604 ! more than pi changes meaningless 00605 IF (rmdihedral(molecule_type) .GT. pi) rmdihedral(molecule_type)=pi 00606 00607 ! clear the counters 00608 move_updates%bias_dihedral%attempts=0 00609 move_updates%bias_dihedral%successes=0 00610 00611 ! write the new displacement to a file 00612 IF(ionode) WRITE(rm,*) nnstep,' rmdihedral = ',& 00613 rmdihedral(molecule_type)/pi*180.0E0_dp,' degrees' 00614 ENDIF 00615 00616 ! update the maximum displacement for molecule translation 00617 IF( move_updates%bias_trans%attempts .GT. 0 ) THEN 00618 00619 ! first account for the extreme cases 00620 IF ( move_updates%bias_trans%successes == 0 ) THEN 00621 rmtrans(molecule_type)=rmtrans(molecule_type)/2.0E0_dp 00622 ELSEIF ( move_updates%bias_trans%successes == & 00623 move_updates%bias_trans%attempts ) THEN 00624 rmtrans(molecule_type)=rmtrans(molecule_type)*2.0E0_dp 00625 ELSE 00626 ! now for the middle case 00627 test_ratio=REAL(move_updates%bias_trans%successes,dp) 00628 /REAL(move_updates%bias_trans%attempts,dp)/0.5E0_dp 00629 IF (test_ratio .GT. 2.0E0_dp) test_ratio=2.0E0_dp 00630 IF (test_ratio .LT. 0.5E0_dp) test_ratio=0.5E0_dp 00631 rmtrans(molecule_type)=rmtrans(molecule_type)*test_ratio 00632 ENDIF 00633 00634 ! make an upper bound...10 a.u. 00635 IF (rmtrans(molecule_type) .GT. 10.0E0_dp) & 00636 rmtrans(molecule_type )= 10.0E0_dp 00637 00638 ! clear the counters 00639 move_updates%bias_trans%attempts=0 00640 move_updates%bias_trans%successes=0 00641 00642 ! write the new displacement to a file 00643 IF(ionode) WRITE(rm,*) nnstep,' rmtrans = ',& 00644 rmtrans(molecule_type)*angstrom,' angstroms' 00645 ENDIF 00646 00647 ! update the maximum displacement for molecule rotation 00648 IF( move_updates%bias_rot%attempts .GT. 0 ) THEN 00649 00650 ! first account for the extreme cases 00651 IF ( move_updates%bias_rot%successes == 0 ) THEN 00652 rmrot=rmrot/2.0E0_dp 00653 00654 IF (rmrot(molecule_type) .GT. pi) rmrot(molecule_type)=pi 00655 00656 ELSEIF ( move_updates%bias_rot%successes == & 00657 move_updates%bias_rot%attempts ) THEN 00658 rmrot(molecule_type)=rmrot(molecule_type)*2.0E0_dp 00659 00660 ! more than pi rotation is meaningless 00661 IF (rmrot(molecule_type) .GT. pi) rmrot(molecule_type)=pi 00662 00663 ELSE 00664 ! now for the middle case 00665 test_ratio=REAL(move_updates%bias_rot%successes,dp) 00666 /REAL(move_updates%bias_rot%attempts,dp)/0.5E0_dp 00667 IF (test_ratio .GT. 2.0E0_dp) test_ratio=2.0E0_dp 00668 IF (test_ratio .LT. 0.5E0_dp) test_ratio=0.5E0_dp 00669 rmrot(molecule_type)=rmrot(molecule_type)*test_ratio 00670 00671 ! more than pi rotation is meaningless 00672 IF (rmrot(molecule_type) .GT. pi) rmrot(molecule_type)=pi 00673 00674 ENDIF 00675 00676 ! clear the counters 00677 move_updates%bias_rot%attempts=0 00678 move_updates%bias_rot%successes=0 00679 00680 ! write the new displacement to a file 00681 IF(ionode) WRITE(rm,*) nnstep,' rmrot = ',& 00682 rmrot(molecule_type)/pi*180.0E0_dp,' degrees' 00683 ENDIF 00684 00685 CASE ("volume") 00686 00687 ! update the maximum displacement for volume displacement 00688 IF ( move_updates%volume%attempts .NE. 0) THEN 00689 00690 ! first account for the extreme cases 00691 IF ( move_updates%volume%successes == 0 ) THEN 00692 rmvolume=rmvolume/2.0E0_dp 00693 00694 ELSEIF ( move_updates%volume%successes == & 00695 move_updates%volume%attempts ) THEN 00696 rmvolume=rmvolume*2.0E0_dp 00697 ELSE 00698 ! now for the middle case 00699 test_ratio=REAL(move_updates%volume%successes,dp)/ 00700 REAL(move_updates%volume%attempts,dp)/0.5E0_dp 00701 IF (test_ratio .GT. 2.0E0_dp) test_ratio=2.0E0_dp 00702 IF (test_ratio .LT. 0.5E0_dp) test_ratio=0.5E0_dp 00703 rmvolume=rmvolume*test_ratio 00704 00705 ENDIF 00706 00707 ! clear the counters 00708 move_updates%volume%attempts=0 00709 move_updates%volume%successes=0 00710 00711 ! write the new displacement to a file 00712 IF(ionode) WRITE(rm,*) nnstep,' rmvolume = ',& 00713 rmvolume*angstrom**3,' angstroms^3' 00714 00715 ENDIF 00716 00717 END SELECT 00718 00719 ! set some of the MC parameters 00720 CALL set_mc_par(mc_par,rmbond=rmbond,rmangle=rmangle,rmrot=rmrot,& 00721 rmtrans=rmtrans,rmvolume=rmvolume,rmdihedral=rmdihedral) 00722 00723 ! end the timing 00724 CALL timestop(handle) 00725 00726 END SUBROUTINE mc_move_update 00727 00728 END MODULE mc_move_control 00729
1.7.3