diff --git a/mrdna/model/nbPot.py b/mrdna/model/nbPot.py
index 6941f479cb05100ec26c207886dc8e9741582683..1ad9022c88a33326b1fe9382e45093183f87a47d 100644
--- a/mrdna/model/nbPot.py
+++ b/mrdna/model/nbPot.py
@@ -7,37 +7,13 @@ from mrdna.model.nonbonded import NonbondedScheme
from mrdna import get_resource_path
from mrdna.config import CACHE_DIR
-## TODO: clean this file up a bit, make things protected
-
-_dxParam = -1
-
-def _maxForce(x,u0,maxForce=40):
- maxForce = maxForce * 0.014393265 # convert pN to kcal_mol/AA
- dx = np.diff(x)
- f = -np.diff(u0)/dx
- f[ np.isnan(f) ] = 0
-
- 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])
- assert( np.isnan(u).sum() == 0 )
- return u
+from scipy.optimize import curve_fit
+from pathlib import Path
-def integrate_potential(potfn,x,y,*parms):
- u = np.zeros(np.shape(x))
- for h in y:
- r = np.sqrt(x**2 + h**2)
- r[r<1] = 1
- u = u + potfn(r,*parms) # * x/r
- return u
+## TODO: clean this file up a bit, make things protected
+## units "291 k K" kcal_mol
+kT = 0.57827709
def _round_bp(bp):
return 10**(np.around(np.log10(bp)*100)/100.0)
@@ -50,13 +26,21 @@ class AbstractNbDnaScheme(NonbondedScheme):
x,y = self.load_pmf()
self.target_x = x
self.target_y = y
- self.target_x_interp = np.linspace(np.min(x),np.max(x),1000) # interpolate potential
+ base = 4
+ self.param_x = np.logspace(np.log(15)/np.log(base),np.log(120)/np.log(base),80,base=base)
+ # print(self.param_x)
+
+ self.target_x_interp = np.linspace(np.min(x),np.max(x),2*len(x)) # interpolate potential
+ self.target_x_interp = np.arange(23.5,38,0.1) # interpolate potential
self.target_y_interp = interp1d(x,y)(self.target_x_interp)
+
+ self.bead_distribution_y = None
""" Initialize caches """
- self._pot_dict = dict() # cache including only what has already been acessed
try:
- self._cached_params = np.load( self.get_package_cache_file() )
+ self.single_bp_pair_pot = np.loadtxt( self.get_package_cache_file() ).T
+ self.param_x = self.single_bp_pair_pot[0,:]
+ assert( len(self.param_x > 2) )
except:
pass
@@ -66,7 +50,7 @@ class AbstractNbDnaScheme(NonbondedScheme):
def load_pmf(self):
raise NotImplementedError()
- def get_cache_file_for_pair(self,key):
+ def load_bead_distributions(self):
raise NotImplementedError()
def potential(self, r, typeA, typeB):
@@ -82,40 +66,233 @@ class AbstractNbDnaScheme(NonbondedScheme):
def load_np(self,filename):
return np.load(filename)['arr_0'].tolist()
- def set_cache_params(self,key):
- if key not in self._pot_dict:
- try:
- self._pot_dict[key] = self._cached_params['{}:{}'.format(*key)]
- except:
- cache_file = self.get_cache_file_for_pair(key)
- try:
- self._pot_dict[key] = self.load_np(cache_file)
- except:
- self._pot_dict[key] = self.fit_potential(*key)
- np.savez(cache_file, self._pot_dict[key])
-
def parametric_potential(self,x,*parms):
- x0 = self.target_x
+ x0 = self.param_x
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 fit_potential(self,bp1,bp2):
- assert(bp2 >= bp1)
+ return u-u[-1]
- def fitFun(x,*parms):
- ## 1-turn iteracting with at least 4 turns (more is not needed)
- y = np.arange(-20,21)*(bp2*3.4)
+ def get_bead_distributions(self, interhelical_distances):
+ if self.bead_distribution_y is None:
+ self.load_bead_distributions()
+ return self.bead_distribution_x, self.bead_distribution_y(interhelical_distances)
- return (10.0/bp1) * integrate_potential(self.parametric_potential, x,y,*parms)
+ def iterative_fit(self):
+ try:
+ import bead_dist_data as bdd
+ x = bdd.bead_centers.x
+ except:
+ ## First iteration
+ x = np.arange(0,200,0.2)
+ f = self.get_target_force(x+2)
+ f[f>0.5] = 0.5
+ u = np.cumsum(f[::-1])[::-1]
+ u = 0.5* u / (10*20)
+ np.savetxt("pair-pot.dat", np.array((x,u)).T)
+ return interp1d(x,u,
+ bounds_error = False, fill_value="extrapolate",
+ assume_sorted=True)(self.param_x).T
+
+ x = x[x<50]
+ data = np.loadtxt("last/pair-pot.dat")
+ u_last = interp1d(data[:,0], data[:,1],
+ bounds_error = False, fill_value="extrapolate",
+ assume_sorted=True)(x).T
+
+ particle_force = np.diff(u_last) / np.diff(x)
+
+ for b in bdd.bead_dists:
+ ids = bdd.bead_centers < 10
+ b[:,ids] = 0
+
+ R = bdd.interhelical_spacing
+ R,y = self.load_rau_pressure()
+
+ P = self.get_pressure(R)
+ P_t = self.get_target_pressure(R)
+ dP = P_t-P
+
+ np.savetxt("raw_pressure.dat", self.get_raw_pressure())
+ np.savetxt("pressure.dat", np.array((R,P)).T)
+ np.savetxt("target_pressure.dat", np.array((R,P_t)).T)
+
+ def fitFun(R,*uvals):
+ df_interp = dfFun(*uvals)
+ dP_local = self.get_virial_pressure(R, df_interp)
+ print( np.std(dP_local-dP) )
+ return dP_local
+
+ def duFun(x,*popt):
+ f = dfFun(*popt)(x)
+ u = -np.cumsum(f) * np.mean( np.diff(x) )
+ u = u - u[-1]
+ return u
+
+ def dfFun(x0,f0,x1_par,x2):
+ x1 = (1-x1_par)*x0 + x1_par*x2
+ t = np.linspace(0,1,100)
+ x = x0*(1-t)**2 + 2*x1*(1-t)*t + x2*t**2
+ f = f0*(1-t)**2
+ return interp1d(x, f,
+ bounds_error = False,
+ fill_value=0,
+ assume_sorted=True)
+
+ dP_sofar = np.zeros(dP.shape)
+ popt_list = []
+ for i in range(3):
+ popt, pcov = opt.curve_fit(fitFun, R, dP-dP_sofar,
+ bounds=((10,-0.1,0.1,35),
+ (34,0.1,0.9,50)),
+ p0=(20,0,0.25,50))
+ popt_list.append(popt)
+ dP_sofar += fitFun(R,*popt)
+
+ du = np.zeros(x.shape)
+ for popt in popt_list:
+ du = du + duFun(x,*popt)
+
+ du = savgol(du,int(2.5/(x[1]-x[0]))*2+1,3)
+ u = u_last+0.75*du
+ u = u-u[-1]
+
+ np.savetxt("pair-pot.dat", np.array((x,u)).T)
+
+ return interp1d(x,u,
+ bounds_error = False, fill_value="extrapolate",
+ assume_sorted=True)(self.param_x).T
+
+
+ def get_virial_pressure(self, R, bead_force):
+ """ Return pressure as fn of interhelical spacing """
+ import bead_dist_data as bdd
+
+ interhelical_spacing = bdd.interhelical_spacing
+ volume = (np.pi*np.array(bdd.radii)**2) * bdd.height
+
+ # units "kcal_mol/AA**3" atm
+ # * 68568.409
+
+ xy = bdd.bead_centers[np.newaxis,:]
+ r = bdd.bead_centers[:,np.newaxis]
+ force = bead_force(r)
+
+ """ Pxy = 1/(4V) Sum_ij (xy_ij)**2 * bead_force(r_ij) / r_ij """
+ ## 2*count because count only includes i < j
+ pressure = [68568.409*np.sum( 2*count * (xy**2) * force / r ) / (4 * vol)
+ for vol,count in zip(volume,bdd.bead_dists)]
+ return pchip_interpolate(interhelical_spacing, pressure, R)
+ """
+ return interp1d(interhelical_spacing, pressure,
+ bounds_error = False,
+ kind='cubic',
+ fill_value="extrapolate",
+ assume_sorted=True)(R)
+
+ ## Fit double exponential
- popt, pcov = opt.curve_fit(fitFun, self.target_x_interp,
- self.target_y_interp,
- p0=self.target_y/10**2 ) # find fit to interpolated potential
- return popt
+ try:
+ def fitfn(x,A,a,b):
+ return A*a**(-b**x)
+ def fit_wrap(x,*popt):
+ print("pressure:",pressure)
+ print(" ",fitfn(interhelical_spacing, *popt) - pressure)
+ return fitfn(x,*popt)
+ popt, pcov = opt.curve_fit(fit_wrap, interhelical_spacing, pressure,
+ bounds=((-np.inf,0,0),
+ (np.inf,np.inf,np.inf)),
+ p0=(0,1,1))
+ return fitfn(R,*popt)
+ except:
+ return interp1d(interhelical_spacing, pressure,
+ bounds_error = False,
+ fill_value="extrapolate",
+ assume_sorted=True)(R)
+ """
+
+
+ def get_interhelical_density(self,R):
+ data = np.loadtxt("last/density.dat")
+ r,rho = data.T
+ rho = rho / r**2
+ rho = rho/np.trapz(rho,r)
+ filter_points = int(2.5/(r[1]-r[0]))*2 + 1
+ rho = savgol(rho, filter_points, 3) # smooth density
+ return interp1d(r,rho,
+ bounds_error = False, fill_value=0,
+ assume_sorted=True)(R).T
+
+ def get_raw_pressure(self):
+ from bead_dist_data import interhelical_spacing, pressure
+ distance = interhelical_spacing
+ assert( len(distance) > 2 )
+ return distance,pressure
+
+ def get_pressure(self,R):
+ distance,pressure = self.get_raw_pressure()
+ # distance,pressure = np.loadtxt("last/pressure.dat")
+ return pchip_interpolate(distance, pressure, R)
+ return interp1d(distance, pressure,
+ bounds_error = False,
+ kind='cubic',
+ fill_value="extrapolate",
+ assume_sorted=True)(R)
+ def fitfn(x,a,b):
+ return a**(-b**x)
+
+ popt, pcov = opt.curve_fit(fitfn, distance, pressure,
+ bounds=(0,np.inf),
+ p0=(1,1))
+ return fitfn(R,*popt)
+
+
+ def get_interhelical_force(self,R):
+ y = pressure = self.get_pressure(R)
+ x = R
+ force = (x*y / np.sqrt(3)) * 34 * 1.4583976e-05 # convert atm AA**2 to kcal/mol AA
+ return force
+
+ def get_target_interhelical_density(self,R):
+ r,u = self.load_pmf()
+ rho = np.exp( -u/kT )
+ rho = rho / np.trapz(rho,r)
+ return interp1d(r,rho,
+ bounds_error = False, fill_value=0,
+ assume_sorted=True)(R).T
+
+ def get_target_force(self,R):
+ x,y = self.load_rau_force()
+ # return pchip_interpolate(x, y, R)
+ # return interp1d(x,y,
+ # bounds_error = False,
+ # # kind='cubic',
+ # fill_value="extrapolate",
+ # assume_sorted=True)(R)
+
+ def fitfn(x,A,a):
+ return A * np.exp(-(x-18)/a)
+ popt,pcov = curve_fit(fitfn, x,y, p0=(10,20))
+ return fitfn(R,*popt)
+
+ def get_target_pressure(self,R):
+ x,y = self.load_rau_pressure()
+ return pchip_interpolate(x, y, R)
+ return interp1d(x,y,
+ bounds_error = False,
+ kind='cubic',
+ fill_value="extrapolate",
+ assume_sorted=True)(R)
+ def fitfn(x,A,a):
+ return A * np.exp(-(x-18)/a)
+ popt,pcov = curve_fit(fitfn, x,y, p0=(10,20))
+ return fitfn(R,*popt)
+
+ def load_rau_force(self):
+ x,y = self.load_rau_pressure()
+ y = (x*y / np.sqrt(3)) * 34 * 1.4583976e-05 # convert atm AA**2 to kcal/mol AA
+ return x,y
def get_rounded_bp(self,bps1,bps2):
larger = np.max([bps1,bps2])
@@ -126,45 +303,11 @@ class AbstractNbDnaScheme(NonbondedScheme):
return smaller,larger,smaller_shift,larger_shift
def nbPot(self, x, bps1, bps2):
- smaller,larger,smaller_shift,larger_shift = self.get_rounded_bp(bps1,bps2)
- key = (smaller_shift,larger_shift)
- self.set_cache_params(key)
-
- u = self.parametric_potential(x,*self._pot_dict[key])
- u = u * ( (smaller/smaller_shift)*(larger/larger_shift) )
-
- # u1 = _maxForce(x,u, 100)
- # u2 = savgol(u1, int(5.0/(self.x0[1]-self.x0[0])), 3)
- u2 = u
- assert(np.sum(np.isnan(u2)) == 0)
- return u2
-
- """ Creating package cache """
-
- def _write_cache(self):
- keys = []
- for i in np.arange(0.4,5.5,0.01):
- ir = _round_bp(i)
- for j in np.arange(0.4,5.5,0.01):
- if j < i: continue
- keys.append( (ir,_round_bp(j)) )
-
- keys = sorted(list(set(keys)))
- print(len(keys))
-
- new_cache_dict = dict()
-
- args = [(self, new_cache_dict, key) for key in keys]
-
- from multiprocessing import Pool
- p = Pool(15)
- print("Starting pool")
- for arg in args:
- _add_parameter_to_cache_dict(arg)
-
- p.map( _add_parameter_to_cache_dict, args )
- print("done with pool")
- np.savez( self.get_package_cache_file(), **new_cache_dict)
+ x0,u0 = self.single_bp_pair_pot
+ u = self.parametric_potential(x,*u0)
+ u = u * bps1 * bps2
+ assert(np.sum(np.isnan(u)) == 0)
+ return u
def _add_parameter_to_cache_dict( args ):
@@ -175,58 +318,133 @@ def _add_parameter_to_cache_dict( args ):
class nbDnaScheme100Mg(AbstractNbDnaScheme):
- def get_package_cache_file(self):
- return get_resource_path("nbPot.jj-nar-dna-pmf-100Mg.npz")
+ def get_package_cache_file(self):
+ return get_resource_path("nbPot.rau_refined_100Mg.dat")
+
+ def load_rau_pressure(self):
+ data = np.loadtxt("rau-data-100Mg.dat")
+ x = data[:,0]
+ y0 = data[:,1]
+ y = (10**y0) * 9.8692327e-07 # convert dyne/cm^2 to atm
+ return x,y
+
def load_pmf(self):
d = np.loadtxt(get_resource_path("jj-nar-dna-pmf-100Mg.dat"))
- x0 = d[:,1] + _dxParam
+ x0 = d[:,1]
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
+ y0 = y0 - np.mean(y0[x0 > 38.5]) # set 0 to tail of potential
+ filter_points = int(2.5/(x0[1]-x0[0]))*2 + 1
+ y0 = savgol(y0, filter_points, 3) # smooth potential
return x0,y0
- def get_cache_file_for_pair(self,key):
- filename = "nbPot-100Mg-{}-{}.npz".format(*key)
- return CACHE_DIR / filename
+ def load_bead_distributions(self):
+ from bead_dist_data import helix_bins, bead_bins
+ from bead_dist_data import helix_centers, bead_centers
+ import bead_dist_data as bdd
+
+ self.bead_distribution_x = x = bdd.bead_centers
+
+ all_x = []
+ all_y = []
+ for i in range(0,len(bdd.helix_centers),1):
+ y = bdd.bead_dists_for_helix_dist[i]
+
+ # if y.sum() == 0: continue
+ if bdd.helix_centers[i] < 20: continue
+ all_x.append( bdd.helix_centers[i] )
+ filter_points = int(5/(x[1]-x[0]))*2 + 1
+ y = savgol( y, filter_points, 3) # smooth and normalize data
+ y[y<0] = 0 # ensure positive
+ y = 10*20 * y/np.trapz(y,x) # normalize
+ all_y.append( y )
+
+ self.bead_distribution_y = interp1d( np.array(all_x), np.array(all_y),
+ axis=0,
+ kind='linear',
+ bounds_error = False, fill_value="extrapolate",
+ assume_sorted=True)
+
class nbDnaScheme20Mg(AbstractNbDnaScheme):
- def get_package_cache_file(self):
- return get_resource_path("nbPot.jj-nar-dna-pmf-20Mg-200Na.npz")
+ def get_package_cache_file(self):
+ return get_resource_path("nbPot.rau_refined_20Mg_200Na.dat")
+
+ def load_rau_pressure(self):
+ data = np.loadtxt("rau-data-20Mg.dat")
+ x = data[:,0]
+ y0 = data[:,1]
+ y = (10**y0) * 9.8692327e-07 # convert dyne/cm^2 to atm
+ return x,y
+
def load_pmf(self):
d = np.loadtxt(get_resource_path("jj-nar-dna-pmf-20Mg-200Na.dat"))
- x0 = d[:,1] + _dxParam
+ x0 = d[:,1]
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
+ y0 = y0 - np.mean(y0[x0 > 38.5]) # set 0 to tail of potential
+ filter_points = int(2.5/(x0[1]-x0[0]))*2 + 1
+ y0 = savgol(y0, filter_points, 3) # smooth potential
return x0,y0
- def get_cache_file_for_pair(self,key):
- filename = "nbPot-20Mg-200Na-{}-{}.npz".format(*key)
- return CACHE_DIR / filename
+ def load_bead_distributions(self):
+ import bead_dist_data as bdd
-nbDnaScheme = nbDnaScheme100Mg()
+ self.bead_distribution_centers = bdd.bead_centers
+ self.bead_distribution_radii = bdd.radii
+ self.bead_distribution_counts = bdd.bead_dists
+ ...
+
+
+ for i in range(0,len(bdd.helix_centers),1):
+ y = bdd.bead_dists_for_helix_dist[i]
+
+ # if y.sum() == 0: continue
+ if bdd.helix_centers[i] < 20: continue
+ all_x.append( bdd.helix_centers[i] )
+ filter_points = int(5/(x[1]-x[0]))*2 + 1
+ y = savgol( y, filter_points, 3) # smooth and normalize data
+ y[y<0] = 0 # ensure positive
+ y = 10*20 * y/np.trapz(y,x) # normalize
+ all_y.append( y )
-def _plot_things():
- np.savetxt("jj-filtered.dat", np.array((x0,y0)).T)
- x = np.arange(0, 50, 0.1)
- ## run some tests
- for i in (0.5,1,1.05,2,3.333,5):
- for j in (1,1.05,1.1,2,3.333, 4,5):
- if i > j: continue
-
- # y = (10.0/i)*integrate_potential(nbPot_full, x, np.arange(-20,20)*(3.4*j) , i, j)
- # np.savetxt("test1-{}-{}.dat".format(i,j), np.array((x,y)).T),
+ self.bead_distribution_y = interp1d( np.array(all_x), np.array(all_y),
+ axis=0,
+ bounds_error = False, fill_value="extrapolate",
+ assume_sorted=True)
- y = (10.0/i)*integrate_potential(nbPot, x, np.arange(-20,20)*(3.4*j) , i, j)
- np.savetxt("test2-{}-{}.dat".format(i,j), np.array((x,y)).T),
+
+# nbDnaScheme100MgObj = nbDnaScheme100Mg()
+nbDnaScheme = nbDnaScheme20MgObj = nbDnaScheme20Mg()
if __name__ == "__main__":
# import cProfile
# cProfile.run("_run()")
# _generate_cache()
+ # tmp = nbDnaScheme20Mg()
+ # tmp._write_cache()
+ # pass
+
tmp = nbDnaScheme20Mg()
- tmp._write_cache()
- pass
+
+ # dR = 1
+ # for R in np.arange(22,35,2):
+ # x,ys = tmp.get_bead_distributions(np.array([R-0.5*dR,R+0.5*dR]))
+ # y1,y2 = ys
+ # d_rho_dR = (y2-y1)/dR
+ # # plt.plot(x,d_rho_dR)
+
+ # x = np.arange(20,50,1)
+ # tmp.nbPot(x,1,1)
+
+
+ # plt.gca().set_xlim((10,40))
+ # plt.savefig("test.pdf")
+
+ # x = np.arange(20,50,1)
+ # tmp.nbPot(x,1,1)
+
+ nb = nbDnaScheme20MgObj
+ x = np.arange(10,50,1)
+ u = nb.nbPot(x,1,1)
diff --git a/mrdna/resources/nbPot.jj-nar-dna-pmf-100Mg.npz b/mrdna/resources/nbPot.jj-nar-dna-pmf-100Mg.npz
deleted file mode 100644
index 3424b86a976637c5ef74a864b87487d1e863ec56..0000000000000000000000000000000000000000
Binary files a/mrdna/resources/nbPot.jj-nar-dna-pmf-100Mg.npz and /dev/null differ
diff --git a/mrdna/resources/nbPot.jj-nar-dna-pmf-20Mg-200Na.npz b/mrdna/resources/nbPot.jj-nar-dna-pmf-20Mg-200Na.npz
deleted file mode 100644
index f596d0268c3ebc161b03c4b4d2c22d6c8dfd3dca..0000000000000000000000000000000000000000
Binary files a/mrdna/resources/nbPot.jj-nar-dna-pmf-20Mg-200Na.npz and /dev/null differ
diff --git a/mrdna/resources/nbPot.rau_refined_20Mg_200Na.dat b/mrdna/resources/nbPot.rau_refined_20Mg_200Na.dat
new file mode 100644
index 0000000000000000000000000000000000000000..0bfb24d0cabf3d1a73b91b8c977d1fd1ff76ae09
--- /dev/null
+++ b/mrdna/resources/nbPot.rau_refined_20Mg_200Na.dat
@@ -0,0 +1,250 @@
+1.000000000000000056e-01 1.482384481649214880e-01
+3.000000000000000444e-01 1.469884481649214869e-01
+5.000000000000000000e-01 1.457384481649215136e-01
+7.000000000000000666e-01 1.444884481649215402e-01
+9.000000000000000222e-01 1.432384481649215668e-01
+1.100000000000000089e+00 1.419884481649215380e-01
+1.300000000000000266e+00 1.407384481649215091e-01
+1.500000000000000000e+00 1.394884481649215358e-01
+1.700000000000000178e+00 1.382384481649215069e-01
+1.899999999999999911e+00 1.369884481649215335e-01
+2.100000000000000089e+00 1.357384481649215602e-01
+2.300000000000000266e+00 1.344884481649215313e-01
+2.500000000000000000e+00 1.332384481649215302e-01
+2.700000000000000178e+00 1.319884481649215013e-01
+2.900000000000000355e+00 1.307384481649215002e-01
+3.100000000000000089e+00 1.294884481649214991e-01
+3.300000000000000266e+00 1.282384481649215535e-01
+3.500000000000000000e+00 1.269884481649214969e-01
+3.700000000000000178e+00 1.257384481649215235e-01
+3.900000000000000355e+00 1.244884481649214808e-01
+4.099999999999999645e+00 1.232384481649214797e-01
+4.300000000000000711e+00 1.219884481649214786e-01
+4.500000000000000000e+00 1.207384481649214775e-01
+4.700000000000001066e+00 1.194884481649214764e-01
+4.900000000000000355e+00 1.182384481649214752e-01
+5.099999999999999645e+00 1.169884481649214741e-01
+5.300000000000000711e+00 1.157384481649214730e-01
+5.500000000000000000e+00 1.144884481649214719e-01
+5.700000000000001066e+00 1.132384481649214708e-01
+5.900000000000000355e+00 1.119884481649214697e-01
+6.099999999999999645e+00 1.107384481649214686e-01
+6.300000000000000711e+00 1.094884481649214675e-01
+6.500000000000000000e+00 1.082384481649214664e-01
+6.700000000000001066e+00 1.069884481649214653e-01
+6.900000000000000355e+00 1.057384481649214641e-01
+7.099999999999999645e+00 1.044884481649214908e-01
+7.300000000000000711e+00 1.032384481649214897e-01
+7.500000000000000000e+00 1.019884481649214886e-01
+7.700000000000001066e+00 1.007384481649214875e-01
+7.900000000000000355e+00 9.948844816492148635e-02
+8.100000000000001421e+00 9.823844816492148524e-02
+8.300000000000000711e+00 9.698844816492148413e-02
+8.500000000000000000e+00 9.573844816492148302e-02
+8.699999999999999289e+00 9.448844816492148191e-02
+8.900000000000000355e+00 9.323844816492148080e-02
+9.100000000000001421e+00 9.198844816492147969e-02
+9.300000000000000711e+00 9.073844816492147858e-02
+9.500000000000000000e+00 8.948844816492147747e-02
+9.700000000000001066e+00 8.823844816492147636e-02
+9.900000000000000355e+00 8.698844816492147525e-02
+1.010000000000000142e+01 8.573844816492147414e-02
+1.030000000000000071e+01 8.448844816492147303e-02
+1.050000000000000000e+01 8.323844816492147192e-02
+1.070000000000000107e+01 8.198844816492147081e-02
+1.090000000000000036e+01 8.073844816492146970e-02
+1.110000000000000142e+01 7.948844816492146859e-02
+1.130000000000000071e+01 7.823844816492146748e-02
+1.150000000000000000e+01 7.698844816492146637e-02
+1.170000000000000107e+01 7.573844816492146526e-02
+1.190000000000000036e+01 7.448844816492146415e-02
+1.210000000000000142e+01 7.323844816492146303e-02
+1.230000000000000071e+01 7.198844816492146192e-02
+1.250000000000000000e+01 7.073844816492146081e-02
+1.270000000000000107e+01 6.948844816492145970e-02
+1.290000000000000036e+01 6.823844816492145859e-02
+1.310000000000000142e+01 6.698844816492145748e-02
+1.330000000000000071e+01 6.573844816492145637e-02
+1.350000000000000000e+01 6.448844816492148302e-02
+1.370000000000000107e+01 6.323844816492148191e-02
+1.390000000000000036e+01 6.198844816492148080e-02
+1.410000000000000142e+01 6.073844816492147969e-02
+1.430000000000000071e+01 5.948844816492148552e-02
+1.450000000000000000e+01 5.823844816492149135e-02
+1.470000000000000107e+01 5.698844816492148330e-02
+1.490000000000000036e+01 5.573844816492148219e-02
+1.510000000000000142e+01 5.448844816492147414e-02
+1.530000000000000071e+01 5.323844816492147997e-02
+1.550000000000000000e+01 5.198844816492148579e-02
+1.570000000000000107e+01 5.073844816492148468e-02
+1.590000000000000036e+01 4.948844816492148357e-02
+1.610000000000000142e+01 4.823844816492146859e-02
+1.630000000000000071e+01 4.698844816492148135e-02
+1.650000000000000000e+01 4.573844816492149412e-02
+1.670000000000000284e+01 4.448844816492147219e-02
+1.689999999999999858e+01 4.323844816492149884e-02
+1.710000000000000142e+01 4.198844816492146997e-02
+1.730000000000000071e+01 4.073844816492148274e-02
+1.750000000000000000e+01 3.948844816492149551e-02
+1.770000000000000284e+01 3.823844816492146664e-02
+1.789999999999999858e+01 3.698844816492149329e-02
+1.810000000000000142e+01 3.573890750804423200e-02
+1.830000000000000071e+01 3.449964148174038125e-02
+1.850000000000000000e+01 3.326543479507449580e-02
+1.870000000000000284e+01 3.203174952939432524e-02
+1.889999999999999858e+01 3.079467276921514507e-02
+1.910000000000000142e+01 2.955087453884154558e-02
+1.930000000000000426e+01 2.829757488742942056e-02
+1.950000000000000000e+01 2.703251082681309980e-02
+1.970000000000000284e+01 2.575391270614731457e-02
+1.989999999999999858e+01 2.446048025400142725e-02
+2.010000000000000142e+01 2.315136325026871320e-02
+2.030000000000000426e+01 2.182635155962668605e-02
+2.050000000000000000e+01 2.051131269836524787e-02
+2.070000000000000284e+01 1.924601393329058649e-02
+2.089999999999999858e+01 1.804203627945345778e-02
+2.110000000000000142e+01 1.689542978922086575e-02
+2.130000000000000426e+01 1.580294468056044413e-02
+2.150000000000000000e+01 1.476199686941700659e-02
+2.170000000000000284e+01 1.377063537642885492e-02
+2.189999999999999858e+01 1.282751146689690633e-02
+2.210000000000000142e+01 1.193185020470221951e-02
+2.230000000000000426e+01 1.108342375715123083e-02
+2.250000000000000000e+01 1.028252618964296122e-02
+2.270000000000000284e+01 9.529949800822945516e-03
+2.289999999999999858e+01 8.826962744718941725e-03
+2.310000000000000142e+01 8.174601725929608181e-03
+2.330000000000000426e+01 7.559802863360493355e-03
+2.350000000000000000e+01 6.980227972423127571e-03
+2.370000000000000284e+01 6.433697866858183699e-03
+2.389999999999999858e+01 5.918183580024932201e-03
+2.410000000000000142e+01 5.431797696638207551e-03
+2.430000000000000426e+01 4.972783058729532239e-03
+2.450000000000000000e+01 4.539507565223987323e-03
+2.470000000000000284e+01 4.130454369144049079e-03
+2.489999999999999858e+01 3.744215735107158695e-03
+2.510000000000000142e+01 3.379484616440264159e-03
+2.530000000000000426e+01 3.034735017073319965e-03
+2.550000000000000000e+01 2.708751698659090778e-03
+2.570000000000000284e+01 2.400396973682830063e-03
+2.589999999999999858e+01 2.108607205669885826e-03
+2.610000000000000142e+01 1.832387574320819192e-03
+2.630000000000000426e+01 1.570807118662753687e-03
+2.650000000000000000e+01 1.322994270503637581e-03
+2.670000000000000284e+01 1.088133472036523827e-03
+2.689999999999999858e+01 8.654613447925279569e-04
+2.710000000000000142e+01 6.542629252750590098e-04
+2.730000000000000426e+01 4.538686452942744076e-04
+2.750000000000000000e+01 2.634321949839485508e-04
+2.770000000000000284e+01 8.247267022967929770e-05
+2.789999999999999858e+01 -8.946541459389451901e-05
+2.810000000000000142e+01 -2.528133885021067941e-04
+2.830000000000000426e+01 -4.091560659871982739e-04
+2.850000000000000000e+01 -5.589405244258737783e-04
+2.870000000000000284e+01 -7.023220036435457673e-04
+2.889999999999999858e+01 -8.395902371019729055e-04
+2.910000000000000142e+01 -9.706573749801771043e-04
+2.930000000000000426e+01 -1.094885234970539205e-03
+2.950000000000000000e+01 -1.210868142349666486e-03
+2.970000000000000284e+01 -1.318294497706409092e-03
+2.989999999999999858e+01 -1.416376082145733798e-03
+3.010000000000000142e+01 -1.507897700854660128e-03
+3.030000000000000426e+01 -1.593664383871955985e-03
+3.050000000000000000e+01 -1.673736318902464002e-03
+3.070000000000000284e+01 -1.740630644905084368e-03
+3.089999999999999858e+01 -1.789005209551238103e-03
+3.110000000000000142e+01 -1.816262369186042325e-03
+3.130000000000000426e+01 -1.819590548340866906e-03
+3.150000000000000000e+01 -1.797513546768186215e-03
+3.170000000000000284e+01 -1.749719398738666274e-03
+3.189999999999999858e+01 -1.676916296299224389e-03
+3.210000000000000142e+01 -1.580713066499942572e-03
+3.229999999999999716e+01 -1.463516874192155072e-03
+3.250000000000000000e+01 -1.328120318387734922e-03
+3.270000000000000284e+01 -1.178210741843272423e-03
+3.290000000000000568e+01 -1.018015332235708715e-03
+3.310000000000000142e+01 -8.521844400039598876e-04
+3.329999999999999716e+01 -6.839823655063989897e-04
+3.350000000000000000e+01 -5.172414360357352860e-04
+3.370000000000000284e+01 -3.557896324776453623e-04
+3.390000000000000568e+01 -2.028567483376484239e-04
+3.410000000000000142e+01 -6.142784627760602327e-05
+3.429999999999999716e+01 6.543236827204639602e-05
+3.450000000000000000e+01 1.739210005797695263e-04
+3.470000000000000284e+01 2.614542459275110971e-04
+3.490000000000000568e+01 3.253756029878940802e-04
+3.510000000000000142e+01 3.686125820956898873e-04
+3.529999999999999716e+01 3.957017214793738958e-04
+3.550000000000000000e+01 4.140565726203760235e-04
+3.570000000000000284e+01 4.243128039879684715e-04
+3.590000000000000568e+01 4.271424833245643826e-04
+3.610000000000000142e+01 4.247443843440323603e-04
+3.629999999999999716e+01 4.180692204939637701e-04
+3.650000000000000000e+01 4.079833121649558856e-04
+3.670000000000000284e+01 3.952583380214028629e-04
+3.690000000000000568e+01 3.805613489470237461e-04
+3.710000000000000142e+01 3.644521186922300622e-04
+3.729999999999999716e+01 3.473838900283582418e-04
+3.750000000000000000e+01 3.297170368332490577e-04
+3.770000000000000284e+01 3.117972700258186502e-04
+3.790000000000000568e+01 2.939477526717323096e-04
+3.810000000000000142e+01 2.764309452821347751e-04
+3.830000000000000426e+01 2.594507763954954655e-04
+3.850000000000000000e+01 2.431533217796447534e-04
+3.870000000000000284e+01 2.276286476681602155e-04
+3.890000000000000568e+01 2.129166931495652846e-04
+3.910000000000000142e+01 1.990201027572733412e-04
+3.930000000000000426e+01 1.859199247042467339e-04
+3.950000000000000000e+01 1.735712963263425282e-04
+3.970000000000000284e+01 1.619215854356025129e-04
+3.990000000000000568e+01 1.509221298071244705e-04
+4.010000000000000142e+01 1.405363274813113260e-04
+4.030000000000000426e+01 1.307243348890880900e-04
+4.050000000000000000e+01 1.214435491656517510e-04
+4.070000000000000284e+01 1.126593927418429196e-04
+4.090000000000000568e+01 1.043511996648476751e-04
+4.110000000000000142e+01 9.650082646819679128e-05
+4.130000000000000426e+01 8.909086209109149572e-05
+4.150000000000000000e+01 8.210450640018338183e-05
+4.170000000000000284e+01 7.552549097972430388e-05
+4.190000000000000568e+01 6.933785642572585942e-05
+4.210000000000000142e+01 6.352631491411889881e-05
+4.230000000000000426e+01 5.807583448498591569e-05
+4.250000000000000000e+01 5.297178910661089711e-05
+4.270000000000000284e+01 4.820003894337416628e-05
+4.290000000000000568e+01 4.374650406156861210e-05
+4.310000000000000142e+01 3.959790209395871660e-05
+4.330000000000000426e+01 3.574061806543175413e-05
+4.350000000000000000e+01 3.216219288570860030e-05
+4.370000000000000284e+01 2.884952289972540562e-05
+4.390000000000000568e+01 2.579078759602929034e-05
+4.410000000000000142e+01 2.297351224541677726e-05
+4.430000000000000426e+01 2.038623200756899386e-05
+4.450000000000000000e+01 1.801711888583558746e-05
+4.470000000000000284e+01 1.585494734997432351e-05
+4.490000000000000568e+01 1.388858690309734434e-05
+4.510000000000000142e+01 1.210713181401948185e-05
+4.530000000000000426e+01 1.049999177710488341e-05
+4.550000000000000000e+01 9.056737155256784215e-06
+4.570000000000000284e+01 7.767064003557306708e-06
+4.590000000000000568e+01 6.621014064012466623e-06
+4.610000000000000142e+01 5.608581617831055883e-06
+4.630000000000000426e+01 4.720264514606751047e-06
+4.650000000000000000e+01 3.946644315468124158e-06
+4.670000000000000284e+01 3.278020415190346346e-06
+4.690000000000000568e+01 2.705258556026685530e-06
+4.710000000000000142e+01 2.218653873663369839e-06
+4.730000000000000426e+01 1.809039669579279237e-06
+4.750000000000000000e+01 1.467178681034118453e-06
+4.770000000000000284e+01 1.189494982933484530e-06
+4.790000000000000568e+01 9.667878669621938840e-07
+4.810000000000000142e+01 7.909738701770091379e-07
+4.830000000000000426e+01 6.540029661753313390e-07
+4.850000000000000000e+01 5.478563620465557584e-07
+4.870000000000000284e+01 4.645444404766585417e-07
+4.890000000000000568e+01 3.961048374423187819e-07
+4.910000000000000142e+01 3.346006465608565721e-07
+4.930000000000000426e+01 2.721187417509974269e-07
+4.950000000000000000e+01 2.007682104092737559e-07
+4.970000000000000284e+01 1.126788898204200922e-07
+4.990000000000000568e+01 0.000000000000000000e+00