Density

LAMMPS LAMMPS Mailing List Mirror
1

Hello,

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…

dimension 3
boundary p p p
units real
atom_style full
pair_style lj/cut/coul/long 11.0
pair_modify tail yes mix arithmetic
kspace_style ewald 1e-4

bond_style harmonic
angle_style harmonic
dihedral_style opls
special_bonds lj/coul 0.0 0.0 0.5 #1,4 interaction scaling by 0,5

read_data data.LAMMPS
neighbor 2.0 bin

thermo 5000
thermo_style custom step press
timestep 1

#equilibration
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
run 5000

fix 4 all shake 0.0001 20 100000 b 5 6 8 9 11 #fix certain bond lengths
run 10000

unfix 1
unfix 2

fix 1 all npt temp 298 298 100 iso 0.98692 0.98692 500 #298 K and 1 bar
run 700000

#production run
run 700000

2

Hi Thomas

If you let me know what the molecule is, I can very easily build it independently in MedeA, assign the OPLS forcefield and run a calculation identical to yours. This will give an independent check of your results. If our results differ, we can compare the LAMMPS files to see if in spite of your care something is wrong with the forcefield parameters or system.

Paul

psaxe@…1795…

Phone: +1.760.495.4924 x 102
Cell: +1.505.603.8182

This email is covered by the Electronic Communications Privacy Act, 18 U.S.C. Section 2510-2521. This message and any attachments hereto may contain confidential information intended only for the use of the individual or entity named above. If you are not the intended recipient(s), or the employee or agent responsible for delivery of this message to the intended recipient(s), you are hereby notified that any dissemination, distribution or copying of this email message is strictly prohibited. If you have received this message in error, please immediately notify the sender and delete this email from your computer. The sender does not waive any privilege in the event this message was inadvertently disseminated.

MedeA® and Materials Design® are registered trademarks of Materials Design, Inc.

3

I would question the error +/-0.002 in the literature. 0.2 % fluctuations in the density seem quite small to me.