[gmx-developers] Update partial charges after every iteration, add a term to the potential energy

[Anton Feenstra]

Dear Miguel,

On 13/12/11 00:39, Miguel Garcia wrote:
> Hi
>
> I am new at modifying modifying the gromacs source code, so I need a place
> to start. What I intend to do, is update the partial charges of a specific
> small molecule after every iteraction, depending on the electric field
> felt by the molecule. My question is basically which c files I need to
> modify? More precisely:

First off, excuse my ignorance, but should you not worry about energy 
conservation? And, as an aside, isn't this some sort of polarizability 
feature? Then you should be aware that several different options to 
implement this are conceivable, and that one (or a few) are under 
construction in various stages (search the mailinglists and/or web).

> a) Which c file gives me the electric field at a particular position?

I'm not sure if that is explicitly available anywhere

> b) Which c file reads the partial charges from the topology? (So I can
> change it to read the charges from an internally updated table)

Only read once, then stored (iirc, in the mdatoms datastructure).
We may have an option for you that does not involve diving hip-deep in 
the Gromacs source.

A post-doc in our group, René Pool, has used the GromPy python interface 
to implement a library interface to the mdrun integrator (mdrunner). 
This library interface exposes the md datastructures, and makes it 
trivial to do whatever coordinate transformations that you might 
require, and then continue the simulation. I cc him here now.

There is an open source repo for GromPy
at gitub: https://github.com/GromPy
Also, a JCC paper is almost out (accepted with minor revisions).

The usecase we developed GromPy for is to 'intervene' in the MD at 
regular intervals - in effect you start a short simulation, get back 
final coordinates, forces, etc. (the whole md datastructure is 
accessible, so also the charges), then do whatever you want (we have 
this embedded in an MC scheme) and start another (short) simulation. But 
because you deal with ready-made topologies, have everything in memory, 
and avoid system calls or shell scripting (e.g. for grompp), there is 
virtually no overhead in the runtime (less than 5%, I think generally 
even below 2%). This is not quite what you do, but may get you close.

What we have been thinking about, but not yet implemented, is to have a 
callback function from Gromacs, at one or several points in the 
integration timestep. That would allow you to do arbitrary manipulations 
at any time during the simulation. GromPy would be a good framework for 
that; the thing you need to implement in the gromacs c code are the 
callbacks.

> An alternative is for me to add a new term to the non bonded potential
> energy. I cannot use the tabulated table feature, since I want this new
> term to be proportional to the Field.
>
> c) So in which c file is the potential function calculated?

There is a howto on the Gromacs FAQ/wiki somewhere that describes all 
the steps needed to add an interaction function.


-- 
Groetjes,

Anton
  _____________ _______________________________________________________
|             |                                                       |
|  _   _  ___,| K. Anton Feenstra                                     |
| / \ / \'| | | IBIVU/Bioinformatics - Free University  Amsterdam     |
|(   |   )| | | De Boelelaan 1083A - 1081 HV Amsterdam - Netherlands  |
| \_/ \_/ | | | Tel +31 20 59 87783 - Fax +31 20 59 87653 - Room P136 |
|             | Feenstra at few.vu.nl - www.few.vu.nl/~feenstra/         |
|             | "Does All This Money Really Have To Go To Charity ?"  |
|             | (Rick)                                                |
|_____________|_______________________________________________________|