Disagreement of potential energies calculated using DSF vs. Ewald

LAMMPS LAMMPS Mailing List Mirror
1

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

2

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.

it really is not an error. absolute potential energies between
different models have no meaning for the accuracy of a model. only
relative energy differences (or gradients/forces) matter. that being
said, most likely the major contribution to the difference you see is
the kspace energy term.

with any damped+cutoff interaction, you exclude a lot of contributions
to the total energy. in the real space part of ewald summation you do
the same thing, but you compute those interactions in the kspace part.
so if you compare epair for the two cases, you should get much closer
than by comparing pe.

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

this is indeed nonsense. the approximation is quite good. also, since
you use tabulation this is irrelevant. the tables are computed by
calling erfc() from the math library. the polynomial approximation is
only used if you turn off tabulation.

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:

3

Dr. Kohlmeyer,

Thanks for your reply. I'd just like to ask for some clarification. First, I don't know if I understand what you mean by 'using tabulation.' Next, I understand that only forces matter for integration, but I interpreted the large difference in potentials to mean that an error is present somewhere. I thought that the potential energies (and therefore total energies), as calculated by LAMMPS has physical meaning and could be a useful tool/sanity check for comparing to other results and real data. If this is not the case, would you say that agreement in radial distribution functions is sufficient to justify that using DSF instead of an Ewald summation is satisfactory for producing an accurate simulation?

You are correct that the kspace term is very different, but so is E_coul. Here are two sample outputs from the final step of the simulation for both simulations:

Using Ewald:
---------------- Step 24000 ----- CPU = 1020.5894 (sec) ----------------
TotEng = -7767.5605 KinEng = 1805.6629 Temp = 303.0328
PotEng = -9573.2234 E_bond = 0.0000 E_angle = 0.0000
E_dihed = 0.0000 E_impro = 0.0000 E_vdwl = 1539.5859
E_coul = 44630.1718 E_long = -55742.9812 Press = 632.5939

Using DSF:
---------------- Step 24000 ----- CPU = 1211.4333 (sec) ----------------
TotEng = -48516.6115 KinEng = 1813.2609 Temp = 304.3079
PotEng = -50329.8724 E_bond = 0.0000 E_angle = 0.0000
E_dihed = 0.0000 E_impro = 0.0000 E_vdwl = 1436.3031
E_coul = -51766.1755 E_long = 0.0000 Press = -15.1859

In Dr. Kast's paper that I linked before, a table of results gives the potential energies for various simulations and the Ewald and DSF energies are within 5% of each other. Does LAMMPS measure it differently or is it normal for the energies to not agree so well?

I just want to make sure that I've made no mistakes so that I have a proper justification for using DSF to save computation time. Thanks for your help.
Nathan

4

Dr. Kohlmeyer,

Thanks for your reply. I'd just like to ask for some clarification. First, I don't know if I understand what you mean by 'using tabulation.' Next, I understand that only forces matter for integration, but I interpreted the large difference in potentials to mean that an error is present somewhere. I thought that the potential energies (and therefore total energies), as calculated by LAMMPS has physical meaning and could be a useful tool/sanity check for comparing to other results and real data. If this is not the case, would you say that agreement in radial distribution functions is sufficient to justify that using DSF instead of an Ewald summation is satisfactory for producing an accurate simulation?

re tabulation: please see the pair_modify table option.

re accuracy: a good RDF by itself is not a good justification. or
rather it is difficult to see differences that are relevant. the rdf
has large variations at short range, but for liquids, its long range
structure is quickly averaged out. so that it will require very long
simulations to make any correlations visible over the noise.
furthermore, it only looks at structure averaged over spherical
shells. whether using approximations to the long-range coulomb create
large errors or not, depends a lot on the kind of system you simulate
and the kind of properties you want to compute. if you have a
homogeneous bulk system with a close to cubic simulation cell, then
the approximate coulomb causes the least errors. do you have
inhomogeneities or even an interface, then you are toast and have huge
errors because your coulomb approximation assumes a spherically
symmetric environment.

You are correct that the kspace term is very different, but so is E_coul. Here are two sample outputs from the final step of the simulation for both simulations:

Using Ewald:
---------------- Step 24000 ----- CPU = 1020.5894 (sec) ----------------
TotEng = -7767.5605 KinEng = 1805.6629 Temp = 303.0328
PotEng = -9573.2234 E_bond = 0.0000 E_angle = 0.0000
E_dihed = 0.0000 E_impro = 0.0000 E_vdwl = 1539.5859
E_coul = 44630.1718 E_long = -55742.9812 Press = 632.5939

Using DSF:
---------------- Step 24000 ----- CPU = 1211.4333 (sec) ----------------
TotEng = -48516.6115 KinEng = 1813.2609 Temp = 304.3079
PotEng = -50329.8724 E_bond = 0.0000 E_angle = 0.0000
E_dihed = 0.0000 E_impro = 0.0000 E_vdwl = 1436.3031
E_coul = -51766.1755 E_long = 0.0000 Press = -15.1859

In Dr. Kast's paper that I linked before, a table of results gives the potential energies for various simulations and the Ewald and DSF energies are within 5% of each other. Does LAMMPS measure it differently or is it normal for the energies to not agree so well?

I just want to make sure that I've made no mistakes so that I have a proper justification for using DSF to save computation time. Thanks for your help.

did you notice that your /dsf calculation is actually *slower* than
the one using ewald summation? and that one can be made even faster
and scale better with system size using pppm as kspace style without
losing any accuracy.

some quick profiling reveals that this is due to using the exponential
function a lot, which is expensive to compute (much more than the
approximated erfc() and even more than the tabulated coulomb). so with
LAMMPS you may actually be better off using coul/long with pppm and
then tune the coulomb cutoff for optimal performance.

if you want coul/dsf to run fast, you need to first implement an
efficient tabulation scheme similar as for coul/long. :wink:
...and then the benefit is probably still not worth the risk.

axel.

5

Dear Axel/Nathan,

Incidentally, a few days ago I carried out a similar comparison (9Dec14) between the Wolf method, DSF method and Ewald summation for SPC/E water and based on what I found I strongly agree that something is wrong with DSF. I just haven’t had time yet to post it myself.

First of all, for sufficiently large cutoff radii, DSF and Wolf essentially reduce to the same method and should therefore give very similar answers. For example, I use alpha = 7.2/22 for Wolf and therefore the energy shift in the Wolf method, erfc(11 * alpha)/11 = 1.e-8, is tiny. The problem we are talking about here is a factor 7 difference in the potential energy (DSF vs. Ewald/Wolf) for my system, while pressure and kinetic energies are well comparable. The RDFs for DFS, Wolf and Ewald are all spot on. The shifted force in the DSF method as compared to Wolf does not seem to have a noticeable effect on the dynamics for my system and the energies should be very similar as well.

In summary,

1.) Wolf and DSF should be very similar (at least for my settings),

2.) Wolf and Ewald give almost the same answer for the potential energy,

3.) There is a factor 7 difference in the potential energy between DSF and Wolf,

4.) The Wolf implementation was suffering from exactly the same problem about a year ago, before Axel managed to fix it. In November 2013, I reported the same problem for Wolf, namely a factor 8 energy difference as compared to Ewald, which turned out to be a bug (http://lammps.sandia.gov/threads/msg42566.html). Fortunately Axel found it soon:

"…

using the “inexact” version of Force::pair_match() will catch all of
those cases (note the “0” instead of “1”):

diff --git a/src/neighbor.cpp b/src/neighbor.cpp
index 10f7ce0…0fe4672 100644
— a/src/neighbor.cpp
+++ b/src/neighbor.cpp
@@ -320,7 +320,9 @@ void Neighbor::init()
special_flag[3] = 1;
else special_flag[3] = 2;

  • if (force->kspace) special_flag[1] = special_flag[2] = special_flag[3] = 2;
  • if (force->kspace || force->pair_match(“coul/wolf”,0)
  • || force->pair_match(“coul/dsf”,0))
  • special_flag[1] = special_flag[2] = special_flag[3] = 2;

// maxwt = max multiplicative factor on atom indices stored in neigh list

however, this kind of change is indeed a bit difficult to track down.
at the very least some kind of comment near the constructor explaining
that an update of neighbor.cpp may be needed for derived classes
should be added. but something similar to ewaldflag, pppmflag or
tip4pflag would probably be a cleaner solution.

…"

Axel, is there a chance that the current problem is related to what you found last time?

Best regards,

Peter

6

dear peter,

Dear Axel/Nathan,

Incidentally, a few days ago I carried out a similar comparison (9Dec14)
between the Wolf method, DSF method and Ewald summation for SPC/E water and
based on what I found I strongly agree that something is wrong with DSF. I
just haven't had time yet to post it myself.

First of all, for sufficiently large cutoff radii, DSF and Wolf essentially
reduce to the same method and should therefore give very similar answers.
For example, I use alpha = 7.2/22 for Wolf and therefore the energy shift in
the Wolf method, erfc(11 * alpha)/11 = 1.e-8, is tiny. The problem we are
talking about here is a factor 7 difference in the potential energy (DSF vs.
Ewald/Wolf) for my system, while pressure and kinetic energies are well
comparable. The RDFs for DFS, Wolf and Ewald are all spot on. The shifted
force in the DSF method as compared to Wolf does not seem to have a
noticeable effect on the dynamics for my system and the energies should be
very similar as well.

In summary,

1.) Wolf and DSF should be very similar (at least for my settings),

2.) Wolf and Ewald give almost the same answer for the potential energy,

3.) There is a factor 7 difference in the potential energy between DSF and
Wolf,

4.) The Wolf implementation was suffering from exactly the same problem
about a year ago, before Axel managed to fix it. In November 2013, I
reported the same problem for Wolf, namely a factor 8 energy difference as
compared to Ewald, which turned out to be a bug
(LAMMPS Molecular Dynamics Simulator). Fortunately Axel found it
soon:

well, as it turns out, it is not exactly the same problem, but rather
the inverse. while the neighbor lists are set up correctly as
indicated by the diff below, the coul/dsf styles are not computing
energies and forces correctly with exclusions. the solution is
straightforward:

diff --git a/src/pair_coul_dsf.cpp b/src/pair_coul_dsf.cpp
index 17b3dfa..51cf676 100644
--- a/src/pair_coul_dsf.cpp
+++ b/src/pair_coul_dsf.cpp
@@ -65,7 +65,7 @@ void PairCoulDSF::compute(int eflag, int vflag)
   int i,j,ii,jj,inum,jnum;
   double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
   double r,rsq,r2inv,forcecoul,factor_coul;
- double prefactor,erfcc,erfcd,e_self,t;
+ double prefactor,erfcc,erfcd,t;
   int *ilist,*jlist,*numneigh,**firstneigh;

   ecoul = 0.0;
@@ -97,7 +97,7 @@ void PairCoulDSF::compute(int eflag, int vflag)
     jnum = numneigh[i];

     if (eflag) {
- e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e;
+ double e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e;
       ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0);
     }

@@ -115,14 +115,15 @@ void PairCoulDSF::compute(int eflag, int vflag)
         r2inv = 1.0/rsq;

         r = sqrt(rsq);
- prefactor = factor_coul * qqrd2e*qtmp*q[j]/r;
+ prefactor = qqrd2e*qtmp*q[j]/r;
         erfcd = exp(-alpha*alpha*rsq);
         t = 1.0 / (1.0 + EWALD_P*alpha*r);
         erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
         forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd +
                                  r*f_shift) * r;

7

Dear Axel,

Thanks for looking into this and resolving the problem so quickly!

Best regards,
Peter