Pressure oscilations in NPT ensemble

LAMMPS LAMMPS General Discussion
1

I am simulating gas hydrate crystals in LAMMPS using fix npt. When I specify lower pressure of 20 atm, the average pressure of the system is initially at the range of 20 atm but later during the simulation it socialites a lot (it even becomes negative sometimes). I even used drag 2.0 but it did not solve my problem. Anybody has any idea how I can damp those oscillations effectively? Note that I do not face this problem for higher pressure of 100 atm.

Below is my input file:

#----- conversion factors and coefficients
variable ps2fs equal 1e3
variable ns2fs equal 1e6
variable ms2fs equal 1e9
variable fs2s  equal 1e-15
variable atm2Pa equal 101325.0
variable A2m equal 1.0e-10
variable rnd equal round(random(0,999,${RANDOM})) # random number for velocity initialization

# Number of unit cell replicas
variable Nx equal 20
variable Ny equal 10
variable Nz equal 10


# time-step
variable dt equal 10.0


#----- Set time duration for each ensemble integration
variable tNVT equal 500 # in ps
variable N_NVT equal floor(${tNVT}*${ps2fs}/${dt})
variable tNPT equal 0.001 # in ms
variable N_NPT equal floor(${tNPT}*${ms2fs}/${dt})
variable tNVT_fixedHydrate equal 0.100
variable N_NVT_fixedHydrate equal floor(${tNVT_fixedHydrate}*${ms2fs}/${dt})

#----- set thermodynamic properties
variable Thyd equal 250
variable Phyd equal 493
variable Tdiss equal 450
variable Tdamp equal $(100*v_dt)
variable Pdamp equal $(1000*v_dt)

#----- path to the forcefield folder, unit-cell folder, and the directory to the property calculation files
variable uc_folder string "/scratch/madibi/MDSim/co2Hydrate/cg_model/unit_cell/"
variable ff_folder string "/filtered_ff_candidates"
variable props_in_folder string "../../../propsCalc/"  # not using it currently

# file names for uc
variable uc_fileName string   s1_unit_cell_mw_co2.data


#-------------------------- Set global settings of the simulations --------------------
# time-step
timestep ${dt}


# set neighboring parameter (skin)
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
neigh_modify one 5000


#---------------------------------- Init Section ----------------------------------
units                   real
dimension               3
boundary                p p p
atom_style              full
pair_style              sw
#---------------------------------- Initial Configuration ----------------------------------

# Create the simulation box or load from previous restart run.
read_data ${uc_folder}${uc_fileName}

replicate ${Nx} ${Ny} ${Nz}



# Initialize the velocity if needed
velocity all create ${Thyd} ${rnd}  mom yes rot yes dist gaussian loop geom

# Conserve the linear and angular momentums
fix     removeMomentum all momentum 1 linear 1 1 1 angular

#---------------------------------- Forcefield parameters Section ----------------------------------

# load forcefield parameters
pair_coeff * * water_co2_base.sw   O C


# write-out the initial configuration of the molecules
write_data initial_config.data


#---------------------------------- Energy Minimization (only if num_rest_curr==0) ----------------------------------

thermo 1
thermo_style custom fmax fnorm

minimize 1.0e-5 1.0e-10 1000000  4000000
reset_timestep 0

write_data AfterMinimization.data

#-------------------------- s,p,d for time-averaged calculations --------------------------
variable s equal 10    # sample interval
variable p equal 100   # correlation length
variable d equal $s*$p # dump interval to calculate thermodynamic properties

variable time equal time
variable Vout equal vol
variable Pout equal press
variable Tout equal temp
variable Hout equal enthalpy
variable PEout equal pe
variable KEout equal ke
variable ETOTout equal etotal
variable massDensity equal density

# global scalar properties
fix GlobalPropsTimeAvg  all ave/time $s  $p  $d  v_Pout v_Tout v_Vout  v_Hout v_PEout v_KEout v_ETOTout v_massDensity  file GlobalPropsTimeAvg.prop

# vector properties
compute myRDF all rdf  50 1 1 2 2 1 2
fix rdfOut  all ave/time $s $p $d  c_myRDF[*] file rdf_all.rdf mode vector

#-------------------------- Simulations --------------------------
# Setup thermo output.
thermo_style     custom step pe vol press temp
thermo           500

##### NVT relaxation
print "****************************************"
print "NVT relaxation at Thyd=${Thyd}"
print "****************************************"

dump    trjNVT     all custom 20000  NVT_Trajectory.lammpstrj    id mol type x y z ix iy iz  vx vy vz

dump_modify trjNVT header yes time yes flush yes

fix NVTequil all nvt temp ${Thyd} ${Thyd} ${Tdamp} drag 0

run ${N_NVT}

unfix NVTequil
undump trjNVT

write_restart restart.AfterNVT
write_data    data.AfterNVT

##### NPT for pressure equilibration
print "****************************************"
print "NPT hydrate at Thyd=${Thyd} and Phyd=${Phyd}"
print "****************************************"

dump    trjNPT     all custom 20000  NPT_Trajectory.lammpstrj    id mol type x y z ix iy iz  vx vy vz

dump_modify trjNPT header yes time yes flush yes

fix NPTSim all npt temp ${Thyd} ${Thyd} ${Tdamp} iso ${Phyd} ${Phyd}  ${Pdamp} drag 0

run ${N_NPT}

unfix NPTSim
undump trjNPT

write_restart restart.AfterNPT
write_data    data.AfterNPT


##### Increase temperature in the water_methane zone (NVT)
print "****************************************"
print "NVT equilibration to partially melt hydrate at T=${Tdiss}"
print "****************************************"

dump    trjNVT_eql     all custom 20000  NVT_eql_Trajectory.lammpstrj    id mol type x y z ix iy iz  vx vy vz
dump_modify trjNVT_eql header yes time yes flush yes
# Set different zones in the simulation box
variable DX equal (xhi-xlo)/2
variable Xi_hyd equal xlo+0.25*(xhi-xlo)
variable Xf_hyd equal ${Xi_hyd}+${DX}


region hydrate_reg block ${Xi_hyd} ${Xf_hyd} INF INF INF INF
region water_methane_left  block INF  ${Xi_hyd} INF INF INF INF
region water_methane_right  block ${Xf_hyd} INF  INF INF INF INF
region water_methane_reg  union 2  water_methane_left water_methane_right


group HydrateAtoms region hydrate_reg
group WaterMethaneAtoms region water_methane_reg

fix NVT_fixedHydrate WaterMethaneAtoms nvt temp ${Tdiss} ${Tdiss} ${Tdamp} drag 0.0

run ${N_NVT_fixedHydrate}

unfix NVT_fixedHydrate

undump trjNVT_eql

write_restart restart.AfterNVTFixedHydrate_${Nx}_${Ny}_${Nz}

write_data data.AfterNVTFixedHydrate_${Nx}_${Ny}_${Nz}
2

It is normal to see fluctuations, as long as the average pressure oscillates around the target value.
Generally, the smaller the sample, the larger the oscillations.
What are the units, timestep, and damping constants for the fix npt?

3

The problem is that the average pressure oscilates a lot. By a lot I mean it goes all the way to -5 atm for the set pressure of 20 atm.

I am using real units, tstep is 10 fs, Tdamp and Pdamp for npt are 1000 and 10000 fs respectively. I just realized the code was not fully copied so I edited it again and you can find all these information in my lammps script.

4

The magnitude of pressure fluctuations depends strongly on compressibility, system size, and choice of time step. Condensed systems are not very compressible and solid ones even less so.

There are plenty of discussions on the subject archived in the forum that may be worth looking through.

A time step of 10 fs seems quite aggressive.

5

Thanks.

Actually I am simulating a three-phase liquid-solid-gas to find the equilibrium configuration using the NPT coexistence method. Do you think large pressure fluctuations affects the attained equilibrium point?

This is a coarse-grained SW potential in which every molecule is simulated with a single bead. We using the time-step size that is suggested by the researchers who developed the force field for this system.

6

I prefer not to speculate with so little tangible information.

Still seems aggressive. A single water bead would have a mass of 18 amu, the potential would have to be extremely soft (and thus not very realistic) to have a stable time integration at 10fs. In coarse grain models, often a larger error in time integration is often accepted since the model has a larger error through coarse graining already. You seem to be subjecting your system to a large external pressure and that would create a denser system and would require reducing the timestep. Has the model been parameterized for these conditions? I would suspect not.

The fact that you need to increase your neighbor one parameter this much is a sign for a very compressed system.

7

If I understand correctly, you are calculating several average values of pressure throughout the dynamics by usign a given amount of configurations in each of these calculations and then comparing these average values with one another.

So, additionally to the things that have been said, I wonder if the problem isnt simply the amount of configurations you are considering to compute the averages of the pressure. Maybe you are using very few and breaking the ergodicity principle. Have you tried computing these averages using more data?