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
Received on Tue Jul 26 2011 - 10:00:03 PDT
