[gmx-developers] grand-canonical MC using gromacs

[René Pool]

Dear all,

We're busy implementing grand-canonical Monte Carlo using gromacs. 
Particle translations accomplished by NVT MD, rendering this a hybrid MC 
MD method. Our approach is doing this by a python shell around gromacs 
in which we can insert/remove particles and perform MD NVT moves. For 
optimal efficiency, we use the ctypes module. This allows us to keep the 
gromacs data structures in memory during the MC run.

To this end we had to make minor modifications/additions to the gromacs 
src (ver 4.0.5 because we started this project with this version). Most 
importantly, we have split up mdrunner() into (1) mdrunner_initialize(), 
(2) mdrunner_integrate() and (3) mdrunner_finalize(). The first is for 
initializing the mecessary gromacs data structures in python and storing 
them in memory, The second performs NVT MD Monte Carlo moves and single 
point energy evaluations of a given state. The third is called after the 
MC loop for neatly finalizing the simulation.

For each number of particles N, we have the "tpr" in memory. In the case 
of a MC insertion or removal move, we copy the current state->x and 
state->v arrays into the trial state->x and state->v arrays. For 
insertion we add one molecule to these arrays in the case of insertion 
(i.e. replace state->x[last-1] and state->v[last-1] by a trial position 
and a velocity vector taken from the Maxwell-Boltzmann distribution). 
For removal we select a molecule at random and exclude this entry from 
the lists before copying it to the trial state. These "trial states" 
serve as input for mdrunner_integrate() that computes the single point 
energies (an md of nsteps=0).

So far so good, everything works, BUT...

The energies computed by mdrunner_integrate() do not correspond to the 
state and therefore screws up the MC sampling.

My question therefore is:
Are we overlooking some variable in the gromacs data structures that 
also needs to be set during insertion/removal, or is there some static 
(global) variable that stands in the way?

Thanks in advance for your help!

By the way, we're testing this on a simple LJ fluid system on a single 
CPU. So for now, no complicated stuff such as parallellization or charges.

Cheers,
Rene

-- 
=====================================================
René Pool
IBIVU/Bioinformatics
Vrije Universiteit Amsterdam
De Boelelaan 1081a
1081HV AMSTERDAM, the Netherlands
Room P120
E: r.pool at few.vu.nl
T: +31 20 59 83714
F: +31 20 59 87653
=====================================================