<jupyter_start><jupyter_text>Calculate radius of gyration and persistence length of PLGA/water simulations<jupyter_code>def seg_org(poly_atoms, n_monomers):
    
    pmecl = []
    count = 0
    for i in range(n_monomers):
        fgrp = []
        count += 1
        n_atm = len(poly_atoms.select_atoms("resid "+str(count)))
        
        if count == 1:
            #print(poly_atoms.select_atoms("resid "+str(count)+" and name O13"))
            pmecl.append(poly_atoms.select_atoms("resid "+str(count)+" and (name O13 O9)"))
            #print(pmecl)
        elif count == 2:
            #print(poly_atoms.select_atoms("name O5 and resid "+str(count)))
            pmecl.append(poly_atoms.select_atoms("(resid "+str(count)+" and name O5) or (resid "+str(count-1)+" and name O9)"))
        
        if count != 2 and count != 1:
            if count == n_monomers:
                pmecl.append(poly_atoms.select_atoms("(resid "+str(count)+" and name O1) or (resid "+str(count-1)+" and name O5)"))
            elif count != n_monomers:
            #print(count)
                pmecl.append(poly_atoms.select_atoms("(resid "+str(count)+" and name O5) or (resid "+str(count-1)+" and name O5)"))

    return pmecl
# write function to average across frames to give ensembled averaged cosine theta values at each N - 1 value

def pers_length(polymer_atoms, n_monomers):
    """ This function takes the polymer atoms and number of monomers and outputs the cosine theta values 
    at each bond length"""
    # N : number of monomers 
    # Initialize a zeros matrix of 3 rows, N - 1 columns to store cosine theta values at each frame 
    # row 1: Bond Length x values 
    # row 2:
    # row 3:
    vec_poly = np.zeros(shape=(3,n_monomers-1), dtype=float)
    
    # Inititalize a evenly spaced values with length of N - 1
    len_vec = np.arange(n_monomers-1)
    
    # Add 1 to all values to begin sequence from 1, instead of 0 
    len_vec += 1

    # Store bond length values on the vec_poly matrix
    vec_poly[1,:] = len_vec

    # Initialize a counter 
    count = 0 
    
    # Initialize a vector to store each cosine theta value, dtype is object because each polymer image is a sample set  
    sv_ply = np.zeros(shape=(n_monomers-1), dtype=object)
    
    
    for i in range(n_monomers):
        
        count += 1
        
        # 
        ds_cor = np.zeros(shape=(n_monomers-count))
        
        for j in range(n_monomers - count):
        
            jh = np.arange(n_monomers - count)
            
            jh += count+1

            n6_mon1 = polymer_atoms.select_atoms("resid "+str(count))

            n6_mon2 = polymer_atoms.select_atoms("resid "+str(jh[j]))

            if j == 0: 
                    
                v1 = n6_mon1.center_of_mass() - n6_mon2.center_of_mass()
                    
                v1_norm = v1/(np.linalg.norm(v1))
            
                ds_cor[j] = v1_norm.dot(v1_norm)

            elif j != 0:
                    
                v2 = n6_mon1.center_of_mass() - n6_mon2.center_of_mass()
                    
                v2_norm = v2/(np.linalg.norm(v2))
                
                ds_cor[j] = np.dot(v1_norm, v2_norm)

            sv_ply[i] = ds_cor
    
    cor_avg = []
 
    for j in range(n_monomers-1):
        lss = []
        for i in sv_ply.flat:
            try:
                lss.append(i[j])
            except IndexError:
                pass
        cor_avg.append(np.mean(lss))
    
    nm = np.array(cor_avg)

    ang_vg = []
    for i in nm.flat:
        if i >= float(0.99):
            ang_vg.append(0)
        elif i <= float(0.99):
            ang_vg.append(math.degrees(math.acos(i)))

    vec_poly[0,:] = nm 
    vec_poly[2,:] = np.array(ang_vg)
    
    return vec_poly
def get_rg_pers_poly(polymer_atoms, universe, start, end, ln_bnd):
    """This function will calculate the average radius of gyration and persistence length of a polymer within a 
    trajectory block."""
    
    n_monomers = len(np.unique(polymer_atoms.resids))
    
    rg_ens = np.zeros(shape=(1,end-start))
    
    corr_v = np.zeros(shape=(n_monomers-1,end-start))
    
    #ln_corr = np.zeros(shape=(n_monomers-1, end-start))
    
    angle_v = np.zeros(shape=(n_monomers-1,end-start))
    
    v_poly = np.zeros(shape=(4,n_monomers-1))
    
    count_rg = 0
    
    universe.trajectory[start]
    
    for ts in universe.trajectory[start:end]:
        
        p_mat = pers_length(polymer_atoms, n_monomers)
        
        #print(p_mat[0])
     
        corr_v[:,count_rg] = p_mat[0]
        
        #ln_corr[:,count_rg] = np.log(p_mat[0,:])
        
        angle_v[:,count_rg] = p_mat[2]
        
        rg_ens[0,count_rg] = polymer_atoms.radius_of_gyration()
        
        count_rg += 1 
        
        universe.trajectory[end]
    #Added the calculation of the std dev for the correlation values 
    v_poly[3,:] = p_mat[1]
    
    for i in range(n_monomers-1):
        
        v_poly[0,i] = np.mean(corr_v[i,:])
        
        #Added the calculation of the std dev for the correlation values
        v_poly[1,i] = np.std(corr_v[i,:])
        
        #v_poly[2,i] = np.mean(ln_corr[i,:])
        
        #v_poly[3,i] = np.std(ln_corr[i,:])
        
        v_poly[2,i] = np.mean(angle_v[i,:])
     
    avg_rg = np.mean(rg_ens)
    
    #def expfunc(x, a):
    #    return np.exp(-x*ln_bnd/2*a)
    # Starting with x = 0 for fits
    #tr_n6 = np.arange(n_monomers-1)
    
    #tr_n6 += 1

    #pers_lp = scipy.optimize.curve_fit(expfunc, tr_n6, v_poly[0,:])[0][0]
    
    return  rg_ens, v_poly, corr_v, avg_rg 
# I want a list of total fraction of contacts where length is determined by no. of blocks and a dictionary 
# of contact groups as keys and list of fractional contacts as values(length of list will be no. of blocks)
def bavg_pers_cnt(no_of_blks, polymer_atoms, universe, len_bnd, begin, final):
 
    n_size = (final - begin)/no_of_blks
    
    ot_dab = {}
    universe.trajectory[begin]
    
    sf_lbl = ["Avg Radius of gyration", "Avg persistence length"]

    blk_nparr = np.zeros(shape=(len(sf_lbl),no_of_blks))
    
    count = 0 
    
    for i in range(no_of_blks):
        
        start = universe.trajectory.frame
        print(start)
    
        end = int(start + n_size)
        print(end)
       
        pp_rgens, cor_pp, per_pp, rg_avgpp = get_rg_pers_poly(polymer_atoms, universe, start, end, len_bnd)
        
        blk_nparr[0,count] = rg_avgpp
        
        blk_nparr[1,count] = per_pp        
        
        universe.trajectory[end]
        
        count += 1
    # Save fractional contacts for each AA group type, each element in the value array corresponds to a block 
    # calculated value
    for i in range(len(sf_lbl)):
        ot_dab[sf_lbl[i]] = blk_nparr[i,:]
    
    return ot_dab<jupyter_output><empty_output><jupyter_text>Water only systems  N = 6 water<jupyter_code># For the right Rg calculation using MD Analysis, use trajactory without pbc 
n6_plga_wat = mda.Universe("n6_plga_50ns/0_xa_soln/confout.gro", "n6_plga_50ns/0_xa_soln/nopbc_n6wat.xtc")
len(n6_plga_wat.trajectory)
#Select the polymer heavy atoms 
poly_n6 = n6_plga_wat.select_atoms("resname sPLG PLG tPLG and not type H")
com_bond = np.zeros(shape=(1,5000))
count = 0
for ts in n6_plga_wat.trajectory[0:5000]:
    n6_mon1 = n6_plga_wat.select_atoms("resid 1")
    n6_mon2 = n6_plga_wat.select_atoms("resid 2")
    oo_len = mda.analysis.distances.distance_array(n6_mon1.center_of_mass(), n6_mon2.center_of_mass(), 
                                                   box=n6_plga_wat.trajectory.ts.dimensions)
    com_bond[0, count] = oo_len
    count += 1
    
com_bond
lb_avg = np.mean(com_bond)
lb_avg
np.std(com_bond)
plt.scatter(np.arange(5000),com_bond)
#plt.scatter(np.arange(5000),lbond,c="y")
lb_avg
n6_blks = bavg_pers_cnt(5, poly_n6, n6_plga_wat, lb_avg, 0 , 5000)
# Fully extended length of 6 monomer PLGA is 3.838 nm 
n6_blks["Avg persistence length"]
n6_blks["Avg Radius of gyration"]
np.mean(n6_blks["Avg persistence length"])
np.std(n6_blks["Avg persistence length"])
n6_rgens, cor_n6, N6ens_cos, rg_avgn6 = get_rg_pers_poly(poly_n6, n6_plga_wat, 0, 5000, lb_avg)
cor_n6
cor_n6[0]
cor_n6[1]
N6ens_cos
N6ens_cos[0]
## Bootstrapping attempt
n_iter = 3
n_size = int(len(N6ens_cos[1])*0.5)
stats_n6 = []
for i in range(n_iter):
    
    #prep train and tests values 
    train_n6 = sklearn.utils.resample(N6ens_cos[2], n_samples=n_size)
    test_n6 = np.array([x for x in N6ens_cos[2] if x.tolist() not in train.tolist()])
lb_avg = np.mean(com_bond)
#lb_avg /= 10
lb_avg
# x values
blen = np.arange(5, dtype=float)
blen *= lb_avg
#nt_tt[0] = 0
blen
# Sample standard deviation of Bond Length x values 
np.sqrt((np.sum(blen**2) - ((np.sum(blen)**2)/5))/4)
# population standard deviation of Bond Length x values 
np.std(blen)
# Sample variance of Bond Length x values 
n6svar = (np.sum(blen**2) - ((np.sum(blen)**2)/5))/4
n6svar
# numpy calculation of population variance 
np.var(blen)
# y values
n6lc = np.log(cor_n6[0])
n6lc
mk = cor_n6[1]/cor_n6[0]
mk
plt.figure(figsize=(7,7))
plt.errorbar(blen, np.log(cor_n6[0]), yerr=mk, color='b', linestyle="None",marker='o',
     capsize=5, capthick=1, ecolor='black')
# From fitting all points, I get best fit 
from sklearn.linear_model import LinearRegression
model_n6 = LinearRegression(fit_intercept=False)
model_n6.fit(blen.reshape(-1,1), n6lc)
# Slope here is in nanometers
print("Model slope:    ", model_n6.coef_)
print("Model intercept:", model_n6.intercept_)
gg = model_n6.predict(blen.reshape(-1,1))
gg
print("Mean Std Error:", sklearn.metrics.mean_squared_error(n6lc, gg))
print("R2 score:", sklearn.metrics.r2_score(n6lc, gg))
# Residuals between the true y data and model y data 
resid_n6 = n6lc - gg
resid_n6
# How to calculate mean squared error
np.sum(resid_n6**2)/len(resid_n6)
# Standard error: Square root of the mean squared error 
np.sqrt(np.sum(resid_n6**2)/len(resid_n6))
# population variance of the residuals 
var_n6p = np.var(resid_n6)
var_n6p
# Sample variance: https://www.statisticshowto.datasciencecentral.com/
# probability-and-statistics/descriptive-statistics/sample-variance/

# sample variance of the residuals 
ntres = resid_n6
ntres -= np.mean(resid_n6)
nhres = ntres**2
np.sum(nhres)/4
# How to calculate Sum((Xi - avg(X))^2): X values are the bond length values 
nt_tt = np.arange(5, dtype=float)
nt_tt *= lb_avg
nt_tt -= np.mean(nt_tt)
nhui = nt_tt**2
np.sum(nhui)
# How to calculate Sum((Xi - avg(X))^2), alternate version, Bond Length X values 
nt_tt = np.arange(5, dtype=float)
nt_tt *= lb_avg
np.sum(nt_tt**2) - ((np.sum(nt_tt)**2)/5)
# t-value with 95 % confidence intervals 
scipy.stats.t.ppf(0.975, 4)
# How to calculate 95% confidence interval for the slope 
flc_n6 = scipy.stats.t.ppf(0.975, 4)*np.sqrt((np.sum(resid_n6**2)/len(resid_n6))/(np.sum(nhui)))
flc_n6
pers_n6 =-1/model_n6.coef_
pers_n6
# Pers length error: error propagation from uncertainty in slope 
flc_n6/((model_n6.coef_)**2)
blen
plt.figure(figsize=(7,7))
plt.errorbar(blen, np.log(cor_n6[0]), yerr=mk, color='b', linestyle="None",marker='o',
     capsize=5, capthick=1, ecolor='black')
plt.plot(blen, gg, color='b')
plt.title(r'Ensemble Averaged ln(Cosine $\theta$)', fontsize=15, y=1.01)
plt.xlabel(r'Bond Length', fontsize=15)
plt.ylabel(r'ln$\left< Cos(\theta)\right >$', fontsize=15)
#plt.ylim(-1.9,0)
font = font_manager.FontProperties(family='Arial', style='normal', size='14')
plt.legend([r'$N_{PLGA}$ = 6: $L_{p}$ = 18.7 $\AA$ ± 3.5 $\AA$'], loc=3, frameon=0, fontsize=14, prop=font)
plt.tick_params(labelsize=14)
plt.text(5,-0.15,r'R$^{2}$ = 0.98', fontsize=15, color='blue')
lb_avg
exfun(nt, pers_n6)
def exfun(x, a):
    return np.exp(-x*lb_avg/a)
nt = np.arange(5)
#nt += 1
plt.scatter(nt, cor_n6[0])
# a-axis is the number of bond lengths away 
plt.plot(nt,exfun(nt, pers_n6))
plt.ylim(0,1.2)<jupyter_output><empty_output><jupyter_text>N = 8 water<jupyter_code>n8_plga_wat = mda.Universe("n8_plga_50ns/0_xa_soln/confout.gro", "n8_plga_50ns/0_xa_soln/nopbc_traj.xtc")
n8_plga_wat
#Select the polymer heavy atoms 
poly_n8 = n8_plga_wat.select_atoms("resname sPLG PLG tPLG and not type H")
len(n8_plga_wat.trajectory)
com_bond_n8 = np.zeros(shape=(1,5000))
count = 0
for ts in n8_plga_wat.trajectory[0:5000]:
    n8_mon1 = n8_plga_wat.select_atoms("resid 1")
    n8_mon2 = n8_plga_wat.select_atoms("resid 2")
    oo_len = mda.analysis.distances.distance_array(n8_mon1.center_of_mass(), n8_mon2.center_of_mass(), 
                                                   box=n8_plga_wat.trajectory.ts.dimensions)
    com_bond_n8[0, count] = oo_len
    count += 1
lb_avg_n8 = np.mean(com_bond_n8)
lb_avg_n8
plt.scatter(np.arange(5000),com_bond_n8)
plt.scatter(np.arange(5000),com_bond)
lb_avg_n8
lb_avg
n8_blks = bavg_pers_cnt(5, poly_n8, n8_plga_wat, lb_avg, 0 , 5000)
n8_blks["Avg persistence length"]
n8_blks["Avg Radius of gyration"]
n8_rgens, cor_n8_mat, N8ens_cos, rg_avgn8 = get_rg_pers_poly(poly_n8, n8_plga_wat, 0, 5000, lb_avg)
cor_n8_mat
cor_n8_mat[0]
cor_n8_mat[1]
# error prop into natural log values 
mk_n8 = cor_n8_mat[1]/cor_n8_mat[0]
#mk_n8 /= 2
mk_n8
cor_n8_mat[1]/cor_n8_mat[0]
# For some reason, fit does not work if lb_avg is in angstroms
lb_avg = np.mean(com_bond)
#lb_avg /= 10
lb_avg
# x values
nt_ttn8 = np.arange(7, dtype=float)
nt_ttn8 *= lb_avg
#nt_ttn8[0] = 0
nt_ttn8
# Sample standard deviation of Bond Length x values 
np.sqrt((np.sum(nt_ttn8**2) - ((np.sum(nt_ttn8)**2)/7))/6)
# Sample variance of Bond Length x values 
n8svar = (np.sum(nt_ttn8**2) - ((np.sum(nt_ttn8)**2)/7))/6
n8svar
# y values
n8lc = np.log(cor_n8_mat[0])
n8lc
plt.figure(figsize=(7,7))
plt.errorbar(blen, np.log(cor_n6[0]), yerr=mk, color='b', linestyle="None",marker='o',
     capsize=5, capthick=1, ecolor='black')
from sklearn.linear_model import LinearRegression
model_n8 = LinearRegression(fit_intercept=False)
model_n8.fit(nt_ttn8[:5].reshape(-1,1), n8lc[:5])
# Slope here is in nanometers
print("Model slope:    ", model_n8.coef_)
print("Model intercept:", model_n8.intercept_)
gg_n8 = model_n8.predict(nt_ttn8.reshape(-1,1))
gg_n8
print("Mean Std Error:", sklearn.metrics.mean_squared_error(n8lc[:5], gg_n8[:5]))
print("R2 score:", sklearn.metrics.r2_score(n8lc[:5], gg_n8[:5]))
# Residuals between the true y data and model y data 
resid_n8 = n8lc[:5] - gg_n8[:5]
resid_n8
# How to calculate mean squared error
np.sum(resid_n8**2)/len(resid_n8)
# Standard error: Square root of the mean squared error 
np.sqrt(np.sum(resid_n8**2)/len(resid_n8))
# sample variance of the residuals 
ntres_n8 = resid_n8
ntres_n8 -= np.mean(resid_n8)
nhres_n8 = ntres_n8**2
np.sum(nhres_n8)/4
# How to calculate Sum((Xi - avg(X))^2): X values are the bond length values 
nt_n8 = np.arange(5, dtype=float)
nt_n8 *= lb_avg
nt_n8 -= np.mean(nt_n8)
nhui = nt_n8**2
np.sum(nhui)
# How to calculate 95% confidence interval for the slope 
flc_n8 = scipy.stats.t.ppf(0.975, 6)*np.sqrt((np.sum(resid_n8**2)/len(resid_n8))/(np.sum(nhui)))
flc_n8
pers_n8 =-1/model_n8.coef_
pers_n8
# Pers length error: error propagation from uncertainty in slope 
flc_n8/((model_n8.coef_)**2)
nt_ttn8
# ln(cos theta) error: error propagation from uncertainty in cos theta values 
mk_n8 = cor_n8_mat[1]/cor_n8_mat[0]
mk_n8
plt.figure(figsize=(7,7))
plt.errorbar(nt_ttn8, np.log(cor_n8_mat[0]), yerr=mk_n8, color='b', linestyle="None",marker='o',
     capsize=5, capthick=1, ecolor='black')
plt.plot(nt_ttn8, gg_n8, color='b')
plt.title(r'Ensemble Averaged ln(Cosine $\theta$)', fontsize=15, y=1.01)
plt.xlabel(r'Bond Length', fontsize=15)
plt.ylabel(r'ln$\left< Cos(\theta)\right >$', fontsize=15)
#plt.ylim(-1.9,0)
font = font_manager.FontProperties(family='Arial', style='normal', size='14')
plt.legend([r'$N_{PLGA}$ = 8: $L_{p}$ = 17.8 $\AA$ ± 1.5 $\AA$'], loc=3, frameon=0, fontsize=14, prop=font)
plt.tick_params(labelsize=14)
plt.text(5,-0.15,r'R$^{2}$ = 0.99', fontsize=15, color='blue')
def exfun(x, a):
    return np.exp(-x*lb_avg_n8/a)

nt = np.arange(7)
plt.scatter(nt, cor_n8_mat[0,:])
plt.plot(nt,exfun(nt, pers_n8))<jupyter_output><empty_output><jupyter_text>N = 10 water<jupyter_code>n10_plga_wat = mda.Universe("n10_plga_50ns/0_xa_soln/N10_sys.gro", "n10_plga_50ns/0_xa_soln/nopbc_10wat.xtc")
n10_plga_wat
len(n10_plga_wat.trajectory)
poly_n10 = n10_plga_wat.select_atoms("resname sPLG PLG tPLG and not type H")
com_bond_n10 = np.zeros(shape=(1,5000))
count = 0
for ts in n10_plga_wat.trajectory[0:5000]:
    n10_mon1 = n10_plga_wat.select_atoms("resid 1")
    n10_mon2 = n10_plga_wat.select_atoms("resid 2")
    oo_len = mda.analysis.distances.distance_array(n10_mon1.center_of_mass(), n10_mon2.center_of_mass(), 
                                                   box=n10_plga_wat.trajectory.ts.dimensions)
    com_bond_n10[0, count] = oo_len
    count += 1
lb_avg_n10 = np.mean(com_bond_n10)
lb_avg_n10
plt.scatter(np.arange(5000),com_bond)
plt.scatter(np.arange(5000),com_bond_n8)
plt.scatter(np.arange(5000),com_bond_n10)
lb_avg
n10_blks = bavg_pers_cnt(5, poly_n10, n10_plga_wat, lb_avg_n10, 0 , 5000)
# Extended Length of 10 monomer length plga: 6.436 nm 
n10_blks["Avg persistence length"]
n10_blks["Avg Radius of gyration"]
np.mean(n10_blks["Avg Radius of gyration"])
np.std(n10_blks["Avg Radius of gyration"])
np.mean(n10_blks["Avg persistence length"])
np.std(n10_blks["Avg persistence length"])
n10_rgens, cor_n10_mat, N10ens_cos, rg_avgn10 = get_rg_pers_poly(poly_n10, n10_plga_wat, 0, 5000, lb_avg)
cor_n10_mat
cor_n10_mat[0]
cor_n10_mat[1]
mk_n10 = cor_n10_mat[1]/cor_n10_mat[0]
#mk_n10 /= 2
mk_n10
cor_n10_mat[1]/cor_n10_mat[0]
lb_avg = np.mean(com_bond)
#lb_avg /= 10
lb_avg
# x values
nt_ttn10 = np.arange(9, dtype=float)
nt_ttn10 *= lb_avg
#nt_ttn10[0] = 0
nt_ttn10
# Sample standard deviation of Bond Length x values 
np.sqrt((np.sum(nt_ttn10**2) - ((np.sum(nt_ttn10)**2)/9))/8)
# Sample variance of Bond Length x values 
n10svar = (np.sum(nt_ttn10**2) - ((np.sum(nt_ttn10)**2)/9))/8
n10svar
# y values
n10lc = np.log(cor_n10_mat[0])
n10lc
plt.figure(figsize=(7,7))
plt.errorbar(nt_ttn10, np.log(cor_n10_mat[0]), yerr=mk_n10, color='b', linestyle="None",marker='o',
     capsize=5, capthick=1, ecolor='black')
nt_ttn10[:5]
from sklearn.linear_model import LinearRegression
model_n10 = LinearRegression(fit_intercept=False)
model_n10.fit(nt_ttn10[:5].reshape(-1,1), n10lc[:5])
# Slope here is in nanometers
print("Model slope:    ", model_n10.coef_)
print("Model intercept:", model_n10.intercept_)
gg_n10 = model_n10.predict(nt_ttn10.reshape(-1,1))
gg_n10
print("Mean Std Error:", sklearn.metrics.mean_squared_error(n10lc[:5], model_n10.predict(nt_ttn10[:5].reshape(-1,1))))
print("R2 score:", sklearn.metrics.r2_score(n10lc[:5], model_n10.predict(nt_ttn10[:5].reshape(-1,1))))
# Residuals between the true y data and model y data 
resid_n10 = n10lc[:5] - gg_n10[:5]
resid_n10
# How to calculate mean squared error
np.sum(resid_n10**2)/len(resid_n10)
# Standard error: Square root of the mean squared error 
np.sqrt(np.sum(resid_n10**2)/len(resid_n10))
# sample variance of the residuals 
ntres_n10 = resid_n10
ntres_n10 -= np.mean(resid_n10)
nhres_n10 = ntres_n10**2
np.sum(nhres_n10)/4
# How to calculate Sum((Xi - avg(X))^2): X values are the bond length values 
nt_n10 = np.arange(5, dtype=float)
nt_n10 *= lb_avg
nt_n10 -= np.mean(nt_n10)
nhui = nt_n10**2
np.sum(nhui)
# How to calculate 95% confidence interval for the slope 
flc_n10 = scipy.stats.t.ppf(0.975, 4)*np.sqrt((np.sum(resid_n10**2)/len(resid_n10))/(np.sum(nhui)))
flc_n10
pers_n10 =-1/model_n10.coef_
pers_n10
# Pers length error: error propagation from uncertainty in slope 
flc_n10/((model_n10.coef_)**2)
mk_n10
nt_ttn10
np.log(cor_n10_mat[0])
plt.figure(figsize=(7,7))
plt.errorbar(nt_ttn10, np.log(cor_n10_mat[0]), yerr=mk_n10, color='b', linestyle="None",marker='o', capsize=5, capthick=1, ecolor='black')
plt.plot(nt_ttn10, gg_n10, color='b')
plt.title(r'Ensemble Averaged ln(Cosine $\theta$)', fontsize=15, y=1.01)
plt.xlabel(r'Bond Length', fontsize=15)
plt.ylabel(r'ln$\left< Cos(\theta)\right >$', fontsize=15)
#plt.ylim(-1.9,0)
font = font_manager.FontProperties(family='Arial', style='normal', size='14')
plt.legend([r'$N_{PLGA}$ = 10: $L_{p}$ = 14.7 $\AA$ ± 3.4 $\AA$'], loc=3, frameon=0, fontsize=14, prop=font)
plt.tick_params(labelsize=14)
plt.text(5,-0.15,r'R$^{2}$ = 0.97', fontsize=15, color='blue')
def exfun(x, a):
    return np.exp(-x*lb_avg/a)

nt = np.arange(9)
plt.scatter(nt, cor_n10_mat[0,:])
plt.plot(nt,exfun(nt, pers_n10))<jupyter_output><empty_output><jupyter_text>N = 20 water<jupyter_code>n20_plga_wat = mda.Universe("n20_plga_150ns/0_xa_soln/nowat_n20.gro", "n20_plga_150ns/0_xa_soln/N20_nopbcpp.xtc")
n20_plga_wat
len(n20_plga_wat.trajectory)
poly_n20 = n20_plga_wat.select_atoms("resname sPLG PLG tPLG and not type H")
com_bond_n20 = np.zeros(shape=(1,15000))
count = 0
for ts in n20_plga_wat.trajectory[0:15000]:
    n20_mon1 = n20_plga_wat.select_atoms("resid 1")
    n20_mon2 = n20_plga_wat.select_atoms("resid 2")
    oo_len = mda.analysis.distances.distance_array(n20_mon1.center_of_mass(), n20_mon2.center_of_mass(), 
                                                   box=n20_plga_wat.trajectory.ts.dimensions)
    com_bond_n20[0, count] = oo_len
    count += 1
lb_avg_n20 = np.mean(com_bond_n20)
lb_avg_n20
plt.scatter(np.arange(5000),com_bond)
plt.scatter(np.arange(5000),com_bond_n8)
plt.scatter(np.arange(5000),com_bond_n10)
plt.scatter(np.arange(15000),com_bond_n20)
n20_blks = bavg_pers_cnt(5, poly_n20, n20_plga_wat, lb_avg, 0 , 15000)
n20_blks["Avg persistence length"]
n20_blks["Avg Radius of gyration"]
np.mean(n20_blks["Avg Radius of gyration"])
np.std(n20_blks["Avg Radius of gyration"])
np.mean(n20_blks["Avg persistence length"])
np.std(n20_blks["Avg persistence length"])
n20_rgens, cor_n20_mat, N20ens_cos, rg_avgn20 = get_rg_pers_poly(poly_n20, n20_plga_wat, 0, 15000, lb_avg)
cor_n20_mat
cor_n20_mat[0]
cor_n20_mat[1]
lb_avg = np.mean(com_bond)
#lb_avg /= 10
lb_avg
# x values
nt_ttn20 = np.arange(19, dtype=float)
nt_ttn20 *= lb_avg
#nt_ttn20[0] = 0
nt_ttn20
# Sample standard deviation of Bond Length x values 
np.sqrt((np.sum(nt_ttn20**2) - ((np.sum(nt_ttn20)**2)/19))/18)
# Sample variance of Bond Length x values 
n20svar = (np.sum(nt_ttn20**2) - ((np.sum(nt_ttn20)**2)/19))/18
n20svar
# y values
n20lc = np.log(cor_n20_mat[0])
n20lc
mk_n20 = cor_n20_mat[1]/cor_n20_mat[0]
mk_n20
plt.figure(figsize=(7,7))
plt.errorbar(nt_ttn20, np.log(cor_n20_mat[0]), yerr=mk_n20, color='b', linestyle="None",marker='o',
     capsize=5, capthick=1, ecolor='black')
nt_ttn20[:8]
from sklearn.linear_model import LinearRegression
model_n20 = LinearRegression(fit_intercept=False)
model_n20.fit(nt_ttn20[:8].reshape(-1,1), n20lc[:8])
# Slope here is in nanometers
print("Model slope:    ", model_n20.coef_)
print("Model intercept:", model_n20.intercept_)
gg_n20 = model_n20.predict(nt_ttn20.reshape(-1,1))
gg_n20
print("Mean Std Error:", sklearn.metrics.mean_squared_error(n20lc[:8], model_n20.predict(nt_ttn20[:8].reshape(-1,1))))
print("R2 score:", sklearn.metrics.r2_score(n20lc[:8], model_n20.predict(nt_ttn20[:8].reshape(-1,1))))
# Residuals between the true y data and model y data 
resid_n20 = n20lc[:8] - gg_n20[:8]
resid_n20
# How to calculate mean squared error
np.sum(resid_n20**2)/len(resid_n20)
# Standard error: Square root of the mean squared error 
np.sqrt(np.sum(resid_n20**2)/len(resid_n20))
# sample variance of the residuals 
ntres_n20 = resid_n20
ntres_n20 -= np.mean(resid_n20)
nhres_n20 = ntres_n20**2
np.sum(nhres_n20)/7
# How to calculate Sum((Xi - avg(X))^2): X values are the bond length values 
nt_n20 = np.arange(8, dtype=float)
nt_n20 *= lb_avg
nt_n20 -= np.mean(nt_n20)
nhui = nt_n20**2
np.sum(nhui)
# How to calculate 95% confidence interval for the slope 
flc_n20 = scipy.stats.t.ppf(0.975, 4)*np.sqrt((np.sum(resid_n20**2)/len(resid_n20))/(np.sum(nhui)))
flc_n20
pers_n20 =-1/model_n20.coef_
pers_n20
# Pers length error: error propagation from uncertainty in slope 
flc_n20/((model_n20.coef_)**2)
mk_n20 = cor_n20_mat[1]/cor_n20_mat[0]
plt.figure(figsize=(7,7))
plt.errorbar(nt_ttn20, np.log(cor_n20_mat[0]), yerr=mk_n20, color='b', linestyle="None",marker='o',capsize=5, capthick=1, ecolor='black')
plt.plot(nt_ttn20, gg_n20, color='b')
plt.title(r'Ensemble Averaged ln(Cosine $\theta$)', fontsize=15, y=1.01)
plt.xlabel(r'Bond Length', fontsize=15)
plt.ylabel(r'ln$\left< Cos(\theta)\right >$', fontsize=15)
#plt.ylim(-1.9,0)
font = font_manager.FontProperties(family='Arial', style='normal', size='14')
plt.legend([r'$N_{PLGA}$ = 20: $L_{p}$ = 19.1 $\AA$ ± 1.0 $\AA$'], loc=3, frameon=0, fontsize=14, prop=font)
plt.tick_params(labelsize=14)
plt.text(5,-0.15,r'R$^{2}$ = 0.99', fontsize=15, color='blue')
def exfun(x, a):
    return np.exp(-x*lb_avg/a)

nt = np.arange(19)
plt.scatter(nt, cor_n20_mat[0,:])
plt.plot(nt,exfun(nt, pers_n10))<jupyter_output><empty_output><jupyter_text>N= 30 water<jupyter_code>n30_plga_wat = mda.Universe("n30_plga_150ns/0_xa_soln/0xa_n30sys.gro", "n30_plga_150ns/0_xa_soln/N30_nopbcpp_0xa.xtc")
n30_plga_wat
len(n30_plga_wat.trajectory)
poly_n30 = n30_plga_wat.select_atoms("resname sPLG PLG tPLG and not type H")
com_bond_n30 = np.zeros(shape=(1,15000))
count = 0
for ts in n30_plga_wat.trajectory[0:15000]:
    n30_mon1 = n30_plga_wat.select_atoms("resid 1")
    n30_mon2 = n30_plga_wat.select_atoms("resid 2")
    oo_len = mda.analysis.distances.distance_array(n30_mon1.center_of_mass(), n30_mon2.center_of_mass(), 
                                                   box=n30_plga_wat.trajectory.ts.dimensions)
    com_bond_n30[0, count] = oo_len
    count += 1
lb_avg_n30 = np.mean(com_bond_n30)
lb_avg_n30
lb_avg
n30_blks = bavg_pers_cnt(5, poly_n30, n30_plga_wat, lb_avg, 0 , 15000)
n30_blks["Avg persistence length"]
n30_blks["Avg Radius of gyration"]
np.mean(n30_blks["Avg persistence length"])
np.std(n30_blks["Avg persistence length"])
n30_rgens, cor_n30_mat, N30ens_cos, rg_avgn30 = get_rg_pers_poly(poly_n30, n30_plga_wat, 0, 15000, lb_avg)
cor_n30_mat
cor_n30_mat[0]
# error propagation into natural log values 
mk_n30 = cor_n30_mat[1]/cor_n30_mat[0]
# For some reason, fit does not work if lb_avg is in angstroms
lb_avg = np.mean(com_bond)
#lb_avg /= 10
lb_avg
# x values
nt_ttn30 = np.arange(29, dtype=float)
nt_ttn30 *= lb_avg
#nt_ttn10[0] = 0
nt_ttn30
nt_ttn30[:20]
# Sample standard deviation of Bond Length x values 
np.sqrt((np.sum(nt_ttn30**2) - ((np.sum(nt_ttn30)**2)/29))/28)
# Sample variance of Bond Length x values 
n30svar = (np.sum(nt_ttn30**2) - ((np.sum(nt_ttn30)**2)/29))/28
n30svar
# y values
n30lc = np.log(cor_n30_mat[0])
n30lc
mk_n30
plt.figure(figsize=(7,7))
plt.errorbar(nt_ttn30, np.log(cor_n30_mat[0]), yerr=mk_n30, color='b', linestyle="None",marker='o',
     capsize=5, capthick=1, ecolor='black')
n30lc[:6]
nt_ttn30[:6]
from sklearn.linear_model import LinearRegression
model_n30 = LinearRegression(fit_intercept=False)
model_n30.fit(nt_ttn30[:6].reshape(-1,1), n30lc[:6])
# Slope here is in angstroms
print("Model slope:    ", model_n30.coef_)
print("Model intercept:", model_n30.intercept_)
gg_n30 = model_n30.predict(nt_ttn30.reshape(-1,1))
gg_n30
print("Mean Std Error:", sklearn.metrics.mean_squared_error(n30lc[:6], model_n30.predict(nt_ttn30[:6].reshape(-1,1))))
print("R2 score:", sklearn.metrics.r2_score(n30lc[:6], model_n30.predict(nt_ttn30[:6].reshape(-1,1))))
# Residuals between the true y data and model y data 
resid_n30 = n30lc[:6] - gg_n30[:6]
resid_n30
# How to calculate mean squared error
np.sum(resid_n30**2)/len(resid_n30)
# Standard error: Square root of the mean squared error 
np.sqrt(np.sum(resid_n30**2)/len(resid_n30))
# sample variance of the residuals 
ntres_n30 = resid_n30
ntres_n30 -= np.mean(resid_n30)
nhres_n30 = ntres_n30**2
# divide by n-1 samples
np.sum(nhres_n30)/5
# How to calculate Sum((Xi - avg(X))^2): X values are the bond length values 
nt_n30 = np.arange(6, dtype=float)
nt_n30 *= lb_avg
nt_n30 -= np.mean(nt_n30)
nhui = nt_n30**2
np.sum(nhui)
# How to calculate 95% confidence interval for the slope 
flc_n30 = scipy.stats.t.ppf(0.975, 4)*np.sqrt((np.sum(resid_n30**2)/len(resid_n30))/(np.sum(nhui)))
flc_n30
pers_n30 =-1/model_n30.coef_
pers_n30
# Pers length error: error propagation from uncertainty in slope 
flc_n30/((model_n30.coef_)**2)
mk_n30 = cor_n30_mat[1]/cor_n30_mat[0]
def slope_plen(x,m):
    return m*x

ttsam = np.arange(12, dtype=float)
ttsam *= lb_avg
#nt_ttn10[0] = 0
ttsam
model_n10.coef_[0]
slope_plen(ttsam, model_n10.coef_[0])
blen
np.log(cor_n6[0])
mk
plt.figure(figsize=(7,7))
plt.errorbar(nt_ttn30, np.log(cor_n30_mat[0]), yerr=mk_n30, color='b', linestyle="None",marker='h',
     capsize=5, capthick=1, ecolor='black')
plt.plot(ttsam, slope_plen(ttsam, model_n30.coef_[0]), color='b')
plt.errorbar(nt_ttn20, np.log(cor_n20_mat[0]), yerr=mk_n20, color='g', linestyle="None"
             ,marker='*',capsize=5, capthick=1, ecolor='black')
plt.plot(ttsam, slope_plen(ttsam, model_n20.coef_[0]), color='g')
plt.errorbar(nt_ttn10, np.log(cor_n10_mat[0]), yerr=mk_n10, color='r', linestyle="None"
             ,marker='v', capsize=5, capthick=1, ecolor='black')
plt.plot(ttsam, slope_plen(ttsam, model_n10.coef_[0]), color='r')
plt.errorbar(nt_ttn8, np.log(cor_n8_mat[0]), yerr=mk_n8, color='k', linestyle="None",marker='s'
             ,capsize=5, capthick=1, ecolor='black')
plt.plot(ttsam, slope_plen(ttsam, model_n8.coef_[0]), color='k')
plt.errorbar(blen, np.log(cor_n6[0]), yerr=mk, color='c', linestyle="None",marker='x'
             ,capsize=5, capthick=1, ecolor='black')
plt.plot(ttsam, slope_plen(ttsam, model_n6.coef_[0]), color='m')
plt.title(r'Ensemble Averaged ln(Cosine $\theta$)', fontsize=15, y=1.01)
plt.xlabel(r'Bond Length', fontsize=15)
plt.ylabel(r'ln$\left< Cos(\theta)\right >$', fontsize=15)
plt.ylim(-6, 1)
plt.xlim(0,60)
font = font_manager.FontProperties(family='Arial', style='normal', size='14')
#plt.legend([r'$N_{PLGA}$ = 30: $L_{p}$ = 18.5 $\AA$ ± 0.15 $\AA$,  R$^{2}$ = 0.99', 
#r'$N_{PLGA}$ = 20: $L_{p}$ = 19.1 $\AA$ ± 0.27 $\AA$,  R$^{2}$ = 0.99'], 
#           loc=3, frameon=0, fontsize=14, prop=font)
#plt.legend([r'$N_{PLGA}$ = 20: $L_{p}$ = 19.1 $\AA$ ± 0.27 $\AA$,  R$^{2}$ = 0.99'], 
#loc=3, frameon=0, fontsize=14, prop=font)
plt.tick_params(labelsize=14)
plt.text(1, -5.9,r'$N_{PLGA}$ = 30: $L_{p}$ = 18.5 $\AA$ ± 1.6 $\AA$, R$^{2}$ = 0.99', fontsize=14, color='blue')
plt.text(1,-5.6,r'$N_{PLGA}$ = 20: $L_{p}$ = 19.1 $\AA$ ± 1.0 $\AA$, R$^{2}$ = 0.99', fontsize=14, color='green')
plt.text(1,-5.3,r'$N_{PLGA}$ = 10: $L_{p}$ = 14.7 $\AA$ ± 3.4 $\AA$, R$^{2}$ = 0.97', fontsize=14, color='red')
plt.text(1,-5.0,r'$N_{PLGA}$ = 8: $L_{p}$ = 17.8 $\AA$ ± 1.5, R$^{2}$ = 0.99', fontsize=14, color='black')
plt.text(1,-4.7,r'$N_{PLGA}$ = 6: $L_{p}$ = 18.7 $\AA$ ± 3.5 $\AA$, R$^{2}$ = 0.98', fontsize=14, color='magenta')
#plt.text(2.5,0,r'R$^{2}$ = 0.96', fontsize=15, color='red')
cor_n8_mat[0]
nt_ttn8
def exfun(x, a):
    return np.exp(-x*lb_avg/a)

nt = np.arange(29)
plt.scatter(nt, cor_n30_mat[0,:])
plt.plot(nt,exfun(nt, pers_n30))<jupyter_output><empty_output><jupyter_text>Acetone systems only  N = 6 acetone<jupyter_code>n6_plga_ace = mda.Universe("n6_plga_50ns/1_xa_soln/confout.gro", "n6_plga_50ns/1_xa_soln/nopbc_n6ace.xtc")
n6_plga_ace
len(n6_plga_ace.trajectory)
poly_n6_ace = n6_plga_ace.select_atoms("resname sPLG PLG tPLG and not type H")
com_bond_n6ace = np.zeros(shape=(1,5000))
count = 0
for ts in n6_plga_ace.trajectory[0:5000]:
    n6_mon1 = n6_plga_ace.select_atoms("resid 1")
    n6_mon2 = n6_plga_ace.select_atoms("resid 2")
    oo_len = mda.analysis.distances.distance_array(n6_mon1.center_of_mass(), n6_mon2.center_of_mass(), 
                                                   box=n6_plga_ace.trajectory.ts.dimensions)
    com_bond_n6ace[0, count] = oo_len
    count += 1
lb_avg_n6ace = np.mean(com_bond_n6ace)
lb_avg_n6ace
plt.scatter(np.arange(5000),com_bond_n6ace)

n6_blks_ace = bavg_pers_cnt(5, poly_n6_ace, n6_plga_ace, lb_avg_n6ace, 0 , 5000)
n6_blks_ace["Avg persistence length"]
n6_blks_ace["Avg Radius of gyration"]
np.mean(n6_blks_ace["Avg persistence length"])
np.std(n6_blks_ace["Avg persistence length"])
np.mean(n6_blks_ace["Avg Radius of gyration"])
np.std(n6_blks_ace["Avg Radius of gyration"])
n6_rgensace, cor_n6ace_mat, N30ens_cos, rg_avgn6ace = get_rg_pers_poly(poly_n6_ace, n6_plga_ace, 0, 5000, 
                                                                       lb_avg_n6ace)
cor_n6ace_mat
cor_n6ace_mat[0]
# x values
ace_ttn6 = np.arange(5, dtype=float)
ace_ttn6 *= lb_avg_n6ace
#nt_ttn20[0] = 0
ace_ttn6
# Sample standard deviation of Bond Length x values 
np.sqrt((np.sum(ace_ttn6**2) - ((np.sum(ace_ttn6)**2)/5))/4)
# Sample variance of Bond Length x values 
n6sace_var = (np.sum(ace_ttn6**2) - ((np.sum(ace_ttn6)**2)/5))/4
n6sace_var
# y values
n6lc_ace = np.log(cor_n6ace_mat[0])
n6lc_ace
mk_n6ace = cor_n6ace_mat[1]/cor_n6ace_mat[0]
plt.figure(figsize=(7,7))
plt.errorbar(ace_ttn6, np.log(cor_n6ace_mat[0]), yerr=mk_n6ace, color='b', linestyle="None",marker='o',
     capsize=5, capthick=1, ecolor='black')
from sklearn.linear_model import LinearRegression
model_n6ace = LinearRegression(fit_intercept=False)
model_n6ace.fit(ace_ttn6.reshape(-1,1), n6lc_ace)
# Slope here is in angstroms
print("Model slope:    ", model_n6ace.coef_)
print("Model intercept:", model_n6ace.intercept_)
gg_n6ace = model_n6ace.predict(ace_ttn6.reshape(-1,1))
gg_n6ace
print("Mean Std Error:", sklearn.metrics.mean_squared_error(n6lc_ace, model_n6ace.predict(ace_ttn6.reshape(-1,1))))
print("R2 score:", sklearn.metrics.r2_score(n6lc_ace, model_n6ace.predict(ace_ttn6.reshape(-1,1))))
# Residuals between the true y data and model y data 
resid_n6ace = n6lc_ace - gg_n6ace
resid_n6ace
# How to calculate mean squared error
np.sum(resid_n6ace**2)/len(resid_n6ace)
# Standard error: Square root of the mean squared error 
np.sqrt(np.sum(resid_n6ace**2)/len(resid_n6ace))
# sample variance of the residuals 
ntr_n6ace = resid_n6ace
ntr_n6ace -= np.mean(resid_n6ace)
nhr_n6ace = ntr_n6ace**2
np.sum(nhr_n6ace)/4
# How to calculate Sum((Xi - avg(X))^2): X values are the bond length values 
nt_n6ace = np.arange(5, dtype=float)
nt_n6ace *= lb_avg_n6ace
nt_n6ace -= np.mean(nt_n6ace)
nhui = nt_n6ace**2
np.sum(nhui)
# How to calculate 95% confidence interval for the slope 
flc_n6ace = scipy.stats.t.ppf(0.975, 4)*np.sqrt((np.sum(resid_n6ace**2)/len(resid_n6ace))/(np.sum(nhui)))
flc_n6ace
pers_n6ace =-1/model_n6ace.coef_
pers_n6ace
# Pers length error: error propagation from uncertainty in slope 
flc_n6ace/((model_n6ace.coef_)**2)
def slope_plen(x,m):
    return m*x

ttsam = np.arange(5, dtype=float)
ttsam *= lb_avg_n6ace
#nt_ttn10[0] = 0
ttsam
plt.figure(figsize=(7,7))
plt.errorbar(ace_ttn6, np.log(cor_n6ace_mat[0]), yerr=mk_n6ace, color='b', linestyle="None",marker='o',capsize=5, capthick=1, ecolor='black')
plt.plot(ttsam, slope_plen(ttsam, model_n6ace.coef_[0]), color='b')
plt.title(r'Ensemble Averaged ln(Cosine $\theta$) in pure acetone', fontsize=15, y=1.01)
plt.xlabel(r'Bond Length', fontsize=15)
plt.ylabel(r'ln$\left< Cos(\theta)\right >$', fontsize=15)
#plt.ylim(-1.9,0)
font = font_manager.FontProperties(family='Arial', style='normal', size='14')
plt.legend([r'$N_{PLGA}$ = 6: $L_{p}$ = 41.6 $\AA$ ± 5.7 $\AA$'], loc=3, frameon=0, fontsize=14, prop=font)
plt.tick_params(labelsize=14)
plt.text(5,0,r'R$^{2}$ = 0.99', fontsize=15, color='blue')
def exfun(x, a):
    return np.exp(-x*lb_avg_n6ace/a)

nt = np.arange(5)
plt.scatter(nt, cor_n6ace_mat[0,:])
plt.plot(nt,exfun(nt, pers_n6ace))<jupyter_output><empty_output><jupyter_text>N = 8 acetone<jupyter_code>n8_plga_ace = mda.Universe("n8_plga_50ns/1_xa_soln/confout.gro", "n8_plga_50ns/1_xa_soln/nopbc_n8ace.xtc")
n8_plga_ace
len(n8_plga_ace.trajectory)
poly_n8_ace = n8_plga_ace.select_atoms("resname sPLG PLG tPLG and not type H")
com_bond_n8ace = np.zeros(shape=(1,5000))
count = 0
for ts in n8_plga_ace.trajectory[0:5000]:
    n8_mon1 = n8_plga_ace.select_atoms("resid 1")
    n8_mon2 = n8_plga_ace.select_atoms("resid 2")
    oo_len = mda.analysis.distances.distance_array(n8_mon1.center_of_mass(), n8_mon2.center_of_mass(), 
                                                   box=n8_plga_ace.trajectory.ts.dimensions)
    com_bond_n8ace[0, count] = oo_len
    count += 1
lb_avg_n8ace = np.mean(com_bond_n8ace)
lb_avg_n8ace
plt.scatter(np.arange(5000),com_bond_n6ace)
plt.scatter(np.arange(5000),com_bond_n8ace)
n8_blks_ace = bavg_pers_cnt(5, poly_n8_ace, n8_plga_ace, lb_avg_n6ace, 0 , 5000)
n8_blks_ace["Avg persistence length"]
n8_blks_ace["Avg Radius of gyration"]
np.mean(n8_blks_ace["Avg persistence length"])
np.std(n8_blks_ace["Avg persistence length"])
np.mean(n8_blks_ace["Avg Radius of gyration"])
np.std(n8_blks_ace["Avg Radius of gyration"])
n8_rgens_ace, cor_n8ace_mat, N8aceens_cos, rg_avgn8ace = get_rg_pers_poly(poly_n8_ace, n8_plga_ace, 0, 5000, 
                                                                        lb_avg_n6ace)
cor_n8ace_mat
cor_n8ace_mat[0]
# x values
ace_ttn8 = np.arange(7, dtype=float)
ace_ttn8 *= lb_avg_n6ace
#nt_ttn20[0] = 0
ace_ttn8
# Sample standard deviation of Bond Length x values 
np.sqrt((np.sum(ace_ttn8**2) - ((np.sum(ace_ttn8)**2)/7))/6)
# Sample variance of Bond Length x values 
n8sace_var = (np.sum(ace_ttn8**2) - ((np.sum(ace_ttn8)**2)/7))/6
n8sace_var
# y values
n8lc_ace = np.log(cor_n8ace_mat[0])
n8lc_ace
mk_n8ace = cor_n8ace_mat[1]/cor_n8ace_mat[0]
plt.figure(figsize=(7,7))
plt.errorbar(ace_ttn8, np.log(cor_n8ace_mat[0]), yerr=mk_n8ace, color='b', linestyle="None",marker='o',
     capsize=5, capthick=1, ecolor='black')
ace_ttn8[1:]
from sklearn.linear_model import LinearRegression
model_n8ace = LinearRegression(fit_intercept=False)
model_n8ace.fit(ace_ttn8.reshape(-1,1), n8lc_ace)
# Slope here is in angstroms
print("Model slope:    ", model_n8ace.coef_)
print("Model intercept:", model_n8ace.intercept_)
gg_n8ace = model_n8ace.predict(ace_ttn8.reshape(-1,1))
gg_n8ace
print("Mean Std Error:", sklearn.metrics.mean_squared_error(n8lc_ace, model_n8ace.predict(ace_ttn8.reshape(-1,1))))
print("R2 score:", sklearn.metrics.r2_score(n8lc_ace, model_n8ace.predict(ace_ttn8.reshape(-1,1))))
# Residuals between the true y data and model y data 
resid_n8ace = n8lc_ace - gg_n8ace
resid_n8ace
# How to calculate mean squared error
np.sum(resid_n8ace**2)/len(resid_n8ace)
# Standard error: Square root of the mean squared error 
np.sqrt(np.sum(resid_n8ace**2)/len(resid_n8ace))
# sample variance of the residuals 
ntr_n8ace = resid_n8ace
ntr_n8ace -= np.mean(resid_n8ace)
nhr_n8ace = ntr_n8ace**2
np.sum(nhr_n8ace)/6
# How to calculate Sum((Xi - avg(X))^2): X values are the bond length values 
nt_n8ace = np.arange(7, dtype=float)
nt_n8ace *= lb_avg_n6ace
nt_n8ace -= np.mean(nt_n8ace)
nhui = nt_n8ace**2
np.sum(nhui)
# How to calculate 95% confidence interval for the slope 
flc_n8ace = scipy.stats.t.ppf(0.975, 4)*np.sqrt((np.sum(resid_n8ace**2)/len(resid_n8ace))/(np.sum(nhui)))
flc_n8ace
pers_n8ace =-1/model_n8ace.coef_
pers_n8ace
# Pers length error: error propagation from uncertainty in slope 
flc_n8ace/((model_n8ace.coef_)**2)
def slope_plen(x,m):
    return m*x

ttsam = np.arange(7, dtype=float)
ttsam *= lb_avg_n6ace
#nt_ttn10[0] = 0
ttsam
plt.figure(figsize=(7,7))
plt.errorbar(ace_ttn8, np.log(cor_n8ace_mat[0]), yerr=mk_n8ace, color='b', linestyle="None",marker='o',capsize=5, capthick=1, ecolor='black')
plt.plot(ttsam, slope_plen(ttsam, model_n8ace.coef_[0]), color='b')
plt.title(r'Ensemble Averaged ln(Cosine $\theta$) in pure acetone', fontsize=15, y=1.01)
plt.xlabel(r'Bond Length', fontsize=15)
plt.ylabel(r'ln$\left< Cos(\theta)\right >$', fontsize=15)
#plt.ylim(-1.9,0)
font = font_manager.FontProperties(family='Arial', style='normal', size='14')
plt.legend([r'$N_{PLGA}$ = 8: $L_{p}$ = 39.3 $\AA$ ± 5.9 $\AA$'], loc=3, frameon=0, fontsize=14, prop=font)
plt.tick_params(labelsize=14)
plt.text(5,0,r'R$^{2}$ = 0.98', fontsize=15, color='blue')
def exfun(x, a):
    return np.exp(-x*lb_avg_n6ace/a)

nt = np.arange(7)
plt.scatter(nt, cor_n8ace_mat[0,:])
plt.plot(nt,exfun(nt, pers_n8ace))<jupyter_output><empty_output><jupyter_text>N = 10 acetone<jupyter_code>n10_plga_ace = mda.Universe("n10_plga_50ns/1_xa_soln/confout.gro", "n10_plga_50ns/1_xa_soln/nopbc_n10ace.xtc")
n10_plga_ace
len(n10_plga_ace.trajectory)
poly_n10_ace = n10_plga_ace.select_atoms("resname sPLG PLG tPLG and not type H")
com_bond_n10ace = np.zeros(shape=(1,5000))
count = 0
for ts in n10_plga_ace.trajectory[0:5000]:
    n10_mon1 = n10_plga_ace.select_atoms("resid 1")
    n10_mon2 = n10_plga_ace.select_atoms("resid 2")
    oo_len = mda.analysis.distances.distance_array(n10_mon1.center_of_mass(), n10_mon2.center_of_mass(), 
                                                   box=n10_plga_ace.trajectory.ts.dimensions)
    com_bond_n10ace[0, count] = oo_len
    count += 1
lb_avg_n10ace = np.mean(com_bond_n10ace)
lb_avg_n10ace
plt.scatter(np.arange(5000),com_bond_n6ace)
plt.scatter(np.arange(5000),com_bond_n8ace)
plt.scatter(np.arange(5000),com_bond_n10ace)
n10_blks_ace = bavg_pers_cnt(5, poly_n10_ace, n10_plga_ace, lb_avg_n10ace, 0 , 5000)
n10_blks_ace["Avg persistence length"]
n10_blks_ace["Avg Radius of gyration"]
np.mean(n10_blks_ace["Avg persistence length"])
np.std(n10_blks_ace["Avg persistence length"])
np.mean(n10_blks_ace["Avg Radius of gyration"])
np.std(n10_blks_ace["Avg Radius of gyration"])
n10_rgens_ace, cor_n10ace_mat, N10aceens_cos, rg_avgn10ace = get_rg_pers_poly(poly_n10_ace, n10_plga_ace, 0, 5000, 
                                                                        lb_avg_n6ace)
cor_n10ace_mat
cor_n10ace_mat[0]
# x values
ace_ttn10 = np.arange(9, dtype=float)
ace_ttn10 *= lb_avg_n6ace
#nt_ttn20[0] = 0
ace_ttn10
# Sample standard deviation of Bond Length x values 
np.sqrt((np.sum(ace_ttn10**2) - ((np.sum(ace_ttn10)**2)/9))/8)
# Sample variance of Bond Length x values 
n10sace_var = (np.sum(ace_ttn10**2) - ((np.sum(ace_ttn10)**2)/9))/8
n10sace_var
# y values
n10lc_ace = np.log(cor_n10ace_mat[0])
n10lc_ace
mk_n10ace = cor_n10ace_mat[1]/cor_n10ace_mat[0]
plt.figure(figsize=(7,7))
plt.errorbar(ace_ttn10, np.log(cor_n10ace_mat[0]), yerr=mk_n10ace, color='b', linestyle="None",marker='o',
     capsize=5, capthick=1, ecolor='black')
from sklearn.linear_model import LinearRegression
model_n10ace = LinearRegression(fit_intercept=False)
model_n10ace.fit(ace_ttn10.reshape(-1,1), n10lc_ace)
# Slope here is in angstroms
print("Model slope:    ", model_n10ace.coef_)
print("Model intercept:", model_n10ace.intercept_)
gg_n10ace = model_n10ace.predict(ace_ttn10.reshape(-1,1))
gg_n10ace
print("Mean Std Error:", sklearn.metrics.mean_squared_error(n10lc_ace, model_n10ace.predict(ace_ttn10.reshape(-1,1))))
print("R2 score:", sklearn.metrics.r2_score(n10lc_ace, model_n10ace.predict(ace_ttn10.reshape(-1,1))))
# Residuals between the true y data and model y data 
resid_n10ace = n10lc_ace - gg_n10ace
resid_n10ace
# How to calculate mean squared error
np.sum(resid_n10ace**2)/len(resid_n10ace)
# Standard error: Square root of the mean squared error 
np.sqrt(np.sum(resid_n10ace**2)/len(resid_n10ace))
# sample variance of the residuals 
ntr_n10ace = resid_n10ace
ntr_n10ace -= np.mean(resid_n10ace)
nhr_n10ace = ntr_n10ace**2
np.sum(nhr_n10ace)/8
# How to calculate Sum((Xi - avg(X))^2): X values are the bond length values 
nt_n10ace = np.arange(9, dtype=float)
nt_n10ace *= lb_avg_n6ace
nt_n10ace -= np.mean(nt_n10ace)
nhui = nt_n10ace**2
np.sum(nhui)
# How to calculate 95% confidence interval for the slope 
flc_n10ace = scipy.stats.t.ppf(0.975, 4)*np.sqrt((np.sum(resid_n10ace**2)/len(resid_n10ace))/(np.sum(nhui)))
flc_n10ace
pers_n10ace =-1/model_n10ace.coef_
pers_n10ace
# Pers length error: error propagation from uncertainty in slope 
flc_n10ace/((model_n10ace.coef_)**2)
def slope_plen(x,m):
    return m*x

ttsam = np.arange(9, dtype=float)
ttsam *= lb_avg_n6ace
#nt_ttn10[0] = 0
ttsam
plt.figure(figsize=(7,7))
plt.errorbar(ace_ttn10, np.log(cor_n10ace_mat[0]), yerr=mk_n10ace, color='b',
             linestyle="None",marker='o',capsize=5, capthick=1, ecolor='black')
plt.plot(ttsam, slope_plen(ttsam, model_n10ace.coef_[0]), color='b')
plt.title(r'Ensemble Averaged ln(Cosine $\theta$) in pure acetone', fontsize=15, y=1.01)
plt.xlabel(r'Bond Length', fontsize=15)
plt.ylabel(r'ln$\left< Cos(\theta)\right >$', fontsize=15)
#plt.ylim(-1.9,0)
font = font_manager.FontProperties(family='Arial', style='normal', size='14')
plt.legend([r'$N_{PLGA}$ = 10: $L_{p}$ = 40.7 $\AA$ ± 2.6 $\AA$'], loc=3, frameon=0, fontsize=14, prop=font)
plt.tick_params(labelsize=14)
plt.text(5,0,r'R$^{2}$ = 0.995', fontsize=15, color='blue')
def exfun(x, a):
    return np.exp(-x*lb_avg_n6ace/a)

nt = np.arange(9)
plt.scatter(nt, cor_n10ace_mat[0,:])
plt.plot(nt,exfun(nt, pers_n10ace))<jupyter_output><empty_output><jupyter_text>N = 20 acetone<jupyter_code>n20_plga_ace = mda.Universe("n20_plga_150ns/1_xa_soln/1xa_n20sys.gro", "n20_plga_150ns/1_xa_soln/N20_nopbcpp_0xa.xtc")
n20_plga_ace
len(n20_plga_ace.trajectory)
poly_n20_ace = n20_plga_ace.select_atoms("resname sPLG PLG tPLG and not type H")
com_bond_n20ace = np.zeros(shape=(1,5000))
count = 0
for ts in n20_plga_ace.trajectory[0:5000]:
    n20_mon1 = n20_plga_ace.select_atoms("resid 1")
    n20_mon2 = n20_plga_ace.select_atoms("resid 2")
    oo_len = mda.analysis.distances.distance_array(n20_mon1.center_of_mass(), n20_mon2.center_of_mass(), 
                                                   box=n20_plga_ace.trajectory.ts.dimensions)
    com_bond_n20ace[0, count] = oo_len
    count += 1
lb_avg_n20ace = np.mean(com_bond_n20ace)
lb_avg_n20ace
plt.scatter(np.arange(5000),com_bond_n6ace)
plt.scatter(np.arange(5000),com_bond_n8ace)
plt.scatter(np.arange(5000),com_bond_n10ace)
plt.scatter(np.arange(5000),com_bond_n20ace)
n20_blks_ace = bavg_pers_cnt(5, poly_n20_ace, n20_plga_ace, lb_avg_n6ace, 0 , 15000)
n20_blks_ace["Avg persistence length"]
n20_blks_ace["Avg Radius of gyration"]
np.mean(n20_blks_ace["Avg persistence length"])
np.std(n20_blks_ace["Avg persistence length"])
np.mean(n20_blks_ace["Avg Radius of gyration"])
np.std(n20_blks_ace["Avg Radius of gyration"])
n20_rgens_ace, cor_n20ace_mat, N20aceens_cos, rg_avgn20ace = get_rg_pers_poly(poly_n20_ace, n20_plga_ace, 0, 15000, 
                                                                        lb_avg_n6ace)
cor_n20ace_mat
cor_n20ace_mat[0]
# x values
ace_ttn20 = np.arange(19, dtype=float)
ace_ttn20 *= lb_avg_n6ace
#nt_ttn20[0] = 0
ace_ttn20
# Sample standard deviation of Bond Length x values 
np.sqrt((np.sum(ace_ttn20**2) - ((np.sum(ace_ttn20)**2)/19))/18)
# Sample variance of Bond Length x values 
n20sace_var = (np.sum(ace_ttn20**2) - ((np.sum(ace_ttn20)**2)/19))/18
n20sace_var
# y values
n20lc_ace = np.log(cor_n20ace_mat[0])
n20lc_ace
mk_n20ace = cor_n20ace_mat[1]/cor_n20ace_mat[0]
plt.figure(figsize=(7,7))
plt.errorbar(ace_ttn20, np.log(cor_n20ace_mat[0]), yerr=mk_n20ace, color='b', linestyle="None",marker='o',
     capsize=5, capthick=1, ecolor='black')
ace_ttn20[:7]
from sklearn.linear_model import LinearRegression
model_n20ace = LinearRegression(fit_intercept=False)
model_n20ace.fit(ace_ttn20[:7].reshape(-1,1), n20lc_ace[:7])
# Slope here is in angstroms
print("Model slope:    ", model_n20ace.coef_)
print("Model intercept:", model_n20ace.intercept_)
gg_n20ace = model_n20ace.predict(ace_ttn20.reshape(-1,1))
gg_n20ace
print("Mean Std Error:", sklearn.metrics.mean_squared_error(n20lc_ace[:7], 
                                                            model_n20ace.predict(ace_ttn20[:7].reshape(-1,1))))
print("R2 score:", sklearn.metrics.r2_score(n20lc_ace[:7],
                                            model_n20ace.predict(ace_ttn20[:7].reshape(-1,1))))
# Residuals between the true y data and model y data 
resid_n20ace = n20lc_ace[:7] - gg_n20ace[:7]
resid_n20ace
# How to calculate mean squared error
np.sum(resid_n20ace**2)/len(resid_n20ace)
# Standard error: Square root of the mean squared error 
np.sqrt(np.sum(resid_n20ace**2)/len(resid_n20ace))
# sample variance of the residuals 
ntr_n20ace = resid_n20ace
ntr_n20ace -= np.mean(resid_n20ace)
nhr_n20ace = ntr_n20ace**2
np.sum(nhr_n20ace)/7
# How to calculate Sum((Xi - avg(X))^2): X values are the bond length values, I only used 7 points to fit data 
nt_n20ace = np.arange(7, dtype=float)
nt_n20ace *= lb_avg_n6ace
nt_n20ace -= np.mean(nt_n20ace)
nhui = nt_n20ace**2
np.sum(nhui)
# How to calculate 95% confidence interval for the slope 
flc_n20ace = scipy.stats.t.ppf(0.975, 4)*np.sqrt((np.sum(resid_n20ace**2)/len(resid_n20ace))/(np.sum(nhui)))
flc_n20ace
pers_n20ace =-1/model_n20ace.coef_
pers_n20ace
# Pers length error: error propagation from uncertainty in slope 
flc_n20ace/((model_n20ace.coef_)**2)
def slope_plen(x,m):
    return m*x

ttsam = np.arange(19, dtype=float)
ttsam *= lb_avg_n6ace
#nt_ttn10[0] = 0
ttsam
plt.figure(figsize=(7,7))
plt.errorbar(ace_ttn20, np.log(cor_n20ace_mat[0]), yerr=mk_n20ace, color='b',
             linestyle="None",marker='o',capsize=5, capthick=1, ecolor='black')
plt.plot(ttsam, slope_plen(ttsam, model_n20ace.coef_[0]), color='b')
plt.title(r'Ensemble Averaged ln(Cosine $\theta$) in pure acetone', fontsize=15, y=1.01)
plt.xlabel(r'Bond Length', fontsize=15)
plt.ylabel(r'ln$\left< Cos(\theta)\right >$', fontsize=15)
#plt.ylim(-1.9,0)
font = font_manager.FontProperties(family='Arial', style='normal', size='14')
plt.legend([r'$N_{PLGA}$ = 20: $L_{p}$ = 45.6 $\AA$ ± 3.6 $\AA$'], loc=3, frameon=0, fontsize=14, prop=font)
plt.tick_params(labelsize=14)
plt.text(5,0,r'R$^{2}$ = 0.994', fontsize=15, color='blue')
def exfun(x, a):
    return np.exp(-x*lb_avg_n6ace/a)

nt = np.arange(19)
plt.scatter(nt, cor_n20ace_mat[0,:])
plt.plot(nt,exfun(nt, pers_n20ace))<jupyter_output><empty_output><jupyter_text>N = 30 acetone<jupyter_code>n30_plga_ace = mda.Universe("n30_plga_150ns/1_xa_soln/1xa_n30sys.gro", "n30_plga_150ns/1_xa_soln/N30_nopbcpp_1xa.xtc")
n30_plga_ace
len(n30_plga_ace.trajectory)
poly_n30_ace = n30_plga_ace.select_atoms("resname sPLG PLG tPLG and not type H")
com_bond_n30ace = np.zeros(shape=(1,5000))
count = 0
for ts in n30_plga_ace.trajectory[0:5000]:
    n30_mon1 = n30_plga_ace.select_atoms("resid 1")
    n30_mon2 = n30_plga_ace.select_atoms("resid 2")
    oo_len = mda.analysis.distances.distance_array(n30_mon1.center_of_mass(), n30_mon2.center_of_mass(), 
                                                   box=n30_plga_ace.trajectory.ts.dimensions)
    com_bond_n30ace[0, count] = oo_len
    count += 1
lb_avg_n30ace = np.mean(com_bond_n30ace)
lb_avg_n30ace
plt.scatter(np.arange(5000),com_bond_n6ace)
plt.scatter(np.arange(5000),com_bond_n8ace)
plt.scatter(np.arange(5000),com_bond_n10ace)
plt.scatter(np.arange(5000),com_bond_n20ace)
plt.scatter(np.arange(5000),com_bond_n30ace)
n30_blks_ace = bavg_pers_cnt(5, poly_n30_ace, n30_plga_ace, lb_avg_n6ace, 0 , 15000)
n30_blks_ace["Avg persistence length"]
n30_blks_ace["Avg Radius of gyration"]
np.mean(n30_blks_ace["Avg persistence length"])
np.std(n30_blks_ace["Avg persistence length"])
np.mean(n30_blks_ace["Avg Radius of gyration"])
np.std(n30_blks_ace["Avg Radius of gyration"])
n30_rgens_ace, cor_n30ace_mat, N30aceens_cos, rg_avgn30ace = get_rg_pers_poly(poly_n30_ace, n30_plga_ace, 0, 15000, 
                                                                        lb_avg_n6ace)
cor_n30ace_mat
cor_n30ace_mat[0]
# x values
ace_ttn30 = np.arange(29, dtype=float)
ace_ttn30 *= lb_avg_n6ace
#nt_ttn20[0] = 0
ace_ttn30
# Sample standard deviation of Bond Length x values 
np.sqrt((np.sum(ace_ttn30**2) - ((np.sum(ace_ttn30)**2)/29))/28)
# Sample variance of Bond Length x values 
n30sace_var = (np.sum(ace_ttn30**2) - ((np.sum(ace_ttn30)**2)/29))/28
n30sace_var
# y values
n30lc_ace = np.log(cor_n30ace_mat[0])
n30lc_ace
mk_n30ace = cor_n30ace_mat[1]/cor_n30ace_mat[0]
plt.figure(figsize=(7,7))
plt.errorbar(ace_ttn30, np.log(cor_n30ace_mat[0]), yerr=mk_n30ace, color='b', linestyle="None",marker='o',
     capsize=5, capthick=1, ecolor='black')
ace_ttn30[:8]
from sklearn.linear_model import LinearRegression
model_n30ace = LinearRegression(fit_intercept=False)
model_n30ace.fit(ace_ttn30[:7].reshape(-1,1), n30lc_ace[:7])
# Slope here is in angstroms
print("Model slope:    ", model_n30ace.coef_)
print("Model intercept:", model_n30ace.intercept_)
gg_n30ace = model_n30ace.predict(ace_ttn30.reshape(-1,1))
gg_n30ace
print("Mean Std Error:", sklearn.metrics.mean_squared_error(n30lc_ace[:7],
                                                            model_n30ace.predict(ace_ttn30[:7].reshape(-1,1))))
print("R2 score:", sklearn.metrics.r2_score(n30lc_ace[:7], 
                                            model_n30ace.predict(ace_ttn30[:7].reshape(-1,1))))
# Residuals between the true y data and model y data 
resid_n30ace = n30lc_ace[:7] - gg_n30ace[:7]
resid_n30ace
# How to calculate mean squared error
np.sum(resid_n30ace**2)/len(resid_n30ace)
# Standard error: Square root of the mean squared error 
np.sqrt(np.sum(resid_n30ace**2)/len(resid_n30ace))
# sample variance of the residuals 
ntr_n30ace = resid_n30ace
ntr_n30ace -= np.mean(resid_n30ace)
nhr_n30ace = ntr_n30ace**2
# divide by n-1 samples used to fit 
np.sum(nhr_n30ace)/6
# How to calculate Sum((Xi - avg(X))^2): X values are the bond length values 
nt_n30ace = np.arange(7, dtype=float)
nt_n30ace *= lb_avg_n6ace
nt_n30ace -= np.mean(nt_n30ace)
nhui = nt_n30ace**2
np.sum(nhui)
# How to calculate 95% confidence interval for the slope 
flc_n30ace = scipy.stats.t.ppf(0.975, 4)*np.sqrt((np.sum(resid_n30ace**2)/len(resid_n30ace))/(np.sum(nhui)))
flc_n30ace
pers_n30ace =-1/model_n30ace.coef_
pers_n30ace
# Pers length error: error propagation from uncertainty in slope 
flc_n30ace/((model_n30ace.coef_)**2)
def slope_plen(x,m):
    return m*x

ttsam = np.arange(19, dtype=float)
ttsam *= lb_avg_n6ace
#nt_ttn10[0] = 0
ttsam
plt.figure(figsize=(7,7))
plt.errorbar(ace_ttn30, np.log(cor_n30ace_mat[0]), yerr=mk_n30ace, color='b',
             linestyle="None",marker='o',capsize=5, capthick=1, ecolor='black')
plt.plot(ttsam, slope_plen(ttsam, model_n30ace.coef_[0]), color='b')
plt.title(r'Ensemble Averaged ln(Cosine $\theta$) in pure acetone', fontsize=15, y=1.01)
plt.xlabel(r'Bond Length', fontsize=15)
plt.ylabel(r'ln$\left< Cos(\theta)\right >$', fontsize=15)
#plt.ylim(-1.9,0)
font = font_manager.FontProperties(family='Arial', style='normal', size='14')
plt.legend([r'$N_{PLGA}$ = 30: $L_{p}$ = 39.4 $\AA$ ± 4.3 $\AA$'], loc=3, frameon=0, fontsize=14, prop=font)
plt.tick_params(labelsize=14)
plt.text(5,0,r'R$^{2}$ = 0.99', fontsize=15, color='blue')
plt.figure(figsize=(7,7))
plt.errorbar(ace_ttn30, np.log(cor_n30ace_mat[0]), yerr=mk_n30ace, color='b', linestyle="None",marker='h',
     capsize=5, capthick=1, ecolor='black')
plt.plot(ttsam, slope_plen(ttsam, model_n30ace.coef_[0]), color='b')
plt.errorbar(ace_ttn20, np.log(cor_n20ace_mat[0]), yerr=mk_n20ace, color='g', linestyle="None"
             ,marker='*',capsize=5, capthick=1, ecolor='black')
plt.plot(ttsam, slope_plen(ttsam, model_n20ace.coef_[0]), color='g')
plt.errorbar(ace_ttn10, np.log(cor_n10ace_mat[0]), yerr=mk_n10ace, color='r', linestyle="None"
             ,marker='v', capsize=5, capthick=1, ecolor='black')
plt.plot(ttsam, slope_plen(ttsam, model_n10ace.coef_[0]), color='r')
plt.errorbar(ace_ttn8, np.log(cor_n8ace_mat[0]), yerr=mk_n8ace, color='k', linestyle="None",marker='s',
     capsize=5, capthick=1, ecolor='black')
plt.plot(ttsam, slope_plen(ttsam, model_n8ace.coef_[0]), color='k')
plt.errorbar(ace_ttn6, np.log(cor_n6ace_mat[0]), yerr=mk_n6ace, color='c', linestyle="None",marker='x',
     capsize=5, capthick=1, ecolor='black')
plt.plot(ttsam, slope_plen(ttsam, model_n6ace.coef_[0]), color='m')
plt.title(r'Ensemble Averaged ln(Cosine $\theta$): Pure acetone', fontsize=15, y=1.01)
plt.xlabel(r'Bond Length', fontsize=15)
plt.ylabel(r'ln$\left< Cos(\theta)\right >$', fontsize=15)
plt.ylim(-5, 1)
plt.xlim(0,70)
font = font_manager.FontProperties(family='Arial', style='normal', size='14')
#plt.legend([r'$N_{PLGA}$ = 30: $L_{p}$ = 18.5 $\AA$ ± 0.15 $\AA$,  R$^{2}$ = 0.99', r'$N_{PLGA}$ = 20: $L_{p}$ = 19.1 $\AA$ ± 0.27 $\AA$,  R$^{2}$ = 0.99'], 
#           loc=3, frameon=0, fontsize=14, prop=font)
#plt.legend([r'$N_{PLGA}$ = 20: $L_{p}$ = 19.1 $\AA$ ± 0.27 $\AA$,  R$^{2}$ = 0.99'], loc=3, frameon=0, fontsize=14, prop=font)
plt.tick_params(labelsize=14)
plt.text(1, -4.9,r'$N_{PLGA}$ = 30: $L_{p}$ = 39.4 $\AA$ ± 4.3 $\AA$, R$^{2}$ = 0.99', fontsize=14, color='blue')
plt.text(1,-4.65,r'$N_{PLGA}$ = 20: $L_{p}$ = 45.6 $\AA$ ± 3.6 $\AA$, R$^{2}$ = 0.99', fontsize=14, color='green')
plt.text(1,-4.38,r'$N_{PLGA}$ = 10: $L_{p}$ = 40.7 $\AA$ ± 2.6 $\AA$, R$^{2}$ = 0.99', fontsize=14, color='red')
plt.text(1,-4.12,r'$N_{PLGA}$ = 8: $L_{p}$ = 39.3 $\AA$ ± 5.9 $\AA$, R$^{2}$ = 0.98', fontsize=14, color='black')
plt.text(1,-3.85,r'$N_{PLGA}$ = 6: $L_{p}$ = 41.6 $\AA$ ± 5.7 $\AA$, R$^{2}$ = 0.99', fontsize=14, color='magenta')
#plt.text(2.5,0,r'R$^{2}$ = 0.96', fontsize=15, color='red')
def exfun(x, a):
    return np.exp(-x*lb_avg_n6ace/a)

nt = np.arange(29)
plt.scatter(nt, cor_n30ace_mat[0,:])
plt.plot(nt,exfun(nt, pers_n30ace))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1553832180535178, 'fraction_numerical': 0.0679277579822539, 'mean_word_length': 4.765659617321249, 'pattern_counts': {'":': 0, '<': 57, '<?xml version=': 0, '>': 57, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 43, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '38026', 'license': 'MIT', 'n_tokens_mistral': 26992, 'n_tokens_neox': 24538, 'n_words': 4324, 'path': 'Prot_poly_analysis/PLGA_analysis_MDAnalysis_07_31_2019.ipynb', 'repo_name': 'UWPRG/Nance_Enzyme_Encap_MD', 'size': '57252'}