[gmx-users] Re: Water molecules cannot be settled, why?
Mark.Abraham at anu.edu.au
Thu Jul 19 14:31:25 CEST 2012
On 19/07/2012 6:52 PM, Ladasky wrote:
> 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.
You can try a vacuum simulation on your solute if you want a further
test of whether its topology and starting configuration is OK.
> 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?
Different phases of the constraint algorithm, at least. Perhaps as you
say - one would have to read the code.
> 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.
When changing major version (4.0 -> 4.5), the most patch version is
always recommended (i.e. 4.5.5) lest you risk wasting time encountering
minor problems that have been fixed since 4.5 was released.
> 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.
Probably. You'd have to inspect the atom types in the resulting
processed topology (e.g. from grompp -pp) to be sure.
> 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.
There is no default force field. You need to make a choice based on what
you wish to observe.
> 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
Not the force field per se, but
http://www.gromacs.org/Documentation/Terminology/Blowing_Up can happen
any time an unstable numerical integration is attempted (ie. from a
configuration that doesn't belong to the ensemble described by the model
physics in the .mdp and .top files). See that link for discussion and
> 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?
Yes, particularly in early stages of equilibration, and definitely if
you are not using constraints. See above link.
> 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.
The usual coping strategies are contained in that link.
More information about the gromacs.org_gmx-users