Source code for mrdna.model.nbPot

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

from ..arbdmodel.interactions import AbstractPotential
from .. import get_resource_path
from ..config import CACHE_DIR

from scipy.optimize import curve_fit
from pathlib import Path

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

[docs] class AbstractNbDnaScheme(AbstractPotential):
[docs] def __init__(self, *args,**kwargs): AbstractPotential.__init__(self,*args,**kwargs) self.debye_length = self.debye_length0 # set by child class x,y = self.load_pmf() self.target_x = x self.target_y = y 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 """ try: 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
[docs] def get_package_cache_file(self): raise NotImplementedError()
[docs] def load_pmf(self): raise NotImplementedError()
[docs] def load_bead_distributions(self): raise NotImplementedError()
[docs] def potential(self, r, types): typeA, typeB = types """ Implements interface for NonbondedScheme """ if typeA.name[0] == "O" or typeB.name[0] == "O": u = np.zeros( r.shape ) else: u = self.nbPot(r, 0.5*typeA.nts, 0.5*typeB.nts) return u
""" Remainder has to do with nbPot """
[docs] def load_np(self,filename): return np.load(filename)['arr_0'].tolist()
[docs] def parametric_potential(self,x,*parms): 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<13] = u[x<13] + 0.5 * (x[x<13]-13)**2 ## Debye correction if self.debye_length != self.debye_length0: ## units "e**2/(4 * pi * 80 epsilon0 AA)" kcal_mol A = 4.1507964 with np.errstate(divide='ignore',invalid='ignore'): du = (A/x) * (np.exp(-x/self.debye_length) - np.exp(-x/self.debye_length0)) du[x < 10] = du[x>=10][0] u = u + du return u-u[-1]
[docs] 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)
[docs] 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
[docs] 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 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) """
[docs] 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
[docs] def get_raw_pressure(self): from bead_dist_data import interhelical_spacing, pressure distance = interhelical_spacing assert( len(distance) > 2 ) return distance,pressure
[docs] 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)
[docs] 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
[docs] 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
[docs] 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)
[docs] 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)
[docs] 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
[docs] def get_rounded_bp(self,bps1,bps2): larger = np.max([bps1,bps2]) smaller = np.min([bps1,bps2]) smaller_shift, larger_shift = [_round_bp(bp) for bp in (smaller,larger)] key = (smaller_shift,larger_shift) return smaller,larger,smaller_shift,larger_shift
[docs] def nbPot(self, x, bps1, bps2): 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 ): obj, dict_, key = args cache_key = '{}:{}'.format(*key) dict_[cache_key] = obj.fit_potential(*key) return True
[docs] class nbDnaScheme100Mg(AbstractNbDnaScheme):
[docs] def __init__(self,*args,**kwargs): ## units "sqrt( 80 epsilon0 295 k K / ((4*100 mM + 200 mM) e**2/particle) )" AA self.debye_length0 = 5.5771297 AbstractNbDnaScheme.__init__(self,*args,**kwargs)
[docs] def get_package_cache_file(self): return get_resource_path("nbPot.rau_refined_100Mg.dat")
[docs] 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
[docs] def load_pmf(self): d = np.loadtxt(get_resource_path("jj-nar-dna-pmf-100Mg.dat")) x0 = d[:,1] y0 = d[:,0] # kcal/mol for 10bp with 10bp, according to JY 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
[docs] 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)
[docs] class nbDnaScheme20Mg(AbstractNbDnaScheme):
[docs] def __init__(self,*args,**kwargs): ## units "sqrt( 80 epsilon0 295 k K / ((4*25 mM + 50 mM) e**2/particle) )" nm ## 25 mM MgCl2, sorry about name of this class/files self.debye_length0 = 11.154259 AbstractNbDnaScheme.__init__(self,*args,**kwargs)
[docs] def get_package_cache_file(self): return get_resource_path("nbPot.rau_refined_20Mg_200Na.dat")
[docs] 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
[docs] def load_pmf(self): d = np.loadtxt(get_resource_path("jj-nar-dna-pmf-20Mg-200Na.dat")) x0 = d[:,1] y0 = d[:,0] # kcal/mol for 10bp with 10bp, according to JY 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
[docs] def load_bead_distributions(self): import bead_dist_data as bdd 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 ) self.bead_distribution_y = interp1d( np.array(all_x), np.array(all_y), axis=0, bounds_error = False, fill_value="extrapolate", assume_sorted=True)
# nbDnaScheme100MgObj = nbDnaScheme100Mg() nbDnaScheme = nbDnaScheme20MgObj = nbDnaScheme20Mg() if __name__ == "__main__": # import cProfile # cProfile.run("_run()") # _generate_cache() # tmp = nbDnaScheme20Mg() # tmp._write_cache() # pass tmp = nbDnaScheme20Mg() # 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)