Skip to content
GitLab
Commit 79ff74f3 authored by cmaffeo2's avatar cmaffeo2
Browse files

Updated arbdmodel to commit 1578c98 from arbdmodel package 

https://gitlab.engr.illinois.edu/tbgl/tools/arbdmodel
parent abb67908
Expand all Hide whitespace changes
Inline Side-by-side
mrdna/arbdmodel/__init__.py
View file @ 79ff74f3
......@@ -9,6 +9,9 @@ from copy import copy, deepcopy
from inspect import ismethod
import os, sys, subprocess
_RESOURCE_DIR = Path(__file__).parent / 'resources'
def get_resource_path(relative_path):
return _RESOURCE_DIR / relative_path
## Abstract classes
class Transformable():
......@@ -69,8 +72,9 @@ class Parent():
self.impropers = []
self.exclusions = []
## TODO: self.cacheInvalid = True # What will be in the cache?
self.rigid = False
## TODO: self.cacheInvalid = True # What will be in the cache?
def add(self,x):
## TODO: check the parent-child tree to make sure there are no cycles
......@@ -82,6 +86,19 @@ class Parent():
x.parent = self
self.children.append(x)
def insert(self,idx,x):
## TODO: check the parent-child tree to make sure there are no cycles
if not isinstance(x,Child):
raise Exception('Attempted to add an object to a group that does not inherit from the "Child" type')
if x.parent is not None and x.parent is not self:
raise Exception("Child {} already belongs to some group".format(x))
x.parent = self
self.children.insert(idx,x)
def index(self, x):
return self.children.index(x)
def clear_all(self, keep_children=False):
if keep_children == False:
for x in self.children:
......@@ -96,7 +113,8 @@ class Parent():
def remove(self,x):
if x in self.children:
self.children.remove(x)
x.parent = None
if x.parent is self:
x.parent = None
def add_bond(self, i,j, bond, exclude=False):
## TODO: how to handle duplicating and cloning bonds
......@@ -235,11 +253,17 @@ class Child():
def __getattr__(self, name):
"""
Try to get attribute from the parent
"""
# if self.parent is not None:
if "parent" not in self.__dict__ or self.__dict__["parent"] is None or name is "children":
raise AttributeError("'{}' object has no attribute '{}'".format(type(self).__name__, name))
excluded_attributes = ['parent','rigid']
if name in excluded_attributes:
raise AttributeError("'{}' object has no attribute '{}' and cannot look it up from the parent".format(type(self).__name__, name))
## TODO: determine if there is a way to avoid __getattr__ if a method is being looked up
try:
ret = getattr(self.parent,name)
......@@ -307,25 +331,35 @@ class ParticleType():
excludedAttributes = ("idx","type_",
"position",
"children",
"parent", "excludedAttributes",
"parent", "excludedAttributes","rigid"
)
def __init__(self, name, charge=0, **kargs):
# parent = None
# if "parent" in kargs: parent = kargs["parent"]
# Child.__init__(self, parent)
def __init__(self, name, charge=0, parent=None, **kargs):
""" Parent type is used to fall back on for nonbonded interactions if this type is not specifically referenced """
self.name = name
self.charge = charge
self.parent = parent
for key in ParticleType.excludedAttributes:
assert( key not in kargs )
for key,val in kargs.items():
self.__dict__[key] = val
def is_same_type(self, other):
assert(isinstance(other,ParticleType))
if self == other:
return True
elif self.parent is not None and self.parent == other:
return True
else:
return False
def __hash_key(self):
l = [self.name,self.charge]
for keyval in sorted(self.__dict__.items()):
if isinstance(keyval[1], list): keyval = (keyval[0],tuple(keyval[1]))
l.extend(keyval)
return tuple(l)
......@@ -352,7 +386,10 @@ class ParticleType():
def __ge__(a,b):
a._equal_check(b)
return a.name >= b.name
def __repr__(self):
return '<{} {}{}>'.format( type(self), self.name, '[{}]'.format(self.parent) if self.parent is not None else '' )
class PointParticle(Transformable, Child):
def __init__(self, type_, position, name="A", **kwargs):
......@@ -378,6 +415,9 @@ class PointParticle(Transformable, Child):
def get_restraints(self):
return [(self,r) for r in self.restraints]
def duplicate(self):
new = deepcopy(self)
return new
def __getattr__(self, name):
"""
......@@ -401,6 +441,10 @@ class PointParticle(Transformable, Child):
segname = p.segname
except:
segname = "A"
try:
chain = p.chain
except:
chain = "A"
try:
resname = p.resname
except:
......@@ -409,6 +453,10 @@ class PointParticle(Transformable, Child):
resid = p.resid
except:
resid = p.idx+1
try:
mass = p.mass
except:
mass = 1
try:
occ = p.occupancy
......@@ -422,12 +470,12 @@ class PointParticle(Transformable, Child):
data = dict(segname = segname,
resname = resname,
name = str(p.name)[:4],
chain = "A",
chain = chain[0],
resid = int(resid),
idx = p.idx+1,
type = p.type_.name[:4],
type = p.type_.name[:7],
charge = p.charge,
mass = p.mass,
mass = mass,
occupancy = occ,
beta = beta
)
......@@ -622,15 +670,21 @@ class PdbModel(Transformable, Parent):
class ArbdModel(PdbModel):
def __init__(self, children, dimensions=(1000,1000,1000), temperature=291, timestep=50e-6, cutoff=50, decompPeriod=10000, pairlistDistance=None, nonbondedResolution=0.1, remove_duplicate_bonded_terms=True):
def __init__(self, children, dimensions=(1000,1000,1000), temperature=291, timestep=50e-6,
particle_integrator = 'Brown',
cutoff=50, decompPeriod=1000, pairlistDistance=None, nonbondedResolution=0.1,
remove_duplicate_bonded_terms=True):
PdbModel.__init__(self, children, dimensions, remove_duplicate_bonded_terms)
self.temperature = temperature
self.timestep = timestep
self.cutoff = cutoff
self.particle_integrator = particle_integrator
if pairlistDistance == None:
pairlistDistance = cutoff+2
pairlistDistance = cutoff+30
self.decompPeriod = decompPeriod
self.pairlistDistance = pairlistDistance
......@@ -661,11 +715,11 @@ class ArbdModel(PdbModel):
if A is None or B is None:
if A is None and B is None:
return s
elif A is None and typeB == B:
elif A is None and typeB.is_same_type(B):
return s
elif B is None and typeA == A:
elif B is None and typeA.is_same_type(A):
return s
elif typeA == A and typeB == B:
elif typeA.is_same_type(A) and typeB.is_same_type(B):
return s
raise Exception("No nonbonded scheme found for %s and %s" % (typeA.name, typeB.name))
......@@ -697,7 +751,7 @@ class ArbdModel(PdbModel):
if typeA != typeB:
self.nbSchemes.append( (nbScheme, typeB, typeA) )
def simulate(self, output_name, output_directory='output', num_steps=100000000, timestep=None, gpu=0, output_period=1e4, arbd=None, directory='.', replicas=1):
def simulate(self, output_name, output_directory='output', num_steps=100000000, timestep=None, gpu=0, output_period=1e4, arbd=None, directory='.', restart_file=None, replicas=1, dry_run = False):
assert(type(gpu) is int)
num_steps = int(num_steps)
......@@ -716,22 +770,26 @@ class ArbdModel(PdbModel):
if output_directory == '': output_directory='.'
if arbd is None:
for path in os.environ["PATH"].split(os.pathsep):
path = path.strip('"')
fname = os.path.join(path, "arbd")
if os.path.isfile(fname) and os.access(fname, os.X_OK):
arbd = fname
break
if dry_run:
if arbd is None: arbd='arbd'
else:
if arbd is None:
for path in os.environ["PATH"].split(os.pathsep):
path = path.strip('"')
fname = os.path.join(path, "arbd")
if os.path.isfile(fname) and os.access(fname, os.X_OK):
arbd = fname
break
if arbd is None: raise Exception("ARBD was not found")
if arbd is None: raise Exception("ARBD was not found")
if not os.path.exists(arbd):
raise Exception("ARBD was not found")
if not os.path.isfile(arbd):
raise Exception("ARBD was not found")
if not os.access(arbd, os.X_OK):
raise Exception("ARBD is not executable")
if not os.path.exists(arbd):
raise Exception("ARBD was not found")
if not os.path.isfile(arbd):
raise Exception("ARBD was not found")
if not os.access(arbd, os.X_OK):
raise Exception("ARBD is not executable")
if not os.path.exists(output_directory):
os.makedirs(output_directory)
......@@ -741,7 +799,7 @@ class ArbdModel(PdbModel):
self.writePdb( output_name + ".pdb" )
self.writePsf( output_name + ".psf" )
self.writeArbdFiles( output_name, numSteps=num_steps, outputPeriod=output_period )
self.writeArbdFiles( output_name, numSteps=num_steps, outputPeriod=output_period, restart_file=restart_file )
os.sync()
## http://stackoverflow.com/questions/18421757/live-output-from-subprocess-command
......@@ -752,11 +810,14 @@ class ArbdModel(PdbModel):
cmd = cmd + ["%s.bd" % output_name, "%s/%s" % (output_directory, output_name)]
cmd = tuple(str(x) for x in cmd)
print("Running ARBD with: %s" % " ".join(cmd))
process = subprocess.Popen(cmd, stdout=subprocess.PIPE, universal_newlines=True)
for line in process.stdout:
sys.stdout.write(line)
sys.stdout.flush()
if dry_run:
print("Run with: %s" % " ".join(cmd))
else:
print("Running ARBD with: %s" % " ".join(cmd))
process = subprocess.Popen(cmd, stdout=subprocess.PIPE, universal_newlines=True)
for line in process.stdout:
sys.stdout.write(line)
sys.stdout.flush()
except:
raise
......@@ -768,7 +829,7 @@ class ArbdModel(PdbModel):
# Methods for printing model #
# -------------------------- #
def writeArbdFiles(self, prefix, numSteps=100000000, outputPeriod=10000):
def writeArbdFiles(self, prefix, numSteps=100000000, outputPeriod=10000, restart_file=None):
## TODO: save and reference directories and prefixes using member data
d = self.potential_directory = "potentials"
if not os.path.exists(d):
......@@ -787,7 +848,7 @@ class ArbdModel(PdbModel):
self._writeArbdDihedralFile()
self._writeArbdExclusionFile()
self._writeArbdPotentialFiles( prefix, directory = d )
self._writeArbdConf( prefix, numSteps=numSteps, outputPeriod=outputPeriod )
self._writeArbdConf( prefix, numSteps=numSteps, outputPeriod=outputPeriod, restart_file=restart_file )
# def _writeArbdCoordFile(self, filename):
# with open(filename,'w') as fh:
......@@ -796,12 +857,49 @@ class ArbdModel(PdbModel):
def _writeArbdParticleFile(self, filename):
with open(filename,'w') as fh:
for p in self.particles:
data = tuple([p.idx,p.type_.name] + [x for x in p.collapsedPosition()])
fh.write("ATOM %d %s %f %f %f\n" % data)
def _writeArbdConf(self, prefix, randomSeed=None, numSteps=100000000, outputPeriod=10000, restartCoordinateFile=None):
if self.particle_integrator == "Brown":
for p in self.particles:
data = tuple([p.idx,p.type_.name] + [x for x in p.collapsedPosition()])
fh.write("ATOM %d %s %f %f %f\n" % data)
else:
for p in self.particles:
data = [p.idx,p.type_.name] + [x for x in p.collapsedPosition()]
try:
data = data + p.momentum
except:
try:
data = data + p.velocity*p.mass
except:
data = data + [0,0,0]
fh.write("ATOM %d %s %f %f %f %f %f %f\n" % tuple(data))
def _write_rigid_group_file(self, filename, groups):
with open(filename,'w') as fh:
for g in groups:
fh.write("#Group\n")
try:
if len(g.trans_damping) != 3: raise
fh.write(" ".join(str(v) for v in g.trans_damping) + " ")
except:
raise ValueError("Group {} lacks 3-value 'trans_damping' attribute")
try:
if len(g.rot_damping) != 3: raise
fh.write(" ".join(str(v) for v in g.rot_damping) + " ")
except:
raise ValueError("Group {} lacks 3-value 'rot_damping' attribute")
fh.write("{}\n".format(len(g)))
particles = [p for p in g]
def chunks(l, n):
"""Yield successive n-sized chunks from l."""
for i in range(0, len(l), n):
yield l[i:i + n]
for c in chunks(particles,8):
fh.write(" ".join(str(p.idx) for p in c) + "\n")
def _writeArbdConf(self, prefix, randomSeed=None, numSteps=100000000, outputPeriod=10000, restart_file=None):
## TODO: raise exception if _writeArbdPotentialFiles has not been called
filename = "%s.bd" % prefix
......@@ -816,10 +914,10 @@ class ArbdModel(PdbModel):
# params['coordinateFile'] = "%s.coord.txt" % prefix
params['particleFile'] = "%s.particles.txt" % prefix
if restartCoordinateFile is None:
if restart_file is None:
params['restartCoordinates'] = ""
else:
params['restartCoordinates'] = "restartCoordinates %s" % restartCoordinateFile
params['restartCoordinates'] = "restartCoordinates %s" % restart_file
params['outputPeriod'] = outputPeriod
for k,v in zip('XYZ', self.dimensions):
......@@ -828,6 +926,26 @@ class ArbdModel(PdbModel):
params['pairlistDistance'] -= params['cutoff']
""" Find rigid groups """
rigid_groups = []
def get_rigid_groups(parent):
ret_list = []
for c in parent.children:
is_rigid = c.rigid if 'rigid' in c.__dict__ else False
if is_rigid:
rigid_groups.append(c)
elif isinstance(c,Group):
get_rigid_groups(c)
get_rigid_groups(self)
if len(rigid_groups) > 0:
self.particle_integrator = 'FusDynamic'
rb_group_filename = "{}.rb-group.txt".format(prefix)
params['particle_integrator'] = """FusDynamics
groupFileName {}
scaleFactor 0.05""".format(rb_group_filename)
self._write_rigid_group_file(rb_group_filename, rigid_groups)
## Actually write the file
with open(filename,'w') as fh:
fh.write("""{randomSeed}
......@@ -838,6 +956,7 @@ numberFluct 0 # deprecated
interparticleForce 1 # other values deprecated
fullLongRange 0 # deprecated
temperature {temperature}
ParticleDynamicType {particle_integrator}
outputPeriod {outputPeriod}
## Energy doesn't actually get printed!
......@@ -858,10 +977,28 @@ systemSize {dimX} {dimY} {dimZ}
## TODO create new particle types if existing has grid
particleParams = pt.__dict__.copy()
particleParams['num'] = num
if self.particle_integrator in ('Brown','Brownian'):
try:
D = pt.diffusivity
except:
""" units "k K/(amu/ns)" "AA**2/ns" """
D = 831447.2 * self.temperature / (pt.mass * pt.damping_coefficient)
particleParams['dynamics'] = 'diffusion {D}'.format(D = D)
elif self.particle_integrator in ('Langevin','FusDynamic'):
try:
gamma = pt.damping_coefficient
except:
""" units "k K/(AA**2/ns)" "amu/ns" """
gamma = 831447.2 * self.temperature / (pt.mass*pt.diffusivity)
particleParams['dynamics'] = """mass {mass}
transDamping {g} {g} {g}
""".format(mass=pt.mass, g=gamma)
else:
raise ValueError("Unrecognized particle integrator '{}'".format(self.particle_integrator))
fh.write("""
particle {name}
num {num}
diffusion {diffusivity}
{dynamics}
""".format(**particleParams))
if 'grid' in particleParams:
if not isinstance(pt.grid, list): pt.grid = [pt.grid]
......@@ -996,7 +1133,7 @@ component "data" value 3
def _writeArbdBondFile( self ):
for b in list( set( [b for i,j,b,ex in self.get_bonds()] ) ):
if type(b) is not str:
if type(b) is not str and not isinstance(b, Path):
b.write_file()
with open(self._bond_filename,'w') as fh:
......@@ -1009,7 +1146,7 @@ component "data" value 3
def _writeArbdAngleFile( self ):
for b in list( set( [b for i,j,k,b in self.get_angles()] ) ):
if type(b) is not str:
if type(b) is not str and not isinstance(b, Path):
b.write_file()
with open(self._angle_filename,'w') as fh:
......@@ -1019,7 +1156,7 @@ component "data" value 3
def _writeArbdDihedralFile( self ):
for b in list( set( [b for i,j,k,l,b in self.get_dihedrals()] ) ):
if type(b) is not str:
if type(b) is not str and not isinstance(b, Path):
b.write_file()
with open(self._dihedral_filename,'w') as fh:
......@@ -1075,7 +1212,7 @@ set prefix {prefix}
set nLast 0; # increment when continueing a simulation
set n [expr $nLast+1]
set out {output_directory}/$prefix-$n
set temperature 300
set temperature {temperature}
structure $prefix.psf
coordinates $prefix.pdb
......
mrdna/arbdmodel/abstract_polymer.py 0 → 100644
View file @ 79ff74f3
This diff is collapsed.
mrdna/arbdmodel/coords.py 0 → 100644
View file @ 79ff74f3
import numpy as np
from scipy.optimize import newton
def minimizeRmsd(coordsB, coordsA, weights=None, maxIter=100):
## Going through many iterations wasn't really needed
tol = 1
count = 0
R = np.eye(3)
comB = np.zeros([3,])
cNext = coordsB
while tol > 1e-6:
q,cB,comA = _minimizeRmsd(cNext,coordsA, weights)
R = R.dot(quaternion_to_matrix(q))
assert( np.all(np.isreal( R )) )
comB += cB
cLast = cNext
cNext = (coordsB-comB).dot(R)
tol = np.sum(((cNext-cLast)**2)[:]) / np.max(np.shape(coordsB))
if count > maxIter:
Exception("Exceeded maxIter (%d)" % maxIter)
count += 1
print("%d iterations",count)
return R, comB, comA
def minimizeRmsd(coordsB, coordsA, weights=None):
q,comA,comB = _minimizeRmsd(coordsB, coordsA, weights)
assert( np.all(np.isreal( q )) )
return quaternion_to_matrix(q),comA,comB
## http://onlinelibrary.wiley.com/doi/10.1002/jcc.21439/full
def _minimizeRmsd(coordsB, coordsA, weights=None):
A = coordsA
B = coordsB
shapeA,shapeB = [np.shape(X) for X in (A,B)]
for s in (shapeA,shapeB): assert( len(s) == 2 )
A,B = [X.T if s[1] > s[0] else X for X,s in zip([A,B],(shapeA,shapeB))] # TODO: print warning
shapeA,shapeB = [np.shape(X) for X in (A,B)]
assert( shapeA == shapeB )
for X,s in zip((A,B),(shapeA,shapeB)):
assert( s[1] == 3 and s[0] >= s[1] )
# if weights is None: weights = np.ones(len(A))
if weights is None:
comA,comB = [np.mean( X, axis=0 ) for X in (A,B)]
else:
assert( len(weights[:]) == len(B) )
W = np.diag(weights)
comA,comB = [np.sum( W.dot(X), axis=0 ) / np.sum(W) for X in (A,B)]
A = np.array( A-comA )
B = np.array( B-comB )
if weights is None:
s = A.T.dot(B)
else:
s = A.T.dot(W.dot(B))
sxx,sxy,sxz = s[0,:]
syx,syy,syz = s[1,:]
szx,szy,szz = s[2,:]
K = [[ sxx+syy+szz, syz-szy, szx-sxz, sxy-syx],
[syz-szy, sxx-syy-szz, sxy+syx, sxz+szx],
[szx-sxz, sxy+syx, -sxx+syy-szz, syz+szy],
[sxy-syx, sxz+szx, syz+szy, -sxx-syy+szz]]
K = np.array(K)
# GA = np.trace( A.T.dot(W.dot(A)) )
# GB = np.trace( B.T.dot(W.dot(B)) )
## Finding GA/GB can be done more quickly
# I = np.eye(4)
# x0 = (GA+GB)*0.5
# vals = newtoon(lambda x: np.det(K-x*I), x0 = x0)
vals, vecs = np.linalg.eig(K)
i = np.argmax(vals)
q = vecs[:,i]
# RMSD = np.sqrt( (GA+GB-2*vals[i]) / len(A) )
# print("CHECK:", K.dot(q)-vals[i]*q )
return q, comB, comA
def quaternion_to_matrix(q):
assert(len(q) == 4)
## It looks like the wikipedia article I used employed a less common convention for q (see below
## http://en.wikipedia.org/wiki/Rotation_formalisms_in_three_dimensions#Rotation_matrix_.E2.86.94_quaternion
# q1,q2,q3,q4 = q
# R = [[1-2*(q2*q2 + q3*q3), 2*(q1*q2 - q3*q4), 2*(q1*q3 + q2*q4)],
# [ 2*(q1*q2 + q3*q4), 1-2*(q1*q1 + q3*q3), 2*(q2*q3 - q1*q4)],
# [ 2*(q1*q3 - q2*q4), 2*(q1*q4 + q2*q3), 1-2*(q2*q2 + q1*q1)]]
q = q / np.linalg.norm(q)
q0,q1,q2,q3 = q
R = [[1-2*(q2*q2 + q3*q3), 2*(q1*q2 - q3*q0), 2*(q1*q3 + q2*q0)],
[ 2*(q1*q2 + q3*q0), 1-2*(q1*q1 + q3*q3), 2*(q2*q3 - q1*q0)],
[ 2*(q1*q3 - q2*q0), 2*(q1*q0 + q2*q3), 1-2*(q2*q2 + q1*q1)]]
return np.array(R)
def quaternion_from_matrix( R ):
e1 = R[0]
e2 = R[1]
e3 = R[2]
# d1 = 0.5 * np.sqrt( 1+R[0,0]+R[1,1]+R[2,2] )
# d2 = 0.5 * np.sqrt( 1+R[0,0]-R[1,1]-R[2,2] )
# d2 = 0.5 * np.sqrt( 1+R[0,0]-R[1,1]-R[2,2] )
# d2 = 0.5 * np.sqrt( 1+R[0,0]-R[1,1]-R[2,2] )
d1 = 1+R[0,0]+R[1,1]+R[2,2]
d2 = 1+R[0,0]-R[1,1]-R[2,2]
d3 = 1-R[0,0]+R[1,1]-R[2,2]
d4 = 1-R[0,0]-R[1,1]+R[2,2]
maxD = max((d1,d2,d3,d4))
d = 0.5 / np.sqrt(maxD)
if d1 == maxD:
return np.array(( 1.0/(4*d),