AMBER Archive (2007)

Subject: Re: AMBER: QMMM to follow reaction involving water

From: M. L. Dodson (mldodson_at_houston.rr.com)
Date: Tue Jul 31 2007 - 10:14:57 CDT

First off, thanks to Gustavo and Ross for some followups with a
lot of scientific content. I appreciate it.

Ross Walker wrote:
> Hi Bud,
>
>> I have been doing QMMM (using a development version of Amber10
>> sander furnished by Ross Walker). The basic design is steered
>> QMMM with the steered coordinate along a proposed scissile C-N
>> bond. I have been having a lot of trouble with this calculation
>> when I incorporate close waters into the QM region. These
>> problems seems to be of two kinds: (1) the QM water(s) diffuse
>> away from the enzyme active site, (2) the QM calculation (SCC-DFTB
>> w/o dispersion corrections) has great difficulty achieving SCF.
>
> To follow up on what Gustavo said the issue is definitely with the waters
> diffusing and your QM region essentially expanding as a bubble until it gets
> close to the bxo size. This occurs because of the way the direct space QM
> calculation is done. For stability all MM atoms that are within the cut off
> distance of ANY QM atom are included in the direct space sum for ALL QM
> atoms. At present using restraints is probably the only real solution
> available. An alternative would be to write some script that ran Amber for
> say 100 steps, then quit, checked the restart file to find the nearest X
> waters to the main QM region, wrote a new input file with a different list
> of iqmatoms and then reran sander for another 100 steps and so on.
>

The starting structure (this is a small protein + DNA) is from
0.2ns of constant vol dynamics (to remove restraints on the
solutes and equilibrate the energy, the box size having been fixed
with an earlier 50ps constant pressure dynamics run) followed by
1ns of free dynamics (NVE). Getting the IDENTITY of the closest
waters turned out to be a pain. So much so that I may try
developing a patch for ptraj to put out the identity of the
closest waters as well as their numbers. The QM region contained
the 2 closest waters (almost equidistant from the site of the
chemistry). The starting structure was well equilibrated to
FF99SB.

> Gustavo and I ultimately plan to automate this process within the code so
> that you can just run a QM simulation and have it recompute the nearest
> waters on the fly. This will obviously play havoc with the energies but
> should be okay for conformational analysis or to look at say proton transfer
> with a water molecule. As soon as we have it implemented I'll let you know
> so that you can test it.
>
> As for the SCF convergence DFTB does sometimes have problems but it
> shouldn't normally be that bad unless you give it a very bad structure - do
> you see the convergence problems right from the beginning or after so much
> simulation? If it is from the beginning then you could try running either
> with qm_ewald and qm_pme turned off for 1000 steps or so (since convergence
> is easier without this) to get you away from the initial bad structure and
> then turn it back on. Alternatively you could reduce the convergence
> criteria slightly to see if that helps - this causes inaccuracies with the
> gradients but again can be useful if you have bad structures. You might also
> want to set nstlim to 1 and then set verbosity=4 in the qmmm namelist and
> take a look at the convergence process in the DFTB calculation to see if it
> ever gets close to converging.
>
> There is also as you say the potential issue with the convergence problems
> being a function of the waters being a long way away. This is likely
> particularly true when running with PME since you start to get weird direct
> and reciprocal charge distributions and the error in the Ewald sum may be
> getting large. E.g. you can get things where you have an MM water wedged in
> between two QM waters.
>
> Is it always after the waters have started to diffuse that you get the
> convergence problems?
>

The SCF problems were immediately apparent, and the SCF problems
did not appear in earlier simulations without the waters. The
initial positions of the waters were several angstroms away from
the rest of the QM region. My idea was that, as the developing
charge on the C of the scissile C-N bond appeared (that bond
length was the steered coordinate), then one of the waters would
be "drawn in" electrostatically. So I think the QM waters started
out too far away, and this was never relieved. No matter what
course I choose, I think I am going to have to handle each of
these two waters in separate simulations.

> As for the restraint force for the time being I would just set it using
> standard harmonic restraints and a very weak force constant - as weak as you
> can get away with (<1.0?). You might want to run a full classical
> equilibration first to make sure the density is good before you add
> restraints. If you don't get a decent sphere of waters picked out though you
> may have issues with unrestrained waters tunnelling in and out of the
> restrained sphere. Ideally you'd probably want to restrain the sphere of
> classical waters around the outside of the chosen QM waters to try to stop
> the QM waters from escaping through the boundary rather than restraining the
> QM waters themselves - this would have less of an effect on the reaction
> pathway I would guess. I'm not sure how well it would work in practice
> though - something to try. What you ideally want is a water cap type
> restraint but keeping all the periodic boundary stuff - this would require
> some major code changes though.
>
> Good luck.
> Ross

I like the idea of restraints on the second shell of waters
relative to the QM region. Thanks for the hint about the
restraint force constant. I had just about come to the opposite
conclusion: a near vertical wall at some maximum allowed distance
with a 0 force constant interior to that. I.e., my thinking was
to let them rattle around in a box. 

-- 
M. L. Dodson
Email:	mldodson-at-houston-dot-rr-dot-com
Phone:	eight_three_two-five_63-386_one
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu