Adding DPD System initialization

flowermd/assets/forcefields/hoomd-dpd-hhp.xml

@@ -0,0 +1,11 @@
+<ForceField version="0.0.1" name="HOOMD-DPD">
+    <AtomTypes>
+       <Type name="A" element="A" mass="1.0" def="A" desc="CG Bead" overrides=""/>
+    </AtomTypes>
+    <HarmonicBondForce>
+        <Bond class1="A" class2="A" length="1.0" k="20000"/>
+    </HarmonicBondForce>
+    <NonbondedForce>
+       <Atom type="A" A="800" gamma="800" r_cut="1.15"/>
+    </NonbondedForce>
+</ForceField>

flowermd/library/__init__.py

@@ -2,13 +2,15 @@
 """Library of predefined molecules, recipes and forcefields."""
 
 from .forcefields import (
+    DPD,
     GAFF,
     OPLS_AA,
     OPLS_AA_BENZENE,
     OPLS_AA_DIMETHYLETHER,
     OPLS_AA_PPS,
     BaseHOOMDForcefield,
     BaseXMLForcefield,
+    Bead_Spring_DPD,
     BeadSpring,
     EllipsoidFF_DPD,
     EllipsoidForcefield,
@@ -29,4 +31,4 @@
 )
 from .simulations.tensile import Tensile
 from .surfaces import Graphene
-from .systems import SingleChainSystem, mbuildSystem
+from .systems import RandomWalk, SingleChainSystem, mbuildSystem

flowermd/library/forcefields.py

@@ -70,6 +70,14 @@ def __init__(self, forcefield_files=f"{FF_DIR}/dimethylether_opls.xml"):
         )
 
 
+class Bead_Spring_DPD(BaseXMLForcefield):
+    """Forcefield class for loading a forcefield from an XML file."""
+
+    def __init__(self, forcefield_files=f"{FF_DIR}/hoomd-dpd-hhp.xml"):
+        super(Bead_Spring_DPD, self).__init__(forcefield_files=forcefield_files)
+        self.description = "DPD forcefield loaded from an XML file."
+
+
 class FF_from_file(BaseXMLForcefield):
     """Forcefield class for loading a forcefield from an XML file."""
 
@@ -807,3 +815,93 @@ def _create_forcefield(self):
             dpd.params[pair].r_cut = 0.0
         forces.append(dpd)
         return forces
+
+
+class DPD(BaseHOOMDForcefield):
+    """A DPD forcefield to use with bead-spring systems.
+
+    Notes
+    -----
+    This is designed to be used with `flowermd.library.polymers.LJChain`
+
+    The set of interactions are:
+    1. `hoomd.md.bond.Harmonic`
+    3. `hoomd.md.pair.DPD`
+
+    Parameters
+    ----------
+    epsilon : float, required
+        energy
+    lpar: float, required
+        Semi-axis length of the ellipsoid along the major axis.
+    lperp : float, required
+        Semi-axis length of the ellipsoid along the minor axis.
+    A : int, required
+        DPD pair-wise drag force coefficient
+    gamma : int, required
+        DPD pair-wise random force coefficient
+    kT : float, required
+        Temperature used in pair-wise drag force
+    r_cut : float, required
+        Cut off radius for pair interactions
+    angle_k : float, required
+        Spring constant in harmonic angle.
+    angle_theta0: float, required
+        Equilibrium angle between 2 consecutive beads.
+    bond_k : float, required
+        Spring constant in harmonic bond.
+    bond_r0: float, required
+        Equilibrium distance between 2 ellipsoid tips.
+    nlist : type, default hoomd.md.nlist.Cell
+        A class (not an instance) of the HOOMD neighbor list
+        to use for the pair force.
+    nlist_buffer : float, default 0.40
+        The buffer value (distance) used by the neighbor list.
+
+    """
+
+    def __init__(
+        self,
+        A,
+        gamma,
+        kT,
+        r_cut,
+        angle_k=None,
+        angle_theta0=None,
+        bond_k=100,
+        bond_r0=1.1,
+        nlist=hoomd.md.nlist.Cell,
+        nlist_buffer=0.40,
+    ):
+        self.gamma = gamma
+        self.A = A
+        self.kT = kT
+        self.r_cut = r_cut
+        self.angle_k = angle_k
+        self.angle_theta0 = angle_theta0
+        self.bond_k = bond_k
+        self.bond_r0 = bond_r0
+        self.nlist = nlist
+        self.nlist_buffer = nlist_buffer
+        hoomd_forces = self._create_forcefield()
+        super(DPD, 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)
+        forces.append(bond)
+        # Angles
+        if all([self.angle_k, self.angle_theta0]):
+            angle = hoomd.md.angle.Harmonic()
+            angle.params["A-A-A"] = dict(k=self.angle_k, t0=self.angle_theta0)
+            forces.append(angle)
+        # DPD Pairs
+        nlist = self.nlist(buffer=self.nlist_buffer, exclusions=["bond"])
+        dpd = hoomd.md.pair.DPD(
+            nlist=nlist, kT=self.kT, default_r_cut=self.r_cut
+        )
+        dpd.params[("A", "A")] = dict(A=self.A, gamma=self.gamma)
+        forces.append(dpd)
+        return forces

flowermd/library/simulations/dpd_init.py

@@ -0,0 +1,62 @@
+"""DPD energy relaxation simulation class."""
+
+import hoomd
+
+from flowermd.base.simulation import Simulation
+from flowermd.utils.dpd_utils import simulation_energy_end
+
+
+class DPDInit(Simulation):
+    """ """
+
+    def __init__(
+        self,
+        initial_state,
+        forcefield,
+        A,
+        r,
+        r_cut,
+        N,
+        sim_steps_incr,
+        box,
+        reference_values=dict(),
+        dt=0.0001,
+        device=hoomd.device.auto_select(),
+        seed=42,
+        gsd_write_freq=1e4,
+        gsd_file_name="trajectory.gsd",
+        log_write_freq=1e3,
+        log_file_name="log.txt",
+    ):
+        self.A = A
+        self.r = r
+        self.r_cut = r_cut
+        self.N = N
+        self.sim_steps_incr = sim_steps_incr
+        self.L = box
+        super(DPDInit, self).__init__(
+            initial_state=initial_state,
+            forcefield=forcefield,
+            reference_values=reference_values,
+            dt=dt,
+            device=device,
+            seed=seed,
+            gsd_write_freq=gsd_write_freq,
+            gsd_file_name=gsd_file_name,
+            log_write_freq=log_write_freq,
+            log_file_name=log_file_name,
+        )
+        print(self.A, self.r, self.r_cut, self.density)
+        self.run_NVE(n_steps=1)
+        while not simulation_energy_end(
+            A=self.A,
+            r=self.r,
+            r_cut=self.r_cut,
+            N=self.N,
+            L=self.L,
+            log_file_name=self.log_file_name,
+        ):
+            self.run_NVE(n_steps=self.sim_steps_incr)
+            for writer in self.operations.writers:
+                if hasattr(writer, "flush"):
+                    writer.flush()