nbPot.py 3.08 KiB
import numpy as np
import scipy.optimize as opt
from scipy.interpolate import interp1d
from scipy.signal import savgol_filter as savgol
from nonbonded import NonbondedScheme
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("resources/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 pot(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 potForBp(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):
return 10 * potForBp(pot,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 * potForBp(pot,x,largerInt,*popt)
else:
y = smaller * larger * pot(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.nbPot(r, typeA.nts*2, typeB.nts*2)
u = nbPot(r, typeA.nts, typeB.nts)
return u
nbDnaScheme = nbDnaScheme()