AMBER Archive (2006)

Subject: Re: AMBER: SHAKE problem after minimization

From: Steven Winfield (saw44_at_cam.ac.uk)
Date: Wed Nov 22 2006 - 15:53:07 CST

Dear all,

Could this be a problem with the SHAKE algorithm itself? I've recently
been involved in coding up a general version of RATTLE for our MD
program, and found that the algorithm failed if two constraints on one
atom happened to have parallel gradients w.r.t. that atom's position,
which seems to be happening here to all the atoms by construction.

Looking at Frenkel and Smit: Understanding Molecular Simulations in the
Constraints section, I believe this could be a specific case of when the
matrix M becomes singular (this matrix is the one involving dot products
between the gradients of the constraints at times t and t+dt), and so
cannot be inverted to solve for the Lagrange multipliers.

Physically, the two constraints on the dummy atom are working against
each other and so as the algorithm progresses the dummy atom could
oscillate about the mid-point of the molecule but might not settle down.

Any comments anyone? I'm going to tackle this issue fairly soon in our
program (called Learn on the Fly if anyone's interested ;-) ), so if I'm
way off the mark please tell me!

Regards,
Steve.

Hayden Eastwood wrote:
> Sorry,
>
> In the previous correspondence when I said:
>
> "If HD is the dummy atom..."
>
> What I really meant was to say
>
> " if DH is the dummy atom..."
>
> i.e DH should replace my reference to HD
>
>
> hayden
>
> -----Original Message-----
> From: owner-amber_at_scripps.edu [mailto:owner-amber_at_scripps.edu] On Behalf Of
> Rachel
> Sent: 22 November 2006 19:37
> To: amber_at_scripps.edu
> Subject: Re: AMBER: SHAKE problem after minimization
>
> Dear Hayden,
>
> Thanks for your reply, yes, i do have some home-made residues, which are
> hydrogen molecules, i used a three-atom model, with two hydrogen atoms (type
> HD) and one dummy atom in-between the two hydrogen atoms (type DH), and the
> charge of both HD are - 0.475, and +0.95 on DH, and the following is the
> .frcmod file i used for this residue:
> ######################################
> # H2 molecule parameters
>
> MASS
> HD 1.008
> DH 1.008
>
> BOND
> HD-DH 150.00 0.742
> HD-HD 150.0 1.484
>
> ANGLE
> HD-DH-HD 500.0 180.0
> HD-HD-DH 0.000 180.0
>
> DIHEDRAL
>
> IMPROPER
>
> NONBON
> HD 0.000 0.000
> DH 2.920 0.2722
> #########################################
>
> As you can seen, I assumed the mass of DH is the same as HD (which is not
> correct). And as i wanted to keep these three atoms in the same line, so i
> added one extra bond between two HD atoms as suggested by Yong. So shall i
> increase the force constant for HD-DH-HD or HD-HD-DH or both?
>
> Thanks very much!
> Best regards,
> Rachel
>
>
>
> On 11/22/06, Hayden Eastwood <s0237717_at_sms.ed.ac.uk> wrote:
>
> Dear Rachel
>
> I have had similar problems to you when running dynamics with
> home-made
> residues. There is a defect in the amber force field with regards to
>
> treating various highly charged groups in bonding proximity. If you
> have a
> home made molecule then I have the following question:
>
> Do you have any Phosphate-oxygen bonding neighbours with very large
> opposite
> charges? Or any other bonding pairs with huge charge differences? If
> you do,
> this could be the cause of the crash. AMBER assumes that the bond
> angle
> barrier is large enough to prevent neighbouring connected atoms from
>
> collapsing into one another. However, in some situations (i.e. with
> oppositely highly-charged atom pairs) this barrier is overcome and
> the atoms
> are pulled down a pairwise electrostatic interaction slope until
> they
> overlap and cause your energy to go crazy.
>
> The cure: artificially increase the bond angle force constants in
> the force
> field parameter file and this should prevent them running into each
> other
> (although the dynamical behaviour of some of the affected groups
> will of
> course be "stiffer").
>
> I don't really understand the nuts and bolts of why the atoms
> descend into
> one another given that there should be van der waals terms that
> prevent
> this...Mike Crowley personally explained it to me once and whilst it
> made
> sense at the time the explanation now seems a little fuzzy!
>
> If it's not a home made residue with the above mentioned
> characteristics
> then I'd suggest you follow Ross's advice.
>
> Best
>
> Hayden
>
>
> -----------------------------------------------------------------------
> The AMBER Mail Reflector
> To post, send mail to amber_at_scripps.edu
> To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu
>
>
>
>
> -----------------------------------------------------------------------
> The AMBER Mail Reflector
> To post, send mail to amber_at_scripps.edu
> To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu