From 3973fd7286e62b53e3f20b76a7026ef0b29c00ce Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Tue, 11 Jul 2023 12:23:49 -0500
Subject: [PATCH] Create RigidBody class
---
arbdmodel/__init__.py | 89 ++++++++++++++++++++++++++++++-----
arbdmodel/interactions.py | 99 ++++++++++++++++++++++-----------------
2 files changed, 133 insertions(+), 55 deletions(-)
diff --git a/arbdmodel/__init__.py b/arbdmodel/__init__.py
index c899cbf..478d7b5 100644
--- a/arbdmodel/__init__.py
+++ b/arbdmodel/__init__.py
@@ -11,6 +11,7 @@ from inspect import ismethod
import os, sys, subprocess
import shutil
+
_RESOURCE_DIR = Path(__file__).parent / 'resources'
def get_resource_path(relative_path):
return _RESOURCE_DIR / relative_path
@@ -25,6 +26,8 @@ def _get_properties_and_dict_keys(obj):
properties = [name for name,type_ in filter(filter_props, inspect.getmembers(obj.__class__))]
return properties + list(obj.__dict__.keys())
+__active_model = None # variable to be set temporarily by models as they are doing something, used by AbstractPotential classes
+
## Abstract classes
class Transformable():
def __init__(self, position, orientation=None):
@@ -88,7 +91,6 @@ class Parent():
self.group_sites = []
self.rigid = False
-
## TODO: self.cacheInvalid = True # What will be in the cache?
def add(self,x):
@@ -398,6 +400,7 @@ class ParticleType():
excludedAttributes = ("idx","type_",
"position",
+ "orientation",
"children",
"parent", "excludedAttributes",
)
@@ -417,9 +420,9 @@ class ParticleType():
self.rigid_body_potential = rigid_body_potential
for key in ParticleType.excludedAttributes:
- assert( key not in kargs )
+ assert( key not in kwargs )
- for key,val in kargs.items():
+ for key,val in kwargs.items():
self.__dict__[key] = val
def is_same_type(self, other):
@@ -501,6 +504,7 @@ class PointParticle(Transformable, Child):
self.name = name
self.counter = 0
self.restraints = []
+ self.rigid = False
for key,val in kwargs.items():
self.__dict__[key] = val
@@ -578,6 +582,61 @@ class PointParticle(Transformable, Child):
)
return data
+class RigidBody(PointParticle):
+
+ def __init__(self, type_, position, orientation, name="A", attached_particles=tuple(), **kwargs):
+ parent = None
+ if 'parent' in kwargs:
+ parent = kwargs['parent']
+ Child.__init__(self, parent=parent)
+ Transformable.__init__(self,position, orientation)
+
+ self.type_ = type_
+ self.idx = None
+ self.name = name
+ self.counter = 0
+ self.restraints = []
+ self.attached_particles = []
+ for p in attached_particles:
+ self.attach_particle(p)
+
+ for key,val in kwargs.items():
+ self.__dict__[key] = val
+
+ def add_restraint(self, restraint):
+ raise NotImplementedError('Harmonic restraints are not yet supported for rigid bodies')
+ ## TODO: how to handle duplicating and cloning bonds
+ self.restraints.append( restraint )
+
+ def get_restraints(self):
+ return [(self,r) for r in self.restraints]
+
+ def duplicate(self):
+ new = deepcopy(self)
+ return new
+
+ def attach_particle(self, particle):
+ """ The particle argument can be a PointParticle or Group (RigidBody children will be ignored). The position/orientation of the attached particle/group is in the RigidBody frame. """
+ if particle.parent is not None:
+ raise ValueError('RigidBody-attached particles are not allowed to have a parent')
+ self.attached_particles.append( particle )
+
+ def __getattr__(self, name):
+ """
+ First try to get attribute from the parent, then type_
+
+ Note that this data structure seems to be fragile, can result in stack overflow
+
+ """
+ # return Child.__getattr__(self,name)
+ try:
+ return Child.__getattr__(self,name)
+ except Exception as e:
+ if 'type_' in self.__dict__:
+ return getattr(self.type_, name)
+ else:
+ raise AttributeError(r"'{type(self).__name__}' object has no attribute '{name}'")
+
class Group(Transformable, Parent, Child):
@@ -821,13 +880,12 @@ class ArbdModel(PdbModel):
def __init__(self, children, origin=None, dimensions=(1000,1000,1000),
- remove_duplicate_bonded_terms=True, nonbonded_resolution=0.1,
+ remove_duplicate_bonded_terms=True,
configuration=None, dummy_types=tuple(), **conf_params):
PdbModel.__init__(self, children, dimensions, remove_duplicate_bonded_terms)
self.origin = origin
- self.nonbonded_resolution = nonbonded_resolution
if configuration is None:
configuration = SimConf(**conf_params)
self.configuration = configuration
@@ -890,7 +948,7 @@ class ArbdModel(PdbModel):
def _updateParticleOrder(self):
## Create ordered list
- self.particles = [p for p in self]
+ self.particles = [p for p in self if not p.rigid]
# self.particles = sorted(particles, key=lambda p: (p.type_, p.idx))
## Update particle indices
@@ -905,12 +963,12 @@ class ArbdModel(PdbModel):
def useNonbondedScheme(self, nbScheme, typeA=None, typeB=None):
""" deprecated """
- self.add_nonbonded_scheme(nbScheme, typeA, typeB)
+ self.add_nonbonded_interaction(nbScheme, typeA, typeB)
- def add_nonbonded_interaction(self, nonbonded_scheme, typeA=None, typeB=None):
- self.nonbonded_interactions.append( (nonbonded_scheme, typeA, typeB) )
+ def add_nonbonded_interaction(self, nonbonded_potential, typeA=None, typeB=None):
+ self.nonbonded_interactions.append( (nonbonded_potential, typeA, typeB) )
if typeA != typeB:
- self.nonbonded_interactions.append( (nonbonded_scheme, typeB, typeA) )
+ self.nonbonded_interactions.append( (nonbonded_potential, typeB, typeA) )
def prepare_for_simulation(self):
...
@@ -1284,11 +1342,16 @@ class ArbdEngine(SimEngine):
if configuration is None:
configuration = self._get_combined_conf(model, **conf_params)
- x = np.arange(0, configuration.cutoff, model.nonbonded_resolution)
+ x = np.arange(0, configuration.cutoff)
for i,j,t1,t2 in model._particleTypePairIter():
f = "%s.%s-%s.dat" % (prefix, t1.name, t2.name)
- scheme = model._get_nonbonded_interaction(t1,t2)
- scheme.write_file(f, t1, t2, rMax = configuration.cutoff)
+ interaction = model._get_nonbonded_interaction(t1,t2)
+ old_range = interaction.range_
+
+ interaction.range_ = [0, configuration.cutoff]
+ interaction.write_file(f, (t1, t2))
+ interaction.range_ = old_range
+
model._nonbonded_interaction_files.append(f)
def _write_restraint_file( self, model, filename ):
diff --git a/arbdmodel/interactions.py b/arbdmodel/interactions.py
index 69609d8..db3a7da 100644
--- a/arbdmodel/interactions.py
+++ b/arbdmodel/interactions.py
@@ -3,20 +3,28 @@ import os, sys
import numpy as np
from shutil import copyfile
-class NonbondedInteraction(metaclass=ABCMeta):
- """ Abstract class for writing nonbonded interactions """
+""" Module providing classes used to describe potentials in ARBD """
+
+
+""" Abstract classes """
+class AbstractPotential(metaclass=ABCMeta):
+ """ Abstract class for writing potentials """
- def __init__(self, typesA=None, typesB=None, resolution=0.1, rMin=0):
- """If typesA is None, and typesB is None, then """
+ def __init__(self, range_=(0,None), resolution=0.1,
+ max_force = None, max_potential = None):
self.resolution = resolution
- self.rMin = rMin
+ self.range_ = range_
- def add_sim_system(self, simSystem):
- self.rMax = simSystem.cutoff
- self.r = np.arange(rMin,rMax,resolution)
+ @property
+ def range_(self):
+ return self.__range_
+ @range_.setter
+ def range_(self,value):
+ assert(len(value) == 2)
+ self.__range_ = tuple(value)
@abstractmethod
- def potential(self, r, typeA, typeB):
+ def potential(self, r, types):
raise NotImplementedError
def __remove_nans(self, u):
@@ -26,16 +34,29 @@ class NonbondedInteraction(metaclass=ABCMeta):
right[np.isnan(u[right])] += 1
u[:-1][left] = u[right]
- def write_file(self, filename, typeA, typeB, rMax):
- r = np.arange(self.rMin, rMax+self.resolution, self.resolution)
+ def write_file(self, filename, types=None):
+ rmin,rmax = self.range_
+ r = np.arange(rmin, rmax+self.resolution, self.resolution)
with np.errstate(divide='ignore',invalid='ignore'):
- u = self.potential(r, typeA, typeB)
+ 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
+
+class BondedPotential(AbstractPotential):
+ """ Abstract class for writing bonded interactions """
-class LennardJones(NonbondedInteraction):
- def potential(self, r, typeA, typeB):
+ # def add_sim_system(self, simSystem):
+ # rmin,rmax = self.range_
+ # self.r = np.arange(rmin,rmax,resolution)
+
+""" Concrete nonbonded pontentials """
+class LennardJones(AbstractPotential):
+ def potential(self, r, types):
+ typeA, typeB = types
epsilon = np.sqrt( typeA.epsilon**2 + typeB.epsilon**2 )
r0 = 0.5 * (typeA.radius + typeB.radius)
r6 = (r0/r)**6
@@ -43,33 +64,31 @@ class LennardJones(NonbondedInteraction):
u = 4 * epsilon * (r12-r6)
# u[0] = u[1] # Remove NaN
return u
-# LennardJones = LennardJones()
-class HalfHarmonic(NonbondedInteraction):
- def potential(self, r, typeA, typeB):
+class HalfHarmonic(AbstractPotential):
+ def potential(self, r, types):
+ typeA, typeB = types
k = 10 # kcal/mol AA**2
r0 = (typeA.radius + typeB.radius)
u = 0.5 * k * (r-r0)**2
u[r > r0] = np.zeros( np.shape(u[r > r0]) )
return u
-# HalfHarmonic = HalfHarmonic()
-class TabulatedNonbonded(NonbondedInteraction):
- def __init__(self, tableFile, typesA=None, typesB=None, resolution=0.1, rMin=0):
- """If typesA is None, and typesB is None, then """
+class TabulatedNonbonded(AbstractPotential):
+ def __init__(self, tableFile, *args, **kwargs):
self.tableFile = tableFile
- # self.resolution = resolution
- # self.rMin = rMin
+ AbstractPotential.__init__(self,*args,**kwargs)
## TODO: check that tableFile exists and is regular file
- def write_file(self, filename, typeA, typeB, rMax):
+ def write_file(self, filename):
if filename != self.tableFile:
copyfile(self.tableFile, filename)
-## Bonded potentials
-class BasePotential():
- def __init__(self, r0, rRange=(0,50), resolution=0.1, maxForce=None, max_potential=None, zero="min", prefix="potentials/"):
+""" 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/"):
+ self.k = k
self.r0 = r0
self.rRange = rRange
self.resolution = resolution
@@ -128,11 +147,11 @@ class BasePotential():
np.savetxt( self.filename(), np.array([r, u]).T, fmt="%f" )
-class HarmonicPotential(BasePotential):
+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
- BasePotential.__init__(self, r0, rRange, resolution, maxForce, max_potential, zero, prefix)
+ 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_,
@@ -150,13 +169,8 @@ class HarmonicPotential(BasePotential):
return False
return True
-# class NonBonded(HarmonicPotential):
-# def _init_hook(self):
-# self.type = "nonbonded"
-# 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/", zero="min", correct_geometry=False, temperature=295):
+ 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
@@ -174,22 +188,23 @@ class HarmonicBond(HarmonicPotential):
class HarmonicAngle(HarmonicPotential):
- def __init__(self, k, r0, rRange=(0,181), resolution=0.1, maxForce=None, max_potential=None, zero="min", prefix="potentials/"):
+ 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
class HarmonicDihedral(HarmonicPotential):
- 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)
+ 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
-class WLCSKBond(BasePotential):
+class WLCSKBond(BaseBondedPotential):
""" ## 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/"):
- BasePotential.__init__(self, 0, rRange, resolution, maxForce, max_potential, zero, prefix)
+ ## 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"
self.d = d # separation
self.lp = lp # persistence length
@@ -224,10 +239,10 @@ class WLCSKBond(BasePotential):
return False
return True
-class WLCSKAngle(BasePotential):
+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/"):
- BasePotential.__init__(self, 180, rRange, resolution, maxForce, max_potential, zero, prefix)
+ BaseBondedPotential.__init__(self, 180, rRange, resolution, maxForce, max_potential, zero, prefix)
self.type_ = "wlcangle"
self.d = d # separation
self.lp = lp # persistence length
--
GitLab