Adding DPD System initialization

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

@@ -0,0 +1,58 @@
+<ForceField version="0.0.1" name="HOOMD-DPD">
+    <FFMetaData combiningRule="lorentz">
+        <Units energy="kJ/mol" distance="nm" mass="amu" time="s" charge="coulomb"/>
+    </FFMetaData>
+    <AtomTypes expression="0">
+       <AtomType name="_A" element="_A" atomclass="_A" mass="1.0" charge="0.0" definition="_A" description="CG Bead">
+           <Parameters>
+           </Parameters>
+       </AtomType>
+    </AtomTypes>
+    <BondTypes expression="0.5 * k * (r-r_eq)**2">
+        <ParametersUnitDef parameter="r_eq" unit="nm"/>
+        <ParametersUnitDef parameter="k" unit="kJ/mol/(nm**2)"/>
+        <BondType name="HarmonicBondPotential" type1="_A" type2="_A">
+            <Parameters>
+                <Parameter name="r_eq" value="1.0"/>
+                <Parameter name="k" value="10"/>
+            </Parameters>
+        </BondType>
+    </BondTypes>
+    <AngleTypes expression="0.5 * k * (theta-theta_eq)**2">
+        <ParametersUnitDef parameter="theta_eq" unit="rad"/>
+        <ParametersUnitDef parameter="k" unit="kJ/mol/(rad**2)"/>
+        <AngleType name="HarmonicAnglePotential" type1="_A" type2="_A" type3="_A">
+            <Parameters>
+                <Parameter name="theta_eq" value="1.0"/>
+                <Parameter name="k" value="10"/>
+            </Parameters>
+        </AngleType>
+    </AngleTypes>
+    <DihedralTypes expression="0.5 * k * (1 + d * cos(n*phi-phi_eq))">
+        <ParametersUnitDef parameter="phi_eq" unit="rad"/>
+        <ParametersUnitDef parameter="k" unit="kJ/mol"/>
+        <ParametersUnitDef parameter="d" unit="dimensionless"/>
+        <ParametersUnitDef parameter="n" unit="dimensionless"/>
+        <DihedralType name="HOOMDPeriodicDihedralPotential" type1="_A" type2="_A" type3="_A" type4="_A">
+            <Parameters>
+                <Parameter name="phi_eq" value="0"/>
+                <Parameter name="k" value="10"/>
+                <Parameter name="d" value="-1"/>
+                <Parameter name="n" value="3"/>
+            </Parameters>
+        </DihedralType>
+    </DihedralTypes>
+
+  <PairPotentialTypes expression="A * (1-(r/r_cut)) - γ">
+    <ParametersUnitDef parameter="A" unit="N"/>
+    <ParametersUnitDef parameter="r_cut" unit="nm"/>
+    <ParametersUnitDef parameter="γ" unit="amu/s"/>
+    <PairPotentialType name="HOOMDDPDForce" type1="_A" type2="_A">
+        <Parameters>
+            <Parameter name="A" value="40.0"/>
+	    <Parameter name="r_cut" value="1.0"/>
+	    <Parameter name="γ" value="8.0"/>
+        </Parameters>
+    </PairPotentialType>
+  </PairPotentialTypes>
+</ForceField>

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

@@ -0,0 +1,51 @@
+<ForceField version="0.0.1" name="HOOMD-DPD">
+    <FFMetaData combiningRule="lorentz">
+        <Units energy="kJ/mol" distance="nm" mass="amu" time="s" charge="coulomb"/>
+    </FFMetaData>
+    <AtomTypes name="HOOMDDPDForce" expression="A*(-r/r_cut + 1) - γ">
+        <ParametersUnitDef parameter="A" unit="N"/>
+        <ParametersUnitDef parameter="r_cut" unit="nm"/>
+        <ParametersUnitDef parameter="γ" unit="amu/s"/>
+        <AtomType name="_A" element="_A" atomclass="_A" mass="1.0" charge="0.0" definition="_A" description="CG Bead">
+            <Parameters>
+                <Parameter name="A" value="40.0"/>
+	            <Parameter name="r_cut" value="1.0"/>
+	            <Parameter name="γ" value="8.0"/>
+            </Parameters>
+        </AtomType>
+    </AtomTypes>
+    <BondTypes expression="0.5 * k * (r-r_eq)**2">
+        <ParametersUnitDef parameter="r_eq" unit="nm"/>
+        <ParametersUnitDef parameter="k" unit="kJ/mol/(nm**2)"/>
+        <BondType name="HarmonicBondPotential" type1="_A" type2="_A">
+            <Parameters>
+                <Parameter name="r_eq" value="1.0"/>
+                <Parameter name="k" value="10"/>
+            </Parameters>
+        </BondType>
+    </BondTypes>
+    <AngleTypes expression="0.5 * k * (theta-theta_eq)**2">
+        <ParametersUnitDef parameter="theta_eq" unit="radian"/>
+        <ParametersUnitDef parameter="k" unit="kJ/mol/(radian**2)"/>
+        <AngleType name="HarmonicAnglePotential" type1="_A" type2="_A" type3="_A">
+            <Parameters>
+                <Parameter name="theta_eq" value="1.0"/>
+                <Parameter name="k" value="10"/>
+            </Parameters>
+        </AngleType>
+    </AngleTypes>
+    <DihedralTypes expression="0.5 * k * (1 + d * cos(n*phi-phi0))">
+        <ParametersUnitDef parameter="phi0" unit="radian"/>
+        <ParametersUnitDef parameter="k" unit="kJ/mol"/>
+        <ParametersUnitDef parameter="d" unit="dimensionless"/>
+        <ParametersUnitDef parameter="n" unit="dimensionless"/>
+        <DihedralType name="HOOMDPeriodicDihedralPotential" type1="_A" type2="_A" type3="_A" type4="_A">
+            <Parameters>
+                <Parameter name="phi0" value="0"/>
+                <Parameter name="k" value="10"/>
+                <Parameter name="d" value="-1"/>
+                <Parameter name="n" value="3"/>
+            </Parameters>
+        </DihedralType>
+    </DihedralTypes>
+</ForceField>

flowermd/base/forcefield.py

@@ -1,19 +1,25 @@
 """Base forcefield classes."""
 
 import forcefield_utilities as ffutils
-import foyer
+from gmso.core.forcefield import ForceField
 
 
-class BaseXMLForcefield(foyer.Forcefield):
+class BaseXMLForcefield:
     """Base XML forcefield class."""
 
-    def __init__(self, forcefield_files=None, name=None):
-        super(BaseXMLForcefield, self).__init__(
-            forcefield_files=forcefield_files, name=name
-        )
-        self.gmso_ff = (
-            ffutils.FoyerFFs().load(forcefield_files or name).to_gmso_ff()
-        )
+    def __init__(self, forcefield_files=None, name=None, gmso_xml=None):
+        self.forcefield_files = forcefield_files
+        self.name = name
+        self.gmso_xml = gmso_xml
+        if all([name, forcefield_files]):
+            raise ValueError("Give only one of `name` or `forcefield_files`.")
+        if self.gmso_xml is True:
+            self.gmso_ff = ForceField(forcefield_files or name)
+
+        else:
+            self.gmso_ff = (
+                ffutils.FoyerFFs().load(forcefield_files or name).to_gmso_ff()
+            )
 
 
 class BaseHOOMDForcefield:

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

@@ -13,7 +13,7 @@
 class GAFF(BaseXMLForcefield):
     """General Amber forcefield class."""
 
-    def __init__(self, forcefield_files=f"{FF_DIR}/gaff.xml"):
+    def __init__(self, forcefield_files=f"{FF_DIR}/gaff.xml", gmso_xml=False):
         super(GAFF, self).__init__(forcefield_files=forcefield_files)
         self.description = (
             "The General Amber Forcefield written in foyer XML format. "
@@ -25,15 +25,17 @@ def __init__(self, forcefield_files=f"{FF_DIR}/gaff.xml"):
 class OPLS_AA(BaseXMLForcefield):
     """OPLS All Atom forcefield class."""
 
-    def __init__(self, name="oplsaa"):
+    def __init__(self, name="oplsaa", gmso_xml=False):
         super(OPLS_AA, self).__init__(name=name)
         self.description = "opls-aa forcefield found in the Foyer package."
 
 
 class OPLS_AA_PPS(BaseXMLForcefield):
     """OPLS All Atom for PPS molecule forcefield class."""
 
-    def __init__(self, forcefield_files=f"{FF_DIR}/pps_opls.xml"):
+    def __init__(
+        self, forcefield_files=f"{FF_DIR}/pps_opls.xml", gmso_xml=False
+    ):
         super(OPLS_AA_PPS, self).__init__(forcefield_files=forcefield_files)
         self.description = (
             "Based on flowermd.forcefields.OPLS_AA. "
@@ -49,7 +51,9 @@ def __init__(self, forcefield_files=f"{FF_DIR}/pps_opls.xml"):
 class OPLS_AA_BENZENE(BaseXMLForcefield):
     """OPLS All Atom for benzene molecule forcefield class."""
 
-    def __init__(self, forcefield_files=f"{FF_DIR}/benzene_opls.xml"):
+    def __init__(
+        self, forcefield_files=f"{FF_DIR}/benzene_opls.xml", gmso_xml=False
+    ):
         super(OPLS_AA_BENZENE, self).__init__(forcefield_files=forcefield_files)
         self.description = (
             "Based on flowermd.forcefields.OPLS_AA. "
@@ -60,7 +64,11 @@ def __init__(self, forcefield_files=f"{FF_DIR}/benzene_opls.xml"):
 class OPLS_AA_DIMETHYLETHER(BaseXMLForcefield):
     """OPLS All Atom for dimethyl ether molecule forcefield class."""
 
-    def __init__(self, forcefield_files=f"{FF_DIR}/dimethylether_opls.xml"):
+    def __init__(
+        self,
+        forcefield_files=f"{FF_DIR}/dimethylether_opls.xml",
+        gmso_xml=False,
+    ):
         super(OPLS_AA_DIMETHYLETHER, self).__init__(
             forcefield_files=forcefield_files
         )
@@ -70,11 +78,23 @@ 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, gmso_xml=True
+        )
+        self.description = "DPD forcefield loaded from an XML file."
+
+
 class FF_from_file(BaseXMLForcefield):
     """Forcefield class for loading a forcefield from an XML file."""
 
-    def __init__(self, forcefield_files):
-        super(FF_from_file, self).__init__(forcefield_files=forcefield_files)
+    def __init__(self, forcefield_files, gmso_xml):
+        super(FF_from_file, self).__init__(
+            forcefield_files=forcefield_files, gmso_xml=gmso_xml
+        )
         self.description = "Forcefield loaded from an XML file. "
 
 
@@ -807,3 +827,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