[gmx-users] Segmentation fault while using external water model

[Irem Altan]

Hi,

I am still having trouble using TIP4P2005. When I start running my MD simulation, it “runs” for a minute or so (although it does not start the simulation itself), and then crashes with a segmentation fault.

I have (modified) copy of the amber99sb forcefields in my working directory, and I have made the following changes:

- added a tip4p2005.itp file with the correct parameters
  — named OW as OW_tip4p2005 and HW as HW_tip4p2005
- added the LJ parameters to ffnonbonded.itp:
OW_tip4p2005 8      16.00    0.0000  A   3.15890e-01  7.74908e-01
        HW_tip4p2005 1       1.008   0.0000  A   0.00000e+00  0.00000e+00
- added the new atom types to atomtypes.atp:
       HW_tip4p2005    1.00080 ; tip4p/2005 HW
       OW_tip4p2005   16.00000 ; tip4p/2005 OW
- added the following line to watermodels.dat
       tip4p2005 TIP4P/2005 tip4p/2005

My system consists of a protein in its unit cell, unstripped of its crystal waters. Below is how my topology file looks like (after solvation, adding ions, and energy minimization):

; 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_Other_chain_A2.itp"
#include "topol_Other_chain_B2.itp"

; Include water topology
#include "./amber99sb.ff/tip4p2005.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 "./tip4p2005cw.itp"

[ system ]
; Name
Protein in water

[ molecules ]
; Compound        #mols
Protein_chain_A     1
Protein_chain_B     1
Other_chain_A2      1
Other_chain_B2      1
HO4                96
HO4                62
SOL         5763
NA               6
CL               10

where tip4p2005cw.itp is the same file, only with [settles] removed, and molname changed to HO4 for labeling purposes.

What could be the problem? (PS, see below for the contents of tip4p2005.itp and tip4p2005cw.itp)

Best,
Irem

=======================================================================================================
tip4p2005.itp (taken from here and slightly modified: https://github.com/wesbarnett/trappeua/blob/master/trappeua.ff/tip4p2005.itp )
————————

 TIP4P/2005 Water Model

; J. L. F. Abascal and C. Vega, J. Chem. Phys 123 234505 (2005)
; http://dx.doi.org/10.1063/1.2121687


[ moleculetype ]
; molname    nrexcl
  SOL        2

[ atoms ]
; nr type resnr residu atom cgnr charge    mass
  1  OW_tip4p2005   1     SOL    OW   1     0.0000   15.994
  2  HW_tip4p2005   1     SOL    HW1  1     0.5564    1.008
  3  HW_tip4p2005   1     SOL    HW2  1     0.5564    1.008
  4  MW             1     SOL    MW   1    -1.1128    0.0

#ifndef FLEXIBLE
[ settles ]
; i  funct   doh       dhh
  1  1       0.09572   0.15139
#else
[ bonds ]
; i     j     funct   length  force.c.
  1     2     1       0.09572 502416.0
  1     3     1       0.09572 502416.0

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

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

; The position of the dummy is computed as follows:
;
;         O
;
;         D
;
;   H           H
;
; a = b = (1/2) * (distance(OD)  / [ cos (angle(DOH)) * distance(OH) ] )
; = (1/2) * (0.01546 nm    / [ cos (52.26 deg)  * 0.09572 nm   ] )
; = 0.13193828
; Dummy pos x4 = x1 + a*(x2-x1) + b*(x3-x1)

[ virtual_sites3 ]
; Vsite from                    funct   a               b
  4     1 2 3 1 0.13193828 0.13193828


==============
tip4p2005cw.itp
————————

; TIP4P/2005 Water Model

; J. L. F. Abascal and C. Vega, J. Chem. Phys 123 234505 (2005)
; http://dx.doi.org/10.1063/1.2121687


[ moleculetype ]
; molname    nrexcl
  HO4        2

[ atoms ]
; nr type resnr residu atom cgnr charge    mass
  1  OW_tip4p2005   1     SOL    OW   1     0.0000   15.994
  2  HW_tip4p2005   1     SOL    HW1  1     0.5564    1.008
  3  HW_tip4p2005   1     SOL    HW2  1     0.5564    1.008
  4  MW             1     SOL    MW   1    -1.1128    0.0

;#ifndef FLEXIBLE
;[ settles ]
; i  funct   doh       dhh
;  1  1       0.09572   0.15139
;#else
[ bonds ]
; i     j     funct   length  force.c.
  1     2     1       0.09572 502416.0
  1     3     1       0.09572 502416.0

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

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

; The position of the dummy is computed as follows:
;
;         O
;
;         D
;
;   H           H
;
; a = b = (1/2) * (distance(OD)  / [ cos (angle(DOH)) * distance(OH) ] )
; = (1/2) * (0.01546 nm    / [ cos (52.26 deg)  * 0.09572 nm   ] )
; = 0.13193828
; Dummy pos x4 = x1 + a*(x2-x1) + b*(x3-x1)

[ virtual_sites3 ]
; Vsite from                    funct   a               b
  4     1 2 3 1 0.13193828 0.13193828