Add classes for simulating ellipsoid polymer systems.

flowermd/library/__init__.py

@@ -15,6 +15,7 @@
     PEEK,
     PEKK,
     PPS,
+    EllipsoidChain,
     LJChain,
     PEKK_meta,
     PEKK_para,

flowermd/library/forcefields.py

@@ -481,3 +481,61 @@ def _check_widths(self):
                         "number of points for table energies and forces."
                     )
         return bond_width, angle_width, dih_width
+
+
+class EllipsoidForcefield(BaseHOOMDForcefield):
+    """A forcefield for modeling anisotropic bead polymers."""
+
+    def __init__(
+        self,
+        epsilon,
+        lper,
+        lpar,
+        bead_length,
+        r_cut,
+        bond_k,
+        bond_r0,
+        angle_k=None,
+        angle_theta0=None,
+    ):
+        self.epsilon = epsilon
+        self.lper = lper
+        self.lpar = lpar
+        self.bead_length = bead_length
+        self.r_cut = r_cut
+        self.bond_k = bond_k
+        self.bond_r0 = bond_r0
+        self.angle_k = angle_k
+        self.angle_theta0 = angle_theta0
+        hoomd_forces = self._create_forcefield()
+        super(TableForcefield, self).__init__(hoomd_forces)
+
+    def _create_forcefield(self):
+        forces = []
+        # Bonds
+        bond = hoomd.md.bond.Harmonic()
+        bond.params["A-A"] = dict(k=self.bond_k, r0=self.bond_r0)
+        bond.params["B-B"] = dict(k=0, r0=self.bead_length / 2)
+        forces.append(bond)
+        # Angles
+        if all([self.angle_k, self.angle_theta0]):
+            angle = hoomd.md.angle.Harmonic()
+            angle.params["B-B-B"] = dict(k=self.angle_k, t0=self.angle_theta0)
+            forces.append(angle)
+        # Gay-Berne Pairs
+        nlist = hoomd.md.nlist.Cell(buffer=0.40)
+        gb = hoomd.md.pair.aniso.GayBerne(nlist=nlist, default_r_cut=self.r_cut)
+        gb.params[("R", "R")] = dict(
+            epsilon=self.epsilon, lperp=self.lperp, lpar=self.lpar
+        )
+        # Add zero pairs
+        for pair in [
+            ("A", "A"),
+            ("B", "B"),
+            ("A", "B"),
+            ("A", "R"),
+            ("B", "R"),
+        ]:
+            gb.params[pair] = dict(epsilon=0.0, lperp=0.0, lpar=0.0)
+        forces.append(gb)
+        return forces

flowermd/library/polymers.py

@@ -3,6 +3,7 @@
 
 import mbuild as mb
 from mbuild.coordinate_transform import z_axis_transform
+from mbuild.lib.recipes import Polymer as mbPolymer
 
 from flowermd import CoPolymer, Polymer
 from flowermd.assets import MON_DIR
@@ -261,3 +262,42 @@ def _build(self, length):
                     chain.add_bond([next_bead, last_bead])
                 last_bead = next_bead
         return chain
+
+
+class EllipsoidChain(Polymer):
+    def __init__(
+        self, lengths, num_mols, bead_length, bead_name, bead_mass, bond_length
+    ):
+        self.bead_name = bead_name
+        self.bead_mass = bead_mass
+        self.bead_bond_length = bond_length
+        self.bead_length = bead_length
+        super(EllipsoidChain, self).__init__(lengths=lengths, num_mols=num_mols)
+
+    def _build(self, length):
+        # Build up ellipsoid bead
+        bead = mb.Compound(name="ellipsoid")
+        head = mb.Compound(
+            pos=(self.bead_length / 2, 0, 0), name="A", mass=self.bead_mass / 2
+        )
+        tail = mb.Compound(
+            pos=(-self.bead_length / 2, 0, 0), name="A", mass=self.bead_mass / 2
+        )
+        head_mid = mb.Compound(
+            pos=(self.bead_length / 4, 0, 0), name="B", mass=0
+        )
+        tail_mid = mb.Compound(
+            pos=(-self.bead_length / 4, 0, 0), name="B", mass=0
+        )
+        bead.add([head, tail, head_mid, tail_mid])
+
+        chain = mbPolymer()
+        chain.add_monomer(
+            bead,
+            indices=[0, 1],
+            orientation=[[1, 0, 0], [-1, 0, 0]],
+            replace=False,
+            separation=self.bead_bond_length,
+        )
+        chain.build(n=length, add_hydrogens=False)
+        return chain