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

[David van der Spoel]

Mark Abraham wrote:
> 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!

C is good enough for this analysis type computations, although doing it 
during mdrun will make it probably ten times slower than normal mdrun. A 
further question is whether to average the field or to dump the 
instantaneous field each 20 fs as suggested by John. The answer to this 
depends on the type of experiment one compares to. It would also be nice 
to compute the electrostatic potential while we're at it. Finally if I 
understand it correctly for this calculation one would need to store a 
3x3 tensor per atom, which means it will give some overhead for parallel 
computation as well for passing around this data.

Finally 2, one will need to do this in two parts: based on the direct 
space and based on the PME grid.

> 
> 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
> _______________________________________________
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the www 
> interface or send it to gmx-developers-request at gromacs.org.


-- 
David van der Spoel, Ph.D.
Molec. Biophys. group, Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205. Fax: +4618511755.
spoel at xray.bmc.uu.se	spoel at gromacs.org   http://folding.bmc.uu.se