Add restraint on Phosphorous atoms of lipid headgroups that were farther than 0.5 nm from protein

User discussions
1

GROMACS version: 2018.8
I am doing membrane simulation, where I am doing umbrella sampling by pulling my substrate through the transporter and generating bins and doing umbrella sampling.

I am trying to add positional restraint on Phosphate atoms of lipid during the pulling part, for which I am using “gmx select” command to select phosphates of lipid whichare farther than 0.5nm from protein, and using the output index file to feed in “gmx genrestr” to generate .ipt file for the positional restraints.

I am getting this error:
ERROR 1 [file lipid_05.itp, line 5]:
Atom index (3687) in position_restraints out of bounds (1-134).
This probably means that you have inserted topology section
“position_restraints”
in a part belonging to a different molecule than you intended to.
In that case move the “position_restraints” section to the right molecule.

My topol.top file:
;;
;; Generated by CHARMM-GUI FF-Converter
;;
;; Correspondance:
;; jul316@lehigh.edu or wonpil@lehigh.edu
;;
;; The main GROMACS topology file
;;

; Include forcefield parameters
#include “toppar/forcefield.itp”
#include “toppar/SOD.itp”
#include “toppar/CLA.itp”
#include “toppar/TIP3.itp”
#include “toppar/CARA.itp”
#include “toppar/PROA.itp”
#include “toppar/POPC.itp”

[ system ]
; Name
Title

[ molecules ]
; Compound #mols
PROA 1
CARA 1
POPC 161
SOD 32
CLA 43
TIP3 12246

My POPC.itp file: (last few lines)
.
.
.
125 128 131 132 9
125 128 131 133 9
125 128 131 134 9
129 128 131 132 9
129 128 131 133 9
129 128 131 134 9
130 128 131 132 9
130 128 131 133 9
130 128 131 134 9

[ dihedrals ]
; ai aj ak al funct q0 cq
31 30 33 32 2
40 39 42 41 2

#ifdef POSRES
[ position_restraints ]
20 1 0.0 0.0 POSRES_FC_LIPID
#endif

#ifdef DIHRES
[ dihedral_restraints ]
25 36 28 30 1 -120.0 2.5 DIHRES_FC
60 63 65 67 1 0.0 0.0 DIHRES_FC
#endif

; Include restraints on lipids which are farther than 0.5nm from protein
#ifdef PORSE_PHOSPHATE
#include “lipid_05.itp”
#endif

here I am defining new restrain file i have generated using “gmx genrestr”

my lipid_05.itp file: (first few lines)
; position restraints for phosphate_notwithin0.5_protein_f0_t0.000 of Title

[ position_restraints ]
; i funct fcx fcy fcz
3687 1 0 0 2000
4759 1 0 0 2000
4893 1 0 0 2000
5027 1 0 0 2000
5563 1 0 0 2000
5697 1 0 0 2000
5831 1 0 0 2000
5965 1 0 0 2000
6099 1 0 0 2000
6233 1 0 0 2000
6501 1 0 0 2000
6769 1 0 0 2000
6903 1 0 0 2000
7037 1 0 0 2000
7171 1 0 0 2000
7305 1 0 0 2000
7439 1 0 0 2000
7573 1 0 0 2000
7707 1 0 0 2000
7841 1 0 0 2000
7975 1 0 0 2000
.
.
.

pulling.mdp file: (first few lines, where i am defining restrains)
title = Umbrella pulling simulation
define = -DPORSE_PULL -DPORSE_PHOSPHATE -DPORSE_CA
; Run parameters
integrator = md
dt = 0.002
tinit = 0
nsteps = 4000000 ; 8000 ps
nstcomm = 10
.
.
.

When i am doing the same thing for Protein C-alpha restrain it is working fine, so why not for lipids?

how I can add lipid restrains for selective lipids only?
Please suggest me a correct way of doing this

2

First of all, the main problem is that position restraints are defined PER molecule. Take a look at some position restraints you already have, like the one you are also reporting here

#ifdef POSRES
[ position_restraints ]
20 1 0.0 0.0 POSRES_FC_LIPID
#endif

You will notice that they are defined just once inside the lipid. This is referring to atom 20 of your lipid, for example.

What you are doing with gmx select is asking for a list of the P atoms of the lipids, and then you generate the restraints on that list. When then you run the simulation, the grompp passage sees POPC (that has something like 134 atoms) and the associated restraint file, which however contains indexes for ALL the POPC P in the system, which have indexed > 134, and gets confused because it expects position restraints defined on the atoms of ONE lipid. In fact, the error is

Atom index (3687) in position_restraints out of bounds (1-134).

which is telling you ‘I see that you want to restrain atom 3687, but your topology has 134 atoms. You are doing something wrong’. So, restraints need to be defined WITHIN the topology and referring ONLY to that atom numbering inside of the topology. If you restrain the position of P inside the topology then this will be applied to all the POPC of the system. You can read further here.

Secondly, you can’t use position restraints the way you want to use them in this case. The position restraint works by keeping the position of something (e.g. P atoms of POPC) restrained with respect to a position of reference that you provide with the -r flag during grompp preprocessing. As such, you can’t dynamically pick those POPC that are within 0.5nm of your protein. What you can do is, for example, see in the starting frame which POPC molecules respect this condition, and use the pull-code to define pull-groups only for those molecules, and keep them – for example – within some cut-off distance from the centre of mass of the protein. Again, this won’t select which lipids are within the cut-off dynamically, won´t prevent other lipids to come within the cut-off, and will be ill-defined if the protein is not reasonably spherical and with its centre of mass corresponding to its geometrical centre. To be more elegant you may want to rely on external tools such as PLUMED or the colvars module, which is now native in gromacs 2024 distributions. I do not have enough experience in these to tell you if you can use these to do what you want.

3

Thank you @obZehn for your prompt answer.
I know there is a mismatch in the lipid index I am referring to the new lipid_05.ipt vs the topology file (where only a representative lipid is taken into account).
But is there any way I can define those lipids using the distance restraints or positional restraints, that we want to apply only on those lipids that are >0.5nm from protein (I want to select such lipids just after the equilibration steps, using the last .gro file).
If I restrain P present in the topology of lipids, it will restrain all Phosphates, not the specific ones.
As we do in Flat bottom potential, we define a radius of some x nm to restrain our substrate, can we define a hollow cylinder that will take care of our problem? I don’t know how to do so.

4

I am not completely sure I am getting the right idea about what you are trying to achieve! Do you want to restrain the position of the POPC P atoms where they are but just of those POPC that are > 0.5 nm away from the protein? In that case, I do not see a simple way via straight plain GROMACS to do that. There are dirty ways to do this. An idea could be to double the topologies of your system, that is, generate POPC of type A and of type B, that are identical except for the name and the restraint on P. But that would mean that you have to modify and reorder both your box (.gro) and your topology (.top), which is very error prone, plus this is not dynamic, and I would argue that it is a very borderline approach. It’s bit hard to justify something like this and keep it under control.

Honestly, I would take a look at the external tools I linked before to see if it is doable. I am pretty sure that there are selection criteria within PLUMED and colvars that let you dynamically keep under control everything and keep this lipid list updated. Then you would have to find a way to apply the restraining potential only to selected atoms, which again should be feasible. I hope I am wrong and someone with more experience can point a way to do this in native GROMACS, but from my (limited) experience you really can’t as both restraints and pull-groups have no dynamical selections in their definitions.

5

Thank you @obZehn , I will look into the external tools that you have mentioned for this purpose.