# Normal Mode Analysis, TIP3P Constraint Error
**URL:** https://gromacs.bioexcel.eu/t/normal-mode-analysis-tip3p-constraint-error/1012
**Category:** User discussions
**Created:** [October 27, 2020, 6:37pm UTC](https://gromacs.bioexcel.eu/t/normal-mode-analysis-tip3p-constraint-error/1012 "2020-10-27T18:37:19Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@antoszewski](https://gromacs.bioexcel.eu/u/antoszewski)
#### Post date: [October 27, 2020, 6:37pm UTC](https://gromacs.bioexcel.eu/t/normal-mode-analysis-tip3p-constraint-error/1012/1 "2020-10-27T18:37:19Z")
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
---
### Author:  [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)
#### Post date: [October 27, 2020, 8:23pm UTC](https://gromacs.bioexcel.eu/t/normal-mode-analysis-tip3p-constraint-error/1012/2 "2020-10-27T20:23:47Z")
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.
---
### Author:  [@antoszewski](https://gromacs.bioexcel.eu/u/antoszewski)
#### Post date: [October 28, 2020, 5:17pm UTC](https://gromacs.bioexcel.eu/t/normal-mode-analysis-tip3p-constraint-error/1012/3 "2020-10-28T17:17:52Z")
> [@jalemkul](#):
>
> ```auto
> #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
>
> ```
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