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

[Irem Altan]

On 3/3/16 5:16 PM, Irem Altan wrote:
Hi,

Thanks for your reply. I have checked step12b.pdb and step12c.pdb, and I didn’t observe a deformation in the shape of the molecule in question. It doesn’t look like there is any overlap or clash either. I have also tried to comment out the #ifndef FLEXIBLE block and only leave the settles block uncommented. It results in the same error. How can I check if it is the combination of water models and protein force fields that cause the problem?


Do people use AMBER + TIP5P often?  Does it run with AMBER + TIP3P?

I have no knowledge of AMBER + TIP5P, but I have used AMBER + TIP3P, TIP4P, SPC, and SPC/E without problems.

Also realize that water not being settled doesn't necessarily mean you need to troubleshoot the water, it's just that it's the first thing that fails when something goes haywire.

What would you recommend me to do in order to troubleshoot the problem?


-Justin

Best,
Irem

On Mar 3, 2016, at 5:01 PM, Justin Lemkul <jalemkul at vt.edu<mailto:jalemkul at vt.edu><mailto:jalemkul at vt.edu>> wrote:



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