[gmx-users] Re: Water molecules cannot be settled, why?

[Ladasky]

Tsjerk Wassenaar wrote
> 
> Hi John,
> 
> Check where the unsettling water molecule is placed. If it's in the
> protein. that may be the cause of the problem. Otherwise, it's some of
> the other stuff you're doing, but rule out the simple things first.
> 
> 

Thank you Tsjerk.

I posted a reply several days ago, and I tried again earlier today.  My
first attempt appeared in the email digest but not on this web page.  My
second attempt hasn't appeared, and I expect that my ISP is responsible
(Yahoo started attaching ads to email signatures, which might set off spam
filters).  I'm now posting directly from Nabble.

I looked at two different simulations which crashed.  The unsettled water
molecules appear to be far from the protein in all cases.  Tentatively,
then, I would conclude that my protein is fine, because the errors are not
occurring in atoms which are proximal to the protein.  Instead, I suspect
there is something wrong with my protocol.  But here are the details, so you
can check my reasoning.

As the details of the protocol in the last message that I posted showed, I'm
using five CPU cores for my simulations.  So the domain decomposition in
mdrun is 5 x 1 x 1.  When the simulations crash, the error handling routine
outputs only the domain which contains the unsettled water molecule.  

In the first simulation that I examined, the domain was all solvent.  No
protein atoms were present.  The crash occurred on step 81,240 of my final
molecular dynamics run, and the atom that was identified as a problem was
number 45191.  So I looked at the two PDB files that were dumped, and
searched for atom 45191.  I think that the first file describes the system
immediately before the crash, and the second file shows the results of the
bad calculations, am I correct?  Here is a selected portion of each file,
showing atom 45191 and its neighbors:

step81240b_n4.pdb

ATOM  41009  OW  SOL     1      74.590  46.654  65.507
ATOM  41010  HW1 SOL     1      74.580  47.272  66.294
ATOM  41011  HW2 SOL     1      73.700  46.205  65.425
ATOM  43679  OW  SOL     1      71.579  49.011  70.894
ATOM  43680  HW1 SOL     1      71.730  49.197  69.924
ATOM  43681  HW2 SOL     1      71.013  49.735  71.288
ATOM  45191  OW  SOL     1      72.656  46.524  71.131
ATOM  45192  HW1 SOL     1      72.614  45.938  70.322
ATOM  45193  HW2 SOL     1      72.389  47.455  70.883
ATOM   6926  OW  SOL     1      75.926  46.701  73.317
ATOM   6927  HW1 SOL     1      75.024  46.674  72.886
ATOM   6928  HW2 SOL     1      76.564  47.194  72.724
ATOM  16016  OW  SOL     1      73.244  51.473  72.961
ATOM  16017  HW1 SOL     1      73.712  51.310  72.092
ATOM  16018  HW2 SOL     1      73.472  52.387  73.296

step81240c_n4.pdb

ATOM  41009  OW  SOL     1      74.582  46.669  65.502
ATOM  41010  HW1 SOL     1      74.582  47.320  66.260
ATOM  41011  HW2 SOL     1      73.689  46.221  65.445
ATOM  43679  OW  SOL     1    -66321836.00084864872.00046148104.000
ATOM  43680  HW1 SOL     1      71.706  49.216  69.937
ATOM  43681  HW2 SOL     1 1085546624.000-1490344960.000-271069120.000
ATOM  45191  OW  SOL     1    -176865.797532401.250-1222093.250
ATOM  45192  HW1 SOL     1    19264928.000-53927836.000125858544.000
ATOM  45193  HW2 SOL     1      72.351  47.488  70.877
ATOM   6926  OW  SOL     1      75.925  46.703  73.302
ATOM   6927  HW1 SOL     1      75.021  46.680  72.875
ATOM   6928  HW2 SOL     1      76.562  47.196  72.710
ATOM  16016  OW  SOL     1      73.240  51.478  72.957
ATOM  16017  HW1 SOL     1      73.728  51.296  72.103
ATOM  16018  HW2 SOL     1      73.442  52.407  73.266

WOW.  Something went seriously wrong with four atoms in that step. Are those
even coordinates?  For some reason, atom 45191 was selected as the one to
trigger the crash report.  I looked through the whole PDB file for other
problems, and I found one more oxygen atom in another water molecule which
had similarly wild coordinates.  It was half-way across the domain, not even
a neighbor of atom 45191.

I had a look at atom 45191 in PyMol, and there was nothing obviously wrong
with its position right before the crash. Stereo image here:

http://www.flickr.com/photos/15579975@N00/7553843726/

In the second crashed simulation, only a single abnormal atom was found in
the PDB file.  The domain containing that water molecule did contain a part
of the protein, but the water molecule was far from the protein.  Here's an
image of that domain (unsettled water molecule marked with an amber sphere):

http://www.flickr.com/photos/15579975@N00/7554770182/

All right, so what could be wrong with my protocols?  Since my problems
began, I have made exactly four changes that I can see:

1) I have upgraded from GROMACS 4.0 to 4.5.

2) I once used the -deuterate option in pdb2gmx, and I am presently trying
to avoid it.  I asked whether -deuterate has any effect on water molecules
in a previous post, or only on the hydrogens in proteins.  So far, I have
not been made to understand that genbox silently deuterates water hydrogens,
if it receives a deuterated protein as input.

3) I have started to write automated scripts to handle my work flow.  I
showed the sequence of commands in an earlier post.  I specify the
GROMOS45a3 force field, which I thought was the default (choice 1) for
pdb2gmx in GROMACS 4.0.  However, I am not certain of this.  GROMOS45a3 is
not the default choice for GROMACS 4.5, the AMBER03 force field is.

4) I got a warning message from GROMACS 4.5 that I did not observe with 4.0. 
"For proper integration of the Nose-Hoover thermostat, tau_t (0.1) should be
at least 20 times larger than nsttcouple*dt (0.01)."   I chose to double
tau_t to 0.2.  This was documented in my md.mdp file.

I understand that some force fields are more appropriate for certain types
of simulations than others, but should any force field result in an outright
crash?

I have never used a step time other than 2 femtoseconds.    I know that
people who want to model longer time spans will use, e.g., MARTINI and
larger step times, but I've never actually read a recommendation to try a
SMALLER step time.  Might this be necessary?

I could try different force fields, I could try -deuterate again, I could
reduce the step time.  But every single test that I want to try amounts to
approximately one day of simulation.  It could take me weeks to find a
problem, if I search blindly.

Thanks once again for any help you can provide.


--
View this message in context: http://gromacs.5086.n6.nabble.com/Re-Water-molecules-cannot-be-settled-why-tp4999302p4999553.html
Sent from the GROMACS Users Forum mailing list archive at Nabble.com.