I am new to molecular dynamics and have a quick question…
I have been running a few simulations in the NPT ensemble, and calculating the density of my liquid system using LAMMPS in order to validate my force field. Using the same simulation conditions as in the literature:
- 400 molecules (this makes for a pretty large system, ~40 angstrom cube)
- 11 angstrom cutoff
- Ewald sum approximation for long range coulombic interactions
- long range tail corrections for lj interactions
- OPLS 1,4 interactions
I have 0.7 ns equilibration and 0.7 production runs.
For the production run, I obtain an average density of 1.116 g/cm3, whereas the published density is 1.123+/-0.002 g/cm3 for the system & force field in question. I have checked multiple times my force field parameters, and I am sure there are no errors.
I know that this is a very small error compared to the published value, but I am unsure if this is a product of using a different simulation platform (I believe that the publication authors used DL_POLY, and not lammps), or if it necessarily indicates something is amiss in my simulation set-up…what do you guys think?
Thanks for your time !
Here is my run file, if it helps…
boundary p p p
pair_style lj/cut/coul/long 11.0
pair_modify tail yes mix arithmetic
kspace_style ewald 1e-4
special_bonds lj/coul 0.0 0.0 0.5 #1,4 interaction scaling by 0,5
neighbor 2.0 bin
thermo_style custom step press
fix 1 all nve/limit 0.1
fix 2 all temp/rescale 10 298 298 0 1
fix 3 all recenter INIT INIT INIT units box
fix 4 all shake 0.0001 20 100000 b 5 6 8 9 11 #fix certain bond lengths
fix 1 all npt temp 298 298 100 iso 0.98692 0.98692 500 #298 K and 1 bar