I am trying to set up the coordinate restraints for my HMMM lipid bilayer system following the CHARMM-GUI HMMM Builder for Membrane Systems Paper (found in supporting information). Here they restrain DCLE molecules to be between -1.053 and +1.053 nm on the Z axis. The terminal carbon on the lipid tails were restrained to be between -0.853 and +0.853 nm on the Z axis. Finally, the water molecules were restrained to be outside the layer of DCLE.
I think this can all be implemented with pull code in the .mdp files however I am unfamiliar with setting this up so any help is very appreciated.
Thanks
Flat bottom position restraints on the lipid tails is the closest built in restraint that I’ve seen. However the goal of the restraint is to force the tails to interact with the DCLE layer. A flat bottom position restraint does not accomplish this as the tails would be able to move freely +/- the set distance starting from the initial position.
Ideally I would like to tell GROMACS “when this atom is above or below a specified Z coordinate apply a force”. The problem with flat bottom is that the reference position that is used is centered on each individual atom and I am not sure if I am able to change that.
That’s what the flag -r does, that is, when you pre-process the binary with grompp you have to pass a reference structure with respect to which the restraints are applied. As such, if you project all the atoms of interest to the same z value (which is a the center of the region where you want the lipids to be confined) and pass that structure to -r, then the flat bottom will confine everything to be inside that slab.
It is a bit confusing at the beginning, you can find some discussions here and there about their implementation though, like here, here, and here.
Maybe you can achieve something similar with the pull module, but I would say that the flat bottomed potential seems much more straightforward, while in the pull module you have to i) define each and every lipid as separate pull groups or ii) rely on their center of mass, which is arguably very degenerate if you want to restraint every lipid within a slab (the COM can still satisfy the distance restraints but individual lipids may be everywhere, as only their COM contributes to the collective variable).
I second @obZehn 's answer.
On a related note, I found the name “restraint.gro” confusing. I found it actually means the reference position for any restraints you applied in your simulation.
For example, applying a harmonic position restraint requires two information:
force constants (usually intuitive enough; there are obvious places to specify these in itp/mdp/top);
reference position (I find this very cryptic as it is currently implemented; these positions are specified in the gro file you give to -r flag in grompp, with atom-to-atom correspondence to what you give to the -c flag.
But then, usually you don’t restrain every single atom in your system, so only the xyz values for the atoms with restraints matter in restraint.gro; the rest are just dummies. The corollary of this is that restraint.gro itself does not have to make physical sense. You can edit it to meet the needs of your restraints. The atoms can clash as long as clashes are against dummies).
I’ve made mistakes due to not knowing how to specify reference positions. I did some tests and this is my understanding.
Yes, you are correct. The -r flag takes a file that contains the same number of atoms which will be read one by one to set the reference position for the applied restraint, e.g., the reference position for a spring potential. In most cases the interest is to keep something ‘there’, where ‘there’ means locally where it is, so the -r file is - geometrically speaking - the same as the one in the -c flag. But for fancier restraints, such a flat bottom potential to keep all atoms restrained in a certain region of space, the file need to be modified so that the reference positions for the restraints for the file in the -c coordinates has meaning, e.g., projecting all atoms into the center of the potential region.