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/modules/ThermoModule.f90 b/src/modules/ThermoModule.f90 index d8036dd..f8e7f63 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,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 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