[gmx-developers] Virtual site out of the plane
hess at kth.se
Fri Sep 29 11:33:30 CEST 2017
On 2017-09-29 00:38, Mohammad Ghahremanpour wrote:
>> On Sep 28, 2017, at 10:51 PM, David van der Spoel <spoel at xray.bmc.uu.se> wrote:
>> On 28/09/17 14:59, Mohammad Ghahremanpour wrote:
>>> Hello developers,
>>> The calc_vsite3out_param() function in vsite_parm.cpp calculates the weights (called a, b, c in the code) to locate the virtual site (i) out of the k-j-l plane later on in vsite.cpp.
>>> I copied some part of the code here:
>>> pijk = std::cos(aijk)*bij;
>>> pijl = std::cos(aijl)*bij;
>>> a = ( pijk + (pijk*std::cos(akjl)-pijl) * std::cos(akjl) / gmx::square(std::sin(akjl)) ) / bjk;
>>> b = ( pijl + (pijl*std::cos(akjl)-pijk) * std::cos(akjl) / gmx::square(std::sin(akjl)) ) / bjl;
>>> c = -std::sqrt( gmx::square(bij) -
>>> ( gmx::square(pijk) - 2*pijk*pijl*std::cos(akjl) + gmx::square(pijl) )
>>> / gmx::square(std::sin(akjl)) )
>>> / ( bjk*bjl*std::sin(akjl) );
>>> There are a lot of projections going on here which make the code confusing without any documentation, unfortunately.
>>> Is there some reference for the way that a, b, and c are calculated here?
>> Hate to suggest it but have you checked the manual?
> In the manual, it has been described how a, b, and c are used in 3out, not how they are calculated. Pages 103 and 104.
> Are you referring to the same pages?
David might have misunderstood your question (I did initially).
The code you pasted implements the inversion of the formulas in the
manual. We do not have those inverted formulas documented, except for
the code itself.
>>> BTW, what if the number under the SQRT is negative? Then c will be NAN. There is no checking to make sure that we are not calculating the SQRT of a negative number.
>> David van der Spoel, Ph.D., Professor of Biology
>> Head of Department, Cell & Molecular Biology, Uppsala University.
>> Box 596, SE-75124 Uppsala, Sweden. Phone: +46184714205.
>> Gromacs Developers mailing list
>> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List before posting!
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers or send a mail to gmx-developers-request at gromacs.org.
More information about the gromacs.org_gmx-developers