mdanalysis
mdanalysis copied to clipboard
MDAnalysis
Nojump transform using the “hybrid unwrap” scheme
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
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
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 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?
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)]
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)]
@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. 😊😊
Naw, I figured it out. @hmacdope , its good to go now, assuming I wrote the test correctly. See #4011