Normal Mode Analysis, TIP3P Constraint Error

User discussions
1

GROMACS version: 2019.4, double precision
GROMACS modification: Patched with PLUMED 2.5.5

Hey all,

I have been running normal mode analysis on a hexameric protein complex whose X-ray crystal structure contains both stabilizing ligands and crystallographic water. Initially, I was running nma of the system stripped of water, but with the stabilizing ligands still present; however, during the (thorough) CG energy minimization, some essential elements of the protein’s secondary structure (namely, six specific alpha helices) would melt without the crystallographic water that surrounded them in the X-ray structure. I think this affected the resultant nma eigenvectors, as the dominant ones involve the motions of these now-melted former helices.

So, I instead wanted to run nma on the full system, including the crystallographic water (TIP3P). The energy minimization went well, with the minimized structure containing stable alpha helices. However, when trying to actually build the hessian matrix, I receive the following error:

Fatal error: Constraints present with Normal Mode Analysis, this combination is not supported

In the .mdp file, I defined -DFLEXIBLE, and explicitly set constraints = none, so I’m not exactly sure where this constraint error is coming from. Does the water inherently have some constraints even with the above flags turned on? Do I need to add something to my .top file or my TIP3P.itp file to make the define = -DFLEXIBLE command work correctly? Is there any way to run nma with a limited amount of water present? The water is the only thing changed from the previous nma runs which worked smoothly, so I imagine that it is the problem here. Reproduced below are my gromacs commands, relevant .mdp files, .top file, and TIP3P.itp file.

Thanks for the help!
Adam Antoszewski
PhD Student, University of Chicago

GROMACS Commands ========================
#Minimize
gmx_mpi_d grompp -f Minimization.mdp -c prot.gro -p Topology.top -o prot_mincg -maxwarn 1
mpirun -np 28 gmx_mpi_d mdrun -v -deffnm prot_mincg -ntomp 1
#NMA
gmx_mpi_d grompp -f hessian.mdp -c prot_mincg.gro -t prot_mincg.trr -p Topology.top -o hes -maxwarn 1
#The following command throws the constraint error
mpirun -np 28 gmx_mpi_d mdrun -v -deffnm hes -ntomp 1
gmx_mpi_d nmeig -f hes.mtx -s hes.tpr -v prot_evec.trr

Minimization.mdp ========================
define = -DFLEXIBLE
integrator = cg
emtol = 0.001
emstep = 0.00005
nsteps = 4000000
nstlist = 10
nstcgsteep = 100
cutoff-scheme = Verlet
rlist = 1.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 1.0
rvdw = 1.2
coulombtype = pme
rcoulomb = 1.2
;
constraints = none

hessian.mdp ========================
define = -DFLEXIBLE
integrator = nm
emtol = 0.001
emstep = 0.00005
nsteps = 1000000
nstlist = 10
cutoff-scheme = Verlet
rlist = 1.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 1.0
rvdw = 1.2
coulombtype = pme
rcoulomb = 1.2
;
constraints = none

TIP3P.itp ========================
[ moleculetype ]
; name nrexcl
TIP3 2

[ atoms ]
; nr type resnr residu atom cgnr charge mass
1 OT 23 TIP3 OH2 1 -0.834 15.9994 ; qtot -0.834
2 HT 23 TIP3 H1 2 0.417 1.0080 ; qtot -0.417
3 HT 23 TIP3 H2 3 0.417 1.0080 ; qtot 0.000

[ settles ]
; i j funct length
1 1 9.572000e-02 1.513900e-01

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

Topology.top ========================
; Include forcefield parameters
#include “toppar/charmm36.itp”
#include “toppar/PROA.itp”
#include “toppar/PROC.itp”
#include “toppar/PROE.itp”
#include “toppar/PROG.itp”
#include “toppar/PROI.itp”
#include “toppar/PROK.itp”
#include “toppar/HETA.itp”
#include “toppar/HETB.itp”
#include “toppar/HETC.itp”
#include “toppar/HETE.itp”
#include “toppar/HETF.itp”
#include “toppar/HETG.itp”
#include “toppar/HETH.itp”
#include “toppar/HETI.itp”
#include “toppar/HETJ.itp”
#include “toppar/HETK.itp”
#include “toppar/POT.itp”
#include “toppar/CLA.itp”
#include “toppar/TIP3.itp”

[ system ]
; Name
Title

[ molecules ]
; Compound #mols
PROA 1
PROC 1
PROE 1
PROG 1
PROI 1
PROK 1
HETA 1
HETB 1
HETC 1
HETE 1
HETF 1
HETG 1
HETH 1
HETI 1
HETJ 1
HETK 1
TIP3 31
TIP3 29
TIP3 24
TIP3 35
TIP3 17
TIP3 29
TIP3 24
TIP3 27
TIP3 32
TIP3 26
TIP3 25
TIP3 32

2

Your use of -DFLEXIBLE does nothing because your topology does not contain an #ifdef FLEXIBLE... block. It is not the standard TIP3P topology, likely because CHARMM-GUI inputs reflect conventional usage and TIP3P should always be treated as rigid. If you want to remove constraints from TIP3P, see the standard tip3p.itp approach that replaces SETTLE with normal bonds.

#ifdef FLEXIBLE

#ifdef ORIGINAL_TIP3P
[ bonds ]
; i j   funct   length  force.c.
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.c.
2    1  3   1   104.52  628.02  104.52  628.02

#else
;CHARMM TIP3p
[ bonds ]
; i j   funct   length  force.c.
1   2   1   0.09572 376560.0 0.09572    376560.0
1   3   1   0.09572 376560.0 0.09572    376560.0

[ angles ]
; i  j  k   funct   angle   force.c.
2    1  3   1   104.52  460.24  104.52  460.24
#endif

I suspect the issue is more related to the force field being inadequate for a desolvated (vacuum) system, in which case you may want to minimize with dihedral or distance restraints to preserve the geometry, then perform NMA.

3

Thanks so much for the quick and thoughtful response, this seems to have worked well. I should have noticed that in order for the flexible flag to work, it should actually do something!

And your suggestion about minimizing with restraints is a good one, I’ll try that as well.

Thanks again,

Adam