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?
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.
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 :)