feat: Add support for reading unit cell dimensions from MOL2 CRYSIN records

package/CHANGELOG

@@ -16,7 +16,7 @@ The rules for this file:
 -------------------------------------------------------------------------------
 ??/??/?? IAlibay, orbeckst, marinegor, tylerjereddy, ljwoods2, marinegor,
          spyke7, talagayev, tanii1125, BradyAJohnston, hejamu, jeremyleung521,
-         harshitgajjela-droid, kunjsinha, aygarwal, jauy123
+         harshitgajjela-droid, kunjsinha, aygarwal, jauy123, Dreamstick9
 
  * 2.11.0
 
@@ -42,6 +42,7 @@ Fixes
    DSSP by porting upstream PyDSSP 0.9.1 fix (Issue #4913)
 
 Enhancements
+ * MOL2Parser now reads unit cell dimensions from @<TRIPOS>CRYSIN records (Issue #3341)
  * Reduces duplication of code in _apply() function (Issue #5247, PR #5294)
  * Added new top-level `MDAnalysis.fetch` module (PR #4943)
  * Added new function `MDAnalysis.fetch.from_PDB` to download structure files from wwPDB

package/MDAnalysis/coordinates/MOL2.py

@@ -69,8 +69,11 @@
 * The MDAnalysis :class:`MOL2Reader` and :class:`MOL2Writer` only handle the
   MOLECULE, SUBSTRUCTURE, ATOM, and BOND record types. Other records are not
   currently read or preserved on writing.
-* As the CRYSIN record type is not parsed / written, MOL2 systems always have
-  dimensions set to ``None`` and dimensionless MOL2 files are written.
+* The MDAnalysis :class:`MOL2Reader` now reads unit cell dimensions from the
+  @<TRIPOS>CRYSIN record. The space group and setting information (the 7th
+  and 8th fields) are currently ignored. The unit cell is interpreted in the
+  common crystallographic convention with box vector **a** parallel to the x-axis,
+  **b** in the xy-plane, and **c** having a positive z-component.
 
 
 MOL2 format notes
@@ -246,6 +249,10 @@ def _read_frame(self, frame):
                 self.ts.data[sect] = sections[sect]
             except KeyError:
                 pass
+        if "crysin" in sections:
+            line = sections["crysin"][0].strip()
+            dims = [float(x) for x in line.split()[:6]]
+            self.ts.dimensions = np.array(dims, dtype=np.float32)
 
         self.ts.positions = np.array(coords, dtype=np.float32)
 

testsuite/MDAnalysisTests/coordinates/test_mol2.py

@@ -21,6 +21,7 @@
 # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
 #
 import pytest
+import numpy as np
 
 import os
 from numpy.testing import (
@@ -29,6 +30,7 @@
     assert_array_almost_equal,
     TestCase,
     assert_almost_equal,
+    assert_allclose,
 )
 
 from MDAnalysisTests.datafiles import (
@@ -39,6 +41,7 @@
     mol2_comments_header,
     mol2_ligand,
     mol2_sodium_ion,
+    mol2_crysin,
 )
 from MDAnalysis import Universe
 import MDAnalysis as mda
@@ -232,3 +235,11 @@ def test_mol2_universe_write(tmpdir):
         assert_almost_equal(u.atoms.positions, u2.atoms.positions)
         # MDA does not current implement @<TRIPOS>CRYSIN reading
         assert u2.dimensions is None
+
+
+def test_mol2_crysin_dimensions():
+    # test that crysin records are read as dimensions
+    u = mda.Universe(mol2_crysin)
+
+    expected = np.array([40.0, 50.0, 60.0, 90.0, 90.0, 90.0], dtype=np.float32)
+    assert_allclose(u.dimensions, expected, atol=1e-3)

testsuite/MDAnalysisTests/data/mol2/test_crysin.mol2

@@ -0,0 +1,10 @@
+@<TRIPOS>MOLECULE
+test_structure
+ 1 0 0 0 0
+SMALL
+USER_CHARGES
+@<TRIPOS>ATOM
+      1 C          0.0000      0.0000      0.0000 C.3     1  ALA        0.0000
+@<TRIPOS>CRYSIN
+   40.0000   50.0000   60.0000   90.0000   90.0000   90.0000    1    1
+@<TRIPOS>END

testsuite/MDAnalysisTests/datafiles.py

@@ -245,6 +245,7 @@
     "mol2_comments_header",
     "mol2_ligand",
     "mol2_sodium_ion",
+    "mol2_crysin",
     "capping_input",
     "capping_output",
     "capping_ace",
@@ -745,6 +746,7 @@
 mol2_zinc = (_data_ref / "mol2/zinc_856218.mol2").as_posix()
 # MOL2 file without bonds
 mol2_sodium_ion = (_data_ref / "mol2/sodium_ion.mol2").as_posix()
+mol2_crysin = (_data_ref / "mol2/test_crysin.mol2").as_posix()
 
 capping_input = (_data_ref / "capping/aaqaa.gro").as_posix()
 capping_output = (_data_ref / "capping/maestro_aaqaa_capped.pdb").as_posix()