added trajectory wrap method

[jamesmkrieger]

Previously there was only a function to wrapAtoms for individual frames. This extends that to the whole trajectory.

[jamesmkrieger]

This should probably have more testing

[AnthonyBogetti]

I generated a trajectory of Na and Cl in a small water box where wrapping occurs. I visualized it and have quantitative support. I'll test the function on this.

[AnthonyBogetti]

@jamesmkrieger, I have a dcd file with wrapping and am trying to run the following code with it:

from prody import *
f = DCDFile("nacl_namd.dcd")
f.hasUnitcell
f.wrap()

The f.hasUnitcell is True, but when I try to wrap I get the following error:

In [5]: f.wrap()
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
File ~/apps/miniforge3/envs/prody-wrap_traj/lib/python3.8/site-packages/ProDy-2.0.2-py3.8-linux-x86_64.egg/prody/measure/transform.py:530, in wrapAtoms(frame, unitcell, center)
    529 try:
--> 530     unitcell = frame.getUnitcell()[:3]
    531 except AttributeError:

AttributeError: 'numpy.ndarray' object has no attribute 'getUnitcell'

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
Cell In[5], line 1
----> 1 f.wrap()

File ~/apps/miniforge3/envs/prody-wrap_traj/lib/python3.8/site-packages/ProDy-2.0.2-py3.8-linux-x86_64.egg/prody/trajectory/trajbase.py:397, in TrajBase.wrap(self, unitcell, center)
    395 LOGGER.progress('Wrapping trajectory ...', len(coordsets))
    396 for i, coordset in enumerate(coordsets):
--> 397     wrapped.addCoordset(wrapAtoms(coordset), unitcell=unitcell, center=center)
    398     LOGGER.update(i)
    400 LOGGER.finish()

File ~/apps/miniforge3/envs/prody-wrap_traj/lib/python3.8/site-packages/ProDy-2.0.2-py3.8-linux-x86_64.egg/prody/measure/transform.py:532, in wrapAtoms(frame, unitcell, center)
    530         unitcell = frame.getUnitcell()[:3]
    531     except AttributeError:
--> 532         raise TypeError('unitcell information must be provided')
    534 half = unitcell / 2
    535 ucmin = center - half

TypeError: unitcell information must be provided

Am I doing something wrong?

[jamesmkrieger]

I don’t really remember these functions so I don’t know. Did you actually set a unit cell or anything like that?

[AnthonyBogetti]

The frames do have a unitcell:

In [17]: f.getFrame(1).getUnitcell()                                                                                                                                         
Out[17]:                                                                                                                                                                     
array([29.05412292, 29.05412292, 29.05412292, 93.09894122, 93.09894122,                                                                                                      
       93.09894122])

[jamesmkrieger]

I wonder where we are ending up with a numpy array. Maybe DCDFile isn’t the right object?

[AnthonyBogetti]

I'll look into this a bit more.

[jamesmkrieger]

Thanks

[jamesmkrieger]

@AnthonyBogetti did you have a chance to look at this further?

Otherwise, can you send me a minimal trajectory please? You could sample it every 20 frames, so I have 5 frames to play with.

[AnthonyBogetti]

Hi @jamesmkrieger, can you merge all recent changes from main into this branch? I will try to get this resolved this week.

[AnthonyBogetti]

Here is a trajectory with some wrapped ions:
ions_unwrapped.zip

[jamesmkrieger]

Hi @jamesmkrieger, can you merge all recent changes from main into this branch? I will try to get this resolved this week.

I just pulled and it merged cleanly. Hopefully all the checks pass now

[jamesmkrieger]

Here is a trajectory with some wrapped ions: ions_unwrapped.zip

great. Would you like me to test something or is it good to go?

[karolamik13]

I can take a look at this later today if you want, James. I have many simulations which requires wrapping.

[jamesmkrieger]

I can take a look at this later today if you want, James. I have many simulations which requires wrapping.

That would be great. Thank you

[karolamik13]

I looked at this and I am unable to wrap my trajectory (see below). I don't know how to provide data with unitcell to wrap all frames at once.

Usually to wrap trajectory programs require xst file with size of the box. Here, we need to provide numpy.ndarray. I tried with collecting the size of the box for each frame and providing them like this dcd.wrap(unitcell=unitCell_data), but that won't work either. Any suggestions?

In [1]: from prody import *

In [2]: dcd = DCDFile('tests_wrap.dcd')

In [3]: dcd.hasUnitcell
Out[3]: <bound method DCDFile.hasUnitcell of <DCDFile: tests_wrap (next 0 of 51 frames; 71107 atoms)>>

In [4]: dcd.wrap()

---

AttributeError Traceback (most recent call last)
File ~/anaconda3/envs/prody_tests/lib/python3.10/site-packages/ProDy-2.5.0-py3.10-linux-x86_64.egg/prody/measure/transform.py:530, in wrapAtoms(frame, unitcell, center)
529 try:
--> 530 unitcell = frame.getUnitcell()[:3]
531 except AttributeError:

AttributeError: 'numpy.ndarray' object has no attribute 'getUnitcell'

During handling of the above exception, another exception occurred:

TypeError Traceback (most recent call last)
Cell In[4], line 1
----> 1 dcd.wrap()

File ~/anaconda3/envs/prody_tests/lib/python3.10/site-packages/ProDy-2.5.0-py3.10-linux-x86_64.egg/prody/trajectory/trajbase.py:397, in TrajBase.wrap(self, unitcell, center)
395 LOGGER.progress('Wrapping trajectory ...', len(coordsets))
396 for i, coordset in enumerate(coordsets):
--> 397 wrapped.addCoordset(wrapAtoms(coordset), unitcell=unitcell, center=center)
398 LOGGER.update(i)
400 LOGGER.finish()

File ~/anaconda3/envs/prody_tests/lib/python3.10/site-packages/ProDy-2.5.0-py3.10-linux-x86_64.egg/prody/measure/transform.py:532, in wrapAtoms(frame, unitcell, center)
530 unitcell = frame.getUnitcell()[:3]
531 except AttributeError:
--> 532 raise TypeError('unitcell information must be provided')
534 half = unitcell / 2
535 ucmin = center - half

TypeError: unitcell information must be provided

[jamesmkrieger]

ok, it looks like you have a similar problem to Anthony so probably there's something that needs fixing here

[jamesmkrieger]

ok, I had a bracket in the wrong place. Now it should work:

In [1]: from prody import *
   ...: f = DCDFile("ions_unwrapped.dcd")

In [2]: g = f.wrap()
                                    
In [3]: f
Out[3]: <DCDFile: ions_unwrapped (next 1 of 10150 frames; 2 atoms)>

In [4]: g
Out[4]: <Ensemble: ions_unwrapped_wrapped (10150 conformations; 2 atoms)>

[jamesmkrieger]

I generated a trajectory of Na and Cl in a small water box where wrapping occurs. I visualized it and have quantitative support. I'll test the function on this.

how did you generate this? Is there something equivalent we can do to compare if the wrapping worked?

[AnthonyBogetti]

The trajectory I uploaded was actually not generated by me (I forget who it came from) but I noticed the wrapping of the ions and thought it would be a good example. To test, calculating the inter-ionic distance and not seeing any jumps would indicate success I think.

[AnthonyBogetti]

Hi @jamesmkrieger, I tested it on my trajectory. It does work (wraps) but I see what you are saying it's hard to know if it is working correctly bc the trajectory I provided is not very simple.

[jamesmkrieger]

Maybe you can use an RMSD?

[AnthonyBogetti]

This is the distance between the two ions:

It's not clear to me if it is working or not... but I don't have any reason to think it is not working.

[AnthonyBogetti]

From visualization it seems to be working.

[jamesmkrieger]

I don't think a system with 2 ions is a very good one to look at. It's not easy to tell what we're seeing at all.

[karolamik13]

I am testing on my system. I will let you know how it goes. In the system I have two protein, ligand and oxygen molecules. Ligand and oxygens should be wrapped.

[karolamik13]

I am wondering..how long should it take? I have ~50 frames, 42 MB file and I am waiting and waiting (already >10 mins).

Also, usually I am using parseDCD instead of DCDFile, so maybe I am doing something wrong, but shouldn't this work differently:
In [3]: from prody import *
In [4]: dcd = DCDFile('tests_wrap.dcd', start=1, stop=5)
In [5]: dcd
Out[5]: <DCDFile: tests_wrap (next 0 of 51 frames; 71107 atoms)>
In [6]: dcd = DCDFile('tests_wrap.dcd', step=10)
In [7]: dcd
Out[7]: <DCDFile: tests_wrap (next 0 of 51 frames; 71107 atoms)>

I wanted to select less number of frames, but I cannot see the difference.

[jamesmkrieger]

I am wondering..how long should it take? I have ~50 frames, 42 MB file and I am waiting and waiting (already >10 mins).

Also, usually I am using parseDCD instead of DCDFile, so maybe I am doing something wrong, but shouldn't this work differently: In [3]: from prody import * In [4]: dcd = DCDFile('tests_wrap.dcd', start=1, stop=5) In [5]: dcd Out[5]: <DCDFile: tests_wrap (next 0 of 51 frames; 71107 atoms)> In [6]: dcd = DCDFile('tests_wrap.dcd', step=10) In [7]: dcd Out[7]: <DCDFile: tests_wrap (next 0 of 51 frames; 71107 atoms)>

I wanted to select less number of frames, but I cannot see the difference.

That's not good. I don't know how long it should take but that's surprising. Have you got lots of atoms?

[jamesmkrieger]

oh, 71k doesn't seem like so many

[AnthonyBogetti]

Hi @karolamik13, were you able to confirm that this functionality works as expected? Is it just slow? @jamesmkrieger do you think we should merge this in for 2.6.0 and improve it later on or should we wait for a future release? I think this may be the last open PR (other than some of the WatFinder improvements and new tools) for the 2.6.0 release.

[karolamik13]

Hi @AnthonyBogetti, I cannot confirm it right now. I have a problem with my Linux machine, where I was testing everything. I will be able to fix it no earlier than Thursday.