Hello,
I have been having some trouble with two different types of simulations producing different results that I don't think should be different. I plotted some of the LAMMPS outputs for simulations of 1000 rigid, TIP3P water molecules (density = 0.998 kg/m^3, T = 300K, NVT simulations) using DSF and Ewald summation for comparison. I uploaded them to imgur here: http://imgur.com/a/2LeMb. As you can see, there is a huge discrepancy between the potential energies. I can't identify the source of this error.
The context is that I'm trying to show that using DSF is a valid approximation to replace the Ewald summation for a particular simulation. Ultimately I'd like to show this is the case for a simulation of tetrolic acid (aka 2-butynoic acid) dissolved in CCl4. However, the potential energies were always off by a factor of 3 (ish) between the two methods. So first I contacted the Computational Chemistry forum and asked what might be wrong and Dr. Stefan Kast replied to me and suggested I try to replicate a simulation performed for this paper that he co-authored: http://pubs.acs.org/doi/abs/10.1021/jp025949h. I attempted to follow his methods and replicate the results in LAMMPS. His work managed to achieve less than 5% error using DSF instead of Ewald summation with a real-space cutoff of 12 angstroms and a Coulomb cutoff of 9. (I used 9 and 9) I got very good-looking RDF plots, as you can see from the imgur album, but the potential energies were still wrong in the same way they were for TTA and CCl4. He suggested there could be a small error that added up for lots of atoms stemming from the use of an approximation for the erfc() instead of the real thing. He modified the appropriate code and I recompiled LAMMPS with that and reran my simulations. The plots for comparison are also in the imgur album, and obviously that is not the issue. This is my input file:
units real
boundary p p p
atom_style full
read_data data.TIP3P
group water type 1 2
# interaction styles: comment or uncomment the following lines to choose a style
pair_style lj/cut/coul/long 9.0 9.0
kspace_style ewald 1.0e-4
# pair_style lj/cut/coul/dsf 0.2 9.0 9.0
bond_style harmonic
angle_style harmonic
dihedral_style opls
improper_style harmonic
pair_modify mix geometric tail no
#field, atom type, atom type, eps(kcal/mol), sigma(angstroms) [http://lammps.sandia.gov/doc/pair_coeff.html]
pair_coeff 1 1 0.0000 0.0000 # H
pair_coeff 2 2 0.1521 3.1500 # O
#field, bond type, kb, r0 [http://lammps.sandia.gov/doc/bond_coeff.html]
bond_coeff 1 450 0.9572 # O-H (this command should be irrelevant for a rigid model)
#field, angle type, kb, theta0 [http://lammps.sandia.gov/doc/bond_coeff.html]
angle_coeff 1 55 104.52 # H-O-H (this command should be irrelevant for a rigid model)
# http://lammps.sandia.gov/doc/Section_howto.html#howto_7
run_style verlet
velocity all create 300.0 239847 dist gaussian mom no rot no
fix 1 water nvt temp 300.0 300.0 100.0
timestep 0.25 #resets timestep to 0.25 fs
run 100
minimize 1.0e-5 1.0e-6 10000 100000
write_restart water-min.restart
fix 2 water shake 0.0001 20 0 b 1 a 1
neigh_modify every 1 delay 1 check yes
thermo 20
thermo_style multi
velocity all scale 300
dump 1 all dcd 20 water-eq.dcd
dump_modify 1 unwrap no
compute rdf_water all rdf 300 1 2 1 1 2 2
fix 3 all ave/time 20 300 8000 c_rdf_water file rdf_water.rdf mode vector
reset_timestep 0
run 24000
write_restart water-eq.restart
Thanks in advance for any suggestions,
Nathan