outputting one component of force, e.g. electric force

Hi folks,

I recently had to figure out how to output the electric force on atoms
in LAMMPS. (It's easy to output the TOTAL force, but I wanted just the
coulombic contribution to the total force.) Here's what I came up


It requires modifying the source code in various places. I don't know
whether there are simpler methods, but this did work for me. I'm
posting it in case anyone else ever needs to do this!

(The instructions could presumably be modified for if you wanted to
output the (fill-in-the-blank) contribution to the total force.)

Why was I doing this? It turns out that the local electric field
experienced by a hydrogen atom in simulated SPC/E water is a
reasonably accurate proxy for the vibrational frequency of its OH
bond. I was using that fact to do an oversimplified simulation of the
infrared nonlinear optical properties of water. (
http://dx.doi.org/10.1016/j.cplett.2011.08.027 )


This is a reasonable way to do this, assuming you
don't want to just post-process the dump files (e.g.
subtract out the vDwl term from the total force).

It is, however, far too invasive to include as an option

If you were willing to run with rRESPA, with loop
settings all = 1 (which is algebraically the same as
Verlet), then the fix_respa.cpp code could have a method
added to allow one level of force to be output.

In your case, that wouldn't help separate the short
range vDwl from Coulombic, but maybe that could
be done with pair hybrid. Have to think about it a bit.


Another thought is to wait until we release
a "rerun" command, hopefully soon, that will
let you take a dump file of snapshots and rerun
them through LAMMPS, and compute new
quantities on each snapshot.

If you did that with only a Coulombic potential
defined, you could re-dump the total force which
would probably be what you want.