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

[Mark Abraham]

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:
>
> 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.

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
> crash?

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 
strategies.

> 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.

Mark