How to output electric forces in LAMMPS

(WRITTEN IN 2010, PROBABLY COMPLETELY OUT OF DATE)

How to make LAMMPS output electric fields. (Similar method should output any given contribution to the force.)

(Written by Steve Byrnes, following advice from Joyce Noah-Vanhoucke.)

For various purposes, you may want to know not just the total force on an atom in an MD simulation, but one of the contributions to the force, for example the electric-force contribution or the Lennard-Jones contribution or whatever. Personally, my goal was to calculate the electric force experienced by each hydrogen atom in SPC/E liquid water, as this can be used to (pretty accurately) infer the OH vibrational frequency of that bond in an isotopically diluted environment (ref). (Obviously, if you know the electric force on a hydrogen atom, then you know the electric field it is in, and vice-versa.) This document says how I did that. I'm writing only the specific thing I did, but with some work it could presumably be adjusted to work for different contributions to the force, different simulation setups, different versions of the software, etc.
I am using...

"cut/coul/long" pair-style
"full" atom-style
"verlet" (not rRESPA) integration
Ewald (not pppm).
15 January 2010 version of the software

I am not a LAMMPS expert. The code passed my tests, but use at your own risk! (It's also possible that there are much easier ways to do this. See thread on LAMMPS email list.)

Good luck!

atom.h

Add into "public" area:

double **fcoul;  // my code

"atom->fcoul" will be the Coulomb component of the force felt by the atom. This is modelled on "atom->f", the total force on the atom, which is stored and accessed in the same way. Note that you could define, store, and extract the other force components (LJ force, SHAKE force, etc.) in a similar way if you wanted to.

atom.cpp

Add a line into the initialization (i.e. Atom::Atom(LAMMPS *lmp) : Pointers(lmp)):

x = v = f = NULL;
fcoul = NULL; // my code

Add a line into the destruction (i.e.Atom::~Atom):

memory->destroy_2d_double_array(f);
memory->destroy_2d_double_array(fcoul); // my code

atom_vec_full.h

Add into "private" area (after "double **x,**v,**f;"):

double **fcoul; // my code

atom_vec_full.cpp

In "void AtomVecFull::grow(int n)", add an extra line:

f = atom->f = memory->grow_2d_double_array(atom->f,nmax,3,"atom:f");
fcoul = atom->fcoul = memory->grow_2d_double_array(atom->fcoul,nmax,3,"atom:fcoul"); // my code

In "double AtomVecFull::memory_usage()", add an extra line:

if (atom->memcheck("f")) bytes += nmax*3 * sizeof(double);
if (atom->memcheck("fcoul")) bytes += nmax*3 * sizeof(double); // my code

In void AtomVecFull::grow_reset(), add an extra line:

x = atom->x; v = atom->v; f = atom->f;
fcoul = atom->fcoul; // my code

There's two more places where "f" is used in atom_vec_full.cpp: pack_reverse and unpack_reverse. These shouldn't get called, as explained below. In fact, to be safe, you can put a little warning at the start of these two functions:

fputs("atom_vec_full.cpp: Warning: E-fields may not be accurate! Did you remember to turn off newton in your script?\n", screen);

verlet.cpp

In void Verlet::force_clear() (which resets the force on atoms to zero), wherever it resets a force, it should also reset fcoul. There should be three places in the program where you modify it to look like:

f[i][0] = 0.0;
f[i][1] = 0.0;
f[i][2] = 0.0;
// my code
atom->fcoul[i][0] = 0.0;
atom->fcoul[i][1] = 0.0;
atom->fcoul[i][2] = 0.0;
// end my code

(f was defined earlier in the program as atom->f).

Note on Coulomb force calculation

If you look at verlet.cpp, the Coulomb force is calculated in two places: The short-range part is calculated in the command pair->compute, and the long-range part in the command kspace->compute. In both cases, we will "intercept" the Coulomb force data as it gets written into atom->f, so that we can also write it into atom->fcoul.

pair_lj_cut_coul_long.cpp

Here is where the short-range part of the Coulomb force is calculated.

There are four main functions here, compute, compute_inner, compute_middle, compute_outer. The last three are for rRESPA and don't get used in Velocity Verlet, so ignore them. You only need to modify compute (more specifically, void PairLJCutCoulLong::compute(int eflag, int vflag)).

You'll see the following code to store the total short-range force:

fpair = (forcecoul + factor_lj*forcelj) * r2inv;

f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {
  f[j][0] -= delx*fpair;
  f[j][1] -= dely*fpair;
  f[j][2] -= delz*fpair;
}

(f was defined earlier in the program as atom->f). There are three things to notice:

The program makes it clear how fpair is the sum of a Coulomb force and an LJ force. Therefore it's easy to just extract the Coulomb part.
Newton's third law is being used: The force on j is minus the force on i. In LAMMPS, Newton's third law is always invoked if j is on the same processor as i (i.e., j<nlocal), and it is also always invoked if the LAMMPS script has the "newton" flag turned on (which is the LAMMPS default).
When forcecoul is calculated (just earlier in the program), it uses factor_coul which comes from special_coul which comes from the LAMMPS special_bonds command. The default for special_bonds is that atoms directly or indirectly bonded do not have their Coulomb forces included. (So if factor_coul=0, the pairwise Coulomb force should be zero, and in fact forcecoul is set so that this short-range force cancels out the contribution to the Ewald-sum long-range force.) If you instead want the total Coulomb force, including Coulomb's law applied to bonded atoms, you need to change the pair_special setting in your LAMMPS script, or else edit the forcecoul calculation just above to have it store extra information. (factor_lj is a similar thing.)

As mentioned, just drop out the LJ part and keep the Coulomb part, and you should be able to store the short-range Coulomb force. So put in the following code:

// my code: fcoul stores the Coulomb contribution to force
atom->fcoul[i][0] += delx * forcecoul * r2inv;
atom->fcoul[i][1] += dely * forcecoul * r2inv;
atom->fcoul[i][2] += delz * forcecoul * r2inv;
if (newton_pair || j < nlocal) {
  atom->fcoul[j][0] -= delx * forcecoul * r2inv;
  atom->fcoul[j][1] -= dely * forcecoul * r2inv;
  atom->fcoul[j][2] -= delz * forcecoul * r2inv;
}
// end my code

Also, in your LAMMPS input script, you must turn off "newton":

newton off

If you want to leave it on, you need to figure out how to modify void Comm::reverse_comm() in comm.cpp, including its subroutines in atom_vec_full.cpp (i.e., AtomVecFull::unpack_reverse and AtomVecFull::pack_reverse). This processor-to-processor communication function corrects for the fact that part of the force on some of the atoms was calculated by the wrong processor when newton is turned on. Therefore if you use newton, the function needs to be modified to also correctly communicate fcoul. But I didn't bother to figure out the details: I'm using serial mode so it doesn't matter for me.

ewald.cpp

In void Ewald::compute(int eflag, int vflag), here are the lines where the Coulomb force is being stored:

f[i][0] += qqrd2e*q[i]*ek[i][0];
f[i][1] += qqrd2e*q[i]*ek[i][1];
f[i][2] += qqrd2e*q[i]*ek[i][2];

(f was defined earlier in the program as atom->f. qqrd2e is a unit conversion factor, also incorporating the dielectric constant, which is the Coulomb energy between two charges q at distance r: 332.06371 in "real" units with dielectric constant 1. [See force.h and force.cpp and update.cpp.]) Modify it to:

f[i][0] += qqrd2e*q[i]*ek[i][0];
f[i][1] += qqrd2e*q[i]*ek[i][1];
f[i][2] += qqrd2e*q[i]*ek[i][2];
// my code: fcoul stores the Coulomb contribution to force
atom->fcoul[i][0] += qqrd2e*q[i]*ek[i][0];
atom->fcoul[i][1] += qqrd2e*q[i]*ek[i][1];
atom->fcoul[i][2] += qqrd2e*q[i]*ek[i][2];
// end my code

Next, if you want to use slab correction (kspace_modify slab command in LAMMPS), you do the same thing in void Ewald::slabcorr(int eflag):

for (int i = 0; i < nlocal; i++) f[i][2] += qqrd2e*q[i]*ffact;
// my code
for (int i = 0; i < nlocal; i++) atom->fcoul[i][2] += qqrd2e*q[i]*ffact;
// end my code

dump_custom.h

I just wanted to output the electric field, so the dump custom command is a good way to do it. We want to add efieldx, efieldy, efieldz to the list of atom attributes that you can dump. In this file, you need to add to the big list of pack (for example, after void pack_fz(int);):

//my code
void pack_efieldx(int);
void pack_efieldy(int);
void pack_efieldz(int);
//end my code

dump_custom.cpp

In void DumpCustom::parse_fields(int narg, char **arg), there is a big list reading the input and calling the appropriate pack function, and you just add the new dump commands into it:

  ...
} else if (strcmp(arg[iarg],"fz") == 0) {
  pack_choice[i] = &DumpCustom::pack_fz;
  vtype[i] = DOUBLE;

  // my code
} else if (strcmp(arg[iarg],"efieldx") == 0) {
  pack_choice[i] = &DumpCustom::pack_efieldx;
  vtype[i] = DOUBLE;
} else if (strcmp(arg[iarg],"efieldy") == 0) {
  pack_choice[i] = &DumpCustom::pack_efieldy;
  vtype[i] = DOUBLE;
} else if (strcmp(arg[iarg],"efieldz") == 0) {
  pack_choice[i] = &DumpCustom::pack_efieldz;
  vtype[i] = DOUBLE;
  // end my code

Next, you need to add your new pack commands. Here is one of the ones in the original file...

void DumpCustom::pack_z(int n)
{
  double **x = atom->x;
  int nlocal = atom->nlocal;

  for (int i = 0; i < nlocal; i++)
    if (choose[i]) {
      buf[n] = x[i][2];
      n += size_one;
    }
}

We basically copy this format. So the new ones to add are:

// my code

void DumpCustom::pack_efieldx(int n)
{
  double **fcoul = atom->fcoul;
  int nlocal = atom->nlocal;

  for (int i = 0; i < nlocal; i++)
    if (choose[i]) {
      buf[n] = fcoul[i][0] / atom->q[i] / force->qe2f;
      n += size_one;
    }
}

void DumpCustom::pack_efieldy(int n)
{
  double **fcoul = atom->fcoul;
  int nlocal = atom->nlocal;

  for (int i = 0; i < nlocal; i++)
    if (choose[i]) {
      buf[n] = fcoul[i][1] / atom->q[i] / force->qe2f;
      n += size_one;
    }
}

void DumpCustom::pack_efieldz(int n)
{
  double **fcoul = atom->fcoul;
  int nlocal = atom->nlocal;

  for (int i = 0; i < nlocal; i++)
    if (choose[i]) {
      buf[n] = fcoul[i][2] / atom->q[i] / force->qe2f;
      n += size_one;
    }
}

// end my code

Note that force->qe2f is a unit-conversion factor: 23.060549 in "real" units (see update.cpp). It converts from a charge ("q") times an electric field ("e") to ("2") a force ("f"). Therefore this code will output an electric field in whatever unit system you're using -- e.g. volts per angstrom in the case of "real" units. (See full list of unit systems here.)

When compiling:

I found that my modified code had a heisenbug: When I compiled it normally, the program crashed, but when I compiled it in debugging mode, it worked fine. I don't know what the bug is. Instead, I just used the debugging-mode program for my calculations.

In your LAMMPS script file:

The code that goes into the LAMMPS script is something like:

dump StevesFavoriteDump all custom 10 stevedumpfile.xyze id type x y z efieldx efieldy efieldz
#10 means every 10 timesteps
#Note that atoms may be written in a RANDOM ORDER. They need to be sorted in postprocessing.

Again, remember to turn off newton:

newton off

and remember that if you want to include the electric field due to directly-bonded atoms, you need

special_bonds coul 1.0 1.0 1.0

or you need to edit pair_lj_cut_coul_long.cpp differently. Note that the special_bonds command can be overwritten by other pair commands if you're not careful, always check your work independently.