Skip to content
Snippets Groups Projects
Commit 9635a94e authored by cmaffeo2's avatar cmaffeo2
Browse files

Added Kim-Hummer model from R. Best's paper 10.1371/journal.pcbi.1005941

parent b6253e16
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
arbdmodel/kh_polymer_model.py 0 → 100644
+ 313
− 0
View file @ 9635a94e
# -*- coding: utf-8 -*-
## Test with `python -m arbdmodel.kh_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
from .kh_polymer_model_pair_epsilon import epsilon_mj
"""Define particle types"""
_types = dict(
A = ParticleType("ALA",
mass = 71.08,
charge = 0,
sigma = 5.04,
),
R = ParticleType("ARG",
mass = 156.2,
charge = 1,
sigma = 6.56,
),
N = ParticleType("ASN",
mass = 114.1,
charge = 0,
sigma = 5.68,
),
D = ParticleType("ASP",
mass = 115.1,
charge = -1,
sigma = 5.58,
),
C = ParticleType("CYS",
mass = 103.1,
charge = 0,
sigma = 5.48,
),
Q = ParticleType("GLN",
mass = 128.1,
charge = 0,
sigma = 6.02,
),
E = ParticleType("GLU",
mass = 129.1,
charge = -1,
sigma = 5.92,
),
G = ParticleType("GLY",
mass = 57.05,
charge = 0,
sigma = 4.5,
),
H = ParticleType("HIS",
mass = 137.1,
charge = 0.5,
sigma = 6.08,
),
I = ParticleType("ILE",
mass = 113.2,
charge = 0,
sigma = 6.18,
),
L = ParticleType("LUE",
mass = 113.2,
charge = 0,
sigma = 6.18,
),
K = ParticleType("LYS",
mass = 128.2,
charge = 1,
sigma = 6.36,
),
M = ParticleType("MET",
mass = 131.2,
charge = 0,
sigma = 6.18,
),
F = ParticleType("PHE",
mass = 147.2,
charge = 0,
sigma = 6.36,
),
P = ParticleType("PRO",
mass = 97.12,
charge = 0,
sigma = 5.56,
),
S = ParticleType("SER",
mass = 87.08,
charge = 0,
sigma = 5.18,
),
T = ParticleType("THR",
mass = 101.1,
charge = 0,
sigma = 5.62,
),
W = ParticleType("TRP",
mass = 186.2,
charge = 0,
sigma = 6.78,
),
Y = ParticleType("TYR",
mass = 163.2,
charge = 0,
sigma = 6.46,
),
V = ParticleType("VAL",
mass = 99.07,
charge = 0,
sigma = 5.86,
)
)
for k,t in _types.items():
t.resname = t.name
class KhNonbonded(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
""" KH scale model """
# alpha = 0.228
# epsilon0 = -1.0
# e_mj = epsilon_mj[(typeA.name,typeB.name)]
# epsilon = alpha * np.abs( e_mj - epsilon0 )
epsilon = epsilon_mj[(typeA.name,typeB.name)]
lambda_ = -1 if epsilon > 0 else 1
epsilon = np.abs(epsilon)
sigma = 0.5 * (typeA.sigma + typeB.sigma)
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 KhBeadsFromPolymer(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
# pdb.set_trace()
# continue
## 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 = 2.3900574,
r0 = 3.8,
rRange = (0,500),
resolution = 0.01,
maxForce = 10)
self.add_bond( i=b1, j=b2, bond = bond, exclude=True )
class KhModel(ArbdModel):
def __init__(self, polymers,
sequences = None,
debye_length = 10,
damping_coefficient = 10,
DEBUG=False,
**kwargs):
"""
[debye_length]: angstroms
[damping_coefficient]: ns
"""
kwargs['timestep'] = 10e-6
kwargs['cutoff'] = max(4*debye_length,20)
if 'decompPeriod' not in kwargs:
kwargs['decompPeriod'] = 1000
""" Assign sequences """
if sequences is None:
raise NotImplementedError("KhModel must be provided a sequences argument")
self.polymer_group = AbstractPolymerGroup(polymers)
self.sequences = sequences
ArbdModel.__init__(self, [], **kwargs)
""" Update type diffusion coefficients """
self.types = all_types = [t for key,t in _types.items()]
self.set_damping_coefficient( damping_coefficient )
""" Set up nonbonded interactions """
nonbonded = KhNonbonded(debye_length)
for i in range(len(all_types)):
t1 = all_types[i]
for j in range(i,len(all_types)):
t2 = all_types[j]
self.useNonbondedScheme( nonbonded, typeA=t1, typeB=t2 )
""" Generate beads """
self.generate_beads()
def update_splines(self, coords):
i = 0
for p in self.polymer_group.polymers:
n = p.num_monomers
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 = [KhBeadsFromPolymer(p,s)
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 self.types:
t.damping_coefficient = damping_coefficient
# t.diffusivity = 831447.2 * temperature / (t.mass * damping_coefficient)
if __name__ == "__main__":
from matplotlib import pyplot as plt
nt = len(_types)
# print("TYPES")
# for n,t in _types.items():
# print("{}\t{}\t{}\t{}\t{}".format(t.name, t.mass, t.charge, t.sigma, t.lambda_))
type_string = 'WYFMLIVAPGCQNTSEDKHR'
d = np.zeros([nt,nt])
for i in range(nt):
n1 = type_string[i]
t1 = _types[n1]
for j in range(nt):
n2 = type_string[j]
t2 = _types[n2]
d[nt-i-1,j] = epsilon_mj[(t1.name,t2.name)]
plt.imshow(d.T)
plt.show()
arbdmodel/kh_polymer_model_pair_epsilon.py 0 → 100644
+ 229
− 0
View file @ 9635a94e
_txtdata = """CYS CYS -0.509719
CYS MET -0.448877
MET MET -0.512423
CYS PHE -0.558393
MET PHE -0.661148
PHE PHE -0.755790
CYS ILE -0.517831
MET ILE -0.588137
PHE ILE -0.699005
ILE ILE -0.658443
CYS LEU -0.562449
MET LEU -0.640867
PHE LEU -0.758494
ILE LEU -0.726045
LEU LEU -0.770663
CYS VAL -0.444821
MET VAL -0.493495
PHE VAL -0.624642
ILE VAL -0.592194
LEU VAL -0.650331
VAL VAL -0.520535
CYS TRP -0.443469
MET TRP -0.524592
PHE TRP -0.607066
ILE TRP -0.555688
LEU TRP -0.604362
VAL TRP -0.474566
TRP TRP -0.458342
CYS TYR -0.336658
MET TYR -0.438061
PHE TYR -0.539464
ILE TYR -0.484030
LEU TYR -0.540816
VAL TYR -0.398852
TRP TYR -0.404260
TYR TYR -0.338010
CYS ALA -0.256888
MET ALA -0.306913
PHE ALA -0.424541
ILE ALA -0.393444
LEU ALA -0.438061
VAL ALA -0.320433
TRP ALA -0.290689
TYR ALA -0.228495
ALA ALA -0.141964
CYS GLY -0.201454
MET GLY -0.232551
PHE GLY -0.332602
ILE GLY -0.285280
LEU GLY -0.336658
VAL GLY -0.231199
TRP GLY -0.236607
TYR GLY -0.181173
ALA GLY -0.086531
GLY GLY -0.077066
CYS THR -0.194694
MET THR -0.248775
PHE THR -0.352882
ILE THR -0.319081
LEU THR -0.360995
VAL THR -0.242015
TRP THR -0.209566
TYR THR -0.181173
ALA THR -0.087883
GLY THR -0.055434
THR THR -0.060842
CYS SER -0.160893
MET SER -0.183877
PHE SER -0.317729
ILE SER -0.250127
LEU SER -0.304209
VAL SER -0.186582
TRP SER -0.178469
TYR SER -0.150076
ALA SER -0.045969
GLY SER -0.020281
THR SER -0.039209
SER SER 0.000000
CYS ASN -0.124388
MET ASN -0.173061
PHE ASN -0.281224
ILE ASN -0.212270
LEU ASN -0.279872
VAL ASN -0.156837
TRP ASN -0.189286
TYR ASN -0.147372
ALA ASN -0.022985
GLY ASN -0.009464
THR ASN -0.028393
SER ASN 0.012168
ASN ASN -0.001352
CYS GLN -0.159541
MET GLN -0.220383
PHE GLN -0.328546
ILE GLN -0.270408
LEU GLN -0.320433
VAL GLN -0.189286
TRP GLN -0.194694
TYR GLN -0.175765
ALA GLN -0.029745
GLY GLN 0.001352
THR GLN -0.031097
SER GLN 0.024337
ASN GLN -0.005408
GLN GLN 0.017577
CYS ASP -0.100051
MET ASP -0.121684
PHE ASP -0.244719
ILE ASP -0.202806
LEU ASP -0.233903
VAL ASP -0.109515
TRP ASP -0.158189
TYR ASP -0.147372
ALA ASP -0.004056
GLY ASP 0.010816
THR ASP -0.017577
SER ASP 0.005408
ASN ASP -0.001352
GLN ASP 0.028393
ASP ASP 0.062194
CYS GLU -0.081122
MET GLU -0.164949
PHE GLU -0.255536
ILE GLU -0.216326
LEU GLU -0.259592
VAL GLU -0.135204
TRP GLU -0.178469
TYR GLU -0.151428
ALA GLU 0.021633
GLY GLU 0.060842
THR GLU -0.009464
SER GLU 0.025689
ASN GLU 0.021633
GLN GLU 0.033801
ASP GLU 0.087883
GLU GLU 0.102755
CYS HIS -0.260944
MET HIS -0.312321
PHE HIS -0.419132
ILE HIS -0.333954
LEU HIS -0.388035
VAL HIS -0.258240
TRP HIS -0.312321
TYR HIS -0.250127
ALA HIS -0.100051
GLY HIS -0.064898
THR HIS -0.101403
SER HIS -0.059490
ASN HIS -0.055434
GLN HIS -0.041913
ASP HIS -0.087883
GLU HIS -0.064898
HIS HIS -0.186582
CYS ARG -0.121684
MET ARG -0.196046
PHE ARG -0.312321
ILE ARG -0.265000
LEU ARG -0.319081
VAL ARG -0.189286
TRP ARG -0.235255
TYR ARG -0.201454
ALA ARG -0.021633
GLY ARG -0.006760
THR ARG -0.031097
SER ARG 0.006760
ASN ARG 0.004056
GLN ARG -0.017577
ASP ARG -0.083826
GLU ARG -0.081122
HIS ARG -0.066250
ARG ARG 0.016224
CYS LYS -0.037857
MET LYS -0.109515
PHE LYS -0.228495
ILE LYS -0.181173
LEU LYS -0.229847
VAL LYS -0.110867
TRP LYS -0.137908
TYR LYS -0.125740
ALA LYS 0.048673
GLY LYS 0.070306
THR LYS 0.048673
SER LYS 0.083826
ASN LYS 0.062194
GLN LYS 0.051378
ASP LYS -0.001352
GLU LYS -0.017577
HIS LYS 0.043265
ARG LYS 0.146020
LYS LYS 0.209566
CYS PRO -0.189286
MET PRO -0.240663
PHE PRO -0.348826
ILE PRO -0.282576
LEU PRO -0.342066
VAL PRO -0.223087
TRP PRO -0.278520
TYR PRO -0.205510
ALA PRO -0.048673
GLY PRO -0.027041
THR PRO -0.031097
SER PRO 0.013520
ASN PRO 0.018929
GLN PRO -0.008112
ASP PRO 0.045969
GLU PRO 0.055434
HIS PRO -0.078418
ARG PRO -0.004056
LYS PRO 0.094643
PRO PRO -0.010816"""
epsilon_mj = dict()
def __add_value(key,v):
if key not in epsilon_mj:
epsilon_mj[key] = v
else:
assert(epsilon_mj[key] == v)
for line in _txtdata.split('\n'):
n1,n2,v = line.split()
v = float(v)
__add_value((n1,n2), v)
__add_value((n2,n1), v)
if __name__ == '__main__':
print(epsilon_mj)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment