[gmx-developers] Writing the electric field (and its gradient) on particular atoms during an MD simulation
Mark.Abraham at anu.edu.au
Tue Apr 15 09:18:44 CEST 2008
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
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.
More information about the gromacs.org_gmx-developers