nbPot.py

import numpy as np
import scipy.optimize as opt
from scipy.interpolate import interp1d
from scipy.signal import savgol_filter as savgol

from mrdna.model.nonbonded import NonbondedScheme
from mrdna import get_resource_path

dxParam = -1

def maxForce(x,u,maxForce=40):
    maxForce = maxForce * 0.014393265 # convert pN to kcal_mol/AA
    dx = np.diff(x)
    f = -np.diff(u)/dx

    ids = np.where(f > maxForce)[0]
    f[ids] = maxForce
    ids = np.where(f < -maxForce)[0]
    f[ids] = -maxForce
   
    u = np.hstack( ((0), np.cumsum(-f*dx)) )
    assert( np.any( x > 40 ) )
    ids = np.where(x > 40)[0]
    
    u = u - np.mean(u[ids])
    return u

# d = np.loadtxt("resources/jj-nar-dna-pmf-20Mg-200Na.dat")
d = np.loadtxt(get_resource_path("jj-nar-dna-pmf-100Mg.dat"))

x0 = d[:,1] + dxParam
y0 = d[:,0]                     # kcal/mol for 10bp with 10bp, according to JY
y0 = y0 - np.mean(y0[x0>(38.5 + dxParam)])    # set 0 to tail of potential
y0 = savgol(y0, int(5.0/(x0[1]-x0[0])), 3) # smooth potential

def parametric_potential(x,*parms):
    x = np.array(x)
    u = interp1d(x0, np.array(parms),
                 bounds_error = False, fill_value="extrapolate",
                 assume_sorted=True)(x).T
    u[x>np.max(x0)] = 0
    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)

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 nbPot(x,bps1,bps2):
    larger  = np.max([bps1,bps2])
    smaller = np.min([bps1,bps2])
    if larger >= 2:
        largerInt = int(0.5*larger)
        smaller *= larger/(2*largerInt)
        y = smaller * integrated_potential(parametric_potential,x,largerInt,*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

x = np.linspace(0,100,1001)

# np.savetxt("potentials/nb.orig.dat",np.array([x0,y0 ]).T)
# np.savetxt("potentials/nb.orig-interp.dat",np.array([xinterp,yinterp ]).T)
# np.savetxt("potentials/nb.fit.dat",np.array([x,fitFun(x,*popt)]).T)
# np.savetxt("potentials/nb.fit2.dat",np.array([x,nbPot(x,10,40)]).T)


def _writeNonbondedParameterFiles(self, prefix):
    x = np.arange(0, 50, 0.1)
    for i,j,t1,t2 in self._particleTypePairIter():
        f = "%s.%s-%s.dat" % (prefix, t1, t2)

        if t1 == "O" or t2 == "O":
            y = np.zeros(np.shape(x))
        else:
            bps1,bps2 = [float( t[1:] )/10 for t in (t1,t2)]
            y = nbPot(x, bps1, bps2)

        np.savetxt( f, np.array([x, y]).T )
        self._nbParamFiles.append(f)

class nbDnaScheme(NonbondedScheme):
    def potential(self, r, typeA, typeB):
        if typeA.name[0] == "O" or typeB.name[0] == "O":
            u = np.zeros(np.shape(r))
        else:
            u = nbPot(r, 0.5*typeA.nts, 0.5*typeB.nts)
        return u
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),