[gmx-users] TIP5P - water molecules cannot be settled

[Justin Lemkul]

On 3/3/16 12:51 AM, Irem Altan wrote:
> Hi,
>
> I am using gromacs 5.1.2, and I’m trying to simulate a system that I’ve previously studied with TIP4P and TIP3P, this time using TIP5P. After I solvate the protein and add ions, I try to do energy minimization, but I get the error
>
> step 12: Water molecule starting at atom 5241 can not be settled.
>
> I have added to aminoacids.rtp the following:
>
> ; tip5p
> [ HO5 ]
>   [ atoms ]
>      OW   OW_tip5p      0.00     0
>     HW1   HW            0.24     0
>     HW2   HW            0.24     0
>     LP1   MW           -0.24     0
>     LP2   MW           -0.24     0
>   [ bonds ]
>      OW   HW1
>      OW   HW2
>
> and the following to aminoacids.hdb:
>
> HO5     1
> 4 11 HW OW
>
> My .top file includes the following:
>
> ; Include forcefield parameters
> #include "./amber99sb.ff/forcefield.itp"
>
> ; Include chain topologies
> #include "topol_Protein_chain_A.itp"
> #include "topol_Protein_chain_B.itp"
> #include "topol_Ion_chain_A2.itp"
> #include "topol_Ion_chain_B2.itp"
>
> ; Include water topology
> #include "./amber99sb.ff/tip5p.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 "./amber99sb.ff/ions.itp"
> #include "./tip5pcw.itp"
>
> [ system ]
> ; Name
> Protein in water
>
> [ molecules ]
> ; Compound        #mols
> Protein_chain_A     1
> Protein_chain_B     1
> Ion_chain_A2        1
> Ion_chain_B2        1
> HO5                96
> HO5                62
> SOL         5681
> NA               6
> CL               10
>
> the included file tip5pcw.itp allows for labeling crystal waters (HO5) and solvation waters (SOL) separately, and its contents are:
>
>
> [ moleculetype ]
> ; molname nrexcl
> HO5             2
>
> [ atoms ]
> ; id  at type     res nr  res name  at name  cg nr  charge    mass
>    1   OW_tip5p    1   HO5       OW       1       0        16.00000
>    2   HW_tip5p    1   HO5       HW1      1       0.241     1.00800
>    3   HW_tip5p    1   HO5       HW2      1       0.241     1.00800
>    4   MW          1   HO5       LP1      1      -0.241     0.00000
>    5   MW          1   HO5       LP2      1      -0.241     0.00000
>
>
> [ 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
>

This means that moleculetype HO5 is flexible, which is probably not stable. 
Most (if not all) water models are supposed to be rigid during MD.

Beyond that, a whole host of problems could lead to a crash, which may nor may 
not be related to the water model at all.  Combining TIP5P with protein force 
fields has caused problems before: http://redmine.gromacs.org/issues/1348

-Justin

>
>
> [ virtual_sites3 ]
> ; Vsite from                    funct   a               b               c
> 4      1       2       3       4        -0.344908262    -0.34490826     -6.4437903493
> 5      1       2       3       4        -0.344908262    -0.34490826     6.4437903493
>
>
> [ exclusions ]
> 1 2 3 4 5
> 2 1 3 4 5
> 3 1 2 4 5
> 4 1 2 3 5
> 5 1 2 3 4
>
> What could I be overlooking?
>
> Best,
> Irem
>

-- 
==================================================

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==================================================