nbPot.py 3.42 KB
Newer Older
1
2
3
4
5
import numpy as np
import scipy.optimize as opt
from scipy.interpolate import interp1d
from scipy.signal import savgol_filter as savgol

6
7
from mrdna.model.nonbonded import NonbondedScheme
from mrdna import get_resource_path
8

9
10
dxParam = -1

11
12
13
14
15
16
17
18
19
20
21
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)) )
cmaffeo2's avatar
cmaffeo2 committed
22
    assert( np.any( x > 40 ) )
23
24
25
26
27
28
    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")
29
d = np.loadtxt(get_resource_path("jj-nar-dna-pmf-100Mg.dat"))
30

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

36
def parametric_potential(x,*parms):
37
38
39
40
41
42
43
    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

44
def integrated_potential(potfn,x,numBp,*parms):
45
46
47
48
49
50
51
52
    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):
53
54
    ## 1-turn iteracting with 4 turns (more is not needed)
    return 10 * integrated_potential(parametric_potential,x,20,*parms)
55
56
57
58
59
60

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):
cmaffeo2's avatar
cmaffeo2 committed
61
62
63
64
65
    larger  = np.max([bps1,bps2])
    smaller = np.min([bps1,bps2])
    if larger >= 2:
        largerInt = int(0.5*larger)
        smaller *= larger/(2*largerInt)
66
        y = smaller * integrated_potential(parametric_potential,x,largerInt,*popt)
cmaffeo2's avatar
cmaffeo2 committed
67
    else:
68
        y = smaller * larger * parametric_potential(x,*popt)
cmaffeo2's avatar
cmaffeo2 committed
69
    y = maxForce(x,y, 100)
70
71
72
    y = savgol(y, int(5.0/(x0[1]-x0[0])), 3)    
    return y

cmaffeo2's avatar
cmaffeo2 committed
73
x = np.linspace(0,100,1001)
74
75

# np.savetxt("potentials/nb.orig.dat",np.array([x0,y0 ]).T)
cmaffeo2's avatar
cmaffeo2 committed
76
# np.savetxt("potentials/nb.orig-interp.dat",np.array([xinterp,yinterp ]).T)
77
# np.savetxt("potentials/nb.fit.dat",np.array([x,fitFun(x,*popt)]).T)
cmaffeo2's avatar
cmaffeo2 committed
78
79
# np.savetxt("potentials/nb.fit2.dat",np.array([x,nbPot(x,10,40)]).T)

80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99

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:
100
            u = nbPot(r, 0.5*typeA.nts, 0.5*typeB.nts)
101
102
        return u
nbDnaScheme = nbDnaScheme()
103
104
105
106
107
108
109

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