CP2K 2.4 (Revision 12889)

csvr_system_utils.f90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------------!
00002 !   CP2K: A general program to perform molecular dynamics simulations         !
00003 !   Copyright (C) 2000 - 2013  CP2K developers group                          !
00004 !-----------------------------------------------------------------------------!
00005 
00006 ! *****************************************************************************
00011 MODULE csvr_system_utils
00012 
00013   USE f77_blas
00014   USE kinds,                           ONLY: dp
00015   USE parallel_rng_types,              ONLY: next_random_number,&
00016                                              rng_stream_type
00017 #include "cp_common_uses.h"
00018 
00019   IMPLICIT NONE
00020  
00021   PRIVATE
00022 
00023   LOGICAL, PARAMETER                   :: debug_this_module=.FALSE.
00024   PUBLIC                               :: rescaling_factor
00025   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'csvr_system_utils'
00026 
00027 CONTAINS
00028 
00029 ! *****************************************************************************
00074   FUNCTION rescaling_factor(kk, sigma, ndeg, taut, rng_stream, error) RESULT(my_res)
00075     REAL(KIND=dp), INTENT(IN)                :: kk, sigma
00076     INTEGER, INTENT(IN)                      :: ndeg
00077     REAL(KIND=dp), INTENT(IN)                :: taut
00078     TYPE(rng_stream_type), POINTER           :: rng_stream
00079     TYPE(cp_error_type), INTENT(INOUT)       :: error
00080     REAL(KIND=dp)                            :: my_res
00081 
00082     CHARACTER(len=*), PARAMETER :: routineN = 'rescaling_factor', 
00083       routineP = moduleN//':'//routineN
00084 
00085     REAL(KIND=dp)                            :: factor, resample, reverse, rr
00086 
00087     my_res=0.0_dp
00088     IF (kk > 0.0_dp) THEN
00089        IF(taut>0.1_dp) THEN
00090           factor=EXP(-1.0_dp/taut)
00091        ELSE
00092           factor=0.0_dp
00093        END IF
00094        rr = next_random_number(rng_stream,error=error)
00095        reverse=1.0_dp
00096        ! reverse of momentum is implemented to have the correct limit to Langevin dynamics for ndeg=1
00097        IF(rr<-SQRT(ndeg*kk*factor/(sigma*(1.0_dp-factor)))) reverse=-1.0_dp
00098        ! for ndeg/=1, the reverse of momentum is not necessary. in principles, it should be there.
00099        ! in practice, it is better to skip it to avoid unnecessary slowing down of the dynamics in the small taut regime
00100        ! anyway, this should not affect the final ensemble
00101        IF(ndeg/=1) reverse=1.0_dp
00102        resample = kk + (1.0_dp-factor)* (sigma*(sumnoises(ndeg-1, rng_stream, error)+rr**2)/REAL(ndeg,KIND=dp)-kk) 
00103                      +  2.0_dp*rr*SQRT(kk*sigma/ndeg*(1.0_dp-factor)*factor)
00104 
00105        resample=MAX(0.0_dp,resample)
00106        my_res = reverse*SQRT(resample/kk)
00107     END IF
00108   END FUNCTION rescaling_factor
00109 
00110 ! *****************************************************************************
00116   FUNCTION sumnoises(nn, rng_stream, error) RESULT(sum_gauss)
00117     INTEGER, INTENT(IN)                      :: nn
00118     TYPE(rng_stream_type), POINTER           :: rng_stream
00119     TYPE(cp_error_type), INTENT(INOUT)       :: error
00120     REAL(KIND=dp)                            :: sum_gauss
00121 
00122     CHARACTER(len=*), PARAMETER :: routineN = 'sumnoises', 
00123       routineP = moduleN//':'//routineN
00124 
00125     INTEGER                                  :: i
00126 
00127     sum_gauss = 0.0_dp
00128     DO i = 1, nn
00129        sum_gauss = sum_gauss + next_random_number(rng_stream,error=error)**2
00130     END DO
00131 
00132   END FUNCTION sumnoises
00133 
00134 END MODULE csvr_system_utils