pyPRISM Basics

In this notebook, we’ll go though the various pieces of pyPRISM and how we go about setting up a “problem” in the pyPRISM environment. We start with a classic physical system: a hard-sphere fluid.

Concepts

• pyPRISM basics
• posing problems in pyPRISM

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:

import pyPRISM
import numpy as np
import matplotlib.pyplot as plt
import holoviews as hv

hv.extension('bokeh')

Create System

To start, we define a “System” object which will contain all of the information that will specify the structure and interactions (chemistry) of a given system. In this case, our goal is to set up the calculation for a simple, hard-sphere monomer solution. We define the system to have 1 site-type called “monomer” and we define the reduced temperature to be \(k_{B}T=1.0\).

sys = pyPRISM.System(['monomer'],kT=1.0)

Define Domain

Next, we need to define the numerical solution domain. In other words, we must define the Real- and Fourier- space grid that we will solve the equations on. We do this by specifying the grid spacing in Real-space (dr) and the number of points in the grid (length). Specifying these values fixes both the Real- and Fourier- space grid, and we could have alternatively specified the Fourier-space grid spacing (dk). Note that it is best for efficiency if the length is a power of 2.

sys.domain = pyPRISM.Domain(dr=0.005,length=32768)
print('r =',sys.domain.r)
print('k = ',sys.domain.k)

r = [5.00000e-03 1.00000e-02 1.50000e-02 ... 1.63830e+02 1.63835e+02
 1.63840e+02]
k =  [1.91747598e-02 3.83495197e-02 5.75242795e-02 ... 6.28280181e+02
 6.28299356e+02 6.28318531e+02]

Define Interactions

Now we’ll define the interactions of the monomer solution. This effectively represents the chemistry of the system by defining how the sites attract or repel one another. For this example, we’ll be examining a simple, hard-sphere fluid.

sys.diameter['monomer'] = 1.0
sys.potential['monomer','monomer'] = pyPRISM.potential.HardSphere()

Define Composition

The composition of a system in PRISM theory is defined via site number densities. Since this is a one component system, we only need to specify one density.

sys.density['monomer'] = 0.8
sys.density

<Density total:0.80>

Define Molecular Structure

In general, the molecular structure is specified via \(\hat{\omega}(k)\) in PRISM calculations. These functions are essentially form-factors which describe the intra-molecular correlations within a molecule. It is through these functions that the molecular structure of a system is described and arbitrarily complicated \(\hat{\omega}(k)\) are possible. In this case of a simple hard-sphere fluid, \(\hat{\omega}(k)=1\) for all \(k\) by definition.

sys.omega['monomer','monomer'] = pyPRISM.omega.SingleSite()

Define Closure Approximation

In order to solve the PRISM equations numerically, an extra closure-approximation must be provided for each pair of sites. A detailed explanation of various closures is slightly beyond the scope of this introductory tutorial. A few guidelines are provided below

• The PercusYevick (PY) closure is a good starting point for most systems with divergent (hard-core, Weeks-Chandler-Andersen, Lennard-Jones) interactions.
• The HypernettedChain (HNC) closure is more accurate than PY for pairs of sites that have very different diameters. This is very useful knowledge when studying composite systems.

sys.closure['monomer','monomer'] = pyPRISM.closure.PercusYevick()

Create PRISM Object and Solve

Finally, we take the System object and use it to spawn a PRISM object. We have created the “System” and “PRISM” classes because they have separate, specific roles. The System class defines the structure and interactions of a given system while the PRISM class contains the machinery for setting up, solving, and storing the results of the mathematical PRISM formalism. After solve() is called, the result of the calculation is stored in the PRISM object.

PRISM = sys.createPRISM()
PRISM.solve()

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

     fun: array([ 4.15803093e-09,  1.25704004e-08,  2.08222783e-08, ...,
        2.85459257e-11, -2.23242033e-11,  4.31762904e-12])
 message: 'A solution was found at the specified tolerance.'
     nit: 7
  status: 1
 success: True
       x: array([ 7.31923234e-02,  2.18136517e-01,  3.61390574e-01, ...,
       -2.84237873e-11,  2.22021204e-11, -4.19551144e-12])

Plotting the Results

To start, we’ll plot the pair correlation function of this system. We calculate this value by passing the solved PRISM object to a calculation function.

%opts Curve Scatter Area [width=500,height=400]
%opts Scatter (size=10,alpha=0.5)

x = sys.domain.r
y = pyPRISM.calculate.pair_correlation(PRISM)['monomer','monomer']
rdf1 = hv.Curve((x,y),extents=(0,0,6,4),label='BaseRDF').redim.label(x='r',y='g(r)')
rdf1

We could instead calculate the structure factor of this system

%opts Curve Scatter Area [width=500,height=400]
%opts Scatter (size=10,alpha=0.5)

x = sys.domain.k
y = pyPRISM.calculate.structure_factor(PRISM)['monomer','monomer']
hv.Curve((x,y),extents=(0,0,20,None)).redim.label(x='k',y='S(k)')

Or maybe the second virial coefficient

B2 = pyPRISM.calculate.second_virial(PRISM)['monomer','monomer']
print('The monomer-monomer second-virial coefficient for this system is',B2)

The monomer-monomer second-virial coefficient for this system is 0.604288942117634

Putting It All Together

The above steps were long and detailed, but it’s important to realize that we only executed a few lines of code, shown below. Try changing the density, closures, and interactions and see how these changes affect the pair-correlation function. A few suggestions are below:

• change the closure to
  • sys.closure[‘monomer’,’monomer’] = pyPRISM.closure.HypernettedChain()
  • sys.closure[‘monomer’,’monomer’] = pyPRISM.closure.MSA(apply_hard_core=True)
• change the site diameter
  • diameter = 0.75
• change the density
  • density = 0.6

%opts Curve Scatter Area [width=500,height=400]
%opts Scatter (size=10,alpha=0.5)

sys = pyPRISM.System(['monomer'],kT=1.0)
sys.domain = pyPRISM.Domain(dr=0.005,length=32768)
sys.diameter['monomer'] = 1.0
sys.potential['monomer','monomer'] = pyPRISM.potential.HardSphere()
sys.density['monomer'] = 0.9
sys.omega['monomer','monomer'] = pyPRISM.omega.SingleSite()
sys.closure['monomer','monomer'] = pyPRISM.closure.HyperNettedChain(apply_hard_core=True)
PRISM = sys.createPRISM()
PRISM.solve()

x = sys.domain.r
y = pyPRISM.calculate.pair_correlation(PRISM)['monomer','monomer']
rdf2 = hv.Curve((x,y),extents=(0,0,6,6),label='ModifiedRDF').redim.label(x='r',y='g(r)')

hv.Overlay([rdf1,rdf2])

0:  |F(x)| = 3.97246; step 0.23001; tol 0.561713
1:  |F(x)| = 3.06534; step 0.382596; tol 0.535897
2:  |F(x)| = 2.67269; step 0.145932; tol 0.684197
3:  |F(x)| = 2.41102; step 0.121401; tol 0.732398
4:  |F(x)| = 2.13474; step 0.21196; tol 0.705555
5:  |F(x)| = 1.7113; step 0.409184; tol 0.578369
6:  |F(x)| = 1.41496; step 0.239333; tol 0.615285
7:  |F(x)| = 1.15999; step 0.232213; tol 0.604872
8:  |F(x)| = 0.902274; step 0.332558; tol 0.544518
9:  |F(x)| = 0.646089; step 0.496285; tol 0.461477
10:  |F(x)| = 0.371822; step 1; tol 0.298076
11:  |F(x)| = 0.042139; step 1; tol 0.0115596
12:  |F(x)| = 0.000206088; step 1; tol 2.15269e-05
13:  |F(x)| = 2.27462e-08; step 1; tol 1.09636e-08

Getting Extra Help

Besides going to the documentation website (see Github page) or compiling the documentation yourself, Jupyter provides a few ways to access the documentation in the notebook. Try running the cell below.

pyPRISM.potential.WeeksChandlerAndersen?

Alternatively, if you are inside the parenthesis of a function call, typing <Shift-Tab> will bring up a small,floating documentation popup. Try this below:

pyPRISM.potential.WeeksChandlerAndersen(epsilon=1.0,sigma=1.0)

<Potential: WeeksChandlerAndersen>

Summary

In this notebook, we stepped through the process of constructing a PRISM calculation in pyPRISM for a standard system in liquid state theory: the hard sphere fluid. We went over each step of the pyPRISM calculation in detail and then calculated several thermodynamic and structural quantities from the converged equations. In the following notebooks, we will move on to covering several real word systems and attempt to reproduce published data from scientific literature.

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