Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Nojump transform using the “hybrid unwrap” scheme #3703

Closed
orbeckst opened this issue Jun 4, 2022 · 10 comments · Fixed by #4031
Closed

Nojump transform using the “hybrid unwrap” scheme #3703

orbeckst opened this issue Jun 4, 2022 · 10 comments · Fixed by #4031
Labels

Comments

@orbeckst
Copy link
Member

orbeckst commented Jun 4, 2022

Is your feature request related to a problem?

MDAnalysis does not have a transformation to create “unwrapped” no-jump trajectories where particles are never put back in the primary unit cell. This is a problem for MSD calculations #3179 and analysis/visualization of oligomeric assemblies.

Describe the solution you'd like

The preprint https://arxiv.org/abs/2111.12052 A Reversible Unwrapping Algorithm for Constant Pressure Molecular Dynamics Simulations by Martin Kulke, Josh V Vermaas shows a simple algorithm for orthorhombic boxes that accurately unwrap long NPT trajectories. It improves over work by von Bülow et al by also maintaining correct bond distances in all instances. It is implemented in various PBC packages for VMD.

The authors do not provide an algorithm for triclinic boxes in the paper but maybe @jvermaas has some ideas that didn’t make it into the paper…?

Describe alternatives you've considered

Do not do anything and have users use native tools such as gmx trajectory -nojump or VMD with pbctools/fastpbc/qwrap

@jvermaas
Copy link
Contributor

jvermaas commented Jun 4, 2022

Martin and I have been talking, and we think the generalization for triclinic boxes is something like this. Written out this way, it looks like it ought to be performant, but we haven't written tests or an implementation yet.
alternative.pdf

@hmacdope
Copy link
Member

hmacdope commented Jun 5, 2022

Would love to solve the lack of a "nojump" transform, as has been bothering me for a while, just havn't found the time to give it the proper attention.

@orbeckst
Copy link
Member Author

orbeckst commented Jun 6, 2022

Just for completeness: PR #2982 (stalled) made an attempt at nojump.

@jvermaas jvermaas mentioned this issue Feb 19, 2023
4 tasks
@jvermaas
Copy link
Contributor

@orbeckst and @hmacdope , I've got a draft that I've tested to work on the trajectories I had on hand with a orthogonal box. I'm pretty sure it would work for non-orthogonal boxes, but if I wanted to write a test that uses existing demo topologies and trajectories, are there any available?

@orbeckst
Copy link
Member Author

orbeckst commented Feb 19, 2023

The standard TPR/TRR trajectory is in a truncated dodecahedron or octahedron IIRC —

>>> import MDAnalysis as mda
>>> from MDAnalysis.tests import datafiles as data
>>> u = mda.Universe(data.TPR, data.TRR)
>>> [u.trajectory.ts.dimensions.copy() for ts in u.trajectory]
[array([80.017006, 80.017006, 80.017006, 60.      , 60.      , 90.      ],
       dtype=float32),
 array([80.13008, 80.13008, 80.13008, 60.     , 60.     , 90.     ],
       dtype=float32),
 array([79.98359, 79.98359, 79.98358, 60.     , 60.     , 90.     ],
       dtype=float32),
 array([80.00516, 80.00516, 80.00516, 60.     , 60.     , 90.     ],
       dtype=float32),
 array([80.13544, 80.13544, 80.13544, 60.     , 60.     , 90.     ],
       dtype=float32),
 array([80.07649 , 80.07649 , 80.076485, 59.999996, 59.999996, 90.      ],
       dtype=float32),
 array([79.93229, 79.93229, 79.93228, 60.     , 60.     , 90.     ],
       dtype=float32),
 array([79.95813, 79.95813, 79.95813, 60.     , 60.     , 90.     ],
       dtype=float32),
 array([80.11252, 80.11252, 80.11252, 60.     , 60.     , 90.     ],
       dtype=float32),
 array([80.08523, 80.08523, 80.08523, 60.     , 60.     , 90.     ],
       dtype=float32)]

@orbeckst
Copy link
Member Author

Other examples: truncated octahedron data.PRM7, data.NCDFtruncoct

>>> u = mda.Universe(data.PRM7, data.NCDFtruncoct)
>>> [u.trajectory.ts.dimensions.copy() for ts in u.trajectory]
[array([ 42.43885,  42.43885,  42.43885, 109.47122, 109.47122, 109.47122],
       dtype=float32),
 array([ 42.42722,  42.42722,  42.42722, 109.47122, 109.47122, 109.47122],
       dtype=float32),
 array([ 42.436703,  42.436703,  42.436703, 109.47122 , 109.47122 ,
        109.47122 ], dtype=float32),
 array([ 42.426884,  42.426884,  42.426884, 109.47122 , 109.47122 ,
        109.47122 ], dtype=float32),
 array([ 42.428925,  42.428925,  42.428925, 109.47122 , 109.47122 ,
        109.47122 ], dtype=float32),
 array([ 42.428345,  42.428345,  42.428345, 109.47122 , 109.47122 ,
        109.47122 ], dtype=float32),
 array([ 42.42935,  42.42935,  42.42935, 109.47122, 109.47122, 109.47122],
       dtype=float32),
 array([ 42.437664,  42.437664,  42.437664, 109.47122 , 109.47122 ,
        109.47122 ], dtype=float32),
 array([ 42.448086,  42.448086,  42.448086, 109.47122 , 109.47122 ,
        109.47122 ], dtype=float32),
 array([ 42.428432,  42.428432,  42.428432, 109.47122 , 109.47122 ,
        109.47122 ], dtype=float32)]

General triclinic data.PSF_TRICLINIC, data.DCD_TRICLINIC

>>> u = mda.Universe(data.PSF_TRICLINIC, data.DCD_TRICLINIC)
>>> [u.trajectory.ts.dimensions.copy() for ts in u.trajectory]
[array([35.446037, 35.06156 , 34.158504, 91.328026, 61.735207, 44.40703 ],
       dtype=float32),
 array([34.659573, 34.226887, 33.098965, 90.562065, 61.791916, 44.14549 ],
       dtype=float32),
 array([34.527714, 34.66422 , 33.538815, 90.55859 , 63.112278, 40.140434],
       dtype=float32),
 array([34.437492, 33.38432 , 34.021328, 88.82457 , 64.980576, 36.773968],
       dtype=float32),
 array([33.73129 , 32.477516, 34.189613, 89.88101 , 65.89031 , 36.10921 ],
       dtype=float32),
 array([33.787033, 31.903166, 34.98833 , 90.03092 , 66.12876 , 35.07141 ],
       dtype=float32),
 array([33.247074, 31.182705, 34.965397, 93.11122 , 68.17744 , 35.736427],
       dtype=float32),
 array([32.925987, 30.313934, 34.991974, 93.89051 , 69.3799  , 33.489452],
       dtype=float32),
 array([32.15295 , 30.430563, 34.961567, 96.01416 , 71.50115 , 32.561108],
       dtype=float32),
 array([31.997482, 30.215181, 35.24292 , 95.858215, 71.08429 , 31.85939 ],
       dtype=float32)]

@hmacdope
Copy link
Member

@jvermaas nojump would be a very welcome addition! Just a note, as far as I am aware you do need seperate code paths for for different box types. However we can restrict to orthorhombic for now. 😊😊

@jvermaas
Copy link
Contributor

Naw, I figured it out. @hmacdope , its good to go now, assuming I wrote the test correctly. See #4011

@orbeckst
Copy link
Member Author

Did you mean PR #4031, @jvermaas ?

@jvermaas
Copy link
Contributor

Yes, I did. Thanks @orbeckst. I've gotten the citations reformatted.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants