Re: [AMBER] Unwrapping in MSD calculations with pytraj/cpptraj

From: Daniel Roe via AMBER <amber.ambermd.org>
Date: Sun, 31 Jul 2022 10:54:14 -0400

Hi,

I've fixed the behavior (I think) as of GitHub version 6.13.0
(https://github.com/Amber-MD/cpptraj/pull/984).

Let me know if you get a chance to try it out.

-Dan

On Fri, Jul 29, 2022 at 7:57 PM Cianna Calia <ccalia.ucsd.edu> wrote:
>
> Interesting, thanks. Right, the diffusion constants are pretty close; I had just been looking at these numbers since it's specifically the MSD that I needed for something.
>
> Thanks so much for the explanation; that helps a lot. Enjoy your vacation.
>
> Cianna
>
>
>
> On Fri, Jul 29, 2022 at 4:20 PM Daniel Roe via AMBER <amber.ambermd.org> wrote:
>>
>> Yeah, I'm pretty sure I know what's going on here. The TL;DR
>> explanation is that the current way I correct for imaging in diffusion
>> will gradually diverge from the way 'unwrap' corrects for imaging.
>> Neither is wrong per se, they just have different strategies to
>> correct for imaging. Arguably the 'unwrap' way is more "correct"
>> because it's not using a summation under the hood. You showed the raw
>> numbers but I'm guessing the diffusion constants calculated from those
>> files are probably really close right?
>>
>> Anyway, the fix (to make 'diffusion' correct imaging the way 'unwrap'
>> does) is simple but will require some testing on my part, and I'm
>> about to head out on vacation. Once I'm back I'll get the fix in, but
>> until then your best bet is probably to "unwrap" first and then
>> calculate diffusion.
>>
>> Thanks for the report.
>>
>> -Dan
>>
>> On Fri, Jul 29, 2022 at 6:56 PM Cianna Calia <ccalia.ucsd.edu> wrote:
>> >
>> > I'm using an orthorhombic box (roughly cubic since I'd used "iso" with solvatebox in tleap).
>> > cpptraj with AmberTools20
>> > pytraj version 2.0.5
>> >
>> > When I compare cpptraj's MSD data before vs. after unwrapping with byatom, the differences are ones to tens of Å^2. The maximum differences in x, y, z, and r datasets are:
>> >
>> > max X err: 7.819200000000137
>> > max Y err: 13.486300000000028
>> > max Z err: 18.626699999999914
>> > max R err: 24.975800000000163
>> >
>> > It's off by more than tenths or hundredths.
>> >
>> > Tailing the two cpptraj data files shows the difference - End of the datafile from before unwrapping:
>> >
>> > 1995 2774.8698 1317.9164 699.3513 757.6022 52.6770
>> > 1996 2590.9952 1182.8005 617.1339 791.0609 50.9018
>> > 1997 2672.6255 1248.8877 684.0725 739.6653 51.6974
>> > 1998 2666.3472 1224.3471 674.2470 767.7531 51.6367
>> > 1999 2700.3229 1240.0459 721.4465 738.8306 51.9646
>> >
>> > End of the datafile from after unwrapping:
>> >
>> > 1995 2785.2869 1312.0513 705.1366 768.0990 52.7758
>> > 1996 2602.1418 1177.2445 622.5693 802.3280 51.0112
>> > 1997 2684.3411 1243.1784 689.7944 751.3682 51.8106
>> > 1998 2677.9833 1218.6942 679.9278 779.3613 51.7492
>> > 1999 2714.8978 1234.3569 727.3224 753.2185 52.1047
>> >
>> > Thanks for giving this your attention,
>> > Cianna
>> >
>> >
>> > On Fri, Jul 29, 2022 at 12:51 PM Daniel Roe via AMBER <amber.ambermd.org> wrote:
>> >>
>> >> Just one more question - just how different are the results? Are they
>> >> pretty close (within a few tenths/hundreths of an Ang.) or are they
>> >> way off?
>> >>
>> >> -Dan
>> >>
>> >> On Fri, Jul 29, 2022 at 3:49 PM Daniel Roe <daniel.r.roe.gmail.com> wrote:
>> >> >
>> >> > Hi,
>> >> >
>> >> > I think I see the behavior you are seeing. I need to look into it more
>> >> > and I'll update you with what I find.
>> >> >
>> >> > What kind of unit cell are you using? Also, what versions of cpptraj
>> >> > and pytraj are you using?
>> >> >
>> >> > Thanks!
>> >> >
>> >> > -Dan
>> >> >
>> >> > On Mon, Jul 25, 2022 at 8:58 AM Cianna Calia via AMBER
>> >> > <amber.ambermd.org> wrote:
>> >> > >
>> >> > > Hello AMBER community!
>> >> > >
>> >> > > I am wondering if anyone can help me understand how the diffusion functions
>> >> > > in pytraj/cpptraj handle unwrapping. When I try using pytraj to calculate
>> >> > > an atom's MSD from a trajectory unwrapped in cpptraj (with "unwrap
>> >> > > byatom"), I am unsure why it deviates from pytraj's result for the original
>> >> > > wrapped trajectory, since I thought pytraj would deal with unwrapping
>> >> > > automatically in a way that would give the same result as using the already
>> >> > > unwrapped trajectory. But the X, Y, and Z datasets differ starting where
>> >> > > unwrapping becomes relevant for each dimension.
>> >> > >
>> >> > > More detail - I have the original wrapped trajectory
>> >> > > 1hg7_WT_cubic_100ns_1_last2ns.nc and an unwrapped trajectory
>> >> > > unwrap_2ns_test1.nc made in cpptraj like this:
>> >> > >
>> >> > > parm 1hg7_WT_cubic.prmtop
>> >> > > trajin 1hg7_WT_cubic_100ns_1_last2ns.nc
>> >> > > unwrap byatom
>> >> > > trajout unwrap_2ns_test1.nc
>> >> > >
>> >> > > I tried the following diffusion calculations with pytraj:
>> >> > >
>> >> > > import pytraj as pt
>> >> > >
>> >> > > original_traj = pt.load('1hg7_WT_cubic_100ns_1_last2ns.nc',
>> >> > > '1hg7_WT_cubic.prmtop')
>> >> > > pt_diffusion_orig_traj = pt.diffusion(original_traj, '.1730', dtype='dict',
>> >> > > frame_indices=[i for i in range(len(original_traj))])
>> >> > >
>> >> > > unwrapped_traj = pt.load('unwrap_2ns_test1.nc', '1hg7_WT_cubic.prmtop')
>> >> > > pt_diffusion_unw_traj = pt.diffusion(unwrapped_traj, '.1730', dtype='dict',
>> >> > > frame_indices=[i for i in range(len(unwrapped_traj))])
>> >> > >
>> >> > > I am hoping for clarification on why pt_diffusion_orig_traj and
>> >> > > pt_diffusion_unw_traj are not the same.
>> >> > >
>> >> > > I observe the same situation when I do analogous calculations with cpptraj
>> >> > > as opposed to pytraj. I also tried seeing how the noimage option with the
>> >> > > cpptraj diffusion function affects the results, and I am similarly confused
>> >> > > about why the two diffusion calculations below give different results
>> >> > > (again, differing starting from where there's unwrapping):
>> >> > >
>> >> > > parm 1hg7_WT_cubic.prmtop
>> >> > > trajin 1hg7_WT_cubic_100ns_1_last2ns.nc
>> >> > > diffusion out orig_traj_1730_diff_cpptraj.dat time 1.0 "@1730" nocalc
>> >> > > unwrap byatom
>> >> > > diffusion out unwrapped_traj_1730_diff_cpptraj.dat time 1.0 noimage "@1730"
>> >> > > nocalc
>> >> > >
>> >> > > Taking out the "byatom" option makes no noticeable difference, nor does
>> >> > > removing the "noimage" option from the diffusion command acting on the
>> >> > > already unwrapped trajectory. The difference between the diffusion data
>> >> > > before and after unwrapping here is not the same as the difference I
>> >> > > observe when I use the diffusion command on the original wrapped trajectory
>> >> > > with vs. without the "noimage" option.
>> >> > >
>> >> > > Thanks in advance for any help with understanding this!
>> >> > > Cianna Calia
>> >> > > _______________________________________________
>> >> > > AMBER mailing list
>> >> > > AMBER.ambermd.org
>> >> > > https://urldefense.com/v3/__http://lists.ambermd.org/mailman/listinfo/amber__;!!Mih3wA!A3PyHn9pZqfoG08rMdV6_T2d61zGEtAPiOSJA7y1FkFtqo18BxS-4b4O4gvuDJdUuUQZn-r_1fI32A$
>> >>
>> >> _______________________________________________
>> >> AMBER mailing list
>> >> AMBER.ambermd.org
>> >> https://urldefense.com/v3/__http://lists.ambermd.org/mailman/listinfo/amber__;!!Mih3wA!A3PyHn9pZqfoG08rMdV6_T2d61zGEtAPiOSJA7y1FkFtqo18BxS-4b4O4gvuDJdUuUQZn-r_1fI32A$
>>
>> _______________________________________________
>> AMBER mailing list
>> AMBER.ambermd.org
>> https://urldefense.com/v3/__http://lists.ambermd.org/mailman/listinfo/amber__;!!Mih3wA!CpPOy-HcJwxbvZgAUfP1-CShCHpDMRwrECMAb5eqw5v2kZoT6_22TXWtiS6fgXbCq_LXhlyI8E-cYw$

_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Aug 04 2022 - 13:37:48 PDT
Custom Search