search

Amber-MD / pytraj

Python interface of cpptraj
https://amber-md.github.io/pytraj
170 stars 38 forks source link

pt.compute 'distance' calcuates nan or zeros with pt.Trajectory #1517

Closed pindakaas42 closed 4 years ago

pindakaas42 commented 4 years ago

If I use pt.compute to calculate the distance of a trajectory loaded with pt.Trajectory pt.compute('distance :1 :2', traj=tr) I get this as a result: OrderedDict([(u'Dis_00000', array([0., 0., 0., ..., 0., 0., 0.]))]) If I use the pt.distance function pt.distance(mask=':1 :2', traj=tr) It returns: array([4.49308083, 4.86024722, 4.24035067, ..., 5.33574095, 4.32270206, 5.9003888 ]) This does not happen when using,TrajectoryIterator in both combinations.

Version: pytraj-2.0.5-py2.7-linux-x86_64.egg

hainm commented 4 years ago

thanks @pindakaas42 for keeping report the issues and I really appreciate.

Would you mind posting toy files so I can try to reproduce (pytraj has very extensive testing for those kind of stuff but it might miss something like you are pointing out).

hainm commented 4 years ago 
In [3]: tr = pt.datafiles.load_tz2_ortho()                                                                                                                                                                         

In [4]: pt.compute('distance :1 :2', traj=tr)                                                                                                                                                                      
Out[4]: 
OrderedDict([('Dis_00000',
              array([4.38414694, 4.57478295, 4.55499558, 4.70996977, 4.6515257 ,
                     4.75688718, 4.81763691, 4.72835557, 4.82638619, 5.30096546]))])

In [5]: pt.distance(mask=':1 :2', traj=tr)                                                                                                                                                                         
Out[5]: 
array([4.38414694, 4.57478295, 4.55499558, 4.70996977, 4.6515257 ,
       4.75688718, 4.81763691, 4.72835557, 4.82638619, 5.30096546])
pindakaas42 commented 4 years ago

Hey, thanks for the responses. Your load_tz2_ortho() uses TrajectoryIterator, so the bug would not show up anyway, but even if I use tr = Trajectory( pt.datafiles.load_tz2_ortho()) I can't recreate the problem. This works though to recreate the problem with the testdata:

top = pt.load_topology('AMBERHOME/lib/python2.7/site-packages/pytraj-2.0.5-py2.7-linux-x86_64.egg/pytraj/datafiles/tz2/tz2.ortho.parm7')
tr = pt.Trajectory('AMBERHOME/lib/python2.7/site-packages/pytraj-2.0.5-py2.7-linux-x86_64.egg/pytraj/datafiles/tz2/tz2.ortho.nc', top)
pt.compute('distance :1 :2', tr)

outputs: 

OrderedDict([(u'Dis_00000',
              array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]))])

while pt.distance(mask=':1 :2', traj=tr) returns 

array([4.38414694, 4.57478295, 4.55499558, 4.70996977, 4.6515257 ,
       4.75688718, 4.81763691, 4.72835557, 4.82638619, 5.30096546])

If it helps I can also give you some sample files but I hope this also helps. The way the files are loaded seems to matter.

hainm commented 4 years ago

thanks. I can reproduce your issue if creating Trajectory directly like you pointed out. It is ok if the traj comes from pytraj.load

In [1]: import pytraj as pt                                                                                                                                                                                        

In [2]: top = pt.load_topology('./tz2.ortho.parm7')                                                                                                                                                                

In [3]: tr = pt.Trajectory('tz2.ortho.nc', top)                                                                                                                                                                    

In [4]: pt.compute('distance :1 :2', traj=tr)                                                                                                                                                                      
Out[4]: 
OrderedDict([('Dis_00000',
              array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]))])

In [5]: tr2 = pt.load('tz2.ortho.nc', 'tz2.ortho.parm7')                                                                                                                                                           

In [6]: pt.compute('distance :1 :2', traj=tr2)                                                                                                                                                                     
Out[6]: 
OrderedDict([('Dis_00000',
              array([4.38414694, 4.57478295, 4.55499558, 4.70996977, 4.6515257 ,
                     4.75688718, 4.81763691, 4.72835557, 4.82638619, 5.30096546]))])
hainm commented 4 years ago

The difference here is that the unitcells is missing in the "nan" case.

In [11]: print(tr.unitcells)                                                                                                                                                                                       
None

In [12]: print(tr2.unitcells)                                                                                                                                                                                      
[[35.26277966 41.84554768 36.16862953 90.         90.         90.        ]
 [35.2563193  41.83788131 36.16200321 90.         90.         90.        ]
 [35.25450718 41.83573091 36.16014454 90.         90.         90.        ]
 ...
 [35.24903289 41.82923469 36.15452962 90.         90.         90.        ]
 [35.26186109 41.84445763 36.16768736 90.         90.         90.        ]
 [35.27351853 41.85829125 36.17964426 90.         90.         90.        ]]
pindakaas42 commented 4 years ago

So that's a bug right?

hainm commented 4 years ago

Yes, that’s a bug.

hainm commented 4 years ago

fixed in https://github.com/Amber-MD/pytraj/pull/1551

In [4]: tr = pt.Trajectory('tz2.ortho.nc', top='tz2.ortho.parm7')                                                                                   

In [5]: pt.compute('distance :1 :2', traj=tr)                                                                                                       
Out[5]: 
OrderedDict([('Dis_00000',
              array([4.38414694, 4.57478295, 4.55499558, 4.70996977, 4.6515257 ,
                     4.75688718, 4.81763691, 4.72835557, 4.82638619, 5.30096546]))])

Thanks @pindakaas42.