Initial commit of bondAngle (one angle, two bonds)

mrdna/arbdmodel/__init__.py

@@ -71,6 +71,7 @@ class Parent():
         self.dihedrals = []
         self.impropers = []
         self.exclusions = []
+        self.bond_angles = []
 
         self.rigid = False
 
@@ -148,6 +149,13 @@ class Parent():
         # for b in (i,j): assert(b in beads)
         self.exclusions.append( (i,j) )
 
+    def add_bond_angle(self, i,j,k, bond_angle, exclude=False):
+        assert( len(set((i,j,k))) == 3 )
+        ## TODO: how to handle duplicating and cloning bonds
+        # beads = [b for b in self]
+        # for b in (i,j): assert(b in beads)
+        self.bond_angles.append( (i,j,k, bond_angle) )
+
     def get_restraints(self):
         ret = []
         for c in self.children:
@@ -200,6 +208,27 @@ class Parent():
         else:
             return ret
 
+    def get_bond_angles(self):
+        ret = self.bond_angles
+        for c in self.children:
+            if isinstance(c,Parent): ret.extend( c.get_bond_angles() )
+        if self.remove_duplicate_bonded_terms:
+            return list(set(ret))
+        else:
+            return ret
+
+    def _get_bond_potentials(self):
+        bonds =  [b for i,j,b,ex in self.get_bonds()]
+        bondangles1 = [b[1] for i,j,k,b in self.get_bond_angles()]
+        bondangles2 = [b[2] for i,j,k,b in self.get_bond_angles()]
+        return list(set( bonds+bondangles1+bondangles2 ))
+
+    def _get_angle_potentials(self):
+        angles =  [b for i,j,k,b in self.get_angles()]
+        bondangles = [b[0] for i,j,k,b in self.get_bond_angles()]
+        return list(set( angles+bondangles ))
+
+
     ## Removed because prohibitively slow
     # def remove_duplicate_terms(self):
     #     for key in "bonds angles dihedrals impropers exclusions".split():
@@ -863,6 +892,7 @@ class ArbdModel(PdbModel):
         self._angle_filename = "%s/%s.angles.txt" % (d, prefix)
         self._dihedral_filename = "%s/%s.dihedrals.txt" % (d, prefix)
         self._exclusion_filename = "%s/%s.exculsions.txt" % (d, prefix)
+        self._bond_angle_filename = "%s/%s.bondangles.txt" % (d, prefix)
         
         # self._writeArbdCoordFile( prefix + ".coord.txt" )
         self._writeArbdParticleFile( prefix + ".particles.txt" )
@@ -871,6 +901,7 @@ class ArbdModel(PdbModel):
         self._writeArbdAngleFile()
         self._writeArbdDihedralFile()
         self._writeArbdExclusionFile()
+        self._writeArbdBondAngleFile()
         self._writeArbdPotentialFiles( prefix, directory = d )
         self._writeArbdConf( prefix, numSteps=numSteps, outputPeriod=outputPeriod, restart_file=restart_file )
         
@@ -1066,13 +1097,13 @@ tabulatedPotential  1
             angles = self.get_angles()
             dihedrals = self.get_dihedrals()
             exclusions = self.get_exclusions()
-
+            bond_angles = self.get_bond_angles()
             if len(bonds) > 0:
-                for b in list(set([b for i,j,b,ex in bonds])):
+                for b in self._get_bond_potentials():
                     fh.write("tabulatedBondFile %s\n" % b)
 
             if len(angles) > 0:
-                for b in list(set([b for i,j,k,b in angles])):
+                for b in self._get_angle_potentials():
                     fh.write("tabulatedAngleFile %s\n" % b)
 
             if len(dihedrals) > 0:
@@ -1089,6 +1120,8 @@ tabulatedPotential  1
                 fh.write("inputDihedrals %s\n" % self._dihedral_filename)
             if len(exclusions) > 0:
                 fh.write("inputExcludes %s\n" % self._exclusion_filename)
+            if len(bond_angles) > 0:
+                fh.write("inputBondAngles %s\n" % self._bond_angle_filename)
      
         write_null_dx = False
         for pt,num in self.getParticleTypesAndCounts():
@@ -1164,7 +1197,7 @@ component "data" value 3
                 fh.write("RESTRAINT %d %f %f %f %f\n" % tuple(item))
 
     def _writeArbdBondFile( self ):
-        for b in list( set( [b for i,j,b,ex in self.get_bonds()] ) ):
+        for b in self._get_bond_potentials():
             if type(b) is not str and not isinstance(b, Path):
                 b.write_file()
 
@@ -1177,7 +1210,7 @@ component "data" value 3
                     fh.write("BOND ADD %d %d %s\n" % item)
 
     def _writeArbdAngleFile( self ):
-        for b in list( set( [b for i,j,k,b in self.get_angles()] ) ):
+        for b in self._get_angle_potentials():
             if type(b) is not str and not isinstance(b, Path):
                 b.write_file()
 
@@ -1202,6 +1235,13 @@ component "data" value 3
                 item = tuple(int(p.idx) for p in ex)
                 fh.write("EXCLUDE %d %d\n" % item)
 
+    def _writeArbdBondAngleFile( self ):
+        if len(self.bond_angles) > 0:
+            with open(self._bond_angle_filename,'w') as fh:
+                for b in self.get_bond_angles():
+                    item = tuple([p.idx for p in b[:-1]] + [str(p) for p in b[-1]])
+                    fh.write("BONDANGLE %d %d %d %s %s %s\n" % item)
+
     def dimensions_from_structure( self, padding_factor=1.5, isotropic=False ):
         ## TODO: cache coordinates using numpy arrays for quick min/max
         raise(NotImplementedError)

mrdna/arbdmodel/nonbonded.py

@@ -58,12 +58,13 @@ class TabulatedPotential(NonbondedScheme):
 
 ## Bonded potentials
 class BasePotential():
-    def __init__(self, r0, rRange=(0,50), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/"):
+    def __init__(self, r0, rRange=(0,50), resolution=0.1, maxForce=None, max_potential=None, zero="min", prefix="potentials/"):
         self.r0 = r0
         self.rRange = rRange
         self.resolution = resolution
         self.maxForce = maxForce
         self.prefix = prefix
+        self.zero = zero
         self.periodic = False
         self.type_ = "None"
         self.max_potential = max_potential
@@ -111,15 +112,16 @@ class BasePotential():
             u[0] = 0
             u[1:] = np.cumsum(f*np.diff(r))
 
+        if self.zero == "min":
             u = u - np.min(u)
 
         np.savetxt( self.filename(), np.array([r, u]).T, fmt="%f" )
 
 class HarmonicPotential(BasePotential):
-    def __init__(self, k, r0, rRange, resolution, maxForce, max_potential, prefix, correct_geometry=False):
+    def __init__(self, k, r0, rRange, resolution, maxForce, max_potential, zero="min", correct_geometry=False, prefix='./'):
         self.k = k
         self.kscale_ = None
-        BasePotential.__init__(self, r0, rRange, resolution, maxForce, max_potential, prefix)
+        BasePotential.__init__(self, r0, rRange, resolution, maxForce, max_potential, zero, prefix)
 
     def filename(self):
         return "%s%s-%.3f-%.3f.dat" % (self.prefix, self.type_,
@@ -143,8 +145,8 @@ class HarmonicPotential(BasePotential):
 #         self.kscale_ = 1.0
 
 class HarmonicBond(HarmonicPotential):
-    def __init__(self, k, r0, rRange=(0,50), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/", correct_geometry=False, temperature=295):
-        HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, max_potential, prefix)
+    def __init__(self, k, r0, rRange=(0,50), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/",  zero="min", correct_geometry=False, temperature=295):
+        HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, max_potential, zero=zero, correct_geometry=correct_geometry, prefix=prefix)
         self.type_ = "gbond" if correct_geometry else "bond"
         self.kscale_ = 1.0
         self.correct_geometry = correct_geometry
@@ -161,22 +163,22 @@ class HarmonicBond(HarmonicPotential):
 
 
 class HarmonicAngle(HarmonicPotential):
-    def __init__(self, k, r0, rRange=(0,181), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/"):
-        HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, max_potential, prefix)
+    def __init__(self, k, r0, rRange=(0,181), resolution=0.1, maxForce=None, max_potential=None, zero="min", prefix="potentials/"):
+        HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, max_potential, zero=zero, prefix=prefix)
         self.type_ = "angle"
         self.kscale_ = (180.0/np.pi)**2
 
 class HarmonicDihedral(HarmonicPotential):
-    def __init__(self, k, r0, rRange=(-180,180), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/"):
-        HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, max_potential, prefix)
+    def __init__(self, k, r0, rRange=(-180,180), resolution=0.1, maxForce=None, max_potential=None, zero="min", prefix="potentials/"):
+        HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, max_potential, zero=zero, prefix=prefix)
         self.periodic = True
         self.type_ = "dihedral"
         self.kscale_ = (180.0/np.pi)**2
 
 class WLCSKBond(BasePotential):
     """ ## https://aip.scitation.org/doi/full/10.1063/1.4968020 """
-    def __init__(self, d, lp, kT, rRange=(0,50), resolution=0.02, maxForce=100, max_potential=None, prefix="potentials/"):
-        BasePotential.__init__(self, 0, rRange, resolution, maxForce, max_potential, prefix)
+    def __init__(self, d, lp, kT, rRange=(0,50), resolution=0.02, maxForce=100, max_potential=None, zero="min", prefix="potentials/"):
+        BasePotential.__init__(self, 0, rRange, resolution, maxForce, max_potential, zero, prefix)
         self.type_ = "wlcbond"
         self.d = d          # separation
         self.lp = lp            # persistence length
@@ -213,8 +215,8 @@ class WLCSKBond(BasePotential):
 
 class WLCSKAngle(BasePotential):
     ## https://aip.scitation.org/doi/full/10.1063/1.4968020
-    def __init__(self, d, lp, kT, rRange=(0,181), resolution=0.5, maxForce=None, max_potential=None, prefix="potentials/"):
-        BasePotential.__init__(self, 180, rRange, resolution, maxForce, max_potential, prefix)
+    def __init__(self, d, lp, kT, rRange=(0,181), resolution=0.5, maxForce=None, max_potential=None, zero="min", prefix="potentials/"):
+        BasePotential.__init__(self, 180, rRange, resolution, maxForce, max_potential, zero, prefix)
         self.type_ = "wlcangle"
         self.d = d          # separation
         self.lp = lp            # persistence length