I’m trying to apply a flat-bottom position restraints to a specific subset of atoms in my system. I was able to apply such a restraint to all atoms of a specific type with the following directive in an .itp file:
[ position_restraints ]
; Flat-bottomed potential to avoid large-scale membrane deformations
; id ft geometry r k ; attractive Z-layer
57 2 5 2.25 1000 ; attractive Z-layer
And #inlcudeing that file in the parent [ moleculetype ] directive where 57 is the id of the atom type I want to restrain:
[ moleculetype ]
; name nrexcl
OOPC 3
[ atoms ]
...
57 cA 1 PC C2 57 0.281242 12.0100 ; qtot 0.282
...
# include "my_restr.itp"
However, I can’t figure out how to do apply this restraint to just a specific list of atoms. I tried generating a restraint reference file, restr_C2.gro, which only included a subset of the atoms, but that, of course, fails with
WARNING 1 [file topol.top, line 30]:
The number of atoms in restr_C2.gro (90) does not match the number of
atoms in the topology (247738). Will assume that the first 90 atoms in
the topology and restr_C2.gro match.
Running gmx genrestr with an index file containing the to-restrain atom IDs produces an .itp file:
But that’s a harmonic restraint, not a flat-bottom restraint. Can I just take the file generated by genrestr, and replace all the lines after the i column with 2 5 2.25 1000, as in my previous file? This is also making me wonder if originally I only restrained atom 57 instead of atom type 57, though the results of the simulation seemed different once I added that restraint compared with the initial simulation.
The way you define and include the position restraints is correct, if you only want to restrain that atom of that molecule (you are trying to open a toroidal pore, I guess?).
Then the .gro file has to have the same number of atoms as the box or it will apply the restraint to the first N atoms, which in your case is clearly wrong as you want to restrain one atom per lipid.
The solution is to generate a .gro (there are also other formats supported) file that contains the atoms in the reference position with respect to which you want them restrained. Generally, when one equilibrates a system, the restraints are the starting conformation so that you harmonically restrain the atoms wrt their starting position. In your case, instead, I think you want to force a lipid atom to stay within 2.25 nm from the center of the lipid bilayer. What I used to do in this case was to produce a frame from the end of the energy minimization and then, with an awk script, rewrite all the z values of the lipids to a reference value which is the COM of the membrane. With this as reference frame, the restraint you define should work the way you want, i.e., keep the lipid atom within 2.25 nm of the COM. In principle, you can send all the atoms of the box to the z of the COM of the membrane, to make your life easier. However, beware that if you are also using other restraints then this will break the simulation.
What you suggest is indeed what I’ve done to restrain all lipid head groups. However what I realized is that this actually has some issues with toroidal pore formation, as it prevents the head groups from entering the pore, making the energy of the fully formed pore artificially high. So what I wanted to try was to select a random subset of lipids (arbitrarily chosen at 10% for now) and just restrain those. Hopefully this is enough to prevent membrane oscillation during pore formation, while still allowing the rest of the lipids to freely rearrange around the pore.
This happens if you put a lower wall, so you do not let your lipids reach the center of the membrane. I never had this problem because I always did the opposite, so put an upper-wall for the membrane to not deform too much and this leaves the lipids free to diffuse towards the center of the membrane to open up the pore. Have you tried in this way?
The random selection wont work simply like this, though. The moment you put the restraints in the topology then those are going to be applied to all the lipid of that type. If you want to wall ~10% of them you will have or to write an individual entry pull-group for each of these lipids in the mdp file or split the topology of these lipids in two groups and generate one with a name with restrains and one with another name without restraints. This will require a lot of work on the box set-up!