Gromacs/MFFA is a customised Gromacs version that supports user-defined, implicit, mean-field force approximation (MFFA) boundary potentials. It offers a computationally efficient method to simulate e.g. large vesicular systems in spherical liquid droplets.
Gromacs/MFFA was originally written by Jelger Risselada for Gromacs 3.1.4 and has been subsequently ported e.g. to Gromacs 4 by Martti Louhivuori. The current version is stable (at least in the machines we have used) and suitable for production runs (AFAIK), but since Gromacs/MFFA is still under active development one should always carefully validate all results obtained by it. For instructions on how to use it, see below.
The current maintainer of Gromacs/MFFA is Martti Louhivuori.
If you have any suggestions / bug reports, please send a sufficiently
detailed description to
m.j.louhivuori@rug.nl.
Download
A mean-field force approximation (MFFA) potential is an implicit interaction potential that mimics a bulk solvent phase. It is intended to match the effective, net interaction a particle would feel from the solvent particles replaced by the potential.
Brooks and Karplus showed [1] that such a mean-field potential can be readily calculated from the radial distribution function (RDF) of a pure liquid system. The exact form of the MFFA potential depends on the geometry of the boundary, but e.g. for a Lennard-Jones liquid it has the same kind of attractive and repulsive regimes as the LJ-potential.
The main advantage of replacing part of the explicit solvent with an implicit MFFA potential is, of course, computational efficiency. It is trivial to show that e.g. for a spherical system the number of particles increases cubicly with an increase in the radius of the system. In practice this means that one spends a relatively large proportion of computational resources on the (less important) bulk solvent surrounding the system of interest. An enveloping MFFA potential allows one to use a thinner bulk layer of explicit solvent, and thus improves computational efficiency dramatically. MFFA boundary potentials have been shown [2] to also have sometimes other secondary advantages e.g. through a molding effect to speed up spontaneous formation of lipid vesicles. It is also possible to create a pressure gradient across a vesicle e.g. by coupling two independent MFFA potentials (one inside and one outside) to different ambient pressures.
Set-up and run an equilibrium simulation of a pure liquid system under the same environmental conditions that you wish to use for the MFFA simulation. Of course, you need to use a large enough simulation box and run for long enough to achieve proper statistics for a good quality RDF. Once done, use g_rdf to obtain the actual solvent-solvent RDF:
Set-up the actual simulation system that you want to run with Gromacs/MFFA. E.g. to simulate something in a spherical droplet of water instead of the normal periodic box, you need to surround the (periodic) system with a spherical water shell of required dimensions. This is readily done using editconf, genbox, and a Python-script called clean-sphere.py (available from here).
So, first edit the box to be large enough and center the system of interest:
Note: the box dimensions should be at least couple of nm larger than the actual size of droplet. This is to ensure that there is no "tunneling" through the MFFA potential in combination with the periodic boundaries that will still be turned on for computational efficiency. Since this won't have a negative effect on the performance of the MD engine, you may just as well use a relatively large vacuum to avoid any "tunneling" events.
Then solvate the rest of the box (van der Waals radius is increased for Martini CG water):
And then remove all water that is beyond the MFFA boundary to make a nice, round water droplet:
If using a different solvent than CG water, you need to replace "W" with the correct residue name.
Once you have your system defined, you need to prepare a matching MFFA boundary that will envelope the droplet. This is simple to do, but may require a little bit of trial and error to get the dimensions correct (we'll come back to this in step 5).
The easiest way to prepare a MFFA boundary potential for a Lennard-Jones type solvent (such as Martini CG water), is to use the tool g_mffa included in Gromacs/MFFA. For our 30nm diameter water droplet, this would be done with the command:
where RHO is the number density of the solvent (step 1). If using a different solvent than Martini CG water, it may be necessary to change some of the other parameters as well (such as the van der Waals radius of solvent particles, the LJ parameters, interaction cut-off etc.). The -quench option selects the type of the boundary potential:
Next, you'll need to set up system topology (.top), MD parameters (.mdp) etc. just as you would do for a normal simulation. The only change you need to do is to set the following options in the .mdp file:
Once everything is ready, you should run a short (few ns) simulation to equilibrate (at least roughly) the system with the MFFA boundary. If all goes well, the boundary doesn't need to contract/expand too much, but will only fluctuate around the original radius, and you are set to go for longer simulations.
If you observe a clear shrinking/expanding of the boundary, i.e. the position (=radius) of the potential (the third line in the .sdb file) changed by more than, let's say, a nanometer, then you should go back to step 3 and re-calculate a more proper MFFA potential; i.e. use the final radius of the short run as the -outer radius of the new MFFA potential. Repeat until you have a MFFA potential that corresponds (roughly) to the geometry of the equilibrated droplet.
The main difference from running normal Gromacs is that you need to provide at least one MFFA potential (.sdb) to mdrun using the -sb option.
If you have multiple potentials (e.g. a fully repulsive one for lipids to prevent them sticking to the boundary) you need to supply them like this: -sb "one.sdb#two.sdb" and also provide a boundary-index file (.ndb) created with ndx2ndb (included in Gromacs/MFFA). The order of the potentials is decided by the order in which you picked the groups in ndx2ndb.
Finally, you can run long production simulations to study the next BIG THING™...