How to treat specific water molecules as ligand?

GROMACS version: 2021
GROMACS modification: No
Here post your question

I have two conserved water molecules, which should be considered as ligands interacting with the protein. The [ molecules ] section in the topology file is like this. The first SOL group has 2 molecules, which represents the 2 conserved water; the second SOL group is the solvent water. But these two SOL groups essentially share the same water topology (i.e. tip3p.itp).

[ molecules ]
; Compound        #mols
Protein_chain_A     1
lig                 1
SOL                 2
SOL         7457
NA               14
CL               8

Justin’s protein-ligand tutorial tells us we can create a new index that includes protein and ligand, so as to separate it from the rest in the tc-grps.

But I think I also need to add a new [ moleculetype ] specifically for the two conserved water, and use the same POSRES position restraints as protein

#ifdef POSRES
; Position restraint for each water oxygen
[ position_restraints ]
;  i funct       fcx        fcy        fcz
   1    1       1000       1000       1000
#endif

to constrain the 2 conserved water, so that define = -DPOSRES will constrain not only the protein but also the 2 conserved water. Can I ask what I need to modify to this new [ moleculetype ]? Is there a explanation for the default tip3p.itp (e.g. #ifndef FLEXIBLE, [ settles ]) ?

Hi,
If you want to position restrain the 2 conserved water molecule you can copy the water topology in your working dir, change the name of the itp file and change the molecular name in the file. Update the [molecules] directive in topol.top and include the new water topology in the correct order in the topol.top together with a new POSRES definition (with a different name from the water one). Then add the new POSRES name in mdp file.

But I wonder if you really need this, since (I guess) you will release the restrain in data production and the water molecules usually exchange with another water molecule.

Kind regards
Alessandra

Hi @alevilla, thank you very much. I think I am following your instructions:

I save the default tip3p.itp as tip3p1.itp and put into my working directory, replace all the SOL strings with FIX. But I am not sure how to deal with the [ settles ] and [ exclusions ] sections. Based on this thread, I replace the [ settles ] and [ exclusions ] sections with [ constraints ], is that correct? But that thread is for “flat-bottom restraint”, what if I want to use the normal “position restraint”?


[ moleculetype ]
; molname	nrexcl
FIX		2

[ atoms ]
; id  at type     res nr  res name  at name  cg nr  charge    mass
  1   OW          1       FIX       OW       1      -0.834    16.00000
  2   HW          1       FIX       HW1      1       0.417     1.00800
  3   HW          1       FIX       HW2      1       0.417     1.00800

#ifndef FLEXIBLE
;;[ settles ]
; OW	funct	doh	dhh
;;1       1       0.09572 0.15139
;;
;;[ exclusions ]
;;1	2	3
;;2	1	3
;;3	1	2

[ constraints ]
1 2 1 0.09572
1 3 1 0.09572
2 3 1 0.15139

#else

[ bonds ]
; i     j       funct   length  force_constant
1       2       1       0.09572 502416.0   0.09572        502416.0 
1       3       1       0.09572 502416.0   0.09572        502416.0 
        

[ angles ]
; i     j       k       funct   angle   force_constant
2       1       3       1       104.52  628.02      104.52  628.02  

#endif

Yes, I have also updated the topology file to include the new water molecule. But I use the same POSRES for the conserved water as protein.

And yes, in the nvt.mdp and npt.mdp files, it has

define		= -DPOSRES
tc-grps		= protein_lig_conserve-water  water_ions_without-conserve-water

while in the md.mdp file, it does not have the restraints and has
tc-grps = protein_lig Water_and_ions

leave [ settles ] and [ exclusions ] how they were originally. You do not want to change the description of water, I understood that you want just position restrain the conserved water molecules.

\Alessandra

Sorry I do not think so. If the [ settles ] and [ exclusions ] sections are kept in the conserved water, this error occurs:

Fatal error:
The [molecules] section of your topology specifies more than one block of
a [moleculetype] with a [settles] block. Only one such is allowed. If you
are trying to partition your solvent into different *groups* (e.g. for
freezing, T-coupling, etc.) then you are using the wrong approach. Index
files specify groups. Otherwise, you may wish to change the least-used
block of molecules with SETTLE constraints into 3 normal constraints.

Any time you wish to treat some subset of waters differently, yes, you do need to replace [settles] with [constraints]. The exclusions are the same either way.

Thank you, @jalemkul, do you mean I should only choose the words of [ settles ] to [ constraints ], or should I change the [ settles ] section to [ constraints ] section?

The below is only changing the words. If the section needs to change, what is the content for the [ constraints ] section?

[ moleculetype ]
; molname	nrexcl
FIX		2

[ atoms ]
; id  at type     res nr  res name  at name  cg nr  charge    mass
  1   OW          1       FIX       OW       1      -0.834    16.00000
  2   HW          1       FIX       HW1      1       0.417     1.00800
  3   HW          1       FIX       HW2      1       0.417     1.00800

#ifndef FLEXIBLE

[ constraints ]
; OW	funct	doh	dhh
1       1       0.09572 0.15139

[ exclusions ]
1	2	3
2	1	3
3	1	2

#else

[ bonds ]
; i     j       funct   length  force_constant
1       2       1       0.09572 502416.0   0.09572        502416.0 
1       3       1       0.09572 502416.0   0.09572        502416.0 
        

[ angles ]
; i     j       k       funct   angle   force_constant
2       1       3       1       104.52  628.02      104.52  628.02  

#endif

You need to replace SETTLE with three constraints.

Thank you Justin, would it be possible if you can show me an example for the “three constraints” for “position restraint”? Is it the same for “position restraint” and “flat-bottom restraint” (earlier thread)?

[ constraints ]
; ai aj funct length_A length_B
1 2 1 0.09572
1 3 1 0.09572
2 3 1 0.15139

Replacing SETTLE by three constraints is not dependent upon the type of restraint being used (in fact it’s not technically related at all). What you have listed for constraints is correct.

Thank you very much!