diff --git a/mrdna/arbdmodel/__init__.py b/mrdna/arbdmodel/__init__.py
index 63d03973478890c63bd7e1c84400326a16edd439..51c26bd9d342a05df71d266bf5f90ff58790253e 100644
--- a/mrdna/arbdmodel/__init__.py
+++ b/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)
diff --git a/mrdna/arbdmodel/nonbonded.py b/mrdna/arbdmodel/nonbonded.py
index f7de00d677b87f554cb1e7bb276ec8a7edd099bb..a9339882acabb531042431c72774793dc2822b6d 100644
--- a/mrdna/arbdmodel/nonbonded.py
+++ b/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))
- u = u - np.min(u)
+ 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