From 4e10dd8d2d8e595def63cf4939fcc8bc770af24a Mon Sep 17 00:00:00 2001 From: Tingyao Lu Date: Wed, 28 Jan 2026 22:08:59 +0800 Subject: [PATCH 1/4] Added ThermoModule (NHC thermostat). --- configureAIMS.sh | 20 ++ src/ThermoModule.f90 | 499 ++++++++++++++++++++++++++++++++++ src/modules/thermomodule.mod0 | 0 3 files changed, 519 insertions(+) create mode 100755 configureAIMS.sh create mode 100644 src/ThermoModule.f90 create mode 100644 src/modules/thermomodule.mod0 diff --git a/configureAIMS.sh b/configureAIMS.sh new file mode 100755 index 0000000..10de423 --- /dev/null +++ b/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 + diff --git a/src/ThermoModule.f90 b/src/ThermoModule.f90 new file mode 100644 index 0000000..d899d74 --- /dev/null +++ b/src/ThermoModule.f90 @@ -0,0 +1,499 @@ +Module ThermoModule +implicit none +private +public :: thermo_init, thermo_NHC_MBDist, thermo_NHC_local, thermo_NHC_global, & + thermo_bussi_global, thermo_bussi_local, thermo_NormDist, thermo_MBDist, vcom_project + +contains + +subroutine thermo_init(iseed) +implicit none +integer iseed + +real(8) foo +integer idum + +idum=sign(iseed,-1) !ensure seed is negative +call thermo_ran(foo,idum) !Initialize random number generator + +end subroutine thermo_init + +function thermo_NormDist(sigma) result(x) +!This subroutine generates a normal distribution exp(-x^2/2sigma^2) +implicit none +real(8) :: x +real(8),intent(in) :: sigma + +x=gasdev()*sigma +end function thermo_NormDist + +subroutine thermo_MBDist(ndim,natom,p,mass,beta,zcom) +implicit none +integer ndim,natom +real(8) p(ndim,natom),mass(natom) +real(8) beta +logical zcom + +real(8) mxw_unitf,sigma +integer iatom,idm + +!Sample momenta from Maxwell-Boltzmann distribution +mxw_unitf=1.0d0/dsqrt(beta) +do iatom=1,natom + sigma=mxw_unitf*dsqrt(mass(iatom)) + do idm=1,ndim + p(idm,iatom)=thermo_NormDist(sigma) + enddo +enddo + + +!Project out center of mass velocity if not a degree of freedom +if (.not.zcom) then + call vcom_project(ndim,natom,p,mass) +endif +end subroutine thermo_MBDist + +subroutine thermo_NHC_MBDist(beta,M,pNHC,mNHC,thermE) +!Initialize Nose variables +implicit none +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 + +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 + + +subroutine thermo_bussi_global(ndim,natom,p,mass,beta,tau,thermE,zcom) +implicit none +integer ndim,natom +real(8) p(ndim,natom),mass(natom) +real(8) beta,tau,thermE +logical zcom + +real(8) ekin_old,ekin_new,ekinfac,sigma,signfac,vscale +integer ndeg,iatom,idm + +ekin_old=0.0d0 +do iatom=1,natom + ekinfac=0.5d0/mass(iatom) + do idm=1,ndim + ekin_old=ekin_old+ekinfac*p(idm,iatom)*p(idm,iatom) + enddo +enddo + +if (.not.zcom) then +!Do not include center of mass as degree of freedom + call vcom_project(ndim,natom,p,mass) + ndeg=ndim*natom-ndim !Rotations are not frozen out, so keep them as degrees of freedom +else + ndeg=ndim*natom +endif + +if (tau.gt.1.0d-15) then !only apply thermostat with nonzero relaxation time + sigma=0.5d0*dble(ndeg)/beta + + ekin_new=resamplekin(ekin_old,sigma,ndeg,tau) + + signfac=sign(1.0d0,ekin_new) + ekin_new=dabs(ekin_new) + + vscale=signfac*dsqrt(ekin_new/ekin_old) + + do iatom=1,natom + do idm=1,ndim + p(idm,iatom)=p(idm,iatom)*vscale + enddo + enddo +end if + +ekin_new=0.0d0 +do iatom=1,natom + ekinfac=0.5d0/mass(iatom) + do idm=1,ndim + ekin_new=ekin_new+ekinfac*p(idm,iatom)*p(idm,iatom) + enddo +enddo +!Calculate new thermostat energy +thermE=thermE+ekin_old-ekin_new + +end subroutine thermo_bussi_global + +subroutine thermo_bussi_local(p,mass,beta,tau,thermE) +implicit none +real(8) p,mass +real(8) beta,tau,thermE + +real(8) c1,c2, c2fac, ekin_old, ekin_new, ekinfac + +if (tau.lt.1.0d-15) return !only apply thermostat with nonzero relaxation time + +ekinfac=0.5d0/mass +ekin_old=ekinfac*p*p + +c1=dexp(-1.0d0/tau) +c2=dsqrt(1.0d0-c1*c1) + +c2fac=c2*dsqrt(mass/beta) +p=c1*p+ c2fac*gasdev() + +ekin_new=ekinfac*p*p + +!Calculate new thermostat energy +thermE=thermE+ekin_old-ekin_new + +end subroutine thermo_bussi_local + +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 ::dndeg + dndeg=dble(ndeg) + if(taut>0.1d0) then + factor=dexp(-1.0d0/taut) + else + factor=0.0d0 + end if + rr = gasdev() + resamplekin = kk + (1.0d0-factor)* (sigma*(sumnoises(ndeg-1)+rr**2)/dndeg-kk) & + + 2.0d0*rr*dsqrt(kk*sigma/dndeg*(1.0d0-factor)*factor) + +! Transfer sign to resamplekin + resamplekin=sign(resamplekin,rr+dsqrt(dndeg*factor/(sigma*(1.0d0-factor)))) + +end function resamplekin + +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) + if(nn==0) then + sumnoises=0.0d0 + else if(nn==1) then + sumnoises=gasdev()**2 + else if(modulo(nn,2)==0) then + sumnoises=2.0d0*gamdev(nn/2) + else + sumnoises=2.0d0*gamdev((nn-1)/2) + gasdev()**2 + end if +end function sumnoises + +function gamdev(ia) +! gamma-distributed random number, implemented as described in numerical recipes + +implicit none +real(8) :: gamdev +integer, intent(in) :: ia +integer j +real(8) am,e,s,v1,v2,x,y,z +if(ia.lt.1)pause 'bad argument in gamdev' +if(ia.lt.6)then + x=1. + do 11 j=1,ia + call thermo_ran(z) + x=x*z +11 continue + x=-log(x) +else +1 call thermo_ran(z) + v1=2.0d0*z-1.0d0 + call thermo_ran(z) + v2=2.0d0*z-1.0d0 + if(v1**2+v2**2.gt.1.0d0)goto 1 + y=v2/v1 + am=ia-1 + s=dsqrt(2.0d0*am+1.0d0) + x=s*y+am + if(x.le.0.0d0)goto 1 + e=(1.0d0+y**2)*dexp(am*dlog(x/am)-s*y) + call thermo_ran(z) + if(z.gt.e)goto 1 +endif +gamdev=x +return +end function gamdev + +function gasdev() + +implicit none +real(8) gasdev +real(8) fac,gset,rsq,v1,v2,z +INTEGER iset +SAVE iset,gset +DATA iset/0/ + +if (iset.eq.0) then +1 call thermo_ran(z) + v1=2.0d0*z-1.0d0 + call thermo_ran(z) + v2=2.0d0*z-1.0d0 + rsq=v1**2+v2**2 + if(rsq.ge.1.0d0.or.rsq.eq.0.0d0)goto 1 + fac=sqrt(-2.0d0*log(rsq)/rsq) + gset=v1*fac + gasdev=v2*fac + iset=1 +else + gasdev=gset + iset=0 +endif +return +END FUNCTION GASDEV + +SUBROUTINE THERMO_RAN(rnd,iseed) +! interface to random number generators +implicit none +integer, optional :: iseed +real(8) rnd, foo + +if (present(iseed)) then +! call init_random_seed(iseed) + + foo=ran1(iseed) +endif + +rnd=ran1() +! call RANDOM_NUMBER(rnd) +return +end SUBROUTINE THERMO_RAN + +SUBROUTINE init_random_seed(iseed) +IMPLICIT NONE +INTEGER :: i, n, iseed +INTEGER, DIMENSION(:), ALLOCATABLE :: seed + +CALL RANDOM_SEED(size = n) +ALLOCATE(seed(n)) + +seed = iseed * (/ (i - 1, i = 1, n) /) +CALL RANDOM_SEED(PUT = seed) + +DEALLOCATE(seed) +END SUBROUTINE init_random_seed + + +FUNCTION ran1(iseed) +! random number generator +integer, optional :: iseed +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 +if (present(iseed)) idum=iseed +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 function ran1 + +subroutine vcom_project(ndim,natom,p,mass) +implicit none +integer ndim,natom +real(8) p(ndim,natom) +real(8) mass(natom) + +real(8) vtot(ndim) +real(8) Mtot_inv +integer iatom,idm + +vtot(1:ndim)=0.0d0 +Mtot_inv=1.0d0/sum(mass) + +do iatom=1,natom + do idm=1,ndim + vtot(idm)=vtot(idm)+p(idm,iatom)*Mtot_inv + enddo +enddo +do iatom=1,natom + do idm=1,ndim + p(idm,iatom)=p(idm,iatom)-mass(iatom)*vtot(idm) + enddo +enddo + +end subroutine vcom_project + +end module ThermoModule diff --git a/src/modules/thermomodule.mod0 b/src/modules/thermomodule.mod0 new file mode 100644 index 0000000..e69de29 From 31c358eaff8c4e7cd771db29808c2d274918642c Mon Sep 17 00:00:00 2001 From: Tingyao Lu Date: Wed, 28 Jan 2026 22:10:02 +0800 Subject: [PATCH 2/4] Added ThermoModule (NHC thermostat). --- src/ThermoModule.f90 | 2 +- src/modules/thermomodule.mod0 | 0 2 files changed, 1 insertion(+), 1 deletion(-) delete mode 100644 src/modules/thermomodule.mod0 diff --git a/src/ThermoModule.f90 b/src/ThermoModule.f90 index d899d74..e1f0b78 100644 --- a/src/ThermoModule.f90 +++ b/src/ThermoModule.f90 @@ -352,7 +352,7 @@ function gamdev(ia) integer, intent(in) :: ia integer j real(8) am,e,s,v1,v2,x,y,z -if(ia.lt.1)pause 'bad argument in gamdev' +if (ia < 1) error stop "bad argument in gamdev" if(ia.lt.6)then x=1. do 11 j=1,ia diff --git a/src/modules/thermomodule.mod0 b/src/modules/thermomodule.mod0 deleted file mode 100644 index e69de29..0000000 From dfc569cf9c9f2819925f60dcb1e9d6e33e037bf0 Mon Sep 17 00:00:00 2001 From: Tingyao Lu Date: Thu, 29 Jan 2026 16:56:39 +0800 Subject: [PATCH 3/4] Adding NHC to ThermoModule --- src/ThermoModule.f90 | 499 ----------------------------------- src/modules/ThermoModule.f90 | 182 ++++++++++++- 2 files changed, 181 insertions(+), 500 deletions(-) delete mode 100644 src/ThermoModule.f90 diff --git a/src/ThermoModule.f90 b/src/ThermoModule.f90 deleted file mode 100644 index e1f0b78..0000000 --- a/src/ThermoModule.f90 +++ /dev/null @@ -1,499 +0,0 @@ -Module ThermoModule -implicit none -private -public :: thermo_init, thermo_NHC_MBDist, thermo_NHC_local, thermo_NHC_global, & - thermo_bussi_global, thermo_bussi_local, thermo_NormDist, thermo_MBDist, vcom_project - -contains - -subroutine thermo_init(iseed) -implicit none -integer iseed - -real(8) foo -integer idum - -idum=sign(iseed,-1) !ensure seed is negative -call thermo_ran(foo,idum) !Initialize random number generator - -end subroutine thermo_init - -function thermo_NormDist(sigma) result(x) -!This subroutine generates a normal distribution exp(-x^2/2sigma^2) -implicit none -real(8) :: x -real(8),intent(in) :: sigma - -x=gasdev()*sigma -end function thermo_NormDist - -subroutine thermo_MBDist(ndim,natom,p,mass,beta,zcom) -implicit none -integer ndim,natom -real(8) p(ndim,natom),mass(natom) -real(8) beta -logical zcom - -real(8) mxw_unitf,sigma -integer iatom,idm - -!Sample momenta from Maxwell-Boltzmann distribution -mxw_unitf=1.0d0/dsqrt(beta) -do iatom=1,natom - sigma=mxw_unitf*dsqrt(mass(iatom)) - do idm=1,ndim - p(idm,iatom)=thermo_NormDist(sigma) - enddo -enddo - - -!Project out center of mass velocity if not a degree of freedom -if (.not.zcom) then - call vcom_project(ndim,natom,p,mass) -endif -end subroutine thermo_MBDist - -subroutine thermo_NHC_MBDist(beta,M,pNHC,mNHC,thermE) -!Initialize Nose variables -implicit none -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 - -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 - - -subroutine thermo_bussi_global(ndim,natom,p,mass,beta,tau,thermE,zcom) -implicit none -integer ndim,natom -real(8) p(ndim,natom),mass(natom) -real(8) beta,tau,thermE -logical zcom - -real(8) ekin_old,ekin_new,ekinfac,sigma,signfac,vscale -integer ndeg,iatom,idm - -ekin_old=0.0d0 -do iatom=1,natom - ekinfac=0.5d0/mass(iatom) - do idm=1,ndim - ekin_old=ekin_old+ekinfac*p(idm,iatom)*p(idm,iatom) - enddo -enddo - -if (.not.zcom) then -!Do not include center of mass as degree of freedom - call vcom_project(ndim,natom,p,mass) - ndeg=ndim*natom-ndim !Rotations are not frozen out, so keep them as degrees of freedom -else - ndeg=ndim*natom -endif - -if (tau.gt.1.0d-15) then !only apply thermostat with nonzero relaxation time - sigma=0.5d0*dble(ndeg)/beta - - ekin_new=resamplekin(ekin_old,sigma,ndeg,tau) - - signfac=sign(1.0d0,ekin_new) - ekin_new=dabs(ekin_new) - - vscale=signfac*dsqrt(ekin_new/ekin_old) - - do iatom=1,natom - do idm=1,ndim - p(idm,iatom)=p(idm,iatom)*vscale - enddo - enddo -end if - -ekin_new=0.0d0 -do iatom=1,natom - ekinfac=0.5d0/mass(iatom) - do idm=1,ndim - ekin_new=ekin_new+ekinfac*p(idm,iatom)*p(idm,iatom) - enddo -enddo -!Calculate new thermostat energy -thermE=thermE+ekin_old-ekin_new - -end subroutine thermo_bussi_global - -subroutine thermo_bussi_local(p,mass,beta,tau,thermE) -implicit none -real(8) p,mass -real(8) beta,tau,thermE - -real(8) c1,c2, c2fac, ekin_old, ekin_new, ekinfac - -if (tau.lt.1.0d-15) return !only apply thermostat with nonzero relaxation time - -ekinfac=0.5d0/mass -ekin_old=ekinfac*p*p - -c1=dexp(-1.0d0/tau) -c2=dsqrt(1.0d0-c1*c1) - -c2fac=c2*dsqrt(mass/beta) -p=c1*p+ c2fac*gasdev() - -ekin_new=ekinfac*p*p - -!Calculate new thermostat energy -thermE=thermE+ekin_old-ekin_new - -end subroutine thermo_bussi_local - -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 ::dndeg - dndeg=dble(ndeg) - if(taut>0.1d0) then - factor=dexp(-1.0d0/taut) - else - factor=0.0d0 - end if - rr = gasdev() - resamplekin = kk + (1.0d0-factor)* (sigma*(sumnoises(ndeg-1)+rr**2)/dndeg-kk) & - + 2.0d0*rr*dsqrt(kk*sigma/dndeg*(1.0d0-factor)*factor) - -! Transfer sign to resamplekin - resamplekin=sign(resamplekin,rr+dsqrt(dndeg*factor/(sigma*(1.0d0-factor)))) - -end function resamplekin - -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) - if(nn==0) then - sumnoises=0.0d0 - else if(nn==1) then - sumnoises=gasdev()**2 - else if(modulo(nn,2)==0) then - sumnoises=2.0d0*gamdev(nn/2) - else - sumnoises=2.0d0*gamdev((nn-1)/2) + gasdev()**2 - end if -end function sumnoises - -function gamdev(ia) -! gamma-distributed random number, implemented as described in numerical recipes - -implicit none -real(8) :: gamdev -integer, intent(in) :: ia -integer j -real(8) am,e,s,v1,v2,x,y,z -if (ia < 1) error stop "bad argument in gamdev" -if(ia.lt.6)then - x=1. - do 11 j=1,ia - call thermo_ran(z) - x=x*z -11 continue - x=-log(x) -else -1 call thermo_ran(z) - v1=2.0d0*z-1.0d0 - call thermo_ran(z) - v2=2.0d0*z-1.0d0 - if(v1**2+v2**2.gt.1.0d0)goto 1 - y=v2/v1 - am=ia-1 - s=dsqrt(2.0d0*am+1.0d0) - x=s*y+am - if(x.le.0.0d0)goto 1 - e=(1.0d0+y**2)*dexp(am*dlog(x/am)-s*y) - call thermo_ran(z) - if(z.gt.e)goto 1 -endif -gamdev=x -return -end function gamdev - -function gasdev() - -implicit none -real(8) gasdev -real(8) fac,gset,rsq,v1,v2,z -INTEGER iset -SAVE iset,gset -DATA iset/0/ - -if (iset.eq.0) then -1 call thermo_ran(z) - v1=2.0d0*z-1.0d0 - call thermo_ran(z) - v2=2.0d0*z-1.0d0 - rsq=v1**2+v2**2 - if(rsq.ge.1.0d0.or.rsq.eq.0.0d0)goto 1 - fac=sqrt(-2.0d0*log(rsq)/rsq) - gset=v1*fac - gasdev=v2*fac - iset=1 -else - gasdev=gset - iset=0 -endif -return -END FUNCTION GASDEV - -SUBROUTINE THERMO_RAN(rnd,iseed) -! interface to random number generators -implicit none -integer, optional :: iseed -real(8) rnd, foo - -if (present(iseed)) then -! call init_random_seed(iseed) - - foo=ran1(iseed) -endif - -rnd=ran1() -! call RANDOM_NUMBER(rnd) -return -end SUBROUTINE THERMO_RAN - -SUBROUTINE init_random_seed(iseed) -IMPLICIT NONE -INTEGER :: i, n, iseed -INTEGER, DIMENSION(:), ALLOCATABLE :: seed - -CALL RANDOM_SEED(size = n) -ALLOCATE(seed(n)) - -seed = iseed * (/ (i - 1, i = 1, n) /) -CALL RANDOM_SEED(PUT = seed) - -DEALLOCATE(seed) -END SUBROUTINE init_random_seed - - -FUNCTION ran1(iseed) -! random number generator -integer, optional :: iseed -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 -if (present(iseed)) idum=iseed -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 function ran1 - -subroutine vcom_project(ndim,natom,p,mass) -implicit none -integer ndim,natom -real(8) p(ndim,natom) -real(8) mass(natom) - -real(8) vtot(ndim) -real(8) Mtot_inv -integer iatom,idm - -vtot(1:ndim)=0.0d0 -Mtot_inv=1.0d0/sum(mass) - -do iatom=1,natom - do idm=1,ndim - vtot(idm)=vtot(idm)+p(idm,iatom)*Mtot_inv - enddo -enddo -do iatom=1,natom - do idm=1,ndim - p(idm,iatom)=p(idm,iatom)-mass(iatom)*vtot(idm) - enddo -enddo - -end subroutine vcom_project - -end module ThermoModule diff --git a/src/modules/ThermoModule.f90 b/src/modules/ThermoModule.f90 index d8036dd..43f0c07 100644 --- a/src/modules/ThermoModule.f90 +++ b/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,183 @@ 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 + + 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 From 8ad45ab49329eafbe254ea68a0f89ae705a3dd04 Mon Sep 17 00:00:00 2001 From: Tingyao Lu Date: Sat, 31 Jan 2026 05:29:51 +0800 Subject: [PATCH 4/4] Adding unit tests for ThermoModule NHC functions --- src/modules/ThermoModule.f90 | 3 +- unit_tests/main.F90 | 4 +- unit_tests/test_thermo.F90 | 237 +++++++++++++++++++++++++++++++++++ 3 files changed, 242 insertions(+), 2 deletions(-) create mode 100644 unit_tests/test_thermo.F90 diff --git a/src/modules/ThermoModule.f90 b/src/modules/ThermoModule.f90 index 43f0c07..f8e7f63 100644 --- a/src/modules/ThermoModule.f90 +++ b/src/modules/ThermoModule.f90 @@ -407,7 +407,8 @@ subroutine thermo_NHC_prop(dt,akin,beta,N,M,rNHC,pNHC,mNHC,scale) !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 diff --git a/unit_tests/main.F90 b/unit_tests/main.F90 index a25ddc9..bc55da3 100644 --- a/unit_tests/main.F90 +++ b/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 diff --git a/unit_tests/test_thermo.F90 b/unit_tests/test_thermo.F90 new file mode 100644 index 0000000..68572e4 --- /dev/null +++ b/unit_tests/test_thermo.F90 @@ -0,0 +1,237 @@ +!> Test suite for ThermoModule (only Nose–Hoover chain components) +!! +!! This suite is intentionally RNG-free: it tests +!! - thermo_NHC_local +!! - thermo_NHC_global + +module test_thermo + use testdrive, only: new_unittest, unittest_type, error_type, check + use testutils, only: check_dieerror_called + use GlobalModule, only: DefReal, DefInt + use ParticleModule + use ThermoModule, only: thermo_NHC_local, thermo_NHC_global + implicit none + private + + public :: collect_thermo_suite + +contains + subroutine collect_thermo_suite(tests) + type(unittest_type), allocatable, intent(out) :: tests(:) + + tests = [ & + new_unittest("NHC local: dt=0 leaves state unchanged", test_nhc_local_dt0), & + new_unittest("NHC local: forward/backward reversibility", test_nhc_local_reversible), & + new_unittest("NHC global: dt=0 leaves state unchanged", test_nhc_global_dt0), & + new_unittest("NHC global: forward/backward reversibility", test_nhc_global_reversible), & + new_unittest("NHC global: bypass when any particle mass=0", test_nhc_global_mass0_bypass) & + ] + end subroutine collect_thermo_suite + + + subroutine test_nhc_local_dt0(error) + type(error_type), allocatable, intent(out) :: error + + integer :: M, l + real(8) :: dt, p, p0, mass, beta, thermE, thermE_expected, tol + real(8), allocatable :: mNHC(:), rNHC(:), pNHC(:), r0(:), pN0(:) + + tol = 1.0d-12 + M = 4 + dt = 0.0d0 + beta = 0.7d0 + mass = 2.3d0 + p = 1.234d0 + p0 = p + + allocate(mNHC(M), rNHC(M), pNHC(M), r0(M), pN0(M)) + + ! deterministic initial thermostat state (no RNG) + mNHC = [ 1.0d0, 1.7d0, 2.1d0, 3.3d0 ] + rNHC = [ 0.10d0, -0.05d0, 0.03d0, 0.07d0 ] + pNHC = [ 0.20d0, -0.10d0, 0.40d0, -0.30d0 ] + + r0 = rNHC + pN0 = pNHC + + call thermo_NHC_local(dt, M, p, mass, beta, rNHC, pNHC, mNHC, thermE) + + call check(error, abs(p - p0) <= tol, "p changed with dt=0 in thermo_NHC_local") + call check(error, maxval(abs(rNHC - r0)) <= tol, "rNHC changed with dt=0 in thermo_NHC_local") + call check(error, maxval(abs(pNHC - pN0)) <= tol, "pNHC changed with dt=0 in thermo_NHC_local") + + thermE_expected = 0.0d0 + do l = 1, M + thermE_expected = thermE_expected + 0.5d0*pN0(l)*pN0(l)/mNHC(l) + r0(l)/beta + end do + + call check(error, abs(thermE - thermE_expected) <= 1.0d-12, "thermE formula mismatch (local, dt=0)") + + deallocate(mNHC, rNHC, pNHC, r0, pN0) + end subroutine test_nhc_local_dt0 + + + subroutine test_nhc_local_reversible(error) + type(error_type), allocatable, intent(out) :: error + + integer :: M + real(8) :: dt, p, p_init, mass, beta, thermE, tol + real(8), allocatable :: mNHC(:), rNHC(:), pNHC(:), r_init(:), pN_init(:) + + tol = 5.0d-11 + M = 4 + dt = 5.0d-3 + beta = 0.9d0 + mass = 1.8d0 + p = 0.77d0 + + allocate(mNHC(M), rNHC(M), pNHC(M), r_init(M), pN_init(M)) + + mNHC = [ 0.8d0, 1.1d0, 1.4d0, 2.0d0 ] + rNHC = [ 0.02d0, -0.03d0, 0.01d0, 0.00d0 ] + pNHC = [ 0.05d0, 0.02d0, -0.04d0, 0.01d0 ] + + p_init = p + r_init = rNHC + pN_init = pNHC + + ! Forward step + call thermo_NHC_local( dt, M, p, mass, beta, rNHC, pNHC, mNHC, thermE) + ! Backward step (reversibility test) + call thermo_NHC_local(-dt, M, p, mass, beta, rNHC, pNHC, mNHC, thermE) + + call check(error, abs(p - p_init) <= tol, "local NHC not reversible: p did not return") + call check(error, maxval(abs(rNHC - r_init)) <= tol, "local NHC not reversible: rNHC did not return") + call check(error, maxval(abs(pNHC - pN_init)) <= tol, "local NHC not reversible: pNHC did not return") + + deallocate(mNHC, rNHC, pNHC, r_init, pN_init) + end subroutine test_nhc_local_reversible + + + subroutine test_nhc_global_dt0(error) + type(error_type), allocatable, intent(out) :: error + + integer :: M, ndim, natom, l, NTotDim + real(8) :: dt, beta, thermE, thermE_expected, tol + real(8), allocatable :: mNHC(:), rNHC(:), pNHC(:), r0(:), pN0(:) + real(8) :: p(3,2), p0(3,2), mass(2) + + tol = 1.0d-12 + M = 4 + ndim = 3 + natom = 2 + NTotDim = ndim*natom + + dt = 0.0d0 + beta = 0.8d0 + + mass = [ 12.0d0, 1.0d0 ] + p(:,1) = [ 0.10d0, -0.20d0, 0.30d0 ] + p(:,2) = [ -0.05d0, 0.04d0, -0.02d0 ] + p0 = p + + allocate(mNHC(M), rNHC(M), pNHC(M), r0(M), pN0(M)) + mNHC = [ 1.0d0, 1.3d0, 1.9d0, 2.7d0 ] + rNHC = [ 0.03d0, 0.01d0, -0.02d0, 0.04d0 ] + pNHC = [ -0.10d0, 0.07d0, 0.02d0, -0.05d0 ] + r0 = rNHC + pN0 = pNHC + + call thermo_NHC_global(dt, M, ndim, natom, p, mass, beta, rNHC, pNHC, mNHC, thermE) + + call check(error, maxval(abs(p - p0)) <= tol, "p changed with dt=0 in thermo_NHC_global") + call check(error, maxval(abs(rNHC - r0)) <= tol, "rNHC changed with dt=0 in thermo_NHC_global") + call check(error, maxval(abs(pNHC - pN0)) <= tol, "pNHC changed with dt=0 in thermo_NHC_global") + + thermE_expected = 0.0d0 + thermE_expected = thermE_expected + 0.5d0*pN0(1)*pN0(1)/mNHC(1) + dble(NTotDim)*r0(1)/beta + do l = 2, M + thermE_expected = thermE_expected + 0.5d0*pN0(l)*pN0(l)/mNHC(l) + r0(l)/beta + end do + + call check(error, abs(thermE - thermE_expected) <= 1.0d-12, "thermE formula mismatch (global, dt=0)") + + deallocate(mNHC, rNHC, pNHC, r0, pN0) + end subroutine test_nhc_global_dt0 + + + subroutine test_nhc_global_reversible(error) + type(error_type), allocatable, intent(out) :: error + + integer :: M, ndim, natom + real(8) :: dt, beta, thermE, tol + real(8) :: p(3,2), p_init(3,2), mass(2) + real(8), allocatable :: mNHC(:), rNHC(:), pNHC(:), r_init(:), pN_init(:) + + tol = 1.0d-10 + M = 4 + ndim = 3 + natom = 2 + + dt = 3.0d-3 + beta = 1.1d0 + + mass = [ 14.0d0, 16.0d0 ] + p(:,1) = [ 0.12d0, -0.08d0, 0.05d0 ] + p(:,2) = [ -0.03d0, 0.07d0, -0.09d0 ] + + allocate(mNHC(M), rNHC(M), pNHC(M), r_init(M), pN_init(M)) + mNHC = [ 0.9d0, 1.2d0, 1.5d0, 2.4d0 ] + rNHC = [ 0.01d0, -0.01d0, 0.02d0, 0.00d0 ] + pNHC = [ 0.03d0, -0.02d0, 0.01d0, 0.02d0 ] + + p_init = p + r_init = rNHC + pN_init = pNHC + + call thermo_NHC_global( dt, M, ndim, natom, p, mass, beta, rNHC, pNHC, mNHC, thermE) + call thermo_NHC_global(-dt, M, ndim, natom, p, mass, beta, rNHC, pNHC, mNHC, thermE) + + call check(error, maxval(abs(p - p_init)) <= tol, "global NHC not reversible: p did not return") + call check(error, maxval(abs(rNHC - r_init)) <= tol, "global NHC not reversible: rNHC did not return") + call check(error, maxval(abs(pNHC - pN_init)) <= tol, "global NHC not reversible: pNHC did not return") + + deallocate(mNHC, rNHC, pNHC, r_init, pN_init) + end subroutine test_nhc_global_reversible + + + subroutine test_nhc_global_mass0_bypass(error) + type(error_type), allocatable, intent(out) :: error + + integer :: M, ndim, natom + real(8) :: dt, beta, thermE, thermE0, tol + real(8) :: p(3,2), p0(3,2), mass(2) + real(8), allocatable :: mNHC(:), rNHC(:), pNHC(:) + + tol = 0.0d0 ! should be exactly unchanged due to early RETURN + M = 3 + ndim = 3 + natom = 2 + + dt = 1.0d-2 + beta = 1.0d0 + + mass = [ 12.0d0, 0.0d0 ] ! triggers bypass + p(:,1) = [ 0.40d0, 0.10d0, -0.20d0 ] + p(:,2) = [ 0.30d0, -0.50d0, 0.60d0 ] + p0 = p + + allocate(mNHC(M), rNHC(M), pNHC(M)) + mNHC = [ 1.0d0, 1.0d0, 1.0d0 ] + rNHC = [ 0.0d0, 0.0d0, 0.0d0 ] + pNHC = [ 0.0d0, 0.0d0, 0.0d0 ] + + thermE0 = 123.456d0 + thermE = thermE0 + + call thermo_NHC_global(dt, M, ndim, natom, p, mass, beta, rNHC, pNHC, mNHC, thermE) + + call check(error, maxval(abs(p - p0)) <= tol, "p changed even though mass=0 should bypass thermostat") + call check(error, thermE == thermE0, "thermE changed even though mass=0 bypass should return early") + + deallocate(mNHC, rNHC, pNHC) + end subroutine test_nhc_global_mass0_bypass + + + +end module test_thermo