Improved fit for non-bonded potential to all-atom data

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),
+

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",