Hello LAMMPS users,
I am trying to reproduce experimental densities of methane, ethane and water using LAMMPS in order to validate the models that I am using. I am using parameters exactly as they appear from the OPLS-AA force field (for methane and ethane) and SPC/E water model with simulation boxes of 1000 molecules. Unfortunately, the densities I am obtaining from NPT are consistently low by a lot (70-90%). I have included my input for methane below. Is there something that I am missing in the neighboring or cutoffs that could cause this effect? I ran these models in GROMACS and obtained the correct densities. Any suggestions would be greatly appreciated
Thanks,
Evan
----------------- Init Section -----------------
units real
atom_style full
pair_style lj/charmm/coul/long 11.0 12.0 12.0
bond_style harmonic
angle_style harmonic
dihedral_style opls
pair_modify mix arithmetic
kspace_style pppm 0.0001
boundary p p p
----------------- Atom Definition Section -----------------
read_data “system_test.data”
----------------- Settings Section -----------------
bond_coeff 1 340.0000 1.0900 # CT - HC
angle_coeff 1 33.0000 109.5000 # HC - CT -HC
pair_coeff 1 1 0.06600 3.50000 # opls-CT - opls-CT
pair_coeff 2 2 0.03000 2.50000 # opls-HC - opls-HC
----------------- Run Section -----------------
timestep 0.1
dump 1 all custom 1000 traj_nvt_ngpu.lammpstrj id mol type x y z ix iy iz
thermo_style custom cpuremain temp press density
thermo 100
thermo_modify flush yes
fix 1 all npt temp 300 300 10 iso 10 10 10
restart 100 sim_ngpu.restart sim_ngpu.restart
run 10000000
write_data NPT.data
Hello LAMMPS users,
I am trying to reproduce experimental densities of methane, ethane and water
using LAMMPS in order to validate the models that I am using. I am using
parameters *exactly* as they appear from the OPLS-AA force field (for
methane and ethane) and SPC/E water model with simulation boxes of 1000
molecules. Unfortunately, the densities I am obtaining from NPT are
consistently low by a lot (70-90%). I have included my input for methane
below. Is there something that I am missing in the neighboring or cutoffs
that could cause this effect? I ran these models in GROMACS and obtained the
correct densities. Any suggestions would be greatly appreciated
difficult to say. are you certain, that you are really comparing
apples to apples here? i.e. are the settings exactly the same.
you are using a very coarse convergence for your kspace forces. have
you tested against using, say 1.0e-6?
you are doing isotropic cell deformation, did you do the same for
gromacs. did you use nose-hoover temperature and pressure coupling?
and the same time constants? have you checked the units? partial
charges?
does you system significantly shrink or expand after reading the data file?
axel.
Also check if you need to add (or not) the tail correction
for the LJ term. LAMMPS does that via the pair_modify command.
If you can do it, the simplest test is to take one snapshot
and run it with both codes to determine energy, force/atom,
and pressure. If the last one doesn’t agree for a snapshot,
it likely won’t equilibrate to the same density with NPT.
Steve
One other thought. 0.1 fs is a tiny timestep,
1.0 would be more typical. In which case
using 10 for the tstat and bstat time constants is pretty small.
Values of 100 fs and 1000 fs would be more typical.
Steve
Steve and Axel,
Thank you for your helpful suggestions. The problem was indeed the timestep. Once changed the timestep and equilibrated for a full 10ns I obtained a more accurate density. In the future I will pay closer attention to the barostat timestep since it seems to be a major factor in the equilibration time.
Evan L.
Steve and Axel,
Thank you for your helpful suggestions. The problem was indeed the timestep.
Once changed the timestep and equilibrated for a full 10ns I obtained a more
accurate density. In the future I will pay closer attention to the barostat
timestep since it seems to be a major factor in the equilibration time.
if you know your target density, it is almost always a good idea to
start equilibration with a fix volume.
in fact, this is almost always a good idea. starting with fix npt on a
non-equilibrated configuration will often make equilibration
needlessly time consuming.
initial configuration frequently have a high potential energy, which
will convert (quickly) into kinetic energy, which will lead to an
expansion of your system.
once the system has expanded, it will take a very long time to get
back, since atom/molecules will have to rearrange to get more densely
packed. that is not - unlike expansion - a highly probability process,
so it requires time until the existing local fluctuations "make room"
for it. expansion, on the other hand, is almost always easily
possible, and doesn't require much of a rearrangement. if you look at
grand canonical monte carlo simulations, you have a conceptually
similar problem, especially for molecules: i.e. inserting a molecule
is difficult, since it is unlikely that you by random choice quickly
find a suitable location and orientation; removing a particle doesn't
have that problem.
for all but the simplest systems, i would also recommend a three (or
more) step equilibration, using first a dissipative thermostat (e.g.
fix langevin) to facilitate equipartitioning of kinetic energy. since
fix nvt only looks at the global temperature and thus depends on the
thermal conductivity of the material to redistribute kinetic energy
properly.
HTH,
axel.