[gmx-developers] Virtual site out of the plane

[Berk Hess]

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.

Cheers,

Berk
>
>>> 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.
>>> Best,
>>> Mohammad
>>
>> -- 
>> 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.
>> http://www.icm.uu.se
>> -- 
>> 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.