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

[Mark Abraham]

Gerrit Groenhof wrote:
> Hi John,
> 
> Instead of modifying the code (which is more elegant of course), would 
> it not be possible to do a rerun, with all bonded and Lennard-Jones 
> interactions set to zero in the topology? I think that in that case the 
> forces on the atoms, as printed in the trr file are equal to the 
> electric field
> (f(i)=q(i)*E(i), right?). Or do I overlook something here?
> 
> Depending on the accuracy you want, The gradient (Nable E) might be 
> obtained by finite differences, once you have the electric field 
> associated with each atom position, but I do not know how easy or cheap 
> this is.

Easy to write is a new utility to iterate over all frames in a 
trajectory, iterating over all atom coordinates to perturb the positions 
to create the finite differences and write a new "trajectory". Then 
mdrun -rerun on this new trajectory, updating lists only every 
old-trajectory-frame steps.

There could also be algorithms for finite differences over 3N 
coordinates that don't require order-3N evaluations, I suppose.

Quicker to run would be a customized non-bonded kernel for the field 
gradients. Since the existing routines are already computing/looking up 
1/r and the charge of the interaction partner, there's no overhead for 
the field gradients. There is a cost to storing them so that they can be 
dumped to output *outside* the nonbonded loops (else the loops will 
probably run terribly). Obviously, one should cut one's teeth on the C 
files in src/gmxlib/nonbonded/nb_kernel before touching the assembler!

One should determine whether existing (non-polarizable) force fields 
should be expected to produce a reasonable atom-wise decomposition of 
electric field (or its gradient) - if the force fields can reproduce a 
suitable observable then I suppose the decomposition *might* work.

Mark