Flat-bottom position restraints on water molecules

User discussions
,
1

GROMACS version: 2023.4
GROMACS modification: Yes/No
Here post your question
Hello everyone!
I want to apply position restraints on water molecules that are 0.55nm away from my protein.
Forcefield I’m using: charmm36-jul2022.ff
I’ve extracted just the protein and these water molecules into a separate .gro file. I generated an index file with group “restrain_water_ow” that contains indices of only the oxygen atoms of the water molecules, through which i generated an .itp file using the command:
gmx genrestr -f prot_0.55water.gro -n posres_index.ndx -o posres_water.itp -fc 1000 1000 1000

One thing to note is that, i’ve been using the same topology file as i’ve gotten post npt equilibration (npt.gro was the file from which prot_0.55water.gro was extracted). I just changed the no. of water molecules once i extracted prot_0.55water.gro

In this topol.top, after generated of the itp file, i included:

#ifdef POSRES_WATER
#include "posres_water.itp"
#endif

And in the mdp file i made sure to define = -DPOSRES_WATER

But it did not seem to work, I figured I’d rather include this in the tip3p.itp file since that was the water model that was used and included in the water topology.
This did not seem to work as well, I kept getting errors like index out of bounds- i thought it was because i was only including the oxygens and not the hydrogens and tried this as well, i was still getting the same error.

I changed the format of the .itp to:

[ POSRES_WATER ]
; atom  type  g  r  k
3318   2  1  2  1000

The error i’m encountering now is as follows:

WARNING 3911 [file posres_water.itp, line 3915]:
Too few parameters on line (source file
…/gromacs-2023.4/src/gromacs/gmxpreprocess/toppush.cpp,
line 1428)

multiple warnings for each line in the itp file.

What does this error signify? is it just an indentation error? Or are there not enough parameters?

I’m fairly new to gromacs, I apologize if the information i’ve provided is haphazard. Any help that would point me the right direction will be appreciated.

Thank you
Aarya

2

Dear @aarya

What do you want to do specifically? Restraint these water molecules to stay close to the protein? Or restrain them with respect to something else? I do not thing you can achieve what you want to do with position restraints of this kind.

The error you are getting here is related the the way GROMACS applies restraints. You are correct in the definition with the #ifdef statement at the end of the tip3p topology and the corresponding define = -DPOSRES_WATER in the mdp file. The first problem is the topology entry that should be [ position_restraints ] not [ POSRES_WATER ]. The other is that the positions restraints are applied to the atoms references in the topology. This means that the atom you pass in the [ position_restraints ] is the atom number within the topology of that molecule, not the one of the box. So in your water topology you will have three atoms, and in the restraints section you should specify which one of the three you want to restrain. Here you are saying 3318 2 1 2 1000, so you are telling GROMACS to apply the restraint to the atom 3318 of the topology, but your topology has 3 atoms only (water), and hence you get the out-f-bound index.

This way of setting restraints is very powerful because you can define them once in the topology and they are valid for all that type of molecules, and the same file can be used with different boxes as long as the topology of the molecule is the same. However, now you’ll see the point here, that is, you are applying this restraints to ALL the water molecules, not only to those you want. For that I am afraid you will need a different approach, like some position restraints defined through the pull code. Why do you want to restraint these specific water molecules?