[gmx-developers] Question on SHAKE equations
Anton Feenstra
k.a.feenstra at vu.nl
Tue Jun 11 11:36:17 CEST 2013
On 10/06/13 20:36, "Pablo García Risueño" wrote:
> Dear Anton
>
> Thank you very much for your explanations. They are very important for my
> research.
>
> I regret I still need to ask more. So what is the way of operationg for
> Gromacs when setting h_angles in the constraint block? Is it working with
> P-LINCS, or with the method of your JCC paper?
Constraints are handled by either SHAKE or LINCS for any other bond (or
angle), as chosen by the user. The only exception is SETTLE which is
defined explicitly in the SPC topology file.
Virtual sites (formerly called dummy atoms) are only *used* when
explicitly defined (by the user), or generated by pdb2gmx with the
-vsite option, IIRC. These constructions are quite distinct from
constraints, as v-sites have no mass but are only a point for an
interaction function (charge and/or LJ). (Note that for mass
conservation in a molecule, the mass of the v-site atom should be added
elsewhere.)
You can read all of this and some more, I think, in the paper and the
relevant section(s) in the Gromacs manual.
> Moreover, I cannot see why the C-N-C angle should be constrained when
> constraining the C1-N-H and the C2-N-H angles (see scheme below).
> H
> |
> C1--N--C2
>
> In principle C1-N-H and C2-N-H lie in different planes, they are two
> 'triangles' which share one edge (this edge is the H-N bond), but which
> can rotate around it in 3D, isn't this true?
Yes, but not for the protein backbone, which is planar around the amide
group (imposed by impropers). But also for e.g. a C-CH2-C group,
constriaining all C-C-H angles (4 of them) would effecively also
constrain the C-C-C angle.
> I know of the limitations of LINCS and P-LINCS, but I was wondering if
> there would be any problem if inverting the cosntraint matrix by using
> direct Gaussian elimination (which should not carry any instability
> problem). I don't know if constraining all the H angles a la Gromacs is
> actually leading to some impossible geometry.
For this, I'm bumping back our discussion to the gmx-dev list.
>> Dear Pablo,
>>
>>
>> There is two issues: limitations in LINCS and couplings between H bond
>> angles and other bond angles. My guess is you are referring to the
second,
>> but for sake of being complete I will also go into the first.
>>
>> 1) LINCS solves constraints by using an approximation to the
inversion of
>> the constraints matrix; this relies on the assumption that this
matrix is
>> sparse. For bond constraints only, that is true. But when also (all)
>> angles are constrained, the matix contains to many cross-terms. This
>> happens already in butane (with united atoms for all the - aliphatic -
>> hydrogens; so effecitvely only four atoms in the molecule) with three
>> bonds and two associate bond angles (if they are all constrained). For
>> more details, certainly on the technical and/or mathematical side,
you had
>> best ask Berk.
>>
>> 2) in the case of, for example, a protein backbone, a problem can arise
>> when all H bond angles are constrained. At the N-H, we will
constrain both
>> C-N-H angles, effectively also constraining the C-N-C angle. This
has been
>> show a long time ago to give wrong behaviour of protein dynamics (before
>> 1998 at least, but I do not have a reference).
>>
>> Regarding your question below, no, as far as I know Gromacs does not use
>> v-sites by default for H atom bond angles. However, pdb2gmx has an
option
>> to generate v-sites for H atoms in proteins and in several other
building
>> blocks.
>>
>> Hope this answers your questions.
>>
>>
>> P.s., please call me Anton.
>>
>> Groetjes,
>>
>>
>> Anton.
>>
>>
>>
>> Op 10 jun. 2013 om 16:55 heeft "Pablo García Risueño"
>> <Risueno at physik.hu-berlin.de> het volgende geschreven:
>>> Dear prof Feenstra
>>>
>>> Thank you very much for your reply. Then, I should understand that, in
>>> order to constrain H bond angles, Gromacs is not using P-LINCS, but the
>>> technique explained in J. Comput. Chem. 20 (8), 786-798 by default, is
>>> this correct? Isn't this having a non-negligible impact on observable
>>> quantities?
>>>
>>> I afraid I cannot understand completely the problem of the
>>> incompatibility
>>> between bond angle constraints yet. Why should this happen? I guess
>>> that
>>> we can freeze every H bond angle to a concrete value at the beginning
>>> of
>>> the simulation, and then to keep it throughout the simulation. Looking
>>> at
>>> the protein topologies, I'd swear that no incompatibility should arise,
>>> since a hydrogen atom is always at a extreme (e.g., there is never
>>> C-H-C).
>>> Isn't this true?
>>>
>>> Thank you very much for your help. Best regards
>>>
>>>
>>>
>>>> On 10/06/13 14:20, "Pablo García Risueño" wrote:
>>>>> Dear Berk
>>>>>
>>>>> I am a bit concerned about this point. Let us restrict ourselves to
>>>>> the
>>>>> case of all H bond angles constrained. I was assuming that every H
>>>>> had
>>>>> one
>>>>> single (arbitrary) constrained bond angle. This is, in the scheme
>>>>> below,
>>>>> e.g. the angle C1-C2-H would be constrained, but not the H-C2-C3
>>>>> angle.
>>>>> Isn't this true?
>>>>>
>>>>> H
>>>>> |
>>>>> C1--C2--C3
>>>>>
>>>>> In some literature I read that the vibration period of H bond angles
>>>>> is
>>>>> the same as that of heavy atoms bond lengths (the Gromacs option
>>>>> h-angles
>>>>> for constraints seems to assume this). Isn't this Gromacs option
>>>>> proceeding as stated above?
>>>>
>>>> Contrary to your expectation, all bond angles involving the H atom are
>>>> actually constrained. There is more on that in the paper cited below.
>>>>
>>>>> Can you give me some reference where the procedure of 'virtual sites'
>>>>> to
>>>>> constrain H bond angles is explained?
>>>>
>>>> Hello Pablo,
>>>>
>>>> Berk and I wrote a paper on that a while ago:
>>>> J. Comput. Chem. 20 (8), 786-798
>>>>
>>>>
>>>>> Thank you very much. Best regards
>>>>>
>>>>>
>>>>>> On 6/10/13 12:15 , David van der Spoel wrote:
>>>>>>> On 2013-06-10 10:57, "Pablo García Risueño" wrote:
>>>>>>>> Thank you very much for the answer. It is said that:
>>>>>>>>
>>>>>>>> "A practical issue is that not many people use shake, since LINCS
>>>>>>>> works in
>>>>>>>> parallel too."
>>>>>>>>
>>>>>>>> But I thought that SHAKE was used when bond angles are
>>>>>>>> constrained,
>>>>>>>> this
>>>>>>>> is when bond lengths of heavy atoms are constrained (since bond
>>>>>>>> lengths of
>>>>>>>> heavy atoms and bond angles of hydrogen have the same period).
>>>>>>>> Isn't
>>>>>>>> this
>>>>>>>> correct?
>>>>>>>>
>>>>>>> That's up to the user, but LINCS can be used too for angles when
>>>>>>> the
>>>>>>> number of iterations is increased a bit (although not too many
>>>>>>> people
>>>>>>> use angle constraints either). Nevertheless gromacs supports the
>>>>>>> shake
>>>>>>> algorithm and hence it should work.
>>>>>> But note that constraining all-angles will lead to incompatible
>>>>>> constraints, unless the molecule is very small.
>>>>>> Constraining all angles involving hydrogens has similar issues. We
>>>>>> use
>>>>>> virtual sites to remove H-angle vibrations.
>>>>>> So in practice the only molecule for which angle constraints are
>>>>>> used
>>>>>> is
>>>>>> water, but there SETTLE is algorithm of choice.
>>>>>>
>>>>>> Cheers,
>>>>>>
>>>>>> Berk
>>>>>>>> Thank you very much
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>> On 2013-06-07 19:26, "Pablo García Risueño" wrote:
>>>>>>>>>> Dear Gromacs developers
>>>>>>>>>>
>>>>>>>>>> I am studying the Gromacs implementation of the SHAKE algorithm.
>>>>>>>>>> As
>>>>>>>>>> far
>>>>>>>>>> as
>>>>>>>>>> I know, it corresponds to the cshake procedure, which lies in
>>>>>>>>>> the
>>>>>>>>>> shakef.c
>>>>>>>>>> file, and the equation (5.6) of the SHAKE paper
>>>>>>>>>> (see
>>>>>>>>>>
http://www.sciencedirect.com/science/article/pii/0021999177900985)
>>>>>>>>>> is
>>>>>>>>>> essentially in
>>>>>>>>>>
>>>>>>>>>> acor = omega*diff*m2[ll]/rrpr
>>>>>>>>>> lagr[ll] += acor;
>>>>>>>>>>
>>>>>>>>>> where the relationship between variables is:
>>>>>>>>>>
>>>>>>>>>> SHAKE var Gromacs var
>>>>>>>>>> omega=1
>>>>>>>>>> d^2 - (r')^2 diff
>>>>>>>>>> 1/(2(1/m_i+1/mj)) m2[ll]
>>>>>>>>>> \vec{r}_{ij}\cdot\vec{r}'_{ij} rrpr
>>>>>>>>>> g lagr[ll]
>>>>>>>>>>
>>>>>>>>>> being ll the constraint index. However, I think that the term
>>>>>>>>>> proportional
>>>>>>>>>> to g^2 in the equation (5.6) of SHAKE paper is not in the
>>>>>>>>>> Gromacs
>>>>>>>>>> code.
>>>>>>>>>> This term is relatively small, but I was assuming that it is not
>>>>>>>>>> negligible. Could some developer explain me whether this term is
>>>>>>>>>> neglected, or it is included in the calculations elsewhere?
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>> I don't have the complete overview, so I will assume your
>>>>>>>>> observation is
>>>>>>>>> correct.
>>>>>>>>> Since the algorithm is iterative the second term may not be
>>>>>>>>> needed,
>>>>>>>>> although it might be converge in fewer iterations with it.
>>>>>>>>>
>>>>>>>>> A practical issue is that not many people use shake, since LINCS
>>>>>>>>> works
>>>>>>>>> in parallel too.
>>>>>>>>>
>>>>>>>>>> Thank you very much. Best regards
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> --
>>>>>>>>>>
>>>>>>>>>> Dr. Pablo García Risueño
>>>>>>>>>>
>>>>>>>>>> Institut für Physik und IRIS Adlershof, Humboldt Universität zu
>>>>>>>>>> Berlin,
>>>>>>>>>> Zum Grossen Windkanal 6, 12489 Berlin, Germany
>>>>>>>>>>
>>>>>>>>>> Tel. +49 030 209366369
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> --
>>>>>>>>> David van der Spoel, Ph.D., Professor of Biology
>>>>>>>>> Dept. of Cell & Molec. Biol., Uppsala University.
>>>>>>>>> Box 596, 75124 Uppsala, Sweden. Phone: +46184714205.
>>>>>>>>> spoel at xray.bmc.uu.se http://folding.bmc.uu.se
>>>>>>>>> --
>>>>>>>>> gmx-developers mailing list
>>>>>>>>> gmx-developers at gromacs.org
>>>>>>>>> http://lists.gromacs.org/mailman/listinfo/gmx-developers
>>>>>>>>> Please don't post (un)subscribe requests to the list. Use the
>>>>>>>>> www interface or send it to gmx-developers-request at gromacs.org.
>>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> --
>>>>>>>>
>>>>>>>> Dr. Pablo García Risueño
>>>>>>>>
>>>>>>>> Institut für Physik und IRIS Adlershof, Humboldt Universität zu
>>>>>>>> Berlin,
>>>>>>>> Zum Grossen Windkanal 6, 12489 Berlin, Germany
>>>>>>>>
>>>>>>>> Tel. +49 030 209366369
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>> --
>>>>>> gmx-developers mailing list
>>>>>> gmx-developers at gromacs.org
>>>>>> http://lists.gromacs.org/mailman/listinfo/gmx-developers
>>>>>> Please don't post (un)subscribe requests to the list. Use the
>>>>>> www interface or send it to gmx-developers-request at gromacs.org.
>>>>>>
>>>>>
>>>>>
>>>>> --
>>>>>
>>>>> Dr. Pablo García Risueño
>>>>>
>>>>> Institut für Physik und IRIS Adlershof, Humboldt Universität zu
>>>>> Berlin,
>>>>> Zum Grossen Windkanal 6, 12489 Berlin, Germany
>>>>>
>>>>> Tel. +49 030 209366369
>>>>>
>>>>
>>>>
>>>> --
>>>> Groetjes,
>>>>
>>>> Anton
>>>> _____________ >>>> | |
>>>> | _ _ ___,| K. Anton Feenstra
>>>> | / \ / \'| | | IBIVU/Bioinformatics - Vrije Universiteit Amsterdam
>>>> |( | )| | | De Boelelaan 1081 - 1081 HV Amsterdam - Netherlands
>>>> | \_/ \_/ | | | Tel +31 20 59 87783 - Fax +31 20 59 87653 - Room
>>>> | | Feenstra at few.vu.nl - www.few.vu.nl/~feenstra/
>>>> | | "Is This the Right Room for an Argument ?" (Monty
>>>> | | Python)
>>>> |_____________|________________________________________________
>>>> --
>>>> gmx-developers mailing list
>>>> gmx-developers at gromacs.org
>>>> http://lists.gromacs.org/mailman/listinfo/gmx-developers
>>>> Please don't post (un)subscribe requests to the list. Use the
>>>> www interface or send it to gmx-developers-request at gromacs.org.
>>>>
>>>
>>>
>>> --
>>>
>>> Dr. Pablo García Risueño
>>>
>>> Institut für Physik und IRIS Adlershof, Humboldt Universität zu Berlin,
>>> Zum Grossen Windkanal 6, 12489 Berlin, Germany
>>>
>>> Tel. +49 030 209366369
>>>
>>
>
>
> --
>
> Dr. Pablo García Risueño
>
> Institut für Physik und IRIS Adlershof, Humboldt Universität zu Berlin,
> Zum Grossen Windkanal 6, 12489 Berlin, Germany
>
> Tel. +49 030 209366369
>
--
Groetjes,
Anton
_____________ _______________________________________________________
| | |
| _ _ ___,| K. Anton Feenstra |
| / \ / \'| | | IBIVU/Bioinformatics - Vrije Universiteit Amsterdam |
|( | )| | | De Boelelaan 1081 - 1081 HV Amsterdam - Netherlands |
| \_/ \_/ | | | Tel +31 20 59 87783 - Fax +31 20 59 87653 - Room P136 |
| | Feenstra at few.vu.nl - www.few.vu.nl/~feenstra/ |
| | "Does All This Money Really Have To Go To Charity ?" |
| | (Rick) |
|_____________|_______________________________________________________|
More information about the gromacs.org_gmx-developers
mailing list