Commit 2e1b1332 authored by cmaffeo2's avatar cmaffeo2
Browse files

Merge branch 'dev'

parents bb71bbfb d3c116bf
../LICENSE
\ No newline at end of file
../README.md
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
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),
d * (R[2,1]-R[1,2]),
d * (R[0,2]-R[2,0]),
d * (R[1,0]-R[0,1]) ))
elif d2 == maxD:
return np.array(( d * (R[2,1]-R[1,2]),
1.0/(4*d),
d * (R[0,1]+R[1,0]),
d * (R[0,2]+R[2,0]) ))
elif d3 == maxD:
return np.array(( d * (R[0,2]-R[2,0]),
d * (R[0,1]+R[1,0]),
1.0/(4*d),
d * (R[1,2]+R[2,1]) ))
elif d4 == maxD:
return np.array(( d * (R[1,0]-R[0,1]),
d * (R[0,2]+R[2,0]),
d * (R[1,2]+R[2,1]),
1.0/(4*d) ))
def rotationAboutAxis(axis,angle, normalizeAxis=True):
if normalizeAxis: axis = axis / np.linalg.norm(axis)
angle = angle * 0.5 * np.pi/180
cos = np.cos( angle )
sin = np.sin( angle )
q = [cos] + [sin*x for x in axis]
return quaternion_to_matrix(q)
def readArbdCoords(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):
import MDAnalysis as mda
usel = mda.Universe(psf, dcd)
sel = usel.select_atoms("name D*")
# r0 = ref.xyz[0,ids,:]
ts = usel.trajectory[-1]
r0 = sel.positions
pos = []
for t in range(ts.frame-1,-1,-1):
usel.trajectory[t]
R,comA,comB = minimizeRmsd(sel.positions,r0)
r = np.array( [(r-comA).dot(R)+comB for r in sel.positions] )
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:
break
t0=t+1
print( "Averaging coordinates in %s after frame %d" % (dcd, t0) )
pos = np.mean(pos, axis=0)
return pos
def unit_quat_conversions():
for axis in [[0,0,1],[1,1,1],[1,0,0],[-1,-2,0]]:
for angle in np.linspace(-180,180,10):
R = rotationAboutAxis(axis, angle)
R2 = quaternion_to_matrix( quaternion_from_matrix( R ) )
if not np.all( np.abs(R-R2) < 0.01 ):
import pdb
pdb.set_trace()
quaternion_to_matrix( quaternion_from_matrix( R ) )
if __name__ == "__main__":
unit_quat_conversions()
# -*- coding: utf-8 -*-
## Test with `python -m arbdmodel.hps_polymer_model`
import numpy as np
import sys
## Local imports
from . import ArbdModel, ParticleType, PointParticle, Group, get_resource_path
from .abstract_polymer import PolymerSection, AbstractPolymerGroup
from .interactions import NonbondedScheme, HarmonicBond, HarmonicPotential
from .coords import quaternion_to_matrix
"""Define particle types"""
type_ = ParticleType("X",
damping_coefficient = 40.9,
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)
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
class FjcBeadsFromPolymer(Group):
def __init__(self, polymer, sequence=None,
rest_length = 3.8, spring_constant = 25,
**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)
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]
""" 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 )
class FjcModel(ArbdModel):
def __init__(self, polymers,
sequences = None,
rest_length = 3.8,
spring_constant = 25,
damping_coefficient = 40.9,
DEBUG=False,
**kwargs):
"""
[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))]
self.polymer_group = AbstractPolymerGroup(polymers)
self.sequences = sequences
self.rest_length = rest_length
self.spring_constant = spring_constant
ArbdModel.__init__(self, [], **kwargs)
""" Update type diffusion coefficients """
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()
def set_damping_coefficient(self, damping_coefficient):
for t in [type_]:
t.damping_coefficient = damping_coefficient
# t.diffusivity = 831447.2 * temperature / (t.mass * damping_coefficient)
if __name__ == "__main__":
pass
# -*- coding: utf-8 -*-
## Test with `python -m arbdmodel.hps_polymer_model`
import numpy as np
import sys
## Local imports
from . import ArbdModel, ParticleType, PointParticle, Group, get_resource_path
from .abstract_polymer import PolymerSection, AbstractPolymerGroup
from .interactions import NonbondedScheme, HarmonicBond, HarmonicAngle, HarmonicDihedral
from .coords import quaternion_to_matrix
"""Define particle types"""
_types = dict(
A = ParticleType("ALA",
mass = 71.08,
charge = 0,
sigma = 5.04,
lambda_ = 0.72973,
),
R = ParticleType("ARG",
mass = 156.2,
charge = 1,
sigma = 6.56,
lambda_ = 0.0,
),
N = ParticleType("ASN",
mass = 114.1,
charge = 0,
sigma = 5.68,
lambda_ = 0.432432,
),
D = ParticleType("ASP",
mass = 115.1,
charge = -1,
sigma = 5.58,
lambda_ = 0.378378,
),
C = ParticleType("CYS",
mass = 103.1,
charge = 0,
sigma = 5.48,
lambda_ = 0.594595,
),
Q = ParticleType("GLN",
mass = 128.1,
charge = 0,
sigma = 6.02,
lambda_ = 0.513514,
),
E = ParticleType("GLU",
mass = 129.1,
charge = -1,
sigma = 5.92,
lambda_ = 0.459459,
),
G = ParticleType("GLY",
mass = 57.05,
charge = 0,
sigma = 4.5,
lambda_ = 0.648649,
),
H = ParticleType("HIS",
mass = 137.1,
charge = 0.5,
sigma = 6.08,
lambda_ = 0.513514,
),
I = ParticleType("ILE",
mass = 113.2,
charge = 0,
sigma = 6.18,
lambda_ = 0.972973,
),
L = ParticleType("LEU",
mass = 113.2,
charge = 0,
sigma = 6.18,
lambda_ = 0.972973,
),
K = ParticleType("LYS",
mass = 128.2,
charge = 1,
sigma = 6.36,
lambda_ = 0.513514,
),
M = ParticleType("MET",
mass = 131.2,
charge = 0,
sigma = 6.18,
lambda_ = 0.837838,
),
F = ParticleType("PHE",
mass = 147.2,
charge = 0,
sigma = 6.36,
lambda_ = 1.0,
),
P = ParticleType("PRO",
mass = 97.12,
charge = 0,
sigma = 5.56,
lambda_ = 1.0,
),
S = ParticleType("SER",
mass = 87.08,
charge = 0,
sigma = 5.18,
lambda_ = 0.594595,
),
T = ParticleType("THR",
mass = 101.1,
charge = 0,
sigma = 5.62,
lambda_ = 0.675676,
),
W = ParticleType("TRP",
mass = 186.2,
charge = 0,
sigma = 6.78,
lambda_ = 0.945946,
),
Y = ParticleType("TYR",
mass = 163.2,
charge = 0,
sigma = 6.46,
lambda_ = 0.864865,
),
V = ParticleType("VAL",
mass = 99.07,
charge = 0,
sigma = 5.86,
lambda_ = 0.891892,
)
)
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)
self.debye_length = debye_length
self.maxForce = 50
def potential(self, r, typeA, typeB):
""" Electrostatics """
ld = self.debye_length
q1 = typeA.charge
q2 = typeB.charge
D = 80 # dielectric of water
## units "e**2 / (4 * pi * epsilon0 AA)" kcal_mol
A = 332.06371
u_elec = (A*q1*q2/D)*np.exp(-r/ld) / r
""" Hydrophobicity scale model """
lambda_ = 0.5 * (typeA.lambda_ + typeB.lambda_)
sigma = 0.5 * (typeA.sigma + typeB.sigma)
epsilon = 0.2
r6 = (sigma/r)**6
r12 = r6**2
u_lj = 4 * epsilon * (r12-r6)
u_hps = lambda_ * np.array(u_lj)
s = r<=sigma*2**(1/6)
u_hps[s] = u_lj[s] + (1-lambda_) * epsilon
u = u_elec + u_hps
u[0] = u[1] # Remove NaN
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))
u = u-u[-1]
return u
class HpsBeadsFromPolymer(Group):
# p = PointParticle(_P, (0,0,0), "P")
# b = PointParticle(_B, (3,0,1), "B")
# nt = Group( name = "nt", children = [p,b])
# nt.add_bond( i=p, j=b, bond = get_resource_path('two_bead_model/BPB.dat') )
def __init__(self, polymer, sequence=None, **kwargs):
if sequence is None:
raise NotImplementedError
# ... set random sequence
self.polymer = polymer
self.sequence = sequence
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)
def _clear_beads(self):
...
def _generate_beads(self):
# beads = self.children
for i in range(self.polymer.num_monomers):
c = self.polymer.monomer_index_to_contour(i)
r = self.polymer.contour_to_position(c)
s = self.sequence[i]
bead = PointParticle(_types[s], r,
name = s,
resid = i+1)
self.add(bead)
# import pdb