Dear Lammps-Users,

I am interested in computing the value of the electric field on any given atom in my simulation. I am looking at a simple system of 2 spc/e water molecules located near the center of a cubic box with length = 200 Angstroms.

I compare the following two numbers:

(1) Force and Efield computed from the definition of force and efield due to point charges. (Thus a sum over all charges in the system divided by either r^2 or r, respectively).

(2) Force and Efield computed from the force calculation based on the pair potential. I keep track of the forces due to Lennard-Jones and SHAKE and subtract those from the total force to give me (what I presume to be) the Coulombic part of the force: forceCOUL = forceTOTAL - forceLJ - forceSHAKE

(NOTE: special_bonds is set to special_bonds 1 1 1 1 1 1.)

When I compare values from (1) and (2) using pair lj/cut/coul/cut with cutoff 20, I get the same values. Number (2) is in fact not sensitive to the cutoff length as long as it's larger than about 3.

When I compare the values from (1) and (2) using pair lj/cut/coul/long with cutoff 20, the numbers differ. I assumed this was because of the long-ranged forces, but when I keep track of the long-ranged forces (in ewald.cpp) and get forceCoulomb from forceCoul = forceTOTAL - forceLJ - forceSHAKE - forceLONG, the results still do not match up.

My questions are

1. Is there an additional contribution to the force in lj/cut/coul/long that I am neglecting? What is it, and where is it coming from?

2. Units. Forces in LAMMPS are in kcal/mol-Angstrom. I am interested in the electric field in units Volts/Angstrom = J/C-Angstrom. Since charge is basically a number in LAMMPS (multiples of electron charge), I need a multiplicative factor alpha so that the units work out

Efield = alpha * forceCoul / q

alpha = 0.24 / e Na = 2.49e-6

where e = electron charge in Coulombs, Na = Avogadro's number, q is the charge of the atom of interest. This gives me values that are about 1.0e4 off in magnitude from what I expect.

I've gone through the unit conversions and I don't see my mistake. Any help here would be much appreciated as well!

Thanks in advance for any help!

Joyce