From ea7ba628c32a367719c80ec3c1b5cb0b5ca5a003 Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Thu, 20 Sep 2018 09:34:30 -0500
Subject: [PATCH] Improved fit for non-bonded potential to all-atom data
---
mrdna/model/nbPot.py | 104 ++++++++++++++++++++++++++++++++++++-------
setup.py | 1 +
2 files changed, 90 insertions(+), 15 deletions(-)
diff --git a/mrdna/model/nbPot.py b/mrdna/model/nbPot.py
index 56e46b9..b93bcee 100644
--- a/mrdna/model/nbPot.py
+++ b/mrdna/model/nbPot.py
@@ -8,10 +8,11 @@ from mrdna import get_resource_path
dxParam = -1
-def maxForce(x,u,maxForce=40):
+def maxForce(x,u0,maxForce=40):
maxForce = maxForce * 0.014393265 # convert pN to kcal_mol/AA
dx = np.diff(x)
- f = -np.diff(u)/dx
+ f = -np.diff(u0)/dx
+ f[ np.isnan(f) ] = 0
ids = np.where(f > maxForce)[0]
f[ids] = maxForce
@@ -23,6 +24,7 @@ def maxForce(x,u,maxForce=40):
ids = np.where(x > 40)[0]
u = u - np.mean(u[ids])
+ assert( np.isnan(u).sum() == 0 )
return u
# d = np.loadtxt("resources/jj-nar-dna-pmf-20Mg-200Na.dat")
@@ -41,21 +43,48 @@ def parametric_potential(x,*parms):
u[x>np.max(x0)] = 0
return u
-def integrated_potential(potfn,x,numBp,*parms):
+def integrate_potential(potfn,x,y,*parms):
u = np.zeros(np.shape(x))
- for h in range(-numBp,numBp+1):
- r = np.sqrt(x**2 + (3.4*h)**2)
+ for h in y:
+ r = np.sqrt(x**2 + h**2)
r[r<1] = 1
u = u + potfn(r,*parms) # * x/r
return u
+
+# def integrated_potential(potfn,x,numBp,*parms):
+# u = np.zeros(np.shape(x))
+# for h in range(-numBp,numBp+1):
+# r = np.sqrt(x**2 + (3.4*h)**2)
+# r[r<1] = 1
+# u = u + potfn(r,*parms) # * x/r
+# return u
+
def fitFun(x,*parms):
## 1-turn iteracting with 4 turns (more is not needed)
- return 10 * integrated_potential(parametric_potential,x,20,*parms)
+ y = np.arange(-20,21)*3.4
+ return 10 * integrate_potential(parametric_potential,x,y,*parms)
+
+def fitFun5(x,*parms):
+ ## 1-turn iteracting with 4 turns (more is not needed)
+ y = np.arange(-20,21)*(5*3.4)
+ return 5 * 5 * integrate_potential(parametric_potential,x,y,*parms)
+
xinterp = np.linspace(np.min(x0),np.max(x0),1000) # interpolate potentialb
yinterp = interp1d(x0,y0)(xinterp)
-popt, pcov = opt.curve_fit(fitFun, xinterp, yinterp, p0=y0/10**2 ) # find fit to interpolated potential
+
+def fit_potential(bp1,bp2):
+ assert(bp2 >= bp1)
+
+ def fitFun(x,*parms):
+ ## 1-turn iteracting with at least 4 turns (more is not needed)
+ y = np.arange(-20,21)*(bp2*3.4)
+ return (10.0/bp1) * integrate_potential(parametric_potential,x,y,*parms)
+
+ popt, pcov = opt.curve_fit(fitFun, xinterp, yinterp, p0=y0/10**2 ) # find fit to interpolated potential
+ return lambda x: parametric_potential(x,*popt)
+
def nbPot(x,bps1,bps2):
larger = np.max([bps1,bps2])
@@ -63,12 +92,31 @@ def nbPot(x,bps1,bps2):
if larger >= 2:
largerInt = int(0.5*larger)
smaller *= larger/(2*largerInt)
- y = smaller * integrated_potential(parametric_potential,x,largerInt,*popt)
+ y = np.arange(-largerInt,largerInt+1)*3.4
+ u = smaller * integrate_potential(parametric_potential,x,y,*popt)
else:
- y = smaller * larger * parametric_potential(x,*popt)
- y = maxForce(x,y, 100)
- y = savgol(y, int(5.0/(x0[1]-x0[0])), 3)
- return y
+ u = smaller * larger * parametric_potential(x,*popt)
+ u1 = maxForce(x,u, 100)
+ # u2 = savgol(u1, int(5.0/(x0[1]-x0[0])), 3)
+ u2 = u
+ assert(np.sum(np.isnan(u2)) == 0)
+ return u2
+
+_pot_dict = {}
+
+def nbPot(x,bps1,bps2):
+ larger = np.max([bps1,bps2])
+ smaller = np.min([bps1,bps2])
+ key = (smaller,larger)
+ if key not in _pot_dict:
+ _pot_dict[key] = fit_potential(*key)
+ u = _pot_dict[key](x)
+
+ u1 = maxForce(x,u, 100)
+ # u2 = savgol(u1, int(5.0/(x0[1]-x0[0])), 3)
+ u2 = u
+ assert(np.sum(np.isnan(u2)) == 0)
+ return u2
x = np.linspace(0,100,1001)
@@ -104,6 +152,32 @@ nbDnaScheme = nbDnaScheme()
if __name__ == "__main__":
## run some tests
np.savetxt("jj-filtered.dat", np.array((x0,y0)).T)
- y = fitFun(x,*popt)
- np.savetxt("1-bp-fit.dat", np.array((x,y)).T),
-
+
+ # y = single_bp_pot = fitFun(x,*popt)
+ # np.savetxt("1-bp-fit.dat", np.array((x,y)).T),
+
+ # y = fitFun5(x,*popt5)
+ # np.savetxt("5-bp-fit.dat", np.array((x,y)).T),
+
+
+ # tmp = nbPot(np.linspace(20,45,100), 1,1)
+ # import pdb
+ # pdb.set_trace()
+ y = 10*integrate_potential(nbPot, x, np.arange(-20,20)*3.4 , 1, 1)
+ np.savetxt("nb-1-1.dat", np.array((x,y)).T),
+
+ y = 10*integrate_potential(nbPot, x, np.arange(-20,20)*(3.4*2) , 1, 2)
+ np.savetxt("nb-1-2.dat", np.array((x,y)).T),
+
+ y = 10*integrate_potential(nbPot, x, np.arange(-20,20)*(3.4*5) , 1, 5)
+ np.savetxt("nb-1-5.dat", np.array((x,y)).T),
+
+ y = 5*integrate_potential(nbPot, x, np.arange(-20,20)*(3.4*2) , 2, 2)
+ np.savetxt("nb-2-2.dat", np.array((x,y)).T),
+
+ y = 2*integrate_potential(nbPot, x, np.arange(-20,20)*(3.4*5) , 5, 5)
+ np.savetxt("nb-5-5.dat", np.array((x,y)).T),
+
+ # y = 10*5*integrate_potential(parametric_potential, x, np.arange(-20,20)*(3.4*5) , *popt)
+ # np.savetxt("scale-5-5.dat", np.array((x,y)).T),
+
diff --git a/setup.py b/setup.py
index e339ebf..2d104e7 100644
--- a/setup.py
+++ b/setup.py
@@ -25,6 +25,7 @@ setuptools.setup(
'scipy>=1.1',
'mdanalysis>=0.18',
'cadnano>=2.5.2.1',
+ 'appdirs>=1.4'
),
classifiers=(
"Programming Language :: Python :: 3",
--
GitLab