[gmx-users] Re: Water molecules cannot be settled, why?
blind.watchmaker at yahoo.com
Thu Jul 19 10:52:33 CEST 2012
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:
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
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:
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):
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
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.
More information about the gromacs.org_gmx-users