diff --git a/resamplekin.cpp b/resamplekin.cpp new file mode 100644 index 0000000..f95cd42 --- /dev/null +++ b/resamplekin.cpp @@ -0,0 +1,131 @@ + +#include +double resamplekin_sumnoises(int nn); + +double ran1(); +double gasdev(); +double gamdev(const int ia); + +double resamplekin(double kk,double sigma, int ndeg, double taut){ +/* + kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units) + sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk) + ndeg: number of degrees of freedom of the atoms to be thermalized + taut: relaxation time of the thermostat, in units of 'how often this routine is called' +*/ + double factor,rr; + if(taut>0.1){ + factor=exp(-1.0/taut); + } else{ + factor=0.0; + } + rr = gasdev(); + return kk + (1.0-factor)* (sigma*(resamplekin_sumnoises(ndeg-1)+rr*rr)/ndeg-kk) + + 2.0*rr*sqrt(kk*sigma/ndeg*(1.0-factor)*factor); +} + +double resamplekin_sumnoises(int nn){ +/* + returns the sum of n independent gaussian noises squared + (i.e. equivalent to summing the square of the return values of nn calls to gasdev) +*/ + double rr; + if(nn==0) { + return 0.0; + } else if(nn==1) { + rr=gasdev(); + return rr*rr; + } else if(nn%2==0) { + return 2.0*gamdev(nn/2); + } else { + rr=gasdev(); + return 2.0*gamdev((nn-1)/2) + rr*rr; + } +} + + + +double gamdev(const int ia) +{ + int j; + double am,e,s,v1,v2,x,y; + + if (ia < 1) {}; // FATAL ERROR + if (ia < 6) { + x=1.0; + for (j=1;j<=ia;j++) x *= ran1(); + x = -log(x); + } else { + do { + do { + do { + v1=ran1(); + v2=2.0*ran1()-1.0; + } while (v1*v1+v2*v2 > 1.0); + y=v2/v1; + am=ia-1; + s=sqrt(2.0*am+1.0); + x=s*y+am; + } while (x <= 0.0); + e=(1.0+y*y)*exp(am*log(x/am)-s*y); + } while (ran1() > e); + } + return x; +} + + + +double gasdev() +{ + static int iset=0; + static double gset; + double fac,rsq,v1,v2; + + if (iset == 0) { + do { + v1=2.0*ran1()-1.0; + v2=2.0*ran1()-1.0; + rsq=v1*v1+v2*v2; + } while (rsq >= 1.0 || rsq == 0.0); + fac=sqrt(-2.0*log(rsq)/rsq); + gset=v1*fac; + iset=1; + return v2*fac; + } else { + iset=0; + return gset; + } +} + + +double ran1() +{ + const int IA=16807,IM=2147483647,IQ=127773,IR=2836,NTAB=32; + const int NDIV=(1+(IM-1)/NTAB); + const double EPS=3.0e-16,AM=1.0/IM,RNMX=(1.0-EPS); + static int iy=0; + static int iv[NTAB]; + int j,k; + double temp; + static int idum=0; /* ATTENTION: THE SEED IS HARDCODED */ + + if (idum <= 0 || !iy) { + if (-idum < 1) idum=1; + else idum = -idum; + for (j=NTAB+7;j>=0;j--) { + k=idum/IQ; + idum=IA*(idum-k*IQ)-IR*k; + if (idum < 0) idum += IM; + if (j < NTAB) iv[j] = idum; + } + iy=iv[0]; + } + k=idum/IQ; + idum=IA*(idum-k*IQ)-IR*k; + if (idum < 0) idum += IM; + j=iy/NDIV; + iy=iv[j]; + iv[j] = idum; + if ((temp=AM*iy) > RNMX) return RNMX; + else return temp; +} diff --git a/resamplekin.f90 b/resamplekin.f90 new file mode 100644 index 0000000..d69f548 --- /dev/null +++ b/resamplekin.f90 @@ -0,0 +1,157 @@ + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! Stochastic velocity rescale, as described in +! Bussi, Donadio and Parrinello, J. Chem. Phys. 126, 014101 (2007) +! +! This subroutine implements Eq.(A7) and returns the new value for the kinetic energy, +! which can be used to rescale the velocities. +! The procedure can be applied to all atoms or to smaller groups. +! If it is applied to intersecting groups in sequence, the kinetic energy +! that is given as an input (kk) has to be up-to-date with respect to the previous rescalings. +! +! When applied to the entire system, and when performing standard molecular dynamics (fixed c.o.m. (center of mass)) +! the degrees of freedom of the c.o.m. have to be discarded in the calculation of ndeg, +! and the c.o.m. momentum HAS TO BE SET TO ZERO. +! When applied to subgroups, one can chose to: +! (a) calculate the subgroup kinetic energy in the usual reference frame, and count the c.o.m. in ndeg +! (b) calculate the subgroup kinetic energy with respect to its c.o.m. motion, discard the c.o.m. in ndeg +! and apply the rescale factor with respect to the subgroup c.o.m. velocity. +! They should be almost equivalent. +! If the subgroups are expected to move one respect to the other, the choice (b) should be better. +! +! If a null relaxation time is required (taut=0.0), the procedure reduces to an istantaneous +! randomization of the kinetic energy, as described in paragraph IIA. +! +! HOW TO CALCULATE THE EFFECTIVE-ENERGY DRIFT +! The effective-energy (htilde) drift can be used to check the integrator against discretization errors. +! The easiest recipe is: +! htilde = h + conint +! where h is the total energy (kinetic + potential) +! and conint is a quantity accumulated along the trajectory as minus the sum of all the increments of kinetic +! energy due to the thermostat. +! +function resamplekin(kk,sigma,ndeg,taut) + implicit none + real*8 :: resamplekin + real*8, intent(in) :: kk ! present value of the kinetic energy of the atoms to be thermalized (in arbitrary units) + real*8, intent(in) :: sigma ! target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk) + integer, intent(in) :: ndeg ! number of degrees of freedom of the atoms to be thermalized + real*8, intent(in) :: taut ! relaxation time of the thermostat, in units of 'how often this routine is called' + real*8 :: factor,rr + real*8, external :: gasdev + if(taut>0.1) then + factor=exp(-1.0/taut) + else + factor=0.0 + end if + rr = gasdev() + resamplekin = kk + (1.0-factor)* (sigma*(sumnoises(ndeg-1)+rr**2)/ndeg-kk) & + + 2.0*rr*sqrt(kk*sigma/ndeg*(1.0-factor)*factor) + +contains + +double precision function sumnoises(nn) + implicit none + integer, intent(in) :: nn +! returns the sum of n independent gaussian noises squared +! (i.e. equivalent to summing the square of the return values of nn calls to gasdev) + real*8, external :: gamdev,gasdev + if(nn==0) then + sumnoises=0.0 + else if(nn==1) then + sumnoises=gasdev()**2 + else if(modulo(nn,2)==0) then + sumnoises=2.0*gamdev(nn/2) + else + sumnoises=2.0*gamdev((nn-1)/2) + gasdev()**2 + end if +end function sumnoises + +end function resamplekin + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! THE FOLLOWING ROUTINES ARE TRANSCRIBED FROM NUMERICAL RECIPES + +double precision function gamdev(ia) +! gamma-distributed random number, implemented as described in numerical recipes + +implicit none +integer, intent(in) :: ia +integer j +real*8 am,e,s,v1,v2,x,y +real*8, external :: ran1 +if(ia.lt.1)pause 'bad argument in gamdev' +if(ia.lt.6)then + x=1. + do 11 j=1,ia + x=x*ran1() +11 continue + x=-log(x) +else +1 v1=2.*ran1()-1. + v2=2.*ran1()-1. + if(v1**2+v2**2.gt.1.)goto 1 + y=v2/v1 + am=ia-1 + s=sqrt(2.*am+1.) + x=s*y+am + if(x.le.0.)goto 1 + e=(1.+y**2)*exp(am*log(x/am)-s*y) + if(ran1().gt.e)goto 1 +endif +gamdev=x +end + +double precision function gasdev() +! gaussian-distributed random number, implemented as described in numerical recipes + +implicit none +integer, save :: iset = 0 +real*8, save :: gset +real*8, external :: ran1 +real*8 fac,rsq,v1,v2 +if(iset==0) then +1 v1=2.*ran1()-1.0d0 + v2=2.*ran1()-1.0d0 + rsq=v1**2+v2**2 + if(rsq.ge.1..or.rsq.eq.0.)goto 1 + fac=sqrt(-2.*log(rsq)/rsq) + gset=v1*fac + gasdev=v2*fac + iset=1 +else + gasdev=gset + iset=0 +end if +end + +FUNCTION ran1() +! random number generator +INTEGER IA,IM,IQ,IR,NTAB,NDIV +REAL ran1,AM,EPS,RNMX +PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836, & + NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS) +INTEGER j,k,iv(NTAB),iy +SAVE iv,iy +DATA iv /NTAB*0/, iy /0/ +INTEGER, SAVE :: idum=0 !! ATTENTION: THE SEED IS HARDCODED +if (idum.le.0.or.iy.eq.0) then + idum=max(-idum,1) + do 11 j=NTAB+8,1,-1 + k=idum/IQ + idum=IA*(idum-k*IQ)-IR*k + if (idum.lt.0) idum=idum+IM + if (j.le.NTAB) iv(j)=idum +11 continue + iy=iv(1) +endif +k=idum/IQ +idum=IA*(idum-k*IQ)-IR*k +if (idum.lt.0) idum=idum+IM +j=1+iy/NDIV +iy=iv(j) +iv(j)=idum +ran1=min(AM*iy,RNMX) +return + END +