In [13]:
#this imports all libraries to be used in the code
import numpy as np
import pandas as pd
import scipy.optimize
import scipy.stats as st
from bokeh.models import NumeralTickFormatter
from bokeh.models import FixedTicker
import bokeh.io
import bokeh.plotting
import bokeh
from bokeh.plotting import figure, output_file, show
from bokeh.palettes import d3
from bokeh.palettes import brewer
from sympy.solvers import solve
from scipy.optimize import fsolve
from bokeh.models import LinearAxis, LogAxis, Range1d
import math
from sympy import Symbol, Poly, Eq, Function, exp
import math
import holoviews as hv
import numba
bokeh.io.output_notebook()
hv.extension('bokeh')
In [14]:
##let's define things again here to keep it contained.
##set the concentration range (in mg/ml)
cp = np.arange(0.0509,211.948,0.05) ##doing the concentration range we
##tested experimentally
##################################
#first solve for the depletion layer range (the characteristic polymer size)
#which is equation 3 in the main text
####################################
#hydrodynamic radii of the polymers in nanometers
rh_1M = 41.62 ##1 MDa
rh_100k = 11.1 ##100 kDa
rh_3k = 1.95 ## 3 kDa
##calculate the radius of gyration in
##nanometers for each using Kirkwood-Riseman relation
rg1M = rh_1M/0.665 #1 MDa
rg100k = rh_100k/0.665 #100 kDa
rg3k = rh_3k/0.665 #3 kDa
##below are the overlap concentrations for each in mg/ml
cstar1M = 1.62 #1 MDa
cstar100k = 8.6 #100 kDa
cstar3k = 53 #3 kDa
##now let's get the sizes from the crossover equation
##for the characteristic polymer size
##which is eq 3 in the main text
delta1M = rg1M*((1+2.3*(cp/cstar1M)**1.3))**-0.5 #1 MDa
delta100k = rg100k*((1+2.3*(cp/cstar100k)**1.3))**-0.5 #100 kDa
delta3k = rg3k*((1+2.3*(cp/cstar3k)**1.3))**-0.5 #3 kDa
###################################
####now let's solve for the osmotic pressure####
###################################
k = 1.38064852e-23 ##boltzmann constant (units = m^2*kg*s^-2*K^-1)
T = 298 ##temperature in Kelvins
Navo = 6.022e23 #avogadro's number (mol^-1)
##mws for each polymer in g/mol
mw1M = 1e6 #1 MDa
mw100k = 1e5 #100 kDa
mw3k = 3350 # 3 kDa
##solve for osmotic pressures for each (equation 2 of main text)
osm1M = ((Navo*k*T)/mw1M)*cp*(1+(cp/cstar1M)**1.3) #1 MDa
osm100k = ((Navo*k*T)/mw100k)*cp*(1+(cp/cstar100k)**1.3) #100 kDa
osm3k = ((Navo*k*T)/mw3k)*cp*(1+(cp/cstar3k)**1.3) #3 kDa
##this comes out in units of J/L. To have it in Pascals
##we need it in J/m^3 (Pa = J/m^3).
ltomcubed = 1000 ##liters to m^3
osm1M = osm1M*ltomcubed #osmotic pressure for 1 MDa
osm100k = osm100k*ltomcubed #osmotic pressure for 100 kDa
osm3k = osm3k*ltomcubed #osmotic pressure for 3 kDa
##########################################
#####Now solve for the interparticle potential###
###equation 5 of main text ################
##########################################
##set the range of distance from surface values
#(D) to solve for
D = np.arange(0.01,2*delta1M[0],0.1)
##vectors to store the interparticle potentials for each
#1MDa
fullpotential1M = np.zeros(len(D))
Vdep1M = np.zeros(len(D))
Vmix1M = np.zeros(len(D))
#100kDa
fullpotential100k = np.zeros(len(D))
Vdep100k = np.zeros(len(D))
Vmix100k = np.zeros(len(D))
#3kDa
fullpotential3k = np.zeros(len(D))
Vdep3k = np.zeros(len(D))
Vmix3k = np.zeros(len(D))
##a few variables we need
a = 500 #particle radius in nm
mcubetonmcube = 1e-27 ## m^3 to nm^3 conversion for depletion potential
L = 6.4 ##length of grafted layer on paticles in nm
gamma = 0.48 ##grafting density of polymer on particle surface in nm^-2
xi = 0.45 ##FH interaction parameter for PEG and water (from literature)
v1 = 0.2 ##volume of water molecule in nm^3
avgtheta2a = 0.43 #average volume fraction of grafted layer
##setup a for loop to solve for each MW and conc
##j will index the concentration
Vmin1M = np.zeros(len(cp)) ##storage vector for 1M minima
Vmin100k = np.zeros(len(cp)) ##storage vector for 100k minima
Vmin3k = np.zeros(len(cp)) ##storage vector for 3k minima
##########
for j in np.arange(len(cp)):
###1MDa loop####
##i will index the distance values
for i in np.arange(len(D)):
##first solve for depletion potential for 1MDa (eq 1, main text)
if D[i] < 2*delta1M[j]:
##writing down the depletion potential for 1MDa
temp = -2*np.pi*a*osm1M[j]*(delta1M[j]-D[i]/2)**2
##unit conversion and dividing my kT
Vdep1M[i] = (temp*mcubetonmcube)/(k*T)
else:
Vdep1M[i] = 0
##define the steric potential (eq 4, main text)
if D[i] < 2*L:
temp2 = ((4*np.pi*a*k*T)/v1)*(avgtheta2a**2)*(1/2-xi)*(L-D[i]/2)**2
Vmix1M[i] = temp2/(k*T)
else:
Vmix1M[i] = 0
##add them up
fullpotential1M[i] = Vdep1M[i]+Vmix1M[i]
##find the minimum in the interparticle potential
Vmin1M[j] = np.amin(fullpotential1M)
##############################
######now solve 100k########
##############################
##i will index the distance values
for i in np.arange(len(D)):
##first solve for 100k depletion potential (eq 1, main text)
if D[i] < 2*delta100k[j]:
##writing down the depletion potential for 100kDa
temp = -2*np.pi*a*osm100k[j]*(delta100k[j]-D[i]/2)**2
##unit conversion and dividing my kT
Vdep100k[i] = (temp*mcubetonmcube)/(k*T)
else:
Vdep100k[i] = 0
##define the steric potential (eq 4, main text)
if D[i] < 2*L:
temp2 = ((4*np.pi*a*k*T)/v1)*(avgtheta2a**2)*(1/2-xi)*(L-D[i]/2)**2
Vmix100k[i] = temp2/(k*T)
else:
Vmix100k[i] = 0
##add them up
fullpotential100k[i] = Vdep100k[i]+Vmix100k[i]
##find the minimum in the interparticle potential
Vmin100k[j] = np.amin(fullpotential100k)
##############################
######now solve 3k########
##############################
##i will index the distance values
for i in np.arange(len(D)):
##first solve for 3k depletion potential (eq 1, main text)
if D[i] < 2*delta3k[j]:
##writing down the depletion potential
temp = -2*np.pi*a*osm3k[j]*(delta3k[j]-D[i]/2)**2
##unit conversion and dividing my kT
Vdep3k[i] = (temp*mcubetonmcube)/(k*T)
else:
Vdep3k[i] = 0
##define the steric potential (eq 4, main text)
if D[i] < 2*L:
temp2 = ((4*np.pi*a*k*T)/v1)*(avgtheta2a**2)*(1/2-xi)*(L-D[i]/2)**2
Vmix3k[i] = temp2/(k*T)
else:
Vmix3k[i] = 0
##add them up
fullpotential3k[i] = Vdep3k[i]+Vmix3k[i]
##find the minimum in the interparticle potential
Vmin3k[j] = np.amin(fullpotential3k)
In [15]:
##plot them
p = bokeh.plotting.figure(plot_height=400,
plot_width=500,
title = "Potential minima",
#x_range = [sigma,(sigma+2*rh_1M)],
#y_axis_type = 'log',
x_axis_label='Polymer concentration (mg/ml)',
x_range = (0.05,200),
x_axis_type = 'log',
#y_range = (-50,20),
y_axis_label='V/kT')
# output_backend = 'svg',
# background_fill_color = None,
# border_fill_color = None)
p.circle(cp,Vmin1M,
legend = 'Vmin 1MDa',
color = 'green',
alpha = 0.2)
p.circle(cp,Vmin100k,
legend = 'Vmin 100kDa',
color = 'red',
alpha = 0.2)
p.circle(cp,Vmin3k,
legend = 'Vmin 3kDa',
color = 'blue',
alpha = 0.2)
p.legend.location = 'bottom_left'
bokeh.io.show(p)
Now let's use the same colors we had been using in the figures looking at aggregation in these three solutions and plot the magnitude of the minima. So we can compare them
In [16]:
colorblindpalette = ['#0072B2',
'#E69F00',
'#F0E442',
'#009E73',
'#56B4E9',
'#D55E00',
'#CC79A7',
'#000000']
In [17]:
##plot them
p = bokeh.plotting.figure(plot_height=400,
plot_width=500,
#title = "Potential minima",
#x_range = [sigma,(sigma+2*rh_1M)],
#y_axis_type = 'log',
x_axis_label='Polymer concentration (mg/ml)',
x_range = (0.03,400),
x_axis_type = 'log',
#y_range = (-50,20),
y_axis_label='|V_min|/kT',
output_backend = 'svg',
background_fill_color = None,
border_fill_color = None)
p.line(cp,abs(Vmin1M),
legend = 'PEG 1MDa',
color = colorblindpalette[0],
alpha = 0.7,
line_width = 3)
#line_dash = 'dashed')
p.line(cp,abs(Vmin100k),
legend = 'PEG 100kDa',
color = colorblindpalette[1],
alpha = 0.7,
line_width = 3)
#line_dash = 'dashed')
p.line(cp,abs(Vmin3k),
legend = 'PEG 3kDa',
color = colorblindpalette[2],
alpha = 0.7,
line_width = 3)
#line_dash = 'dashed')
p.legend.location = 'top_left'
p.xaxis[0].formatter = NumeralTickFormatter(format="0.0")
p.extra_y_ranges = {"foo": Range1d(start=0.03, end=400)}
p.add_layout(LogAxis(y_range_name="foo"), 'above')
p.xaxis[0].formatter = NumeralTickFormatter(format="0.0")
bokeh.io.show(p)
In [18]:
##a function for making paper quality figures
def gussyUpFigure(p):
p.title.text_font = "Arial"
p.title.text_font_style = 'bold'
p.title.text_font_size = '10pt'
p.xaxis.axis_label_text_font_size = "10pt"
p.yaxis.axis_label_text_font_size = "10pt"
p.xaxis.axis_label_text_font = "Arial"
p.yaxis.axis_label_text_font = "Arial"
p.xaxis.major_label_text_font_size = '8pt'
p.yaxis.major_label_text_font_size = '8pt'
p.xaxis.major_label_text_font_style = 'bold'
p.yaxis.major_label_text_font_style = 'bold'
#p.yaxis.axis_label_standoff = 3
p.xaxis.axis_label_standoff = 1
p.legend.label_text_font_size = "8pt"
p.legend.label_text_font = "Arial"
p.xgrid.grid_line_color = None
p.ygrid.grid_line_color = None
p.xaxis.axis_label_text_font_style = 'bold'
p.yaxis.axis_label_text_font_style = 'bold'
p.xaxis.axis_label_text_align = 'left'
p.axis.major_tick_in = 5
p.axis.major_tick_out = 0
p.axis.minor_tick_in = 0
p.axis.minor_tick_out = 0
p.xaxis.major_tick_line_width = 3
p.xaxis.minor_tick_line_width = 3
p.yaxis.major_tick_line_width = 3
p.yaxis.minor_tick_line_width = 3
p.outline_line_width = 2
p.outline_line_color = 'black'
#400 pt converts to 6.11 inches
##w = 3.4892 in, h = 2.1979 in
## for illustrator, 400 pix goes to 6.11 in
p.plot_height = 250
p.plot_width = 400
p.legend.border_line_color = None
p.title.align = 'center'
#p.xaxis.major_label_orientation = pi/3
return p
plots = [gussyUpFigure(p)]
plots = [p]
#export_png(p,filename = 'plot.png')
#p.yaxis.ticker = FixedTicker(ticks = [1,5,10,15])
bokeh.io.show(bokeh.layouts.column(plots))