[gmx-developers] ensemble code and restraints

[Mark Abraham]

Robert Fenwick wrote:
> Hi Berk,
> 
> It has been shown a few times now that averaging over a single
> structure is not sensible as it is possible to get strained
> structures, the other alternative is to do as you suggest to average
> over the whole ensemble. This also creates problems and if you average
> over three or more structures then it is possible that one completely
> unfolds while the others compensate. The same work has shown that if
> you average over pairs of structures you get around this problem. In
> theory I can do all of this in CHARMM. What I would like to know is
> how to implement a similar scheme in GROMACS, so that others can use
> it. I am currently lost in the source code and was looking for some
> pointers as to how I could get started.

There are some implementation similarities between what you propose and 
REMD. (I think) that at each(?) MD step you will need to swap certain 
atomic positions between pairs of simulations so that each simulation 
master process can then apply the pair-averaged restraint force without 
further inter-simulation communication. In REMD, at every so many MD 
steps, you swap the potential energy of certain pairs of simulations, 
and then maybe the whole simulation state. You can get a fair idea of 
how you should approach this project by understanding the REMD algorithm.

Probably your most efficient implementation is to use some nonblocking 
MPI sends and receives for the relevant subset of positions before you 
start doing the normal nonbonded kernel loops, so that the communication 
completes by the time you're done with the kernels and are ready to add 
in the restraint forces. I've no idea how well that approach suits 
modifying the existing implementation of intra-simulation restraints.

Mark