diff --git a/arbdmodel/__init__.py b/arbdmodel/__init__.py
index 478d7b5a87960e75d750b2077ea9e58a185f604f..e22e22a59ac9cb2e6484d7ee88789b172735081f 100644
--- a/arbdmodel/__init__.py
+++ b/arbdmodel/__init__.py
@@ -1,4 +1,28 @@
# -*- coding: utf-8 -*-
+
+## Set up loggers
+import logging
+def _get_username():
+ import sys
+ try:
+ return sys.environ['USER']
+ except:
+ return None
+
+logging.basicConfig(format='%(name)s: %(levelname)s: %(message)s')
+logger = logging.getLogger(__name__)
+logger.setLevel(logging.INFO)
+_ch = logging.StreamHandler()
+_ch.setFormatter(logging.Formatter('%(name)s: %(levelname)s: %(message)s'))
+logger.addHandler(_ch)
+logger.propagate = False
+
+devlogger = logging.getLogger(__name__+'.dev')
+# devlogger.setLevel(logging.DEBUG)
+if _get_username() not in ('cmaffeo2',):
+ devlogger.addHandler(logging.NullHandler())
+
+## Import resources
from pathlib import Path
from abc import abstractmethod, ABCMeta
@@ -1138,9 +1162,9 @@ class SimEngine(metaclass=ABCMeta):
cmd = self._generate_command_string(binary, output_name, output_directory, num_procs, gpu, replicas)
if dry_run:
- print("Run with: %s" % " ".join(cmd))
+ logger.info(f'Run with: {" ".join(cmd)}')
else:
- print("Running {} with: {}".format(binary," ".join(cmd)))
+ logger.info(f'Running {self.default_binary} with: {" ".join(cmd)}')
if log_file is None or (hasattr(log_file,'write') and callable(log_file.write)):
fd = sys.stdout if log_file is None else log_file
process = subprocess.Popen(cmd, stdout=subprocess.PIPE, universal_newlines=True)
@@ -1346,6 +1370,7 @@ class ArbdEngine(SimEngine):
for i,j,t1,t2 in model._particleTypePairIter():
f = "%s.%s-%s.dat" % (prefix, t1.name, t2.name)
interaction = model._get_nonbonded_interaction(t1,t2)
+ if interaction is None: continue
old_range = interaction.range_
interaction.range_ = [0, configuration.cutoff]
@@ -1377,7 +1402,11 @@ class ArbdEngine(SimEngine):
with open(self._bond_filename,'w') as fh:
for i,j,b,ex in model.get_bonds():
- item = (i.idx, j.idx, str(b))
+ try:
+ bfile = b.filename()
+ except:
+ bfile = str(b)
+ item = (i.idx, j.idx, bfile)
if ex:
fh.write("BOND REPLACE %d %d %s\n" % item)
else:
@@ -1391,7 +1420,11 @@ class ArbdEngine(SimEngine):
with open(self._angle_filename,'w') as fh:
for b in model.get_angles():
- item = tuple([p.idx for p in b[:-1]] + [str(b[-1])])
+ try:
+ bfile = b[-1].filename()
+ except:
+ bfile = str(b[-1])
+ item = tuple([p.idx for p in b[:-1]] + [bfile])
fh.write("ANGLE %d %d %d %s\n" % item)
def _write_dihedral_file( self, model, filename ):
@@ -1402,7 +1435,11 @@ class ArbdEngine(SimEngine):
with open(self._dihedral_filename,'w') as fh:
for b in model.get_dihedrals():
- item = tuple([p.idx for p in b[:-1]] + [str(b[-1])])
+ try:
+ bfile = b[-1].filename()
+ except:
+ bfile = str(b[-1])
+ item = tuple([p.idx for p in b[:-1]] + [bfile])
fh.write("DIHEDRAL %d %d %d %d %s\n" % item)
def _write_exclusion_file( self, model, filename ):
@@ -1417,7 +1454,14 @@ class ArbdEngine(SimEngine):
if len(model.bond_angles) > 0:
with open(self._bond_angle_filename,'w') as fh:
for b in model.get_bond_angles():
- item = tuple([p.idx for p in b[:-1]] + [str(p) for p in b[-1]])
+ bfiles = []
+ for p in b[-1]:
+ try:
+ bfile = p.filename()
+ except:
+ bfile = str(p)
+ bfiles.append(bfile)
+ item = tuple([p.idx for p in b[:-1]] + bfiles)
fh.write("BONDANGLE %d %d %d %d %s %s %s\n" % item)
def _write_product_potential_file( self, model, filename ):
@@ -1438,8 +1482,12 @@ class ArbdEngine(SimEngine):
b = tb
if type(b) is not str and not isinstance(b, Path):
b.write_file()
+ try:
+ bfile = b.filename()
+ except:
+ bfile = str(b)
line = line+" ".join([str(x.idx) for x in ijk])+" "
- line = line+" ".join([str(x) for x in [type_,b] if x != ""])+" "
+ line = line+" ".join([str(x) for x in [type_,bfile] if x != ""])+" "
fh.write(line)
def _write_group_sites_file( self, model, filename ):
@@ -1638,15 +1686,27 @@ tabulatedPotential 1
if len(bonds) > 0:
for b in model._get_bond_potentials():
- fh.write("tabulatedBondFile %s\n" % b)
+ try:
+ bfile = b.filename()
+ except:
+ bfile = str(b)
+ fh.write("tabulatedBondFile %s\n" % bfile)
if len(angles) > 0:
for b in model._get_angle_potentials():
- fh.write("tabulatedAngleFile %s\n" % b)
+ try:
+ bfile = b.filename()
+ except:
+ bfile = str(b)
+ fh.write("tabulatedAngleFile %s\n" % bfile)
if len(dihedrals) > 0:
for b in list(set([b for i,j,k,l,b in dihedrals])):
- fh.write("tabulatedDihedralFile %s\n" % b)
+ try:
+ bfile = b.filename()
+ except:
+ bfile = str(b)
+ fh.write("tabulatedDihedralFile %s\n" % bfile)
if len(restraints) > 0:
fh.write("inputRestraints %s\n" % self._restraint_filename)
diff --git a/arbdmodel/coords.py b/arbdmodel/coords.py
index 52b8f899186a77c00fce6e810b6571193d2bd66c..34e704685f49b0d3c4b30ed2826f5fc2cea82cbe 100644
--- a/arbdmodel/coords.py
+++ b/arbdmodel/coords.py
@@ -1,3 +1,4 @@
+from . import logger
import numpy as np
from scipy.optimize import newton
@@ -157,13 +158,21 @@ def rotationAboutAxis(axis,angle, normalizeAxis=True):
return quaternion_to_matrix(q)
def readArbdCoords(fname):
+ logger.warning('readArbdCoords is deprecated. Please update your code to use read_arbd_coordinates')
+ return read_arbd_coordinates(fname)
+
+def read_arbd_coordinates(fname):
coords = []
with open(fname) as fh:
for line in fh:
coords.append([float(x) for x in line.split()[1:]])
return np.array(coords)
-def readAvgArbdCoords(psf,pdb,dcd,rmsdThreshold=3.5, first_frame=0, last_frame=-1, stride=1):
+def readAvgArbdCoords(*args,**kwargs):
+ logger.warning('readAvgArbdCoords is deprecated. Please update your code to use read_average_arbd_coordinates')
+ return read_average_arbd_coordinates(*args,**kwargs)
+
+def read_average_arbd_coordinates(psf,pdb,dcd,rmsd_threshold=3.5, first_frame=0, last_frame=-1, stride=1):
import MDAnalysis as mda
usel = mda.Universe(psf, dcd)
@@ -180,11 +189,11 @@ def readAvgArbdCoords(psf,pdb,dcd,rmsdThreshold=3.5, first_frame=0, last_frame=-
rmsd = np.mean( (r-r0)**2 )
r = np.array( [(r-comA).dot(R)+comB for r in usel.atoms.positions] )
pos.append( r )
- if rmsd > rmsdThreshold**2:
+ if rmsd > rmsd_threshold**2:
break
t0=t+1
- print( "Averaging coordinates in %s after frame %d" % (dcd, t0) )
-
+ logger.info("Averaging coordinates in %s after frame %d" % (dcd, t0) )
+
pos = np.mean(pos, axis=0)
return pos
diff --git a/arbdmodel/fjc_polymer_model.py b/arbdmodel/fjc_polymer_model.py
index ee63f7f745380ec55270cfed374ee29a609eb353..abbff5d2d1b08b9d17cf188ede2d9f8a634dc00a 100644
--- a/arbdmodel/fjc_polymer_model.py
+++ b/arbdmodel/fjc_polymer_model.py
@@ -4,11 +4,10 @@
import numpy as np
import sys
-
## Local imports
-from . import ArbdModel, ParticleType, PointParticle, Group, get_resource_path
-from .polymer import PolymerSection, PolymerGroup
-from .interactions import NonbondedScheme, HarmonicBond, HarmonicPotential
+from . import logger, ArbdModel, ParticleType, PointParticle, Group, get_resource_path
+from .polymer import PolymerSection, PolymerGroup, PolymerBeads, PolymerModel
+from .interactions import AbstractPotential, HarmonicBond
from .coords import quaternion_to_matrix
@@ -18,81 +17,48 @@ type_ = ParticleType("X",
mass=120
)
-## Bonded potentials
-class FjcNonbonded(NonbondedScheme):
- def __init__(self, resolution=0.1, rMin=0):
- NonbondedScheme.__init__(self, typesA=None, typesB=None, resolution=resolution, rMin=rMin)
+# ## Bonded potentials
+# class FjcNonbonded(AbstractPotential):
+# """ This potential should apply zero force; however, it is required for the arbd engine """
+# def __init__(self, resolution=0.1, range_=(0,1)):
+# AbstractPotential.__init__(self, resolution=resolution, range_=range_)
- def potential(self, r, typeA, typeB):
- """ Constant force excluded volume """
- force = 10 # kcal_mol/AA
- radius = 6
-
- u = np.zeros(r.shape)
- # s = r < 2*radius
- # u[s] = (2*radius - r[s]) * force
- return u
+# def potential(self, r, types):
+# u = np.zeros(r.shape)
+# return u
-class FjcBeadsFromPolymer(Group):
+class FjcBeadsFromPolymer(PolymerBeads):
def __init__(self, polymer, sequence=None,
- rest_length = 3.8, spring_constant = 25,
+ rest_length = None, spring_constant = 25, monomers_per_bead_group = 1,
**kwargs):
- # if sequence is None:
- # raise NotImplementedError
- # # ... set random sequence
-
- self.polymer = polymer
- self.sequence = sequence
- self.spring_constant = 25
- self.rest_length = 3.8
-
- for prop in ('segname','chain'):
- if prop not in kwargs:
- # import pdb
- # pdb.set_trace()
- try:
- self.__dict__[prop] = polymer.__dict__[prop]
- except:
- pass
-
- # if len(sequence) != polymer.num_monomers:
- # raise ValueError("Length of sequence does not match length of polymer")
- Group.__init__(self, **kwargs)
+ self.spring_constant = spring_constant
+ PolymerBeads.__init__(self, polymer, sequence, monomers_per_bead_group, rest_length, **kwargs)
+ def _generate_ith_bead_group(self, i, r, o):
+ return PointParticle(type_, r,
+ resid = i+1)
- def _clear_beads(self):
- ...
-
- def _generate_beads(self):
- # beads = self.children
- nb = self.polymer.num_monomers
-
- for i in range(nb):
- c = self.polymer.monomer_index_to_contour(i)
- r = self.polymer.contour_to_position(c)
-
- bead = PointParticle(type_, r,
- resid = i+1)
- self.add(bead)
-
- ## Two consecutive nts
- for i in range(len(self.children)-1):
- b1 = self.children[i]
- b2 = self.children[i+1]
+ def _join_adjacent_bead_groups(self, ids):
+ if len(ids) == 2:
+ b1,b2 = [self.children[i] for i in ids]
""" units "10 kJ/N_A" kcal_mol """
bond = HarmonicBond(k = self.spring_constant,
r0 = self.rest_length,
- rRange = (0,500),
+ range_ = (0,500),
resolution = 0.01,
- maxForce = 50)
+ max_force = 50)
+ logger.info(f'Adding bond to {b1} {b2} with k={self.spring_constant} and r0={self.rest_length}')
self.add_bond( i=b1, j=b2, bond = bond, exclude=True )
+ else:
+ pass
-class FjcModel(ArbdModel):
+class FjcModel(PolymerModel):
def __init__(self, polymers,
sequences = None,
- rest_length = 3.8,
+ rest_length = None,
+ monomers_per_bead_group = 1,
spring_constant = 25,
damping_coefficient = 40.9,
DEBUG=False,
@@ -101,55 +67,30 @@ class FjcModel(ArbdModel):
"""
[damping_coefficient]: ns
"""
-
- print("WARNING: diffusion coefficient arbitrarily set to 100 AA**2/ns in FjcModel")
- kwargs['timestep'] = 50e-6
- kwargs['cutoff'] = 10
-
- if 'decompPeriod' not in kwargs:
- kwargs['decompPeriod'] = 100000
-
- """ Assign sequences """
- if sequences is None:
- # raise NotImplementedError("HpsModel must be provided a sequences argument")
- sequences = [None for i in range(len(polymers))]
+ if 'timestep' not in kwargs: kwargs['timestep'] = 50e-6
+ if 'cutoff' not in kwargs: kwargs['cutoff'] = 5 # no non-bonded interactions, so set this small
+ if 'decomp_period' not in kwargs:
+ kwargs['decomp_period'] = 100000 # Make it huge because we never need to compute FJC NB forces
- self.polymer_group = PolymerGroup(polymers)
- self.sequences = sequences
self.rest_length = rest_length
self.spring_constant = spring_constant
- ArbdModel.__init__(self, [], **kwargs)
+ PolymerModel.__init__(self, polymers, sequences, monomers_per_bead_group, **kwargs)
""" Update type diffusion coefficients """
+ logger.warning("Diffusion coefficient arbitrarily set to 100 AA**2/ns in FjcModel")
self.set_damping_coefficient( damping_coefficient )
""" Set up nonbonded interactions """
- nonbonded = FjcNonbonded()
- self.useNonbondedScheme( nonbonded, typeA=type_, typeB=type_ )
-
- """ Generate beads """
- self.generate_beads()
-
- def update_splines(self, coords):
- i = 0
- for p,c in zip(self.polymer_group.polymers,self.children):
- n = len(c.children)
- p.set_splines(np.linspace(0,1,n), coords[i:i+n])
- i += n
-
- self.clear_all()
- self.generate_beads()
- ## TODO Apply restraints, etc
-
- def generate_beads(self):
- self.peptides = [FjcBeadsFromPolymer(p,s, rest_length = self.rest_length,
- spring_constant = self.spring_constant )
- for p,s in zip(self.polymer_group.polymers,self.sequences)]
-
- for s in self.peptides:
- self.add(s)
- s._generate_beads()
+ # nonbonded = FjcNonbonded()
+ # self.add_nonbonded_interaction( nonbonded, typeA=type_, typeB=type_ )
+
+ def _generate_polymer_beads(self, polymer, sequence):
+ return FjcBeadsFromPolymer(polymer, sequence,
+ rest_length = self.rest_length,
+ spring_constant = self.spring_constant,
+ monomers_per_bead_group = self.monomers_per_bead_group,
+ )
def set_damping_coefficient(self, damping_coefficient):
for t in [type_]:
diff --git a/arbdmodel/hps_polymer_model.py b/arbdmodel/hps_polymer_model.py
index 4d4dbfccd03a5b84e0792cc3464212958ccbee75..e3ff157270ef85dd274550d21e42795eceef2c66 100644
--- a/arbdmodel/hps_polymer_model.py
+++ b/arbdmodel/hps_polymer_model.py
@@ -8,7 +8,7 @@ import sys
## Local imports
from . import ArbdModel, ParticleType, PointParticle, Group, get_resource_path
from .polymer import PolymerSection, PolymerGroup
-from .interactions import NonbondedScheme, HarmonicBond, HarmonicAngle, HarmonicDihedral
+from .interactions import AbstractPotential, HarmonicBond, HarmonicAngle, HarmonicDihedral
from .coords import quaternion_to_matrix
"""Define particle types"""
@@ -137,14 +137,15 @@ _types = dict(
for k,t in _types.items():
t.resname = t.name
-class HpsNonbonded(NonbondedScheme):
- def __init__(self, debye_length=10, resolution=0.1, rMin=0):
- NonbondedScheme.__init__(self, typesA=None, typesB=None, resolution=resolution, rMin=rMin)
+class HpsNonbonded(AbstractPotential):
+ def __init__(self, debye_length=10, resolution=0.1, range_=(0,None)):
+ AbstractPotential.__init__(self, resolution=resolution, range_=range_)
self.debye_length = debye_length
self.maxForce = 50
- def potential(self, r, typeA, typeB):
+ def potential(self, r, types):
""" Electrostatics """
+ typeA,typeB = types
ld = self.debye_length
q1 = typeA.charge
q2 = typeB.charge
diff --git a/arbdmodel/interactions.py b/arbdmodel/interactions.py
index db3a7da3586c4d273c8eae7644687fd01f7b3a8b..2017b4d2b5f14ebdfde4ff40c6ea6f0316a3ee73 100644
--- a/arbdmodel/interactions.py
+++ b/arbdmodel/interactions.py
@@ -6,14 +6,17 @@ from shutil import copyfile
""" Module providing classes used to describe potentials in ARBD """
-""" Abstract classes """
+## Abstract classes that others inherit from
class AbstractPotential(metaclass=ABCMeta):
""" Abstract class for writing potentials """
def __init__(self, range_=(0,None), resolution=0.1,
- max_force = None, max_potential = None):
- self.resolution = resolution
+ max_force = None, max_potential = None, zero='last'):
self.range_ = range_
+ self.resolution = resolution
+ self.max_force = max_force
+ self.max_potential = max_potential
+ self.zero = zero
@property
def range_(self):
@@ -23,8 +26,12 @@ class AbstractPotential(metaclass=ABCMeta):
assert(len(value) == 2)
self.__range_ = tuple(value)
+ @property
+ def periodic(self):
+ return False
+
@abstractmethod
- def potential(self, r, types):
+ def potential(self, r, types=None):
raise NotImplementedError
def __remove_nans(self, u):
@@ -34,26 +41,60 @@ class AbstractPotential(metaclass=ABCMeta):
right[np.isnan(u[right])] += 1
u[:-1][left] = u[right]
- def write_file(self, filename, types=None):
+ def filename(self, types=None):
+ raise NotImplementedError('Inherited potential objects should overload this function')
+
+ def write_file(self, filename=None, types=None):
+ if filename is None:
+ filename = self.filename(types)
rmin,rmax = self.range_
r = np.arange(rmin, rmax+self.resolution, self.resolution)
with np.errstate(divide='ignore',invalid='ignore'):
u = self.potential(r, types)
self.__remove_nans(u)
- np.savetxt(filename, np.array([r,u]).T)
-class NonbondedPotential(AbstractPotential):
- """ Abstract class for writing nonbonded interactions """
- ## No specialization needed so far
+ if self.zero == 'min':
+ u = u - np.min(u)
+ elif self.zero == 'first':
+ u = u - u[0]
+ elif self.zero == 'last':
+ u = u - u[-1]
+ else:
+ raise ValueError('Unrecognized option for "zero" argument')
+
+ max_force = self.max_force
+ if max_force is not None:
+ assert(max_force > 0)
+ f = np.diff(u)/np.diff(r)
+ f[f>max_force] = max_force
+ f[f<-max_force] = -max_force
+ u[0] = 0
+ u[1:] = np.cumsum(f*np.diff(r))
+
+ if self.max_potential is not None:
+ f = np.diff(u)/np.diff(r)
+ ids = np.where( 0.5*(u[1:]+u[:-1]) > self.max_potential )[0]
-class BondedPotential(AbstractPotential):
- """ Abstract class for writing bonded interactions """
+ # w = np.sqrt(2*self.max_potential/self.k)
+ # drAvg = 0.5*(np.abs(dr[ids]) + np.abs(dr[ids+1]))
+ # f[ids] = f[ids] * np.exp(-(drAvg-w)/(w))
+ f[ids] = self.max_potential
+
+ u[0] = 0
+ u[1:] = np.cumsum(f*np.diff(r))
- # def add_sim_system(self, simSystem):
- # rmin,rmax = self.range_
- # self.r = np.arange(rmin,rmax,resolution)
+ if self.zero == 'min':
+ u = u - np.min(u)
+ elif self.zero == 'first':
+ u = u - u[0]
+ elif self.zero == 'last':
+ u = u - u[-1]
+ else:
+ raise ValueError('Unrecognized option for "zero" argument')
-""" Concrete nonbonded pontentials """
+ np.savetxt(filename, np.array([r,u]).T)
+
+## Concrete nonbonded pontentials
class LennardJones(AbstractPotential):
def potential(self, r, types):
typeA, typeB = types
@@ -81,139 +122,143 @@ class TabulatedNonbonded(AbstractPotential):
## TODO: check that tableFile exists and is regular file
- def write_file(self, filename):
+ def write_file(self, filename, types):
if filename != self.tableFile:
copyfile(self.tableFile, filename)
-""" Concrete bonded potentials """
-class BaseBondedPotential():
- def __init__(self, k, r0, rRange=(0,50), resolution=0.1, maxForce=None, max_potential=None, zero='min', prefix="potentials/"):
+## Bonded potentials
+class HarmonicBondedPotential(AbstractPotential):
+ def __init__(self, k, r0, filename_prefix='./potentials/', *args, **kwargs):
self.k = k
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
- self.kscale_ = None # only used for
-
- def filename(self):
- raise NotImplementedError("Not implemented")
-
- def potential(self,dr):
- raise NotImplementedError("Not implemented")
-
- def __str__(self):
- return self.filename()
+ self.kscale_ = None
+ self.filename_prefix = filename_prefix
+ if 'zero' not in kwargs: kwargs['zero'] = 'min'
+ AbstractPotential.__init__(self, *args, **kwargs)
- def write_file(self):
- r = np.arange( self.rRange[0],
- self.rRange[1]+self.resolution,
- self.resolution )
+ @property
+ def kscale(self):
+ return 1.0
+
+ @property
+ @abstractmethod
+ def type_(self):
+ return "none"
+
+ def filename(self, types=None):
+ assert(types is None)
+ return "%s%s-%.3f-%.3f.dat" % (self.filename_prefix, self.type_,
+ self.k*self.kscale, self.r0)
+
+ def potential(self, r, types=None):
dr = r-self.r0
-
if self.periodic == True:
- rSpan = self.rRange[1]-self.rRange[0]
- assert(rSpan > 0)
- dr = np.mod( dr+0.5*rSpan, rSpan) - 0.5*rSpan
-
- u = self.potential(dr)
-
- maxForce = self.maxForce
- if maxForce is not None:
- assert(maxForce > 0)
- f = np.diff(u)/np.diff(r)
- f[f>maxForce] = maxForce
- f[f<-maxForce] = -maxForce
- u[0] = 0
- u[1:] = np.cumsum(f*np.diff(r))
-
- if self.max_potential is not None:
- f = np.diff(u)/np.diff(r)
- ids = np.where( 0.5*(u[1:]+u[:-1]) > self.max_potential )[0]
-
- w = np.sqrt(2*self.max_potential/self.k)
- drAvg = 0.5*(np.abs(dr[ids]) + np.abs(dr[ids+1]))
-
- f[ids] = f[ids] * np.exp(-(drAvg-w)/(w))
- 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(BaseBondedPotential):
- def __init__(self, k, r0, rRange, resolution, maxForce, max_potential, zero="min", correct_geometry=False, prefix='./'):
- self.k = k
- self.kscale_ = None
- BaseBondedPotential.__init__(self, k, r0, rRange, resolution, maxForce, max_potential, zero, prefix)
-
- def filename(self):
- return "%s%s-%.3f-%.3f.dat" % (self.prefix, self.type_,
- self.k*self.kscale_, self.r0)
- def potential(self,dr):
+ r_span = self.range_[1]-self.range_[0]
+ assert(r_span > 0)
+ dr = np.mod( dr+0.5*r_span, r_span) - 0.5*r_span
return 0.5*self.k*dr**2
def __hash__(self):
assert(self.type_ != "None")
- return hash((self.type_, self.k, self.r0, self.rRange, self.resolution, self.maxForce, self.max_potential, self.prefix, self.periodic))
+ return hash((self.type_, self.k, self.r0, self.range_, self.resolution, self.max_force, self.max_potential, self.zero, self.filename_prefix, self.periodic)) #
def __eq__(self, other):
- for a in ("type_", "k", "r0", "rRange", "resolution", "maxForce", "max_potential", "prefix", "periodic"):
- if self.__dict__[a] != other.__dict__[a]:
+ for a in ("type_", "k", "r0", "range_", "resolution", "max_force", "max_potential", "filename_prefix", "periodic"):
+ if getattr(self,a) != getattr(other,a):
return False
return True
-
-class HarmonicBond(HarmonicPotential):
- def __init__(self, k, r0, rRange=(0,50), resolution=0.1, maxForce=None, max_potential=None, zero='min', correct_geometry=False, temperature=295, prefix="potentials/"):
- 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
+
+class HarmonicBond(HarmonicBondedPotential):
+ def __init__(self, k, r0, correct_geometry=False, temperature=295, *args, **kwargs):
+ if 'range_' not in kwargs: kwargs['range_'] = (0,50)
+ HarmonicBondedPotential.__init__(self, k, r0, *args, **kwargs)
self.correct_geometry = correct_geometry
self.temperature = temperature
- def potential(self,dr):
- u = HarmonicPotential.potential(self,dr)
+ @property
+ def kscale(self):
+ return 1.0
+
+ @property
+ def type_(self):
+ return 'gbond' if self.correct_geometry else 'bond'
+
+ def potential(self, r, types=None):
+ u = HarmonicBondedPotential.potential(self, r, types)
if self.correct_geometry:
with np.errstate(divide='ignore',invalid='ignore'):
- du = 2*0.58622592*np.log(dr+self.r0) * self.temperature/295
+ du = 2*0.58622592*np.log(r) * self.temperature/295
du[np.logical_not(np.isfinite(du))] = 0
u = u+du
return u
+class HarmonicAngle(HarmonicBondedPotential):
+ def __init__(self, *args, **kwargs):
+ if 'range_' not in kwargs: kwargs['range_'] = (0,181)
+ HarmonicBondedPotential.__init__(self, *args, **kwargs)
+
+ @property
+ def kscale_(self):
+ return (180.0/np.pi)**2
+
+ @property
+ def type_(self):
+ return 'angle'
+
+class HarmonicDihedral(HarmonicBondedPotential):
+ def __init__(self, *args, **kwargs):
+ if 'range_' not in kwargs: kwargs['range_'] = (-180,180)
+ HarmonicBondedPotential.__init__(self, *args, **kwargs)
-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, zero=zero, prefix=prefix)
- self.type_ = "angle"
- self.kscale_ = (180.0/np.pi)**2
+ @property
+ def kscale_(self):
+ return (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, zero, prefix)
- self.periodic = True
- self.type_ = "dihedral"
- self.kscale_ = (180.0/np.pi)**2
+ @property
+ def type_(self):
+ return 'dihedral'
-class WLCSKBond(BaseBondedPotential):
+ @property
+ def periodic(self):
+ return True
+
+class WLCSKPotential(HarmonicBondedPotential, metaclass=ABCMeta):
""" ## 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, zero="min", prefix="potentials/"):
- ## Note, we set k to lp here, but it's not proper and may affect results from max_potential
- BaseBondedPotential.__init__(self, lp, 0, rRange, resolution, maxForce, max_potential, zero, prefix)
- self.type_ = "wlcbond"
+ def __init__(self, d, lp, kT, **kwargs):
+ ## Note, we're leveraging the HarmonicBondedPotential base class and set k to lp here, but it isn't proper
+ HarmonicBondedPotential.__init__(self, d, lp, **kwargs)
self.d = d # separation
self.lp = lp # persistence length
self.kT = kT
- def filename(self):
- return "%s%s-%.3f-%.3f.dat" % (self.prefix, self.type_,
+ def filename(self,types=None):
+ assert(types is None)
+ return "%s%s-%.3f-%.3f.dat" % (self.filename_prefix, self.type_,
self.d, self.lp)
- def potential(self, dr):
+
+ def __hash__(self):
+ return hash((self.d, self.lp, self.kT, HarmonicBondedPotential.__hash__(self)))
+
+ def __eq__(self, other):
+ for a in ("d", "lp", "kT"):
+ if self.__dict__[a] != other.__dict__[a]:
+ return False
+ return HarmonicBondedPotential.__eq__(self,other)
+
+class WLCSKBond(WLCSKPotential):
+ """ ## https://aip.scitation.org/doi/full/10.1063/1.4968020 """
+ def __init__(self, d, lp, kT, **kwargs):
+ if 'range_' not in kwargs: kwargs['range_'] = (0,50)
+ if 'resolution' not in kwargs: kwargs['resolution'] = 0.02
+ if 'max_force' not in kwargs: kwargs['max_force'] = 100
+ WLCSKPotential.__init__(self, d, lp, kT, **kwargs)
+
+ @property
+ def type_(self):
+ return "wlcbond"
+
+ def potential(self, r, types=None):
+ dr = r - self.d
nk = self.d / (2*self.lp)
q2 = (dr / self.d)**2
a1,a2 = 1, -7.0/(2*nk)
@@ -222,48 +267,23 @@ class WLCSKBond(BaseBondedPotential):
a4 = (p0 + p1/(2*nk) + p2*(2*nk)**-2) / (1+p3/(2*nk)+p4*(2*nk)**-2)
with np.errstate(divide='ignore',invalid='ignore'):
u = self.kT * nk * ( a1/(1-q2) - a2*np.log(1-q2) + a3*q2 - 0.5*a4*q2*(q2-2) )
- max_force = np.diff(u[q2<1][-2:]) / np.diff(dr).mean()
- max_u = u[q2<1][-1]
- max_dr = dr[q2<1][-2]
- assert( max_force >= 0 )
- u[q2>=1] = (dr[q2>=1]-max_dr)*max_force + max_u
return u
- def __hash__(self):
- assert(self.type_ != "None")
- return hash((self.type_, self.d, self.lp, self.kT, self.rRange, self.resolution, self.maxForce, self.max_potential, self.prefix, self.periodic))
-
- def __eq__(self, other):
- for a in ("type_", "d", "lp", "kT" "rRange", "resolution", "maxForce", "max_potential", "prefix", "periodic"):
- if self.__dict__[a] != other.__dict__[a]:
- return False
- return True
+class WLCSKAngle(WLCSKPotential):
+ """ ## https://aip.scitation.org/doi/full/10.1063/1.4968020 """
+ def __init__(self, d, lp, kT, **kwargs):
+ if 'range_' not in kwargs: kwargs['range_'] = (0,181)
+ if 'resolution' not in kwargs: kwargs['resolution'] = 0.5
+ WLCSKPotential.__init__(self, d, lp, kT, **kwargs)
-class WLCSKAngle(BaseBondedPotential):
- ## 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, zero="min", prefix="potentials/"):
- BaseBondedPotential.__init__(self, 180, rRange, resolution, maxForce, max_potential, zero, prefix)
- self.type_ = "wlcangle"
- self.d = d # separation
- self.lp = lp # persistence length
- self.kT = kT
+ @property
+ def type_(self):
+ return "wlcangle"
- def filename(self):
- return "%s%s-%.3f-%.3f.dat" % (self.prefix, self.type_,
- self.d, self.lp)
- def potential(self,dr):
+ def potential(self, r, types=None):
+ dr = r - self.d
nk = self.d / (2*self.lp)
p1,p2,p3,p4 = -1.237, 0.8105, -1.0243, 0.4595
C = (1 + p1*(2*nk) + p2*(2*nk)**2) / (2*nk+p3*(2*nk)**2+p4*(2*nk)**3)
u = self.kT * C * (1-np.cos(dr * np.pi / 180))
return u
-
- def __hash__(self):
- assert(self.type_ != "None")
- return hash((self.type_, self.d, self.lp, self.kT, self.rRange, self.resolution, self.maxForce, self.max_potential, self.prefix, self.periodic))
-
- def __eq__(self, other):
- for a in ("type_", "d", "lp", "kT" "rRange", "resolution", "maxForce", "max_potential", "prefix", "periodic"):
- if self.__dict__[a] != other.__dict__[a]:
- return False
- return True
diff --git a/arbdmodel/polymer.py b/arbdmodel/polymer.py
index c121dd24d8db184af6b8bcc46ce2a2c9b679ccb3..c6f0244a499ff4f78df8067288c4fc994d4bc9d6 100644
--- a/arbdmodel/polymer.py
+++ b/arbdmodel/polymer.py
@@ -1,14 +1,16 @@
+from abc import abstractmethod, ABCMeta
+
+from . import logger
+
from pathlib import Path
from copy import copy, deepcopy
import numpy as np
from scipy import interpolate
-from . import ArbdModel
+from . import ArbdModel, Group
from .coords import rotationAboutAxis, quaternion_from_matrix, quaternion_to_matrix
-
-
"""
TODO:
- test for large systems
@@ -158,7 +160,6 @@ class ConnectableElement():
l = B.parent.locations
if B not in l: l.append(B)
-
class PolymerSection(ConnectableElement):
""" Base class that describes a linear section of a polymer """
@@ -401,7 +402,7 @@ class PolymerSection(ConnectableElement):
for c in cl:
yield c
-
+
class PolymerGroup():
def __init__(self, polymers=[],
**kwargs):
@@ -546,3 +547,125 @@ class PolymerGroup():
## iterate through the model polymers
for s in self.polymers:
draw_cylinder(s,radius,"cyan")
+
+class PolymerBeads(Group, metaclass=ABCMeta):
+ """ Abstract class for bead-based models that are built from Polymer objects """
+
+ def __init__(self, polymer, sequence=None, monomers_per_bead_group = 1, rest_length = None, **kwargs):
+ self.polymer = polymer
+ self.sequence = sequence
+ # self.spring_constant = spring_constant
+ self.monomers_per_bead_group = monomers_per_bead_group
+
+ if rest_length is None:
+ rest_length = polymer.monomer_length * self.monomers_per_bead_group
+
+ self.rest_length = rest_length
+
+ for prop in ('segname','chain'):
+ if prop not in kwargs:
+ try:
+ self.__dict__[prop] = polymer.__dict__[prop]
+ except:
+ pass
+
+ Group.__init__(self, **kwargs)
+
+ @property
+ def monomers_per_bead_group(self):
+ return self.__monomers_per_bead_group
+
+ @property
+ def num_bead_groups(self):
+ return self.__num_bead_groups
+
+ @monomers_per_bead_group.setter
+ def monomers_per_bead_group(self,value):
+ if value <= 0:
+ raise ValueError("monomers_per_bead_group must be positive")
+ self.num_bead_groups = int(np.ceil(self.polymer.num_monomers/value))
+
+ @num_bead_groups.setter
+ def num_bead_groups(self,value):
+ if value <= 0:
+ raise ValueError("num_bead_groups must be positive")
+ self.__num_bead_groups = value
+ self.__monomers_per_bead_group = self.polymer.num_monomers/value # get exact ratio
+
+ def _clear_beads(self):
+ ## self.children = []
+ self.clear_all()
+
+ @abstractmethod
+ def _generate_ith_bead_group(self, i, r, o):
+ """ Normally a bead, but sometimes a group of beads """
+ bead = PointParticle(type_, r,
+ resid = i+1)
+ pass
+
+ @abstractmethod
+ def _join_adjacent_bead_groups(self, ids):
+ if len(ids) == 2:
+ b1,b2 = [self.children[i] for i in ids]
+ """ units "10 kJ/N_A" kcal_mol """
+ bond = HarmonicBond(k = self.spring_constant,
+ r0 = self.rest_length,
+ rRange = (0,500),
+ resolution = 0.01,
+ maxForce = 50)
+ self.add_bond( i=b1, j=b2, bond = bond, exclude=True )
+ else:
+ pass
+
+ def _generate_beads(self):
+ for i in range(self.num_bead_groups):
+ c = self.polymer.monomer_index_to_contour(i * self.monomers_per_bead_group)
+ r = self.polymer.contour_to_position(c)
+ o = self.polymer.contour_to_orientation(c)
+
+ obj = self._generate_ith_bead_group(i, r, o)
+ self.add(obj)
+
+ for di in range(1,4):
+ for i in range(len(self.children)-di):
+ ids = tuple(i+_di for _di in range(di+1))
+ self._join_adjacent_bead_groups(ids)
+
+class PolymerModel(ArbdModel, metaclass=ABCMeta):
+ def __init__(self, polymers,
+ sequences = None,
+ monomers_per_bead_group = 1,
+ **kwargs):
+
+ """ Assign sequences """
+ if sequences is None:
+ sequences = [None for i in range(len(polymers))]
+
+ self.polymer_group = PolymerGroup(polymers)
+ self.sequences = sequences
+ self.monomers_per_bead_group = monomers_per_bead_group
+ ArbdModel.__init__(self, [], **kwargs)
+
+ """ Generate beads """
+ self.generate_beads()
+
+ def update_splines(self, coords):
+ i = 0
+ for p,c in zip(self.polymer_group.polymers,self.children):
+ n = len(c.children)
+ p.set_splines(np.linspace(0,1,n), coords[i:i+n])
+ i += n
+
+ @abstractmethod
+ def _generate_polymer_beads(self, polymer, sequence):
+ pass
+
+ def generate_beads(self):
+ self.clear_all()
+ self.peptides = [self._generate_polymer_beads(p,s)
+ for p,s in zip(self.polymer_group.polymers,self.sequences)]
+
+ for s in self.peptides:
+ self.add(s)
+ s._generate_beads()
+