Gromacs not reading POSRES even though it is pointed to in topology

User discussions
, ,
1

GROMACS versions: 2022 and 2021
GROMACS modification: No

Hi gromacs,
I am running a simulation of a protein with water and some other ingredients and I am having trouble restraining the protein heavy atoms, even though I specify
#ifdef POSRES
#include “posre_Protein_chain_A2.itp”
#endif

in the protein topology component itp file and define = -DPOSRES in the mdp file. The posre_* file has all the heavy atoms of the protein listed with restraint strengths. Below is my topology for your reference. The Urea.itp and TMAO.itp files are defined elsewhere and included here, as are the other itp files. I can post them too if needed.

Thank you very much!
Dan

; Include forcefield parameters
#include “/gpfs/loomis/project/hammes_schiffer/as3586/shslocal/apps/gromacs/share/gromacs/top/amber14sb_OL15.ff/forcefield.itp”

[ atomtypes ]
; name mass charge ptype sigma epsilon
CUR 12.0107 0.921 A 0.37700 0.41700
HUR 1.00800 0.285 A 0.15800 0.08800
OUR 15.99940 -0.675 A 0.31000 0.56000
NUR 14.0067 -0.693 A 0.31100 0.50000
OTM 15.99940 -0.910 A 0.32660 0.63790
NTM 14.0067 0.700 A 0.29260 0.83600
CTM 12.01100 -0.260 A 0.36000 0.28260
HTM 1.00800 0.110 A 0.21010 0.07730

[ nonbond_params ]
; i j func sigma epsilon (urea-TMAO epsilon scaled by 1.10)

CUR CTM 1 0.36840 0.37761
CUR HTM 1 0.28144 0.19749
CUR OTM 1 0.35090 0.56733
CUR NTM 1 0.33213 0.64948
HUR CTM 1 0.23849 0.17347
HUR HTM 1 0.18220 0.09072
HUR OTM 1 0.22716 0.26062
HUR NTM 1 0.21501 0.29836
OUR CTM 1 0.33407 0.43760
OUR HTM 1 0.25521 0.22886
OUR OTM 1 0.31819 0.65745
OUR NTM 1 0.30117 0.75264
NUR CTM 1 0.33460 0.41349
NUR HTM 1 0.25562 0.21626
NUR OTM 1 0.31871 0.62123
NUR NTM 1 0.30166 0.71118

; i j func sigma epsilon (TMAO-TMAO epsilon scaled by 0.90)

OTM OTM 1 0.32660 0.57411
OTM NTM 1 0.30913 0.65724
OTM CTM 1 0.34289 0.38212
OTM HTM 1 0.26195 0.19985
NTM NTM 1 0.29260 0.75240
NTM CTM 1 0.32456 0.43745
NTM HTM 1 0.24794 0.22879
CTM CTM 1 0.36000 0.25434
CTM HTM 1 0.27502 0.13302
HTM HTM 1 0.21010 0.06957

#include “/home/dk758/project/ureaTMAO/NetzM/Urea.itp”
#include “/home/dk758/project/ureaTMAO/NetzM/TMAO.itp”

; Include chain topologies
#include “topol_Ion_chain_A.itp”
#include “topol_Protein_chain_A2.itp”

; Include water topology
#include “/gpfs/loomis/project/hammes_schiffer/as3586/shslocal/apps/gromacs/share/gromacs/top/amber14sb_OL15.ff/tip3p.itp”

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

; Include topology for ions
#include “/gpfs/loomis/project/hammes_schiffer/as3586/shslocal/apps/gromacs/share/gromacs/top/amber14sb_OL15.ff/ions.itp”

[ system ]
; Name
Built with Packmol

[ molecules ]
; Compound #mols
Ion_chain_A 1
UR 28
TMO 28
Protein_chain_A2 1
SOL 1351

2

What makes you think the restraints aren’t working?

3

It’s because the heavy atoms move appreciably the traj.trr from the simulation. They clearly are not restrained…they are moving around :-)

4

If you are restraining only the protein, are the amino acid residues in the beginning of the gro file? E.g., if there are 100K atoms (protein + water + ions + any other species) but your posres file has only 10K atoms, grompp will spit a message saying that only the first 10K atoms will be restrained. (you might have missed it if using -maxwarn). If that is the case, try renumbering the protein chain residues to appear at the beginning.

Another point is: what is the force constant being used? Too small fc like 2 or 5 will be just too small to hold the protein.

5

On second thought, I have reason to believe the restraints ARE working, because when I take a system minimized with restrained protein and then try to run dynamics without minimizing the whole system first, the simulation crashes, but after I minimize everything things are fine. I just didn’t expect the atoms to move AT ALL when restrained, but I guess 1000 kj/mol is not that strong of a restraint, and that’s the default. It is probably ok that they vibrate a bit. Is that correct?

6

Yes, a restraint penalizes motion. It doesn’t prevent it.

7

I am coming from mostly working with OpenMM where you often set your own restraints and I used to set them very high, probably a lot higher than is default in GROMACS, that’s why I had this misunderstanding. Thank you for helping me!