Add support for mBuild's periodic polymer bonding (cmelab#148)

* add in mBuild perioidc polymer bond functionality

* [pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci

* update pytest action versions

* add notes to Lattice class

* add check for periodic bond axis param

* use longer molecules in lattice test

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

.github/workflows/pytest.yml

@@ -36,10 +36,10 @@ jobs:
 
     steps:
     - name: Check out repository
-      uses: actions/checkout@v2
+      uses: actions/checkout@v4
 
     - name: Build environment
-      uses: conda-incubator/setup-miniconda@v2
+      uses: conda-incubator/setup-miniconda@v3
       with:
         environment-file: environment-dev.yml
         python-version: ${{ matrix.python-version }}
@@ -56,7 +56,7 @@ jobs:
       run: python -m pytest -rs -v --cov=./ --cov-report=xml
 
     - name: Upload coverage to Codecov
-      uses: codecov/codecov-action@v3
+      uses: codecov/codecov-action@v4
       with:
         token: ${{ secrets.CODECOV_TOKEN }}
         file: ./coverage.xml

flowermd/base/molecule.py

@@ -477,6 +477,11 @@ class Polymer(Molecule):
         The bond length between connected atoms (units: nm)
     bond_orientation: list, default None
         The orientation of the bond between connected atoms.
+    periodic_bond_axis : str, default None
+        Axis which to orient the polymer backbone along.
+        Once the chain is aligned, a periodic bond between
+        head and tail atoms is formed.
+        Options are "x", "y", or "z"
 
     """
 
@@ -490,12 +495,14 @@ def __init__(
         bond_indices=None,
         bond_length=None,
         bond_orientation=None,
+        periodic_bond_axis=None,
         **kwargs,
     ):
         self.lengths = check_return_iterable(lengths)
         self.bond_indices = bond_indices
         self.bond_length = bond_length
         self.bond_orientation = bond_orientation
+        self.periodic_bond_axis = periodic_bond_axis
         num_mols = check_return_iterable(num_mols)
         if len(num_mols) != len(self.lengths):
             raise ValueError("Number of molecules and lengths must be equal.")
@@ -513,14 +520,27 @@ def monomer(self):
         return self._mb_molecule
 
     def _build(self, length):
+        if self.periodic_bond_axis:
+            if not isinstance(
+                self.periodic_bond_axis, str
+            ) or self.periodic_bond_axis.lower() not in ["x", "y", "z"]:
+                raise ValueError(
+                    "Valid choices for a `periodic_bond_axis` are "
+                    "'x', 'y', 'z'"
+                )
+            add_hydrogens = False
+        else:
+            add_hydrogens = True
         chain = mbPolymer()
         chain.add_monomer(
             self.monomer,
             indices=self.bond_indices,
             separation=self.bond_length,
             orientation=self.bond_orientation,
         )
-        chain.build(n=length, sequence="A")
+        chain.build(n=length, sequence="A", add_hydrogens=add_hydrogens)
+        if self.periodic_bond_axis:
+            chain.create_periodic_bond(axis=self.periodic_bond_axis)
         return chain
 
     def _generate(self):

flowermd/base/system.py

@@ -724,6 +724,10 @@ class Lattice(System):
     lattice is made by repeating and translating in the
     x and y directions.
 
+    See the `periodic_bond_axis` paramter in `flowermd.base.Polymer`
+    if you wish to form head-tail bonds across the periodic boundary
+    in the lattice.
+
     """
 
     def __init__(

flowermd/tests/base/test_molecule.py

@@ -306,3 +306,22 @@ def test_align_z_axis(self, polyethylene):
         for mol in pe.molecules:
             backbone = get_backbone_vector(mol.xyz)
             assert np.allclose(np.abs(backbone), np.array([0, 0, 1]), atol=1e-1)
+
+    @pytest.mark.parametrize("axis", ["x", "y", "z"])
+    def test_periodic_bond(self, polyethylene, axis):
+        pe_no_bond = polyethylene(num_mols=1, lengths=20)
+        n_bonds = pe_no_bond.molecules[0].n_bonds
+        n_particles = pe_no_bond.molecules[0].n_particles
+        pe_with_bond = polyethylene(
+            num_mols=1, lengths=20, periodic_bond_axis=axis
+        )
+        n_bonds_with = pe_with_bond.molecules[0].n_bonds
+        n_particles_with = pe_with_bond.molecules[0].n_particles
+        assert n_bonds - n_bonds_with == 1
+        assert n_particles - n_particles_with == 2
+
+    def test_periodic_bond_bad_axis(self, polyethylene):
+        with pytest.raises(ValueError):
+            polyethylene(num_mols=1, lengths=20, periodic_bond_axis=1)
+        with pytest.raises(ValueError):
+            polyethylene(num_mols=1, lengths=20, periodic_bond_axis="a")

flowermd/tests/base/test_system.py

@@ -799,6 +799,33 @@ def test_lattice_polymer(self, polyethylene):
                     np.abs(backbone), np.array([0, 0, 1]), atol=1e-1
                 )
 
+    @pytest.mark.parametrize("axis", ["x", "y", "z"])
+    def test_lattice_periodic_polymer(self, polyethylene, axis):
+        polyethylene = polyethylene(
+            lengths=10, num_mols=32, periodic_bond_axis=axis
+        )
+        system = Lattice(
+            molecules=[polyethylene],
+            x=1,
+            y=1,
+            n=4,
+        )
+        system.apply_forcefield(
+            r_cut=2.5, force_field=[OPLS_AA()], auto_scale=True
+        )
+
+        assert system.n_mol_types == 1
+        assert len(system.all_molecules) == len(polyethylene.molecules)
+        assert len(system.hoomd_forcefield) > 0
+        assert system.n_particles == system.hoomd_snapshot.particles.N
+        assert system.reference_values.keys() == {"energy", "length", "mass"}
+        for mol_class in system._molecules:
+            for mol in mol_class.molecules:
+                backbone = get_backbone_vector(mol.xyz)
+                assert np.allclose(
+                    np.abs(backbone), np.array([0, 0, 1]), atol=1e-1
+                )
+
     def test_lattice_molecule(self, benzene_molecule):
         benzene_mol = benzene_molecule(n_mols=32)
         system = Lattice(