[gmx-developers] Writing the electric field (and its gradient) on particular atoms during an MD simulation

[John Chodera]

Gerrit, Mark, and David:

Thank you all for your insightful comments!

To clarify the matter, I'll respond to a few questions posed by your emails:

Gerrit correctly observes that the electric field can, in fact, be
obtained from the forces if we do not compute all the other
non-electrostatic contributions.  However, since we need the spatial
gradient of the field for ~ 40 atoms, the finite-difference approach
to estimating the spatial gradient would require on the order of this
many electrostatic evaluations for every snapshot we reprocess (20 fs
between snapshots in a 100 ps simulation), which would be too
computationally costly.

Mark suggests looking at the .c files in
src/gmxlib/nonbonded/nb_kernel/, but I admit that all of these look
like gibberish to me.  Is there some sort of guide to the construction
of these nonbonded kernels, and what all these variable names mean?

David notes that the overhead of writing out 3x3 (symmetric) electric
field gradient tensors could be large, but I should note that we only
need the field and field gradients on a small subset of the solute
atoms -- just a few atoms along the backbone of a solvated peptide,
for example.  Even if we have to call some extra routines every 20 fs
to recompute these fields and gradients, it wouldn't be too expensive.
 Even rerunning the PME code every 20 fs wouldn't be all that bad, if
we needed to do so -- though it would obviously be ideal to avoid
this.

As David points out, yes, the electric field should in fact include
the contributions from both the direct-space and reciprocal-space
components when PME is in use.

Thomas notes that we actually need to write out the field and
gradients at these times because he computes two- and four-time
correlation functions from the data, so that simply accumulating
averages would not be sufficient.

In an ideal world, the code would be straightforward enough that I
could just add some extra lines to (1) compute the analytical
contribution to the field gradient from the direct space parts [easy],
and (2) compute the field gradient tensor from the PME spline
interpolation, which I believe should be linear if a cubic spline is
used [also easy!].  But I don't even really know how to begin in the
gromacs codebase.  If any further concrete pointers could be given
(e.g. "modify this function in this file, using the information
contained in these variables"), it would very much be appreciated.

I should note that this would be useful for another project in the lab
-- the calculation of Stark shifts -- as well.

Finally, Mark's point about whether these fields are really meaningful
for fixed-charge forcefields is certainly well-taken.  I am familiar
with the shortcomings of these models and the charge parameterization
philosophies.  For forcefields like AMBER, the field defined in this
way is both (1) consistent with the definition of the potential energy
function, and (2) consistent with the models used to compute
spectroscopic data.  See Thomas's excellent paper on this if you are
curious about how their method works:

Jansen TLC, Knoester J. Nonadiabatic effects in the two-dimensional
infrared spectra of peptides: Application to alanine dipeptide. JPC B
110:22910, 2006.

http://msc.eldoc.ub.rug.nl/FILES/root/KnoesterJ/2006/JPhysChBlaCourJansen/2006JPhysChemBlaCourJansen.pdf

Best,

- John

--
Dr. John D. Chodera <jchodera at gmail.com>      | Mobile    : 415.867.7384
Postdoctoral researcher, Pande lab            | Lab phone : 650.723.1097
Department of Chemistry, Stanford University  | Lab fax   : 650.724.4021
http://www.dillgroup.ucsf.edu/~jchodera