[gmx-developers] grand-canonical MC using gromacs

[René Pool]

Mark Abraham wrote:
> On 8/12/2010 1:18 AM, René Pool wrote:
>> Hi,
>>
>> Thanks Mark, found it!
>> For me the problem is
>> static bool bFirst=TRUE;
>> which is set to FALSE after initialization of the neighbor lists.
>> in function ns() (in force.c).
>>
>> Resetting this variable to TRUE seems to yield the correct average 
>> density in an MC simulation. The downside is that by doing so, a huge 
>> memory leak is created.
> 
> Yep. GROMACS is very bad at managing its "once per simulation" memory, 
> on the theory that it's only doing one simulation. I once had a version 
> of GROMACS based on 4.0.x that had a pile of memory-cleanup routines of 
> the kind you might need. I could make them available if anyone is 
> interested.
> 
> Hopefully the move to C++ will help a bit with some formal destructors.
> 
> Mark

I already made some cleanup routines of my own to solve some memory 
leaks. Maybe I missed some pointer cleanups though, so I'm certainly 
interested in your cleanup routines!

I solved the above neighbor list memory leak by only initializing the 
nlist after a tpr read. This list is then correctly managed by 
mdrunner_integrate().

I'm happy to announce that I now have a working GCMC module for gromacs: 
the Lennard-Jones muVT equation of state at reduced temperature T*=1.5, 
completely agrees with one determined via NVT MD. This temperature is 
well within the fluid region of the LJ phase diagram, so there is no 
phase coexistence in the simulation box.
The module is written in python and uses direct library calls to the 
gromacs libs.

I'm planning to further increase efficiency by exploiting the test 
particle insertion option of mdrun for MC insertion moves and computing 
the energy of a randomly selected molecule for possible removal during a 
consecutive MC removal move. The latter trick will save single point 
energy calculations for MC removal moves that I perform at the moment.

After this I'll share the GCMC module with the gmx community for those 
who are interested :-)

I expect that it will not take too much trouble making the GCMC module 
compatible with the latest version of gromacs. I will explore doing this.

Cheers and thanks for the help,
Rene

> 
>> At least I now know where to find the solution!
>>
>> Cheers,
>> René
>>
>> Mark Abraham wrote:
>>>
>>> On 12/07/10, René Pool <r.pool at vu.nl> wrote:
>>>> 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?
>>>>
>>>
>>> When are you constructing and managing the neighbour lists? These are 
>>> not trivial to manage well.
>>>
>>> Mark
>>>
>>>
>>>>
>>>> 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
>>>> =====================================================
>>>> -- 
>>>> gmx-developers mailing list
>>>> gmx-developers at gromacs.org
>>>> http://lists.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.
>>>>
>>>>
>>>