Add pytests for polyconf

polyconf/polyconf/Monomer.py

@@ -1,4 +1,5 @@
 #!/usr/bin/env python 
+from __future__ import annotations
 import MDAnalysis as mda
 from MDAnalysis import Merge
 
@@ -27,10 +28,12 @@ def __init__(self, monomerName, fromUni = False) -> None:
             self.monomer = mda.Universe(monomerName)
         self.residues = self.monomer.residues
         self.atoms = self.monomer.atoms
+        self.dimensions = self.monomer.dimensions
 
     @classmethod
     def monomer_from_u(cls, universe:mda.Universe) -> 'Monomer':
         """
+        Create Monomer from an MDAnalysis Universe
 
         :param universe: _description_
         :type universe: mda.Universe
@@ -53,4 +56,21 @@ def select_atoms(self, selection:str) -> mda.AtomGroup:
         :return: An MDAnalysis AtomGroup containing only the selected atoms
         :rtype: mda.AtomGroup
         """
-        return self.monomer.select_atoms(selection)
+        return self.monomer.select_atoms(selection)
+
+    def copy(self) -> Monomer:
+        """
+        Use MDAnalysis Universe.copy() to create a new MDAnalysis Universe that
+        is an exact copy of this one containing the polymer. Deepcopy is not
+        used as it can be problematic with open file sockets.
+
+        :return: A new Polymer that is an exact copy of this polymer
+        :rtype: Polymer
+        """
+        new = self.monomer.copy()
+        return Monomer.monomer_from_u(new)
+
+    # def __eq__(self, value):
+    #     new = self.atoms.subtract(value.atoms)
+    #     print(len(new))
+    #     return self.select_atoms("all") == value.select_atoms("all") # and self.residues.equals(value.residues)

polyconf/polyconf/PDB.py

@@ -55,9 +55,9 @@ def save(self, dummyAtoms="X*", fname="polymer", selectionString = None, gmx = F
         """
         if selectionString:
             if gmx:
-                self.select_atoms(f"{selectionString} and not name {dummyAtoms}").atoms._write(f"{fname}.gro")
+                self.select_atoms(f"{selectionString} and not name {dummyAtoms}")._write(f"{fname}.gro")
             else:
-                self.select_atoms(f"{selectionString} and not name {dummyAtoms}").atoms._write(f"{fname}.pdb")
+                self.select_atoms(f"{selectionString} and not name {dummyAtoms}")._write(f"{fname}.pdb")
         else:
             if gmx:
                 self._write(f"not name {dummyAtoms}", f"{fname}.gro")
@@ -72,7 +72,7 @@ def crudesave(self,fname="polymer_crude"):
         :param fname: name of the output file, defaults to "polymer_crude"
         :type fname: str, optional
         """
-        self.atoms._write(f"{fname}.pdb")
+        self._write("all", f"{fname}.pdb")
 
     def select_atoms(self, selection) -> mda.AtomGroup:
         """

polyconf/polyconf/Polymer.py

@@ -45,6 +45,7 @@ def __init__(self, firstMonomer: Monomer, keep_resids=False) -> None:
         self.first = firstMonomer
         self.polymer = self.first
         self.atoms = self.polymer.atoms
+        self.dimensions = self.polymer.dimensions
 
     def select_atoms(self, selection) -> mda.AtomGroup:
         """

tests/data/CODI_bonds.pdb

@@ -0,0 +1,59 @@
+HEADER    UNCLASSIFIED                            22-Feb-24
+TITLE     UNITED ATOM STRUCTURE FOR MOLECULE UNL                                
+AUTHOR    AUTOMATED TOPOLOGY BUILDER (ATB) REVISION 2023-06-14 20:38:16
+AUTHOR   2  https://atb.uq.edu.au
+HETATM    1  O1B CODI    0      -4.879   3.567  -0.008  1.00  0.00           O
+HETATM    2  C1C CODI    0      -3.936   2.789   0.287  1.00  0.00           C
+HETATM    3  O1A CODI    0      -3.278   2.748   1.364  1.00  0.00           O
+HETATM    4  C1B CODI    0      -3.516   1.771  -0.816  1.00  0.00           C
+HETATM    5  C1A CODI    0      -2.746   0.598  -0.238  1.00  0.00           C
+HETATM    6  OG1 CODI    0      -2.378  -0.288  -1.342  1.00  0.00           O
+HETATM    7   CB CODI    0      -1.416  -1.221  -1.313  1.00  0.00           C
+HETATM    8  OG2 CODI    0      -1.267  -1.923  -2.298  1.00  0.00           O
+HETATM    9   CA CODI    0      -0.494  -1.414  -0.098  1.00  0.00           C
+HETATM   10   HB CODI    0       0.197  -2.204  -0.405  1.00  0.00           H
+HETATM   11    C CODI    0      -1.125  -1.829   1.248  1.00  0.00           C
+HETATM   12   CP CODI    0      -1.858  -3.187   1.213  1.00  0.00           C
+HETATM   13   CQ CODI    0      -3.239  -3.189   0.542  1.00  0.00           C
+HETATM   14   S1 CODI    0       0.580   0.090   0.045  1.00  0.00           S
+HETATM   15  CIA CODI    0       2.224  -0.537   0.246  1.00  0.00           C
+HETATM   16   S2 CODI    0       2.625  -2.120   0.529  1.00  0.00           S
+HETATM   17  CIB CODI    0       3.226   0.554   0.162  1.00  0.00           C
+HETATM   18 CIC2 CODI    0       3.092   1.591  -0.782  1.00  0.00           C
+HETATM   19 HIC2 CODI    0       2.256   1.588  -1.472  1.00  0.00           H
+HETATM   20 CID2 CODI    0       4.054   2.596  -0.870  1.00  0.00           C
+HETATM   21 HID2 CODI    0       3.946   3.376  -1.618  1.00  0.00           H
+HETATM   22 CIE1 CODI    0       5.152   2.593  -0.008  1.00  0.00           C
+HETATM   23 HIE1 CODI    0       5.898   3.379  -0.074  1.00  0.00           H
+HETATM   24 CID1 CODI    0       5.291   1.573   0.940  1.00  0.00           C
+HETATM   25 HID1 CODI    0       6.140   1.572   1.617  1.00  0.00           H
+HETATM   26 CIC1 CODI    0       4.344   0.556   1.020  1.00  0.00           C
+HETATM   27 HIC1 CODI    0       4.448  -0.233   1.757  1.00  0.00           H
+CONECT    1    2
+CONECT    2    1    3    4
+CONECT    3    2
+CONECT    4    2    5
+CONECT    5    4    6
+CONECT    6    5    7
+CONECT    7    6    8    9
+CONECT    8    7
+CONECT    9    7   10   11   14
+CONECT   10    9
+CONECT   11    9   12
+CONECT   12   11   13
+CONECT   13   12
+CONECT   14    9   15
+CONECT   15   14   16   17
+CONECT   16   15
+CONECT   17   15   18   26
+CONECT   18   17   19   20
+CONECT   19   18
+CONECT   20   18   21   22
+CONECT   21   20
+CONECT   22   20   23   24
+CONECT   23   22
+CONECT   24   22   25   26
+CONECT   25   24
+CONECT   26   17   24   27
+CONECT   27   26
+END

tests/data/PEI_end.pdb

@@ -0,0 +1,15 @@
+HETATM    1  N1  ELPR    0      -1.757   2.627   0.115  1.00  0.00           N
+HETATM    2  H1  ELPR    0      -1.791   3.280  -0.670  1.00  0.00           H
+HETATM    3  C2  ELPR    0      -2.966   1.803   0.070  1.00  0.00           C
+HETATM    4  C1  ELPR    0      -4.221   2.675   0.069  1.00  0.00           C
+HETATM    5  NX  ELPR    0      -5.436   1.859   0.055  1.00  0.00           N
+HETATM    6  H2  ELPR    0      -1.732   3.143   0.983  1.00  0.00           H
+HETATM    7  CZ  ELPR    0      -0.469   1.795  -0.026  1.00  0.00           C
+CONECT    4    3    5
+CONECT    3    4    1
+CONECT    2    1
+CONECT    6    1
+CONECT    1    3    6    7    2
+CONECT    5    4
+CONECT    7    1
+END

tests/data/PEI_monomer.pdb

@@ -0,0 +1,15 @@
+HETATM    1  N1  ELNR    0      -1.757   2.627   0.115  1.00  0.00           N
+HETATM    2  H1  ELNR    0      -1.791   3.280  -0.670  1.00  0.00           H
+HETATM    3  C2  ELNR    0      -2.966   1.803   0.070  1.00  0.00           C
+HETATM    4  C1  ELNR    0      -4.221   2.675   0.069  1.00  0.00           C
+HETATM    5  NX  ELNR    0      -5.436   1.859   0.055  1.00  0.00           N
+HETATM    6  H2  ELNR    0      -1.732   3.143   0.983  1.00  0.00           H
+HETATM    7  CX  ELNR    0      -0.469   1.795  -0.026  1.00  0.00           C
+CONECT    4    3    5
+CONECT    3    4    1
+CONECT    2    1
+CONECT    6    1
+CONECT    1    3    6    7    2
+CONECT    5    4
+CONECT    7    1
+END

tests/data/PEI_start.pdb

@@ -0,0 +1,15 @@
+HETATM    1  N1  ELPR    0      -1.757   2.627   0.115  1.00  0.00           N
+HETATM    2  H1  ELPR    0      -1.791   3.280  -0.670  1.00  0.00           H
+HETATM    3  C2  ELPR    0      -2.966   1.803   0.070  1.00  0.00           C
+HETATM    4  C1  ELPR    0      -4.221   2.675   0.069  1.00  0.00           C
+HETATM    5  C0  ELPR    0      -5.436   1.859   0.055  1.00  0.00           N
+HETATM    6  H2  ELPR    0      -1.732   3.143   0.983  1.00  0.00           H
+HETATM    7  CX  ELPR    0      -0.469   1.795  -0.026  1.00  0.00           C
+CONECT    4    3    5
+CONECT    3    4    1
+CONECT    2    1
+CONECT    6    1
+CONECT    1    3    6    7    2
+CONECT    5    4
+CONECT    7    1
+END

tests/polyconf/test_PDB.py

@@ -0,0 +1,49 @@
+#!/usr/bin/env python 
+import random
+import pytest
+from polyconf.polyconf.PDB import PDB
+from polyconf.polyconf.Polymer import Polymer
+from polyconf.polyconf.Monomer import Monomer
+import MDAnalysis as mda
+from .test_data import test_pdbs
+
+@pytest.mark.parametrize("pdb_file", test_pdbs)
+def test_pdb_initialization(data_dir, pdb_file):
+    polymer = Polymer(Monomer(f"{data_dir}/{pdb_file}"))
+    pdb = PDB(polymer)
+    assert pdb.polymer == polymer
+    assert len(pdb.atoms) == len(polymer.atoms)
+
+@pytest.mark.parametrize("pdb_file", test_pdbs)
+def test_cleanup(data_dir, pdb_file):
+    polymer = Polymer(Monomer(f"{data_dir}/{pdb_file}"))
+    pdb = PDB(polymer)
+    pdb.cleanup()
+    box = (pdb.atoms.positions).max(axis=0) - (pdb.atoms.positions).min(axis=0) + [10, 10, 10]
+    assert pdb.dimensions[:3] == pytest.approx(box, rel=1e-1), "Box is not at least 10 Angstroms larger than the polymer"
+    assert pdb.dimensions[3:] == [90, 90, 90], "Box angles are not 90 degrees"
+    assert (pdb.atoms.positions >= -box/2).all(), "Some atoms are located outside the box"
+
+@pytest.mark.parametrize("pdb_file", test_pdbs)
+def test_crude_save(data_dir, pdb_file, output_dir):
+    polymer = Polymer(Monomer(f"{data_dir}/{pdb_file}"))
+    pdb = PDB(polymer)
+    output_file = output_dir / "output.pdb"
+    pdb.crudesave(fname=str(output_dir / "output"))
+    assert output_file.exists()
+    with open(output_file, 'r') as f:
+        content = f.read()
+    for atom in pdb.polymer.atoms: # Test for every atom in the polymer
+        assert f" {atom.name} " in content, f"Atom {atom.name} not found in output file"
+
+@pytest.mark.parametrize("pdb_file", test_pdbs)
+def test_save(data_dir, pdb_file, output_dir):
+    polymer = Polymer(Monomer(f"{data_dir}/{pdb_file}"))
+    pdb = PDB(polymer)
+    output_file = output_dir / "output.gro"
+    pdb.save(fname=str(output_dir / "output"), gmx=True)
+    assert output_file.exists()
+    with open(output_file, 'r') as f:
+        content = f.read()
+    for atom in pdb.polymer.atoms: # Test for every atom in the polymer
+        assert f" {atom.name} " in content, f"Atom {atom.name} not found in output file"

tests/polyconf/test_data.py

@@ -0,0 +1,16 @@
+import pytest
+
+test_pdbs = ['CODI_bonds.pdb']
+
+@pytest.mark.parametrize("pdb_filename", test_pdbs)
+def test_data_files_exist(data_dir, pdb_filename):
+    assert (data_dir / pdb_filename).exists(), f"{pdb_filename} does not exist in {data_dir}"
+
+PEI_pdbs = ['PEI_start.pdb', 'PEI_end.pdb', 'PEI_monomer.pdb']
+PEI_start = ['PEI_start.pdb']
+PEI_monomer = ['PEI_monomer.pdb']
+PEI_end = ['PEI_end.pdb']
+
+@pytest.mark.parametrize("PEI_filename", PEI_pdbs)
+def test_PEI_files_exist(data_dir, PEI_filename):
+    assert (data_dir / PEI_filename).exists(), f"{PEI_filename} does not exist in {data_dir}"

tests/polyconf/test_import.py

@@ -0,0 +1,5 @@
+def test_import_polyconf():
+    try:
+        import polyconf
+    except ImportError:
+        assert False, "Failed to import PolyConf"

tests/polyconf/test_monomer.py

@@ -0,0 +1,39 @@
+#!/usr/bin/env python 
+import random
+import pytest
+from polyconf.polyconf.Monomer import Monomer
+import MDAnalysis as mda
+from .test_data import test_pdbs # This is a list of test pdb files that we tested for existence in the test_data.py file
+
+
+@pytest.mark.parametrize("pdb_filename", test_pdbs) # this parameterization will run this test for every element in test_pbs
+def test_monomer_initialization(data_dir, pdb_filename: str):
+    """
+    Test that we can initialize a Monomer object from a pdb file
+    """
+    monomer = Monomer(data_dir / pdb_filename)
+    assert isinstance(monomer.monomer, mda.Universe), "Monomer.monomer should be an MDAnalysis Universe object"
+    # count the lines in the pdb file that begin with HETATM
+    with open(data_dir / pdb_filename) as f:
+        hetatm_lines = [line for line in f if line.startswith("HETATM")]
+    assert len(monomer.atoms) == len(hetatm_lines), "The number of atoms in the Monomer should match the number of HETATM lines in the pdb file"
+
+@pytest.mark.parametrize("pdb_filename", test_pdbs)
+def test_select_atoms(data_dir, pdb_filename: str):
+    """
+    Test that we can find an atom by name in the monomer using the select_atoms method
+    """
+    monomer = Monomer(data_dir / pdb_filename)
+    random_atom = random.choice(monomer.atoms)
+    selected_atoms = monomer.select_atoms(f"name {random_atom.name}")
+    assert len(selected_atoms) == 1, f"There should only be one {random_atom.name} atom in {pdb_filename}"
+    assert selected_atoms[0].name == random_atom.name, f"The selected atom should be named {random_atom.name}"
+
+@pytest.mark.parametrize("pdb_filename", test_pdbs)
+def test_monomer_from_u(data_dir, pdb_filename: str):
+    """
+    Test that the monomer_from_u Monomer Class method works
+    """
+    universe = mda.Universe(data_dir / pdb_filename)
+    monomer = Monomer.monomer_from_u(universe)
+    assert len(monomer.atoms) == len(universe.atoms)