scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

This commit resolves #128 by implementing the stochastic-velocity rescaling (SVR) thermostat into the dynamics module of PW. 

Additions have been made to PW/src/dynamics_module.f90, and to PW/src/input.f90 to allow for options svr (or SVR, or Svr) as ion_temperature (IONS card).
Furthermore, 2 functions were added to Modules/random_numbers.f90: the first function to calculate \sum_i R^2, where R is drawn from the normal distribution,
the second to draw a gamma-distributed random number.
No previous code was touched in this commit, only new functions and cases added.
I updated changelog and documentation.

TODO: future issue and commit, if needed, to implementent SVR to fcp.f90 and vcsmd.f90. Both should be fairly straightforward!
This commit is contained in:
Leonid Kahle 2019-09-10 13:42:53 +02:00 committed by Samuel Ponce
parent 3655553c9c
commit c3a34ea7dc
5 changed files with 186 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
Doc/release-notes
View File

@ -2,7 +2,8 @@ New in development branch:
* turbo_eels code of TDDFPT module now works with ultrasoft pseudopotentials
and spin-orbit coupling together (Oleksandr Motornyi, Andrea Dal Corso,
Nathalie Vast). lr_sm1_psi.f90 of LR_Modules is rewritten and simplified.
* Stochastic-velocity rescaling as a new thermostat for constant-cell MD as
implemented in dynamics_module (Leonid Kahle, Ngoc Linh Nguyen)
Problems fixed in development branch :
* at2celldm wasn't properly converting vectors into celldm parameters
in the ibtrav=91 case (Tone)

82
Modules/random_numbers.f90
View File

@ -204,4 +204,86 @@ MODULE random_numbers
!
END FUNCTION gauss_dist_vect
!
!-----------------------------------------------------------------------
FUNCTION gamma_dist (ialpha)
!-----------------------------------------------------------------------
!
! gamma-distributed random number, implemented as described in
! Numerical recipes (Press, Teukolsky, Vetterling, Flannery)
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: ialpha
REAL(DP) gamma_dist
INTEGER j
REAL(DP) am,e,s,v1,v2,x,y
REAL(DP), external :: ran1
!
IF ( ialpha < 1 ) CALL errore('gamma_dist', 'bad alpha in gamma_dist', 1)
!
! For small alpha, it is more efficient to calculate Gamma as the waiting time
! to the alpha-th event oin a Poisson random process of unit mean.
! Define alpha as small for 0 < alpha < 6:
IF ( ialpha < 6 ) THEN
!
x = 1.0D0
DO j=1,ialpha
x = x * randy()
ENDDO
x = -LOG(x)
ELSE
DO
v1 = 2.0D0*randy()-1.0D0
v2 = 2.0D0*randy()-1.0D0
!
! need to get this condition met:
IF ( v1**2+v2**2 > 1.0D0) CYCLE
!
y = v2 / v1
am = ialpha - 1
s = sqrt(2.0D0 * am + 1.0D0)
x = s * y + am
!
IF ( x <= 0.) CYCLE
!
e = (1.0D0+y**2)* exp( am * log( x / am ) - s * y)
!
IF (randy() > e) THEN
CYCLE
ELSE
EXIT
ENDIF
ENDDO
ENDIF
!
gamma_dist=x
!
ENDFUNCTION gamma_dist
!
!-----------------------------------------------------------------------
FUNCTION sum_of_gaussians2(inum_gaussians)
!-----------------------------------------------------------------------
! returns the sum of inum independent gaussian noises squared, i.e. the result
! is equivalent to summing the square of the return values of inum calls
! to gauss_dist.
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: inum_gaussians
!
REAL(DP) sum_of_gaussians2
!
IF ( inum_gaussians < 0 ) THEN
CALL errore('sum_of_gaussians2', 'negative number of gaussians', 1)
ELSEIF ( inum_gaussians == 0 ) THEN
sum_of_gaussians2 = 0.0D0
ELSEIF ( inum_gaussians == 1 ) THEN
sum_of_gaussians2 = gauss_dist( 0.0D0, 1.0D0 )**2
ELSEIF( MODULO(inum_gaussians,2) == 0 ) THEN
sum_of_gaussians2 = 2.0 * gamma_dist( inum_gaussians/2 )
ELSE
sum_of_gaussians2 = 2.0 * gamma_dist((inum_gaussians-1)/2) + &
gauss_dist( 0.0D0, 1.0D0 )**2
ENDIF
!
ENDFUNCTION sum_of_gaussians2
!
END MODULE random_numbers

9
PW/Doc/INPUT_PW.def
View File

@ -2345,6 +2345,11 @@ input_description -distribution {Quantum Espresso} -package PWscf -program pw.x
control ionic temperature using Andersen thermostat
see parameters @ref tempw and @ref nraise
}
opt -val 'svr' {
control ionic temperature using stochastic-velocity rescaling
(Donadio, Bussi, Parrinello, J. Chem. Phys. 126, 014101, 2007),
with parameters @ref tempw and @ref nraise.
}
opt -val 'initial' {
initialize ion velocities to temperature @ref tempw
and leave uncontrolled further on
@ -2415,6 +2420,10 @@ input_description -distribution {Quantum Espresso} -package PWscf -program pw.x
if @ref ion_temperature == 'andersen' :
the "collision frequency" parameter is given as nu=1/tau
defined above, so nu*dt = 1/nraise
if @ref ion_temperature == 'svr' :
the "characteristic time" of the thermostat is set to
tau = nraise*dt
}
}

87
PW/src/dynamics_module.f90
View File

@ -475,6 +475,13 @@ CONTAINS
& "Characteristic time =",i3,"*timestep")') &
nraise
!
CASE( 'svr', 'Svr', 'SVR' )
!
WRITE( UNIT = stdout, &
FMT = '(/,5X,"temperature is controlled by ", &
& "stochastic velocity rescaling",/,5x,&
& "Characteristic time =",i3,"*timestep")') &
nraise
CASE( 'initial', 'Initial' )
!
WRITE( UNIT = stdout, &
@ -612,6 +619,13 @@ CONTAINS
!
CALL thermalize( nraise, temp_new, temperature )
!
CASE( 'svr', 'Svr', 'SVR' )
!
WRITE( UNIT = stdout, &
FMT = '(/,5X,"Canonical sampling velocity rescaling")' )
!
CALL thermalize_resamp_vscaling( nraise, temp_new, temperature )
!
CASE( 'andersen', 'Andersen' )
!
kt = temperature / ry_to_kelvin
@ -1602,5 +1616,76 @@ CONTAINS
!
END SUBROUTINE smart_MC
!
!
!-----------------------------------------------------------------------
SUBROUTINE thermalize_resamp_vscaling (nraise, system_temp, required_temp)
!-----------------------------------------------------------------------
!
! Sample velocities using stochastic velocity rescaling, based on
! Bussi, Donadio, Parrinello, J. Chem. Phys. 126, 014101 (2007),
! doi: 10.1063/1.2408420
!
! Implemented (2019) by Leonid Kahle and Ngoc Linh Nguyen,
! Theory and Simulations of Materials Laboratory, EPFL.
!
USE ions_base, ONLY : nat, if_pos
USE constraints_module, ONLY : nconstr
USE cell_base, ONLY : alat
USE random_numbers, ONLY : gauss_dist, sum_of_gaussians2
!
IMPLICIT NONE
!
REAL(DP), INTENT(in) :: system_temp, required_temp
INTEGER, INTENT(in) :: nraise
!
INTEGER :: i, ndof
REAL(DP) :: factor, rr
REAL(DP) :: aux, aux2
real(DP), external :: gasdev, sumnoises
INTEGER :: na
!
! ... the number of degrees of freedom
!
IF ( ANY( if_pos(:,:) == 0 ) ) THEN
!
ndof = 3*nat - count( if_pos(:,:) == 0 ) - nconstr
!
ELSE
!
ndof = 3*nat - 3 - nconstr
!
ENDIF
!
IF ( nraise > 0 ) THEN
!
! ... the "rise time" is tau=nraise*dt so dt/tau=1/nraise
! ... Equivalent to traditional rescaling if nraise=1
!
factor = exp(-1.0/nraise)
ELSE
!
factor = 0.0
!
ENDIF
!
IF ( system_temp > 0.D0 .and. required_temp > 0.D0 ) THEN
!
! Applying Eq. (A7) from J. Chem. Phys. 126, 014101 (2007)
!
rr = gauss_dist(0.0D0, 1.0D0)
aux2 = factor + (1.0D0-factor)*( sum_of_gaussians2(ndof-1) +rr**2) &
* required_temp/(ndof*system_temp) &
+ 2*rr*sqrt((factor*(1.0D0-factor)*required_temp)/(ndof*system_temp))
!
aux = sqrt(aux2)
ELSE
!
aux = 0.d0
!
ENDIF
!
! Global rescaling applied to velocities
vel(:,:) = vel(:,:) * aux
!
END SUBROUTINE thermalize_resamp_vscaling
END MODULE dynamics_module

7
PW/src/input.f90
View File

@ -1008,6 +1008,13 @@ SUBROUTINE iosys()
temperature = tempw
nraise_ = nraise
!
CASE( 'svr', 'Svr', 'SVR' )
!
control_temp = .true.
thermostat = trim( ion_temperature )
temperature = tempw
nraise_ = nraise
!
CASE( 'andersen', 'Andersen' )
!
control_temp = .true.