Good luck!
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.
Add a line into the initialization (i.e. Atom::Atom(LAMMPS *lmp) : Pointers(lmp)):
x = v = f = NULL; fcoul = NULL; // my codeAdd a line into the destruction (i.e.Atom::~Atom):
memory->destroy_2d_double_array(f); memory->destroy_2d_double_array(fcoul); // my code
Add into "private" area (after "double **x,**v,**f;"):
double **fcoul; // my code
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);
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).
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:
// 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 codeAlso, in your LAMMPS input script, you must turn off "newton":
newton offIf 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.
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 codeNext, 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
//my code void pack_efieldx(int); void pack_efieldy(int); void pack_efieldz(int); //end my code
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 codeNext, 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 codeNote 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.)
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 offand 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.0or 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.