From 9635a94ec0a267bc22c6c55e72cb34246dcf08c9 Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Mon, 8 Jul 2019 15:31:45 -0500
Subject: [PATCH] Added Kim-Hummer model from R. Best's paper
10.1371/journal.pcbi.1005941
---
arbdmodel/kh_polymer_model.py | 313 +++++++++++++++++++++
arbdmodel/kh_polymer_model_pair_epsilon.py | 229 +++++++++++++++
2 files changed, 542 insertions(+)
create mode 100644 arbdmodel/kh_polymer_model.py
create mode 100644 arbdmodel/kh_polymer_model_pair_epsilon.py
diff --git a/arbdmodel/kh_polymer_model.py b/arbdmodel/kh_polymer_model.py
new file mode 100644
index 0000000..ffb7432
--- /dev/null
+++ b/arbdmodel/kh_polymer_model.py
@@ -0,0 +1,313 @@
+# -*- 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()
diff --git a/arbdmodel/kh_polymer_model_pair_epsilon.py b/arbdmodel/kh_polymer_model_pair_epsilon.py
new file mode 100644
index 0000000..e696b2c
--- /dev/null
+++ b/arbdmodel/kh_polymer_model_pair_epsilon.py
@@ -0,0 +1,229 @@
+
+_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)
--
GitLab