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