[gmx-users] OPLS/AA + TIP5P, anybody?

[gigo]

Dear GMXers,
Since I am interested in interactions of lone electron pairs of water 
oxygen within the active site of an enzyme that I work on, I decided to 
give TIP5P a shot. I use OPLSAA. I run into troubles very fast trying to 
minimize freshly solvated system. I found on the gmx-users 
(http://lists.gromacs.org/pipermail/gmx-users/2008-March/032732.html) 
that cg and constraints don't go together when TIP5P is to be used - 
thats OK. It turned out, however, that I was not able to minimize my 
protein even with steepest descent. The system minimizes with TIP4P 
pretty well (emtol=1.0). In the meantime I tried to minimize short 
peptide - 10aa, did not work as well. What happens? The LP of water used 
to get too close to positively charged hydrogens (without VDW radius) on 
arginine. It looks like this:

Step=  579, Dmax= 8.0e-03 nm, Epot= -1.40714e+05 Fmax= 1.20925e+04, 
atom= 171
Step=  580, Dmax= 9.6e-03 nm, Epot= -1.41193e+05 Fmax= 8.13923e+04, 
atom= 171
Step=  581, Dmax= 1.1e-02 nm, Epot= -1.43034e+05 Fmax= 1.03648e+06, 
atom= 11181
Step=  585, Dmax= 1.7e-03 nm, Epot= -1.46878e+05 Fmax= 4.23958e+06, 
atom= 11181
Step=  587, Dmax= 1.0e-03 nm, Epot= -1.49565e+05 Fmax= 9.43285e+06, 
atom= 11181
Step=  589, Dmax= 6.2e-04 nm, Epot= -1.59042e+05 Fmax= 3.55920e+07, 
atom= 11181
Step=  591, Dmax= 3.7e-04 nm, Epot= -1.69054e+05 Fmax= 7.79944e+07, 
atom= 11181
Step=  593, Dmax= 2.2e-04 nm, Epot= -1.85575e+05 Fmax= 2.27640e+08, 
atom= 11181
Step=  595, Dmax= 1.3e-04 nm, Epot= -2.35034e+05 Fmax= 5.88938e+08, 
atom= 17181
Step=  597, Dmax= 8.0e-05 nm, Epot= -2.39154e+05 Fmax= 1.22615e+09, 
atom= 11181
Step=  598, Dmax= 9.6e-05 nm, Epot= -2.67157e+05 Fmax= 1.96782e+09, 
atom= 11181
Step=  600, Dmax= 5.8e-05 nm, Epot= -4.37260e+05 Fmax= 1.08988e+10, 
atom= 11181
Step=  602, Dmax= 3.5e-05 nm, Epot= -4.65654e+05 Fmax= 1.29609e+10, 
atom= 11181
Step=  604, Dmax= 2.1e-05 nm, Epot= -1.17945e+06 Fmax= 1.31028e+11, 
atom= 11181
Step=  607, Dmax= 6.3e-06 nm, Epot= -3.07551e+06 Fmax= 6.04297e+11, 
atom= 11181
Step=  610, Dmax= 1.9e-06 nm, Epot= -4.26709e+06 Fmax= 1.61390e+12, 
atom= 11181
Step=  611, Dmax= 2.3e-06 nm, Epot= -4.39724e+06 Fmax= 2.14416e+12, 
atom= 11181
Step=  613, Dmax= 1.4e-06 nm, Epot= -1.27489e+07 Fmax= 1.03223e+13, 
atom= 17181
Step=  614, Dmax= 1.6e-06 nm, Epot= -5.23118e+06 Fmax= 3.18465e+12, 
atom= 11181
Energy minimization has stopped, but the forces havenot converged to the
(...)

In this example atom 171 is HH21 of ARG, and 11181 is oxygen of water 
that got close to this ARG. Sometimes the epot turns nan at the end. If 
you would like to reproduce, I put the peptide.pdb, the mdp file and the 
running script at http://shroom.ibb.waw.pl/tip5p . If anybody have any 
suggestions how to minimize (deep) with OPLSAA + TIP5P in gromacs (4.6.3 
preferably...) without constraining bond lengths (which is also 
problematic), I will be very very grateful.
Best,

Grzegorz Wieczorek