image0

Case Studies: Polymer Melts

In this example, we use pyPRISM to study and understand two simple but important systems: A dense, linear homopolymer melt of varying density and a ring homopolymer melt of varying molecular mass. In the linear melt system, we take advantage of the fact that PRISM theory doesn’t have a simulation ‘box’ and study very long (N=16,000) gaussian polymer chains. Note that studing chains of this length using Molecular Dynamics or Monte Carlo techniques would be difficult if not impossible. In the ring melt system, we leverage PRISM’s speed to quickly span a parameter space in molecular mass that would be challenging/time consuming to replicate synthetically.

Concepts

  • 1 component PRISM
  • Phase space “hopping”
  • Gaussian chains

Tools

  • pyPRISM.calculate.structure_factor
  • pyPRISM.calculate.pair_correlation

References

  1. Curro, J.G.; Schweizer K.S.; Integral Equation Theory of Homopolymer Melts, Molecular Crystals and Liquid Crystals Incorporating Nonlinear Optics, 1990, 180:1, 77-89, doi: 10.1080/00268949008025790
  2. Curro, J.G.; Schweizer K.S.; Theory of Polymer Melts: An Integral Equation Approach, Macromolecules, 1987, 20, 1928-1934, doi: 10.1021/ma00174a040

Notebook Setup

To begin, please run Kernel-> Restart & Clear Output from the menu at the top of the notebook. It is a good idea to run this before starting any notebook so that the notebook is fresh for the user. Next, run the cell below (via the top menu-bar or <Shift-Enter>. If the cell throws an import error, there is likely something wrong with your environment.

If successful, you should see a set of logos appear below the cell. Which logos appear depend on what is inside the hv.extension() command at the bottom of the cell. If no logos appear and the cell throws an error, there is likely something wrong with your environment.

Troubleshooting:

  • Did you activate the correct conda environment before starting the jupyter notebook?
  • If not using anaconda, did you install all dependencies before starting the jupyter notebook?
  • Is pyPRISM installed in your current environment on your PYTHONPATH?

Holoviews + Bokeh Logos: Logos

[1]:

import pyPRISM
import numpy as np
import holoviews as hv
hv.extension('bokeh')

Linear Homopolymer Melt: Comparison to Hard Sphere Fluid

First, we solve the hard-sphere fluid (HSF) system from the previous notebook. Note that while variables are “global” within a Jupyter notebook, separate notebooks do not share variables so we have to re-do this calculation here.

[2]:

sys = pyPRISM.System(['A'],kT=1.0)
sys.domain = pyPRISM.Domain(dr=0.005,length=32768)
sys.diameter['A']      = 1.0
sys.potential['A','A'] = pyPRISM.potential.HardSphere()
sys.density['A']       = 0.8
sys.omega['A','A']     = pyPRISM.omega.SingleSite()
sys.closure['A','A']   = pyPRISM.closure.PY(apply_hard_core=True)
PRISM_HSF = sys.createPRISM()
PRISM_HSF.solve()

x_HSF = sys.domain.r
y_HSF = pyPRISM.calculate.pair_correlation(PRISM_HSF)['A','A']

print('Done!')

0:  |F(x)| = 2.62378; step 1; tol 0.266006
1:  |F(x)| = 1.7493; step 1; tol 0.400055
2:  |F(x)| = 0.639398; step 1; tol 0.144039
3:  |F(x)| = 0.171479; step 1; tol 0.0647323
4:  |F(x)| = 0.00617309; step 1; tol 0.00116634
5:  |F(x)| = 3.63074e-06; step 1; tol 3.11335e-07
Done!

Then, we move onto the polmer melt case. Note that the only difference is the intra-molecular structure function (\(\hat{\omega}(k)\)).

[3]:

sys.omega['A','A'] = pyPRISM.omega.Gaussian(length=16000,sigma=sys.diameter.sigma['A','A'])
PRISM_Melt = sys.createPRISM()
PRISM_Melt.solve()

x_Melt = sys.domain.r
y_Melt = pyPRISM.calculate.pair_correlation(PRISM_Melt)['A','A']

print('Done!')

0:  |F(x)| = 49.5243; step 0.339595; tol 0.387728
1:  |F(x)| = 1.78198; step 1; tol 0.135299
2:  |F(x)| = 0.179788; step 1; tol 0.0091614
3:  |F(x)| = 0.00469113; step 1; tol 0.000612737
4:  |F(x)| = 1.06398e-06; step 1; tol 4.62974e-08
Done!

[4]:

%%opts Overlay [width=500,height=400,legend_position='top_right']

HSF = hv.Curve((x_HSF,y_HSF),label='Hard-Sphere Fluid',extents=(0,0,25,None))
Melt = hv.Curve((x_Melt,y_Melt),label='Gaussian Melt',extents=(0,0,25,None))

hv.Overlay([HSF,Melt])

[4]:

Now, let’s explore the effect of \(\hat{\omega}(k)\) on the pair correlation function. Insert the following \(\hat{\omega}(k)\) into the cell below.

Note that the length is purposefully changed for the last two :math:`omega`. - sys.omega[‘A’,’A’] = pyPRISM.omega.GaussianRing(length=16000,sigma=diameter) - sys.omega[‘A’,’A’] = pyPRISM.omega.FreelyJointedChain(length=16000,l=diameter) - sys.omega[‘A’,’A’] = pyPRISM.omega.NFJC(length=50,l=diameter) - sys.omega[‘A’,’A’] = pyPRISM.omega.DiscreteKoyama(sigma=1.0,l=diameter,length=50,lp=4/3)

[5]:

%%opts Overlay [width=500,height=400,legend_position='top_right']

# Change setting below
# ----------------------------------
sys.omega['A','A'] = pyPRISM.omega.GaussianRing(length=16000,sigma=sys.diameter['A','A'])
# ----------------------------------

PRISM_Melt = sys.createPRISM()
PRISM_Melt.solve()

x_Melt2 = sys.domain.r
y_Melt2 = pyPRISM.calculate.pair_correlation(PRISM_Melt)['A','A']

Melt2 = hv.Curve((x_Melt2,y_Melt2),label='Melt',extents=(0,0,25,None))
hv.Overlay([HSF,Melt,Melt2])

0:  |F(x)| = 74.9546; step 0.207941; tol 0.564017
1:  |F(x)| = 66.9276; step 0.108737; tol 0.717558
2:  |F(x)| = 66.6832; step 0.00446985; tol 0.893439
3:  |F(x)| = 66.5865; step 0.00180916; tol 0.89739
4:  |F(x)| = 65.8112; step 0.0146023; tol 0.879164
5:  |F(x)| = 60.6825; step 0.0998973; tol 0.765192
6:  |F(x)| = 60.5201; step 0.0361758; tol 0.895188
7:  |F(x)| = 60.3859; step 0.0353109; tol 0.896015
8:  |F(x)| = 60.2734; step 0.0340261; tol 0.896649
9:  |F(x)| = 60.2609; step 0.00431724; tol 0.899626
10:  |F(x)| = 60.2069; step 0.00769838; tol 0.898388
11:  |F(x)| = 60.1458; step 0.0222937; tol 0.898174
12:  |F(x)| = 60.1201; step 0.0100317; tol 0.899231
13:  |F(x)| = 60.0508; step 0.0274146; tol 0.897928
14:  |F(x)| = 59.9013; step 0.03615; tol 0.895523
15:  |F(x)| = 59.8455; step 0.0269553; tol 0.898326
16:  |F(x)| = 59.7906; step 0.0278476; tol 0.898347
17:  |F(x)| = 59.6937; step 0.0311638; tol 0.897087
18:  |F(x)| = 59.6934; step 0.000105054; tol 0.899989
19:  |F(x)| = 59.5362; step 0.0359007; tol 0.895266
20:  |F(x)| = 59.5225; step 0.00873502; tol 0.899586
21:  |F(x)| = 59.4992; step 0.0149687; tol 0.899296
22:  |F(x)| = 59.4764; step 0.0148359; tol 0.899311
23:  |F(x)| = 59.4029; step 0.0288751; tol 0.897776
24:  |F(x)| = 59.3802; step 0.0151269; tol 0.899313
25:  |F(x)| = 55.5415; step 1; tol 0.787398
26:  |F(x)| = 55.4249; step 1; tol 0.896224
27:  |F(x)| = 55.4169; step 1; tol 0.899741
28:  |F(x)| = 55.4167; step 1; tol 0.899995
29:  |F(x)| = 55.4167; step 1; tol 0.9
30:  |F(x)| = 55.4167; step 1; tol 0.9
31:  |F(x)| = 6.64034; step 1; tol 0.729
32:  |F(x)| = 6.11109; step 1; tol 0.762252
33:  |F(x)| = 3.97474; step 1; tol 0.522925
34:  |F(x)| = 1.54704; step 1; tol 0.246105
35:  |F(x)| = 0.0512046; step 1; tol 0.000985955
36:  |F(x)| = 0.000268186; step 1; tol 2.46887e-05
37:  |F(x)| = 1.11475e-08; step 1; tol 1.55498e-09

[5]:

Linear Homopolymer Melt: Comparison to Literature

Now that we’ve solved a simple, homopolymer system, let’s compare to data from the literature. This will also give us a chance to demonstrate some detailed scripting that we can do with typyPRISM.

First, we’ll load in the data extracted from the literature [Ref. 1] that we’ll be comparing against.

[6]:

gr_compare = []
sk_compare = []
for density in [0.6,0.8,1.0]:
    fname = '../data/LinearMelt-Gr-rho{}.csv'.format(density)
    x,y = np.loadtxt(fname,delimiter=',').T
    gr_compare.append([density,x,y])

    fname = '../data/LinearMelt-Sk-rho{}.csv'.format(density)
    x,y = np.loadtxt(fname,delimiter=',').T
    sk_compare.append([density,x,y])

Next, we define a bunch of variables to define how our plots will look later in this section.

[7]:

%opts Curve Scatter [width=500,height=400] Layout [shared_axes=False] Scatter (size=10,alpha=0.5)
%opts Curve Scatter [fontsize={'xlabel':14,'legend':14,'ylabel':14,'ticks':14}]
%opts Overlay [legend_position='bottom_left']
%opts Layout [shared_axes=False]


colors = {} # empty dictionary
colors[1.0] = 'blue'
colors[0.8] = 'red'
colors[0.6] = 'green'

ls = {} # empty dictionary
ls[1.0] = 'solid'
ls[0.8] = 'dashed'
ls[0.6] = 'dotted'

markers = {} # empty dictionary
markers[1.0] = 'o'
markers[0.8] = '^'
markers[0.6] = 'd'

Next, we’ll use pyPRISM to calculate data at three densities for a linear polymer melt of of gaussian chains of length 16000, as described in Reference 1. Notice how we store the resulting “guess” from the previous density and use it as the starting guess for the next density. After solving, note how many iterations each density took.

[8]:

sys = pyPRISM.System(['polymer'],kT=1.0)
sys.domain = pyPRISM.Domain(dr=0.005,length=32768)
sys.diameter['polymer'] = 1.0
sys.closure['polymer','polymer'] = pyPRISM.closure.PercusYevick()
sys.potential['polymer','polymer'] = pyPRISM.potential.HardSphere()
sys.omega['polymer','polymer'] = pyPRISM.omega.Gaussian(sigma=sys.diameter['polymer','polymer'],length=16000)

gr_results = []
sk_results = []
guess = np.zeros_like(sys.domain.r)
for density in [0.6,0.8,1.0]:
    print('==> Solving for density {}'.format(density))
    sys.density['polymer'] = density
    PRISM = sys.createPRISM()
    result = PRISM.solve(guess)
    guess = np.copy(PRISM.x)

    y = pyPRISM.calculate.structure_factor(PRISM)['polymer','polymer']
    x = sys.domain.k
    sk_results.append([density,x,y])

    x = sys.domain.r
    y = pyPRISM.calculate.pair_correlation(PRISM)['polymer','polymer']
    gr_results.append([density,x,y])
    print('')


print('Done!')

==> Solving for density 0.6
0:  |F(x)| = 81.9231; step 0.186349; tol 0.593894
1:  |F(x)| = 81.8849; step 0.00049023; tol 0.89916
2:  |F(x)| = 81.6474; step 0.003047; tol 0.894787
3:  |F(x)| = 81.6472; step 2.06434e-06; tol 0.899996
4:  |F(x)| = 81.6443; step 3.76541e-05; tol 0.899936
5:  |F(x)| = 81.5191; step 0.00161753; tol 0.897241
6:  |F(x)| = 52.1202; step 0.15956; tol 0.724537
7:  |F(x)| = 46.0427; step 1; tol 0.702348
8:  |F(x)| = 42.036; step 1; tol 0.750178
9:  |F(x)| = 41.0872; step 1; tol 0.859828
10:  |F(x)| = 40.8693; step 1; tol 0.89048
11:  |F(x)| = 40.8594; step 1; tol 0.899568
12:  |F(x)| = 40.8594; step 1; tol 0.899996
13:  |F(x)| = 40.8593; step 1; tol 0.899999
14:  |F(x)| = 40.8593; step 1; tol 0.9
15:  |F(x)| = 40.8593; step 1; tol 0.9
16:  |F(x)| = 48.9979; step 1; tol 0.9999
17:  |F(x)| = 40.8604; step 1; tol 0.89982
18:  |F(x)| = 40.8585; step 0.0185991; tol 0.899916
19:  |F(x)| = 48.2255; step 1; tol 0.9999
20:  |F(x)| = 39.8334; step 1; tol 0.89982
21:  |F(x)| = 39.4845; step 0.494518; tol 0.884299
22:  |F(x)| = 34.0081; step 1; tol 0.703786
23:  |F(x)| = 46.7361; step 1; tol 0.9999
24:  |F(x)| = 41.5151; step 0.431742; tol 0.89982
25:  |F(x)| = 46.705; step 1; tol 0.9999
26:  |F(x)| = 41.7892; step 0.407615; tol 0.89982
27:  |F(x)| = 37.32; step 1; tol 0.728708
28:  |F(x)| = 45.2683; step 1; tol 0.9999
29:  |F(x)| = 41.8015; step 1; tol 0.89982
30:  |F(x)| = 35.7199; step 1; tol 0.728708
31:  |F(x)| = 34.5162; step 0.419989; tol 0.840364
32:  |F(x)| = 47.4471; step 1; tol 0.9999
33:  |F(x)| = 36.6828; step 1; tol 0.89982
34:  |F(x)| = 32.6796; step 1; tol 0.728708
35:  |F(x)| = 32.3298; step 0.0507434; tol 0.880832
36:  |F(x)| = 38.8325; step 1; tol 0.9999
37:  |F(x)| = 32.9971; step 1; tol 0.89982
38:  |F(x)| = 32.4921; step 0.476009; tol 0.872663
39:  |F(x)| = 39.0104; step 1; tol 0.9999
40:  |F(x)| = 33.1882; step 1; tol 0.89982
41:  |F(x)| = 33.0107; step 0.450081; tol 0.890398
42:  |F(x)| = 36.9105; step 1; tol 0.9999
43:  |F(x)| = 32.3099; step 1; tol 0.89982
44:  |F(x)| = 30.6935; step 0.0565517; tol 0.812199
45:  |F(x)| = 27.2482; step 0.282094; tol 0.709296
46:  |F(x)| = 27.1626; step 0.00388862; tol 0.89435
47:  |F(x)| = 24.9831; step 1; tol 0.761364
48:  |F(x)| = 23.4707; step 0.28068; tol 0.794334
49:  |F(x)| = 20.5345; step 1; tol 0.688905
50:  |F(x)| = 18.9605; step 0.17719; tol 0.767312
51:  |F(x)| = 14.6473; step 1; tol 0.537103
52:  |F(x)| = 14.3981; step 0.0218796; tol 0.869634
53:  |F(x)| = 11.786; step 1; tol 0.680637
54:  |F(x)| = 10.4035; step 0.160641; tol 0.701247
55:  |F(x)| = 8.94065; step 0.300383; tol 0.664692
56:  |F(x)| = 7.92109; step 1; tol 0.706438
57:  |F(x)| = 4.49779; step 1; tol 0.449149
58:  |F(x)| = 2.77687; step 1; tol 0.343047
59:  |F(x)| = 0.622838; step 1; tol 0.105913
60:  |F(x)| = 0.0143531; step 1; tol 0.000477952
61:  |F(x)| = 2.56211e-06; step 1; tol 2.86778e-08

==> Solving for density 0.8
0:  |F(x)| = 0.706714; step 1; tol 0.000751233
1:  |F(x)| = 0.0284976; step 1; tol 0.00146343
2:  |F(x)| = 7.05967e-05; step 1; tol 5.52325e-06

==> Solving for density 1.0
0:  |F(x)| = 0.413714; step 1; tol 0.000691077
1:  |F(x)| = 0.018618; step 1; tol 0.00182267
2:  |F(x)| = 5.4341e-05; step 1; tol 7.6671e-06

Done!

We start by reproducing Figure 1 of Reference 1.

[9]:

%%opts Overlay [legend_position='bottom_right']
from math import sqrt

extents = (0,0,0.6,1.0)

gr_plots = []
for rho,x,y in gr_results:
    Rg = sqrt(16000/6.0)
    label = 'rho={} (pyPRISM)'.format(rho)
    style = {'line_dash':ls[rho],'color':colors[rho]}
    c1 = hv.Curve((x/Rg,y),label=label,extents=extents)(style=style)
    gr_plots.append(c1)

for rho,x,y in gr_compare:
    label = 'rho={} (Ref 1)'.format(rho)
    style = {'marker':markers[rho],'color':colors[rho]}
    c1 = hv.Scatter((x,y),label=label,extents=extents)(style=style)
    gr_plots.append(c1)


hv.Overlay(gr_plots).redim.label(x='r',y='g(r)')

[9]:

Next, we reproduce Figure 3:

[10]:

sk_plots = []
for rho,x,y in sk_results:
    yd = y - 1.0
    label = 'rho={} (pyPRISM)'.format(rho)
    style = {'line_dash':ls[rho],'color':colors[rho]}
    c1 = hv.Curve((x,yd/yd[0]),label=label,extents=(0,0,3,None))(style=style)
    sk_plots.append(c1)

for rho,x,y in sk_compare:
    style = {'marker':markers[rho],'color':colors[rho]}
    c1 = hv.Scatter((x,y),label='rho={} (Ref 1)'.format(rho),extents=(0.0,0,3,None))(style=style)
    sk_plots.append(c1)


hv.Overlay(sk_plots).redim.label(x='k',y='(S(k)-1)/(S(0)-1)')

[10]:

Finally, we reproduce the inset of Figure 3:

[11]:

sk0_plots = []
for rho,x,y in sk_results:
    yd = y - 1
    label = 'rho={} (pyPRISM)'.format(rho)
    style = {'line_dash':ls[rho],'color':colors[rho]}
    c1 = hv.Scatter((rho,yd[0]),label=label,extents=(0.4,0,1.2,11))(style=style)
    sk0_plots.append(c1)

hv.Overlay(sk0_plots).redim.label(x='rho',y='S(0)')

[11]:

Ring Homopolymer Melt: Comparison to Literature

As a comparison, we’ll now focus on a ring-homopolymer melt as described in [Ref. 2]. This paper is actually the first full paper in which PRISM theory was described.

Again, to start we define the aesthetics of the plots.

[12]:

%opts Curve Scatter [width=500,height=400] Layout [shared_axes=False] Scatter (size=10,alpha=0.5)
%opts Curve Scatter [fontsize={'xlabel':14,'xlabel':14,'ylabel':14,'ticks':12}]
%opts Overlay [legend_position='bottom_left']
%opts Layout [shared_axes=False]


colors = {}
colors[2000] = 'blue'
colors[500] = 'red'
colors[100] = 'green'

ls = {}
ls[2000] = 'solid'
ls[500] = 'dashed'
ls[100] = 'dotted'

markers = {}
markers[2000] = 'o'
markers[500] = '^'
markers[100] = 'd'

[13]:

gr_compare = []
for density in [2000,500,100]:
    fname = '../data/RingMelt-Gr-N{}.csv'.format(density)
    x,y = np.loadtxt(fname,delimiter=',').T
    gr_compare.append([density,x,y])

Now we do the PRISM calculation…

[14]:

sys = pyPRISM.System(['polymer'],kT=1.0)
sys.domain = pyPRISM.Domain(dr=0.05,length=4096)
sys.diameter['polymer'] = 1.0
sys.density['polymer'] = 0.9
sys.closure['polymer','polymer'] = pyPRISM.closure.PercusYevick()
sys.potential['polymer','polymer'] = pyPRISM.potential.HardSphere()

gr_results = []
guess = np.zeros_like(sys.domain.r)
for length in [2000,500,100]:
    print('==> Solving for gaussian ring of length {}'.format(length))
    sys.omega['polymer','polymer'] = pyPRISM.omega.GaussianRing(sigma=sys.diameter['polymer','polymer'],length=length)
    PRISM = sys.createPRISM()
    result = PRISM.solve(guess)
    guess = np.copy(PRISM.x)

    x = sys.domain.r
    y = pyPRISM.calculate.pair_correlation(PRISM)['polymer','polymer']
    gr_results.append(['ring',length,x,y])

print('Done!')

==> Solving for gaussian ring of length 2000
0:  |F(x)| = 13.564; step 0.116955; tol 0.701605
1:  |F(x)| = 12.8877; step 0.0499993; tol 0.812488
2:  |F(x)| = 12.6488; step 0.0187059; tol 0.866947
3:  |F(x)| = 11.7434; step 0.0724084; tol 0.775771
4:  |F(x)| = 11.7428; step 5.40487e-05; tol 0.899907
5:  |F(x)| = 11.7419; step 7.80313e-05; tol 0.899865
6:  |F(x)| = 11.7406; step 0.000120235; tol 0.899793
7:  |F(x)| = 11.7383; step 0.000203317; tol 0.899649
8:  |F(x)| = 11.7339; step 0.000393242; tol 0.899322
9:  |F(x)| = 11.7233; step 0.000937021; tol 0.898386
10:  |F(x)| = 11.6878; step 0.00316625; tol 0.894554
11:  |F(x)| = 11.4713; step 0.0193873; tol 0.866966
12:  |F(x)| = 11.4571; step 0.00132076; tol 0.897765
13:  |F(x)| = 11.4557; step 0.000126257; tol 0.899787
14:  |F(x)| = 11.453; step 0.0002504; tol 0.899577
15:  |F(x)| = 11.4463; step 0.000621288; tol 0.89895
16:  |F(x)| = 11.4221; step 0.00225248; tol 0.896198
17:  |F(x)| = 11.2524; step 0.0158588; tol 0.873447
18:  |F(x)| = 11.2523; step 6.45276e-06; tol 0.899989
19:  |F(x)| = 11.2522; step 8.00376e-06; tol 0.899987
20:  |F(x)| = 11.2521; step 1.01827e-05; tol 0.899983
21:  |F(x)| = 11.252; step 1.32855e-05; tol 0.899978
22:  |F(x)| = 11.2518; step 1.79953e-05; tol 0.89997
23:  |F(x)| = 11.2515; step 2.55292e-05; tol 0.899958
24:  |F(x)| = 11.2511; step 3.83886e-05; tol 0.899936
25:  |F(x)| = 11.2505; step 6.28229e-05; tol 0.899896
26:  |F(x)| = 11.2493; step 0.000115591; tol 0.899808
27:  |F(x)| = 11.2467; step 0.000254364; tol 0.899578
28:  |F(x)| = 11.2389; step 0.000750361; tol 0.898757
29:  |F(x)| = 11.2005; step 0.00371077; tol 0.893866
30:  |F(x)| = 10.7864; step 0.0403192; tol 0.834678
31:  |F(x)| = 10.4959; step 0.0314757; tol 0.85217
32:  |F(x)| = 10.4731; step 0.00275572; tol 0.896107
33:  |F(x)| = 10.1771; step 0.0362332; tol 0.849843
34:  |F(x)| = 10.1755; step 0.000245052; tol 0.899704
35:  |F(x)| = 10.1582; step 0.00253582; tol 0.896942
36:  |F(x)| = 9.98704; step 0.025387; tol 0.869934
37:  |F(x)| = 9.97755; step 0.0016413; tol 0.898291
38:  |F(x)| = 9.80201; step 0.0306349; tol 0.868609
39:  |F(x)| = 9.80165; step 7.641e-05; tol 0.899935
40:  |F(x)| = 9.8015; step 3.34378e-05; tol 0.899972
41:  |F(x)| = 9.80082; step 0.000147167; tol 0.899875
42:  |F(x)| = 9.79355; step 0.00157967; tol 0.898665
43:  |F(x)| = 9.69962; step 0.0206006; tol 0.882819
44:  |F(x)| = 9.67173; step 0.00710773; tol 0.894833
45:  |F(x)| = 9.66218; step 0.00256072; tol 0.898222
46:  |F(x)| = 9.64319; step 0.0051766; tol 0.896467
47:  |F(x)| = 9.58713; step 0.0158394; tol 0.889566
48:  |F(x)| = 9.54215; step 0.0142497; tol 0.891575
49:  |F(x)| = 9.49661; step 0.0159618; tol 0.89143
50:  |F(x)| = 9.48357; step 0.00511646; tol 0.897529
51:  |F(x)| = 9.48304; step 0.000214479; tol 0.8999
52:  |F(x)| = 9.47626; step 0.00275702; tol 0.898713
53:  |F(x)| = 9.44666; step 0.0122378; tol 0.894388
54:  |F(x)| = 9.41731; step 0.0131669; tol 0.894416
55:  |F(x)| = 9.38668; step 0.0149674; tol 0.894154
56:  |F(x)| = 9.37647; step 0.00548425; tol 0.898045
57:  |F(x)| = 9.35142; step 0.00696148; tol 0.895196
58:  |F(x)| = 9.34054; step 0.00456984; tol 0.897907
59:  |F(x)| = 9.33118; step 0.00387752; tol 0.898198
60:  |F(x)| = 9.33091; step 9.42877e-05; tol 0.899948
61:  |F(x)| = 9.30692; step 0.00906009; tol 0.895377
62:  |F(x)| = 9.26985; step 0.0154426; tol 0.892846
63:  |F(x)| = 9.26204; step 0.00362434; tol 0.898483
64:  |F(x)| = 9.23088; step 0.0160041; tol 0.893954
65:  |F(x)| = 9.22582; step 0.00330641; tol 0.899015
66:  |F(x)| = 9.2188; step 0.00389694; tol 0.898631
67:  |F(x)| = 9.20601; step 0.00567855; tol 0.897503
68:  |F(x)| = 9.18625; step 0.0109425; tol 0.896142
69:  |F(x)| = 9.14046; step 0.0194129; tol 0.891049
70:  |F(x)| = 9.11915; step 0.00957558; tol 0.895808
71:  |F(x)| = 9.1081; step 0.00603236; tol 0.89782
72:  |F(x)| = 9.08375; step 0.0136708; tol 0.895195
73:  |F(x)| = 9.07434; step 0.00476936; tol 0.898136
74:  |F(x)| = 9.07327; step 0.000515113; tol 0.899787
75:  |F(x)| = 9.06892; step 0.0020917; tol 0.899138
76:  |F(x)| = 9.05313; step 0.00790108; tol 0.896869
77:  |F(x)| = 9.0377; step 0.00804661; tol 0.896935
78:  |F(x)| = 9.00401; step 0.0164117; tol 0.893303
79:  |F(x)| = 8.97103; step 0.0167142; tol 0.893419
80:  |F(x)| = 8.97027; step 0.00040741; tol 0.899847
81:  |F(x)| = 8.95825; step 0.00655805; tol 0.897591
82:  |F(x)| = 8.94265; step 0.00873797; tol 0.896866
83:  |F(x)| = 8.93898; step 0.00177432; tol 0.899263
84:  |F(x)| = 8.93109; step 0.00431515; tol 0.898412
85:  |F(x)| = 8.93102; step 3.81577e-05; tol 0.899986
86:  |F(x)| = 8.92359; step 0.0043836; tol 0.898502
87:  |F(x)| = 8.89848; step 0.0127696; tol 0.894942
88:  |F(x)| = 8.89697; step 0.000852819; tol 0.899694
89:  |F(x)| = 8.87556; step 0.0122856; tol 0.895674
90:  |F(x)| = 8.85198; step 0.0132539; tol 0.895225
91:  |F(x)| = 8.84007; step 0.00638331; tol 0.89758
92:  |F(x)| = 8.82702; step 0.00722105; tol 0.897344
93:  |F(x)| = 8.82701; step 3.24097e-06; tol 0.899999
94:  |F(x)| = 7.27724; step 1; tol 0.728998
95:  |F(x)| = 1.87778; step 1; tol 0.478294
96:  |F(x)| = 1.39712; step 0.417561; tol 0.498219
97:  |F(x)| = 0.892148; step 1; tol 0.366988
98:  |F(x)| = 0.19859; step 1; tol 0.121212
99:  |F(x)| = 0.00803493; step 1; tol 0.0014733
100:  |F(x)| = 4.27395e-05; step 1; tol 2.54646e-05
101:  |F(x)| = 8.29099e-10; step 1; tol 3.38685e-10
==> Solving for gaussian ring of length 500
0:  |F(x)| = 0.0018377; step 1; tol 5.01925e-08
1:  |F(x)| = 1.16134e-06; step 1; tol 3.59432e-07
==> Solving for gaussian ring of length 100
0:  |F(x)| = 0.0146501; step 1; tol 5.13412e-06
1:  |F(x)| = 7.48487e-05; step 1; tol 2.34925e-05
2:  |F(x)| = 2.44558e-09; step 1; tol 9.60811e-10
Done!

Here we reproduce Figure 1 of Reference 1.

[15]:

%%opts Overlay [legend_position='bottom_right']
from math import sqrt

extents = (0,0,13,1.0)

gr_plots = []
for name,N,x,y in gr_results:
    label = 'N={} (PRISM)'.format(N)
    style = {'line_dash':ls[N],'color':colors[N]}
    c1 = hv.Curve((x,y),label=label,extents=extents)(style=style)
    gr_plots.append(c1)

for N,x,y in gr_compare:
    label = 'N={} (Ref 2)'.format(N)
    style = {'marker':markers[N],'color':colors[N]}
    c1 = hv.Scatter((x,y),label=label,extents=extents)(style=style)
    gr_plots.append(c1)


hv.Overlay(gr_plots).redim.label(x='r',y='g(r)')


[15]:

Summary

With both the previous and this notebook finished, we have covered a basic monomer fluid system and two “real-world” polymer systems. We have verified that the pyPRISM implementation matches previous implementations by reproducing data from the literature. Finally, we have begun to touch on two more advanced topics: solving difficult to converge systems and using loops to scan parameter spaces. We will go into more detail on both of these topics in the following notebooks.

image0

NB0.Introduction \(\cdot\) NB1.PythonBasics \(\cdot\) NB2.Theory.General \(\cdot\) NB3.Theory.PRISM \(\cdot\) NB4.pyPRISM.Overview \(\cdot\) NB5.CaseStudies.PolymerMelts \(\cdot\) NB6.CaseStudies.Nanocomposites \(\cdot\) NB7.CaseStudies.Copolymers \(\cdot\) NB8.pyPRISM.Internals \(\cdot\) NB9.pyPRISM.Advanced

[ ]: