Thermostat nhc

configureAIMS.sh

@@ -0,0 +1,20 @@
+#!/bin/bash -l
+#SBATCH --job-name=ConfigAIMS
+#SBATCH --partition=rooster,chem          
+#SBATCH --nodes=1
+#SBATCH --ntasks=1                  
+#SBATCH --time=04:00:00
+#SBATCH --output=ConfigAIMS.%j.o
+#SBATCH --error=ConfigAIMS.%j.e
+#SBATCH --gres=gpu:1
+#SBATCH --cpus-per-task=2
+#SBATCH --mail-type=REQUEUE
+##SBATCH --nodelist=gpu193
+#SBATCH --mem=200GB
+
+
+./configure.sh
+make -j
+make unittest
+make test 
+

src/modules/ThermoModule.f90

@@ -3,7 +3,8 @@ module ThermoModule
    use GlobalModule
    implicit none
    private
-   public :: thermo_init, thermo, thermo_bussi_global, thermo_bussi_local, thermo_NormDist, thermo_MBDist
+   public :: thermo_init, thermo, thermo_bussi_global, thermo_bussi_local, thermo_NormDist, &
+           thermo_MBDist, thermo_NHC_MBDist, thermo_NHC_local, thermo_NHC_global
    character(len=8), save :: therm = "" !which thermostat to use in simple interface
    double precision, save :: beta_s = 0.0 !inverse temperature in simple interface
    double precision, save :: thermtime = 0.0 !Thermostat relaxation time (in au) for simple interface
@@ -268,4 +269,184 @@ subroutine vcom_project(ndim, natom, p, mass)
 
    end subroutine vcom_project
 
+   subroutine thermo_NHC_MBDist(beta,M,pNHC,mNHC,thermE)
+   !Initialise Nose cariables
+   integer M
+   real(8) beta,thermE
+   real(8),dimension(M) :: mNHC,pNHC
+   real(8) sigma
+   integer i
+
+   thermE=0.0d0
+   do i=1,M
+     sigma=dsqrt(mNHC(i)/beta) !Maxwell-Boltzmann distribution of Nose momenta
+     pNHC(i)=thermo_NormDist(sigma)
+     thermE=thermE+0.5d0*pNHC(i)*pNHC(i)/mNHC(i)
+   enddo
+   end subroutine thermo_NHC_MBDist
+
+   subroutine thermo_NHC_local(dt,M,p,mass,beta,rNHC,pNHC,mNHC,thermE)
+   implicit none
+   integer M
+   real(8) dt,p,mass,beta,thermE
+   real(8),dimension(M) :: mNHC,rNHC,pNHC
+
+   real(8) akin,scale
+   integer l
+
+   akin=p*p/mass !Note factor of 2. Grr Glenn....
+   call thermo_NHC_prop(dt,akin,beta,1,M,rNHC,pNHC,mNHC,scale)
+
+   !scale momenta
+   p=p*scale
+
+   !Calculate thermostat energy
+   thermE=0.0d0
+   do l=1,M
+     thermE=thermE+0.5d0*pNHC(l)*pNHC(l)/mNHC(l) + rNHC(l)/beta
+   enddo
+
+   end subroutine thermo_NHC_local
+
+   subroutine thermo_NHC_global(dt,M,ndim,natom,p,mass,beta,rNHC,pNHC,mNHC,thermE)
+   implicit none
+   integer M,ndim,natom
+   real(8),dimension(ndim,natom) :: p
+   real(8),dimension(natom) :: mass
+   real(8) dt,beta,thermE
+   real(8),dimension(M) :: mNHC,rNHC,pNHC
+
+   real(8) akin,scale
+   integer iatom,idm,l,NTotDim
+
+   if(any(mass.lt.1.0d-15)) return !Do not thermostat if mass is zero
+   akin=0.0d0
+   do iatom=1,natom
+     do idm=1,ndim
+       akin=akin+p(idm,iatom)*p(idm,iatom)/mass(iatom)
+     enddo
+   enddo
+   NTotDim=ndim*natom
+   call thermo_NHC_prop(dt,akin,beta,NTotDim,M,rNHC,pNHC,mNHC,scale)
+   !scale momenta
+   p=p*scale
+   !Calculate thermostat energy
+   thermE=0.0d0
+   !1st chain
+   thermE=thermE+0.5d0*pNHC(1)*pNHC(1)/mNHC(1) + dble(NTotDim)*rNHC(1)/beta
+   do l=2,M
+     thermE=thermE+0.5d0*pNHC(l)*pNHC(l)/mNHC(l) + rNHC(l)/beta
+   enddo
+   end subroutine thermo_NHC_global
+
+
+   subroutine thermo_NHC_prop(dt,akin,beta,N,M,rNHC,pNHC,mNHC,scale)
+   !Propagates NHC variables from t to dt and scales momenta
+   implicit none
+   integer N     !Number of degrees of freedom
+   integer M     !Length of chain
+   real(8) dt    !Time step for this update (usually half of system timestep)
+   real(8) akin  !twice system kinetic energy to which chain 1 is coupled
+   real(8) beta  !Inverse temperature
+   real(8) scale !Scaling factor for system momenta
+   real(8),dimension(M) :: mNHC & ! NHC masses
+                          ,rNHC &! NHC positions
+                          ,pNHC ! NHC momenta
+
+   !hard coded parameters. These have been tested and appear to be optimal, atleast
+   !for 100 bead harmonic oscillator with 4 NHC variables per mode
+   integer,parameter :: nc=10 ! Number of Trotter Factorizations
+   integer,parameter :: nysmax=7
+   integer,parameter :: nys=5 !Order of Yoshida-Suzuki Integration
+
+   !local variables
+   real(8),dimension(M) :: G,mNHC_inv
+   real(8),dimension(nysmax) :: w
+   real(8) :: beta_inv,dt_nc,dts,AA,dNbeta_inv
+   integer j,k,l
+
+   scale=1.0d0
+   !If any NHC mass is zero, we disable NHC propagation
+   if (any(mNHC.lt.1.0d-15)) then
+     write(6,*)'mNHC is zero'
+     return
+   endif
+
+   !Set up higher-order integration weights
+   select case(nys)
+   case (1)
+     w(1)=1.0d0
+   case (3)
+     w(1)=1.0d0/(2.0d0-2.0d0**(1.0d0/3.0d0))
+     w(2)=1.0d0-2.0d0*w(1)
+     w(3)=w(1)
+   case (5)
+     w(1)=1.0d0/(4.0d0-4.0d0**(1.0d0/3.0d0))
+     w(2)=w(1)
+     w(3)=1.0d0-4.0d0*w(1)
+     w(4)=w(1)
+     w(5)=w(1)
+   case (7)
+     w(1)=-0.117767998417887D1
+     w(2)=0.235573213359357D0
+     w(3)=0.784513610477560D0
+     w(4)=1.0D0 - 2.0D0*(w(1)+w(2)+w(3))
+     w(5)=w(3)
+     w(6)=w(2)
+     w(7)=w(1)
+   case default
+     write(*,*) 'Error, nys=',nys,' not coded'
+   end select
+
+   !Some useful variables
+   beta_inv=1.0d0/beta
+   dt_nc=dt/dble(nc)
+   mNHC_inv(1:M)=1.0d0/mNHC(1:M)
+   dNbeta_inv=beta_inv*dble(N)
+
+   !Update NHC forces
+   G(1)=akin-dNbeta_inv
+   G(2:M)=mNHC_inv(1:M-1)*pNHC(1:M-1)*pNHC(1:M-1)-beta_inv
+
+   !Eqs.35 in Martyna et al. 1996 for the NHC propagator
+   do k=1,nc
+     do j=1,nys
+       dts=w(j)*dt_nc
+       !Update NHC momenta
+       pNHC(M)=pNHC(M)+0.5d0*dts*G(M)
+       do l=M-1,1,-1
+         AA=dexp(-0.25d0*dts*pNHC(l+1)*mNHC_inv(l+1))
+         pNHC(l)=pNHC(l)*AA*AA+0.5d0*dts*G(l)*AA
+       enddo
+       !Update particle momenta
+       AA=dexp(-dts*pNHC(1)*mNHC_inv(1))
+       scale=scale*AA
+       !Update NHC forces
+       G(1)=scale*scale*akin-dNbeta_inv
+       !Update NHC positions
+       do l=1,M
+         rNHC(l)=rNHC(l)+pNHC(l)*dts*mNHC_inv(l)
+       enddo
+       !Update NHC momenta
+       do l=1,M-1
+         AA=dexp(-0.25d0*dts*pNHC(l+1)*mNHC_inv(l+1))
+         pNHC(l)=pNHC(l)*AA*AA+0.5d0*dts*G(l)*AA
+         G(l+1)=mNHC_inv(l)*pNHC(l)*pNHC(l)-beta_inv
+       enddo
+       pNHC(M)=pNHC(M)+G(M)*0.5d0*dts
+     enddo
+   enddo
+   return
+
+   end subroutine thermo_NHC_prop
+
+
+
+
+
+
+
+
+
+
 end module ThermoModule

unit_tests/main.F90

@@ -6,6 +6,7 @@ program tester
    use test_particle, only: collect_particle_suite
    use test_trajectory, only: collect_trajectory_suite
    use test_bundle, only: collect_bundle_suite
+   use test_thermo, only: collect_thermo_suite
    use GlobalModule, only: set_error_handler
 #ifndef __GNUC__
 ! Needed for isatty intrinsic
@@ -29,7 +30,8 @@ program tester
    testsuites = [ &
                 new_testsuite("ParticleModule", collect_particle_suite), &
                 new_testsuite("TrajectoryModule", collect_trajectory_suite), &
-                new_testsuite("BundleModule", collect_bundle_suite) &
+                new_testsuite("BundleModule", collect_bundle_suite), &
+                new_testsuite("ThermoModule",collect_thermo_suite) &
                 ]
 
    ! Swap the default FMS_DieError handler for a unit-test friendly