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

[Irem Altan]

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



[ 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