From 0d46571b2c430a28debe50827bfd8db03fb1f387 Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Thu, 6 Aug 2020 18:12:40 -0500
Subject: [PATCH] Added support to arbdmodel for product potentials
---
mrdna/arbdmodel/__init__.py | 46 +++++++++++++++++++++++++++++++++++++
1 file changed, 46 insertions(+)
diff --git a/mrdna/arbdmodel/__init__.py b/mrdna/arbdmodel/__init__.py
index 0f97148..b4189f6 100644
--- a/mrdna/arbdmodel/__init__.py
+++ b/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)
--
GitLab