Re: [AMBER] Cpptraj - Stat bug?

From: Daniel Roe <daniel.r.roe.gmail.com>
Date: Thu, 28 Jul 2016 16:03:33 -0600

This bug is confirmed. The issue is that this line:

if (dval < 0) dval += 360;

should actually be:

if (dval < 30) dval += 360;

in order to ensure that values below 30 end up in bin 5 (c) instead of
bin 0 (g+). This bug appears to have been present since ptraj. I'll
start working on a fix. Thanks for the report!

-Dan


On Thu, Jul 28, 2016 at 3:05 PM, Daniel Roe <daniel.r.roe.gmail.com> wrote:
> There does seem to be a problem. I'll look into it further and get
> back to you. Thanks for the report.
>
> -Dan
>
> On Fri, Jul 22, 2016 at 6:22 AM, Markus Schneider
> <markus.schneider2.mytum.de> wrote:
>> Hi,
>>
>> during analysis of "multidihedral" generated torsion datasets with the
>> "stat" command, I encountered an unexpected result. Compared to my own
>> calculations, the "c" and "g+" bin occupancies are different. I expected
>> that the "c" bin would contain angles from -30 to +30 deg and "g+" would
>> contain the range from +30 to +90 deg.
>>
>> The raw multidihedral output data is in the range [-180,+180]. When I
>> bin the raw data myself, I can reproduce the result of cpptraj's stat
>> analysis exactly if I set the c bin to [-30,0] and the g+ bin to [0,90].
>> When I set the bins to the expected values above, I get more data points
>> in the c bin, while fewer are placed in the g+ bin. I would like to ask
>> whether there might be an error with how the stat command bins the data.
>>
>> The relevant code piece in Analysis_Statistics.cpp seems to be:
>>
>> double value = ds.Dval( i );
>> double dval = value;
>> if (dval < 0) dval += 360;
>> *curbin = (int) (dval - 30.0) / 60;*
>>
>> I suspect that the last line might be responsible for the observation.
>> The first three lines, as I infer, provide that dval is transformed to
>> the range [0,360]. I created a test in C for the last line:
>>
>> --------------
>>
>> int bin(double angle)
>> {
>> return (int) (angle - 30.0) / 60;
>> }
>>
>> int main()
>> {
>> double angles[] = {0.0,29.0,30.0,89.0,90.0,149.0,
>> 150.0,209.0,210.0,269.0,270.0,329.0,330.0,359.0};
>> for(int i=0;i<14;i++){
>> printf("angle: %.1f bin: %d\n",angles[i],bin(angles[i]));
>> }
>>
>> return 0;
>> }
>>
>> --------------
>>
>> which creates the following output:
>>
>> angle: 0.0 bin: 0
>> angle: 29.0 bin: 0
>> angle: 30.0 bin: 0
>> angle: 89.0 bin: 0
>> angle: 90.0 bin: 1
>> angle: 149.0 bin: 1
>> angle: 150.0 bin: 2
>> angle: 209.0 bin: 2
>> angle: 210.0 bin: 3
>> angle: 269.0 bin: 3
>> angle: 270.0 bin: 4
>> angle: 329.0 bin: 4
>> angle: 330.0 bin: 5
>> angle: 359.0 bin: 5
>>
>> According to bin definitions above, I expect that the angles 0 and 29
>> deg should be part of bin 5. Instead, it seems that stat's g+ bin (0)
>> ranges from 0 to 90, while the c bin (5) ranges from 330 to 0. This is
>> consistent with my manual calculations. What do you think about this?
>>
>> This was tested with a fresh build of AmberTools 16.
>>
>> Best,
>>
>> Markus
>>
>>
>>
>> _______________________________________________
>> AMBER mailing list
>> AMBER.ambermd.org
>> http://lists.ambermd.org/mailman/listinfo/amber
>
>
>
> --
> -------------------------
> Daniel R. Roe, PhD
> Department of Medicinal Chemistry
> University of Utah
> 30 South 2000 East, Room 307
> Salt Lake City, UT 84112-5820
> http://home.chpc.utah.edu/~cheatham/
> (801) 587-9652
> (801) 585-6208 (Fax)



-- 
-------------------------
Daniel R. Roe, PhD
Department of Medicinal Chemistry
University of Utah
30 South 2000 East, Room 307
Salt Lake City, UT 84112-5820
http://home.chpc.utah.edu/~cheatham/
(801) 587-9652
(801) 585-6208 (Fax)
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Jul 28 2016 - 15:30:02 PDT
Custom Search