Inverted Flat-bottomed position restraint

User discussions
1

Hi everyone,

I need to use an inverted flat-bottomed position restraint on water molecules for a lipid bilayer system. According to gromacs manual, I added following section to the tip3p.itp file:

[ position_restraints ]
;i function g (type) r (nm) k(kJ/mol)
1 2 5 -2.0 1000

But it causes an unstable system.

Is it correct the way I use it for an inverted restraint?

Best,
Mehry

1 Like
2

Here you are applying the restraint on atom 1 of your system only. Regardless of the fact that restraining water molecules might have unexpected effects on your system, you need to select all the water molecules you are restraining in the position_restrains list.

Aside from that, what stage of your simulation you are applying this restraint to? Let’s say, in an NPT ensemble, you don’t want to restrain the water molecules, as they are basically the means to adjust the pressure in your system.

[ position_restraints ]
;i function g (type) r (nm) k(kJ/mol)
1 2 5 2.0 1000

Lastly, the fourth value here refers to the distance from the center, where the force acts on, in an inverted flat-bottomed restraint. I suspect that it might not need to be a negative value. See here. Check on that please, as I am not quite sure.

3

Hi Ali,

Thank you for your reply.

I use it for an NVT equilibration. Actually, I already used position restraints and I am almost familiar with using them. My question is just about an inverted one.

As I know, when you use position_restrains in your topology file, you don’t need to have all indices. It means that when I choose “1”, it would consider all oxygen atoms in all water molecules. It’s just much easier than having a big list of indices :)

Best,
Mehry

4

Hi Mehr,
Could you resolve the problem? I am using the flat-bottommed restraints on water molecules to restrict its movements through the membrane. I am using it in the following form in the topology file and using -DPOSRES_FB in the mdp file.

#ifdef POSRES_FB
[ position_restraints ]
;i function g (type) r (nm) k(kJ/mol)
1 2 5 -1 10000
#endif

Is it the correct form? Am I missing something?
It doesn’t seem to restrain the molecules after the simulation and also resulting in an unstable system. Any help would be appreciated.

5

You would need to provide a reference coordinate file to grompp with all water oxygens at the center of the membrane for the approach to work. If you use the initial coordinates of the water molecule as reference, I would expect you system to explode as you are pushing away all water molecules from the initial positions using a large force.

6

Hello,

As Hess explained, you need to define a new gro file as a reference. I applied the restraints only on the oxygen atoms of water molecules, and r=-2.5 nm, k=2.5 kJ/mol gave me the best results.

7

Thanks Hess and Mehr.

If I have understood correctly, I have to manually change the coordinates of the water oxygens and use it as the reference coordinate file. Hope it works this way.

8

Actually, you need to change all z coordinates of water oxygen atoms to the value of COM of the membrane in the reference gro file and apply the flat-bottomed restraint on them. I recommend using a simple bash script to do that.