Added support to arbdmodel for product potentials

mrdna/arbdmodel/__init__.py

@@ -72,6 +72,7 @@ class Parent():
         self.impropers = []
         self.exclusions = []
         self.bond_angles = []
+        self.product_potentials = []
 
         self.rigid = False
 
@@ -156,6 +157,18 @@ class Parent():
         # for b in (i,j): assert(b in beads)
         self.bond_angles.append( (i,j,k,l, bond_angle) )
 
+    def add_product_potential_angle(self, potential_list):
+        """ potential_list: list of tuples of form (particle_i, particle_j,..., TabulatedPotential) """
+        if len(terms) < 2: raise ValueError("Too few potentials")
+        for elem in potential_list:
+            beads = elem[:-1]
+            pot = elem[-1]
+            if len(beads) < 2: raise ValueError("Too few particles specified in product_potential")
+            if len(beads) > 4: raise ValueError("Too many particles specified in product_potential")
+
+        self.product_potentials.append(potential_list)
+        ## TODO: how to handle duplicating and cloning bonds
+
     def get_restraints(self):
         ret = []
         for c in self.children:
@@ -217,6 +230,15 @@ class Parent():
         else:
             return ret
 
+    def get_product_potentials(self):
+        ret = self.product_potentials
+        for c in self.children:
+            if isinstance(c,Parent): ret.extend( c.get_product_potentials() )
+        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,l,b in self.get_bond_angles()]
@@ -893,6 +915,7 @@ class ArbdModel(PdbModel):
         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._product_potential_filename = "%s/%s.prodpot.txt" % (d, prefix)
         
         # self._writeArbdCoordFile( prefix + ".coord.txt" )
         self._writeArbdParticleFile( prefix + ".particles.txt" )
@@ -902,6 +925,7 @@ class ArbdModel(PdbModel):
         self._writeArbdDihedralFile()
         self._writeArbdExclusionFile()
         self._writeArbdBondAngleFile()
+        self._writeArbdProductPotentialFile()
         self._writeArbdPotentialFiles( prefix, directory = d )
         self._writeArbdConf( prefix, numSteps=numSteps, outputPeriod=outputPeriod, restart_file=restart_file )
         
@@ -1098,6 +1122,7 @@ tabulatedPotential  1
             dihedrals = self.get_dihedrals()
             exclusions = self.get_exclusions()
             bond_angles = self.get_bond_angles()
+            prod_pots = self.get_product_potentials()
             if len(bonds) > 0:
                 for b in self._get_bond_potentials():
                     fh.write("tabulatedBondFile %s\n" % b)
@@ -1122,6 +1147,8 @@ tabulatedPotential  1
                 fh.write("inputExcludes %s\n" % self._exclusion_filename)
             if len(bond_angles) > 0:
                 fh.write("inputBondAngles %s\n" % self._bond_angle_filename)
+            if len(prod_pots) > 0:
+                fh.write("inputProductPotentials %s\n" % self._product_potential_filename)
      
         write_null_dx = False
         for pt,num in self.getParticleTypesAndCounts():
@@ -1242,6 +1269,25 @@ component "data" value 3
                     item = tuple([p.idx for p in b[:-1]] + [str(p) for p in b[-1]])
                     fh.write("BONDANGLE %d %d %d %d %s %s %s\n" % item)
 
+    def _writeArbdProductPotentialFile( self ):
+        if len(self.product_potentials) > 0:
+            with open(self._product_potential_filename,'w') as fh:
+                for ijk_tb in self.get_product_potentials():
+                    ijk = ijk_b[:-1]
+                    tb = ijk_tb[-1]
+                    if not (type(tb) is tuple or type(tb) is list):
+                        if len(tb) != 2: raise ValueError("Invalid product potential")
+                        type_,b = tb
+                        if type(type_) is not str: raise ValueError("Invalid product potential: unrecognized specification of potential type")
+                    else:
+                        type_ = ""
+                        b = tb
+
+                    if type(b) is not str and not isinstance(b, Path):
+                        b.write_file()
+                    fh.write("PRODUCTPOTENTIAL "+" ".join([str(x) for x in ijk+[type_,b] if x != ""]))
+
+
     def dimensions_from_structure( self, padding_factor=1.5, isotropic=False ):
         ## TODO: cache coordinates using numpy arrays for quick min/max
         raise(NotImplementedError)