Help understanding DSF calculated potential

Hello,

I am attempting to do some analysis on outputs predicted by lammps using the DSF potential. The e_coul values predicted in lammps looked a little strange to me, so I came up with a toy case where I could calculate the energy by hand. I simply placed two atoms 4 angstroms apart with no bonded interactions or vdw interactions. I made the charge of one atom 1e, and the other 2e. Using a box of length 1000.0 angstroms, a cutoff of 15.0 angstroms, and a damping parameter of 0.22 angstoms^(-1), and remembering to multiply by the appropriate 332 to convert to kcal/mol, I hand calculate an energy of e_coul = 35.1. However, lammps seems to be calculating -170.6715. Can anyone explain to me this? Is e_coul modified here anywhere besides pair_lj_cut_coul_dsf.cpp?

Numbers for the hand DSF calculation:

term1 = erfc(alphar)/r = 0.053
term2 = -erfc(alpha
rcut)/rcut = -2.03810^(-7)
term3 = erfc(alpha
rcut)/rcut/rcut = 1.358910^(-8)
term4 = 2
alpha/sqrt(pi) * exp(-alphaalpharcutrcut)/rcut = 3.08510^(-7)

e_coul = 33212*(term1+term2+(term3+term4)*(4-15)) = 35.1

Here is the data file

LAMMPS data file. CGCMM style. atom_style full generated by VMD/TopoTools v1.2 on Wed Jul 01 22:13:52 EDT 2015
2 atoms
0 bonds
0 angles
0 dihedrals
0 impropers
1 atom types
0 bond types
0 angle types
0 dihedral types
0 improper types
0.0000 1000.00 xlo xhi
0.0000 1000.00 ylo yhi
0.0000 1000.00 zlo zhi

Pair Coeffs

This is because DSF contains the self-energy terms, which are absent in your manual calculation.

You can look at the source code (pair_lj_cut_coul_dsf.cpp) lines 115-116 for the self energy term calculation. Commenting out line 116 (where ev_tally() is called to accumulate e_self into ecoul) and rebuilding LAMMPS give you the result without the self-energy terms (i.e. 35.414).

-Trung