Re: [AMBER] selective qmshake?

From: Jason Swails <jason.swails.gmail.com>
Date: Tue, 26 Jul 2011 19:57:29 -0400

On Jul 26, 2011, at 2:29 PM, Brian Radak <radak004.umn.edu> wrote:

> Thanks Ross!
>
> It would appear that "noshakemask" DOES work when applied to QM atoms. For
> some reason I did not expect this behavior. Is this a bug that has become a
> feature? :P

Amber only has the latter. Though I would not have expected it to work either (but only out of implementation convenience). Glad it does though.

>
> My reason for using SHAKE is not to increase the time step. Of course that
> could be a good idea,

But only if you SHAKE all H- bonds. It only takes one to ruin the integrator :).

> but I'm not really sure what the constrained C-H bond
> lengths should be (although I'd imagine most reasonable guesses would give
> rather similar results). In any event I am trying to keep apples to apples
> comparison with my previous simulations and I would think that new
> constraints or changes in the time step should give one pause in that
> regard.

A perfectly understandable reason for doing what you're doing.

All the best,
Jason

>
> Brian
>
> On Tue, Jul 26, 2011 at 12:40 PM, Ross Walker <ross.rosswalker.co.uk> wrote:
>
>> Hi Brian,
>>
>> I have not looked at the code in a while but there is the option for
>> noshakemask in the &ctrl namelist which allows you to specify atoms to NOT
>> be shaken. This applies to ALL atoms in a MM simulation. I am not sure if
>> this will work for QM or not although if you have things setup already it
>> would be quick to test. Just add the QM atoms you do not want shaken to
>> this
>> mask and then run a short MD simulation and take a look at the bond
>> lengths,
>> are they changing or not?
>>
>> Essentially what happens in the code is that all the bonds in the QM
>> region,
>> between QM-QM atoms are deleted during two calls to setbon during setup.
>> One
>> for bonds without hydrogen and one with. You can find this code in set.f
>> inside sander. If the above does not work you could very easily go in and
>> hack this code for your specific system, i.e. inside the loop just have an
>> 'if atom number = blah then ...' That way you can have it delete the bond,
>> instead of preserving the bond if shake is in use, and I think you will
>> then
>> not be shaking that atom.
>>
>> Of all the suggestions tweaking the code is probably going to be the
>> easiest
>> if noshakemask does not work.
>>
>> Note, you don't gain from having atoms unshaken. E.g. if you have just a
>> single hydrogen that is not shaken you will need to run with 1fs timestep
>> or
>> less. The main reason for wanting to keep shake on some of the atoms is
>> because things like TIP3P were parameterized to be used ONLY with shake.
>>
>> Good luck,
>>
>> All the best
>> Ross
>>
>>> -----Original Message-----
>>> From: Brian Radak [mailto:radak004.umn.edu]
>>> Sent: Tuesday, July 26, 2011 8:57 AM
>>> To: AMBER Mailing List
>>> Subject: [AMBER] selective qmshake?
>>>
>>> I would like to run QM/MM MD where only *part* of the QM region (waters
>>> in
>>> TIP3P geometry) is constrained by SHAKE. That is, there are several C-
>>> H
>>> bonds which I would like to leave unconstrained if I can help it. The
>>> reason being, and perhaps people have opinions on this as well, that I
>>> would
>>> like to compare to other simulations with no QM waters in which I did
>>> not
>>> use SHAKE.
>>>
>>> My understanding from the manual is that qmshake is all or none wrt
>>> bonds
>>> with hydrogen. I can think of two doable, but arguably outlandish,
>>> ways to
>>> accomplish this. Will either of these work? Is there a better way to
>>> do
>>> this?
>>>
>>> 1.) Manually go into the parm7 and change the BONDS_INC_HYDROGEN and
>>> BONDS_WITHOUT_HYDROGEN tables... yuck! and will probably fail at first
>>>
>>> 2.) Create a new atom type for hydrogens bound to carbons (maybe I can
>>> steal
>>> one of the halide symbols?) so that leap does not recognize them as
>>> bonds
>>> with hydrogen. Of course this new atom type would have the same mass
>>> as
>>> hydrogen. Is that enough to trick leap?
>>>
>>> Thanks,
>>> Brian
>>>
>>>
>>> --
>>> ================================ Current Address
>>> =======================
>>> Brian Radak : BioMaPS
>>> Institute for Quantitative Biology
>>> PhD candidate - York Research Group : Rutgers, The State
>>> University of New Jersey
>>> University of Minnesota - Twin Cities : Wright-Rieman Hall
>>> 101
>>> Graduate Program in Chemical Physics : 610 Taylor Road,
>>> Department of Chemistry : Piscataway, NJ
>>> 08854-8066
>>> radak004.umn.edu :
>>> radakb.biomaps.rutgers.edu
>>> ====================================================================
>>> Sorry for the multiple e-mail addresses, just use the institute
>>> appropriate
>>> address.
>>> _______________________________________________
>>> AMBER mailing list
>>> AMBER.ambermd.org
>>> http://lists.ambermd.org/mailman/listinfo/amber
>>
>>
>> _______________________________________________
>> AMBER mailing list
>> AMBER.ambermd.org
>> http://lists.ambermd.org/mailman/listinfo/amber
>>
>
>
>
> --
> ================================ Current Address =======================
> Brian Radak : BioMaPS
> Institute for Quantitative Biology
> PhD candidate - York Research Group : Rutgers, The State
> University of New Jersey
> University of Minnesota - Twin Cities : Wright-Rieman Hall 101
> Graduate Program in Chemical Physics : 610 Taylor Road,
> Department of Chemistry : Piscataway, NJ
> 08854-8066
> radak004.umn.edu :
> radakb.biomaps.rutgers.edu
> ====================================================================
> Sorry for the multiple e-mail addresses, just use the institute appropriate
> address.
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber

_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Jul 26 2011 - 17:00:04 PDT
Custom Search