add ghost particles for bonds

flowermd/library/forcefields.py

@@ -613,7 +613,8 @@ def _create_forcefield(self):
         forces = []
         # Bonds
         bond = hoomd.md.bond.Harmonic()
-        bond.params["A-X"] = dict(k=50, r0=1.0)
+        bond.params["T-T"] = dict(k=50, r0=0.01)
+        bond.params["A-X"] = dict(k=0, r0=0)
         forces.append(bond)
         # Angles
         if all([self.angle_k, self.angle_theta0]):
@@ -630,10 +631,14 @@ def _create_forcefield(self):
         # Add zero pairs
         for pair in [
             ("R", "R"),
+            ("T", "T"),
+            ("T", "R"),
             ("A", "A"),
             ("A", "X"),
+            ("A", "T"),
             ("A", "R"),
             ("X", "R"),
+            ("X", "T"),
         ]:
             gb.params[pair] = dict(epsilon=0.0, lperp=0.0, lpar=0.0)
             gb.params[pair].r_cut = 0.0

flowermd/library/polymers.py

@@ -303,32 +303,40 @@ class EllipsoidChain(Polymer):
         The mass of the ellipsoid bead.
     """
 
-    def __init__(self, lengths, num_mols, lpar, bead_mass):
+    def __init__(self, lengths, num_mols, lpar, bead_mass, bond_L=0.1):
         self.bead_mass = bead_mass
         self.lpar = lpar
+        self.bond_L = bond_L
         # get the indices of the particles in a rigid body
         self.bead_constituents_types = ["X", "A"]
         super(EllipsoidChain, self).__init__(lengths=lengths, num_mols=num_mols)
 
     def _build(self, length):
         # Build up ellipsoid bead
         bead = mb.Compound(name="ellipsoid")
-        center = mb.Compound(pos=(0, 0, 0), name="X", mass=self.bead_mass / 2)
+        center = mb.Compound(pos=(0, 0, 0), name="X", mass=self.bead_mass / 4)
         head = mb.Compound(
-            pos=(self.lpar, 0, 0), name="A", mass=self.bead_mass / 2
+            pos=(self.lpar, 0, 0), name="A", mass=self.bead_mass / 4
         )
-        bead.add([center, head])
+        tether_head = mb.Compound(
+            pos=(self.lpar, 0, 0), name="T", mass=self.bead_mass / 4
+        )
+        tether_tail = mb.Compound(
+            pos=(-self.lpar, 0, 0), name="T", mass=self.bead_mass / 4
+        )
+        bead.add([center, head, tether_head, tether_tail])
         bead.add_bond([center, head])
 
         chain = mb.Compound()
         last_bead = None
         for i in range(length):
-            translate_by = np.array([i * self.lpar * 2, 0, 0])
+            translate_by = np.array([(i * self.lpar * 2) + self.bond_L, 0, 0])
             this_bead = mb.clone(bead)
             this_bead.translate(by=translate_by)
             chain.add(this_bead)
             if last_bead:
                 chain.add_bond([this_bead.children[0], last_bead.children[1]])
+                chain.add_bond([this_bead.children[3], last_bead.children[2]])
             last_bead = this_bead
 
         return chain

flowermd/utils/constraints.py

@@ -125,7 +125,7 @@ def create_rigid_ellipsoid_chain(
     """
     # find typeid sequence of the constituent particles types in a rigid bead
     # X and A
-    bead_len = 2
+    bead_len = 4
     # find indices that matches the constituent particle types
     typeids = snapshot.particles.typeid.reshape(-1, bead_len)
     matches = np.where((typeids == typeids))
@@ -212,7 +212,7 @@ def create_rigid_ellipsoid_chain(
 
     rigid_constrain = hoomd.md.constrain.Rigid()
     rigid_constrain.body["R"] = {
-        "constituent_types": ["X", "A"],
+        "constituent_types": ["X", "A", "T", "T"],
         "positions": local_coords,
         "orientations": [initial_orientation] * len(local_coords),
     }