pyPRISM Internals¶
Here, we highlight some implementation details of the pyPRISM. This notebook is only necessary for those wishing to extend pyPRISM or better understand how the code functions. As such, this notebook will assume a higher familiarity Python and code development than the others.
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: 
[1]:
import numpy as np
import pyPRISM
import holoviews as hv
hv.extension('bokeh')
MatrixArray Data Structure¶
One of the challenges of using PRISM theory is the nature of the data itself: for every space-point in the solution domain, the user must store a \(n \times n\) matrix of correlation values. The theory requires that all of these matrices be used in mathematical operations e.g., matrix multiplcation and matrix inversion. Furthermore, as these operations are being carried out on many (>1024) matrices, these operations need to be optimized. Complicating this further, in order to transfer
between Fourier- and Real-space, the correlation values across matrices must be extracted as correlation curves. The schematic above describes the primary data structure that is used to carry out these tasks in pyPRISM: the MatrixArray.
Below, we create a 3x3 matrix array with 1024 matrices. MatrixArrays are type-aware so we also provide the site types. For this example, we will consider a polymer blend in a solvent. Finally, MatrixArrays keep track of whether they represent Real- or Fourier- space data, so we specify the “space” of the array.
[2]:
MA1 = pyPRISM.MatrixArray(length=1024,rank=3,types=['polymer1','polymer2','solvent'],space=pyPRISM.Space.Real)
MA1
[2]:
<MatrixArray rank:3 length:1024>
We can view the underlying Numpy array that stores the data and see that it’s intialized to all zeros by default. Note that you should not access this array directly and instead should rely on the defined methods as shown below.
[3]:
MA1.data
[3]:
array([[[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]],
[[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]],
[[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]],
...,
[[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]],
[[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]],
[[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]]])
We can assign and extract curves from the MatrixArray using the square-bracket operators. We’ll fill this MatrixArray with intra-molecular correlation functions. To do this, we’ll need to employ a few other features from pyPRISM.
Note that matrix symmetry is automatically applied during assignment so assigning MA[‘A’,’B’] automatically assigns M[‘B’,’A’].
[4]:
domain = pyPRISM.Domain(length=1024,dr=0.25)
MA1['polymer1','polymer1'] = pyPRISM.omega.Gaussian(length=1000,sigma=1.0).calculate(domain.k)
MA1['polymer2','polymer2'] = pyPRISM.omega.Gaussian(length=10,sigma=1.0).calculate(domain.k)
MA1['solvent','solvent'] = pyPRISM.omega.SingleSite().calculate(domain.k)
print('Full MatrixArray')
print(MA1.data)
Full MatrixArray
[[[991.68567821 0. 0. ]
[ 0. 9.99917183 0. ]
[ 0. 0. 1. ]]
[[967.35716814 0. 0. ]
[ 0. 9.99668766 0. ]
[ 0. 0. 1. ]]
[[928.76831425 0. 0. ]
[ 0. 9.9925496 0. ]
[ 0. 0. 1. ]]
...
[[ 1. 0. 0. ]
[ 0. 1. 0. ]
[ 0. 0. 1. ]]
[[ 1. 0. 0. ]
[ 0. 1. 0. ]
[ 0. 0. 1. ]]
[[ 1. 0. 0. ]
[ 0. 1. 0. ]
[ 0. 0. 1. ]]]
pyPRISM also provides methods to loop over the curves so they can be operated on in sequence.
[5]:
for i,t,curve in MA1.itercurve():
print('{}-{} Curve'.format(t[0],t[1]))
print(curve,'\n')
polymer1-polymer1 Curve
[991.68567821 967.35716814 928.76831425 ... 1. 1.
1. ]
polymer1-polymer2 Curve
[0. 0. 0. ... 0. 0. 0.]
polymer1-solvent Curve
[0. 0. 0. ... 0. 0. 0.]
polymer2-polymer2 Curve
[9.99917183 9.99668766 9.9925496 ... 1. 1. 1. ]
polymer2-solvent Curve
[0. 0. 0. ... 0. 0. 0.]
solvent-solvent Curve
[1. 1. 1. ... 1. 1. 1.]
We’ll create another MatrixArray which will start out as an array of identity matrices before we modify the values.
[6]:
# This creates a matrix Array of identity matrices
MA2 = pyPRISM.IdentityMatrixArray(length=1024,rank=3,types=['polymer1','polymer2','solvent'],space=pyPRISM.Space.Real)
print('Identity MatrixArray')
print(MA2.data,'\n')
MA2 *= 2.0 #multply all values by 2
MA2['solvent','solvent'] = 1.2 # set the entire 'solvent-solvent' curve equal to 1.2
MA2 += 3.0 # add three to all values
print('Modified MatrixArray')
print(MA2.data,'\n')
Identity MatrixArray
[[[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
[[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
[[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
...
[[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
[[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
[[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]]
Modified MatrixArray
[[[5. 3. 3. ]
[3. 5. 3. ]
[3. 3. 4.2]]
[[5. 3. 3. ]
[3. 5. 3. ]
[3. 3. 4.2]]
[[5. 3. 3. ]
[3. 5. 3. ]
[3. 3. 4.2]]
...
[[5. 3. 3. ]
[3. 5. 3. ]
[3. 3. 4.2]]
[[5. 3. 3. ]
[3. 5. 3. ]
[3. 3. 4.2]]
[[5. 3. 3. ]
[3. 5. 3. ]
[3. 3. 4.2]]]
An important step in the PRISM solution process involves inverting the matrices at all space points. The MatrixArray has a reasonably fast vectorized version of this operation built in.
[7]:
result = MA1.invert()
print('Inverted MatrixArray')
print(result.data,'\n')
Inverted MatrixArray
[[[0.00100838 0. 0. ]
[0. 0.10000828 0. ]
[0. 0. 1. ]]
[[0.00103374 0. 0. ]
[0. 0.10003313 0. ]
[0. 0. 1. ]]
[[0.00107669 0. 0. ]
[0. 0.10007456 0. ]
[0. 0. 1. ]]
...
[[1. 0. 0. ]
[0. 1. 0. ]
[0. 0. 1. ]]
[[1. 0. 0. ]
[0. 1. 0. ]
[0. 0. 1. ]]
[[1. 0. 0. ]
[0. 1. 0. ]
[0. 0. 1. ]]]
Finally, the MatrixArray provides a number of different operations between multiple arrays. Note that * represents in-place multiplication while the .dot() method is matrix mulitplication.
[8]:
result = MA1 * MA2
print('In-place Multiplication')
print(result.data,'\n')
result = MA1.dot(MA2)
print('Matrix Multiplication')
print(result.data,'\n')
In-place Multiplication
[[[4.95842839e+03 0.00000000e+00 0.00000000e+00]
[0.00000000e+00 4.99958591e+01 0.00000000e+00]
[0.00000000e+00 0.00000000e+00 4.20000000e+00]]
[[4.83678584e+03 0.00000000e+00 0.00000000e+00]
[0.00000000e+00 4.99834383e+01 0.00000000e+00]
[0.00000000e+00 0.00000000e+00 4.20000000e+00]]
[[4.64384157e+03 0.00000000e+00 0.00000000e+00]
[0.00000000e+00 4.99627480e+01 0.00000000e+00]
[0.00000000e+00 0.00000000e+00 4.20000000e+00]]
...
[[5.00000000e+00 0.00000000e+00 0.00000000e+00]
[0.00000000e+00 5.00000000e+00 0.00000000e+00]
[0.00000000e+00 0.00000000e+00 4.20000000e+00]]
[[5.00000000e+00 0.00000000e+00 0.00000000e+00]
[0.00000000e+00 5.00000000e+00 0.00000000e+00]
[0.00000000e+00 0.00000000e+00 4.20000000e+00]]
[[5.00000000e+00 0.00000000e+00 0.00000000e+00]
[0.00000000e+00 5.00000000e+00 0.00000000e+00]
[0.00000000e+00 0.00000000e+00 4.20000000e+00]]]
Matrix Multiplication
[[[4.95842839e+03 2.97505703e+03 2.97505703e+03]
[2.99975155e+01 4.99958591e+01 2.99975155e+01]
[3.00000000e+00 3.00000000e+00 4.20000000e+00]]
[[4.83678584e+03 2.90207150e+03 2.90207150e+03]
[2.99900630e+01 4.99834383e+01 2.99900630e+01]
[3.00000000e+00 3.00000000e+00 4.20000000e+00]]
[[4.64384157e+03 2.78630494e+03 2.78630494e+03]
[2.99776488e+01 4.99627480e+01 2.99776488e+01]
[3.00000000e+00 3.00000000e+00 4.20000000e+00]]
...
[[5.00000000e+00 3.00000000e+00 3.00000000e+00]
[3.00000000e+00 5.00000000e+00 3.00000000e+00]
[3.00000000e+00 3.00000000e+00 4.20000000e+00]]
[[5.00000000e+00 3.00000000e+00 3.00000000e+00]
[3.00000000e+00 5.00000000e+00 3.00000000e+00]
[3.00000000e+00 3.00000000e+00 4.20000000e+00]]
[[5.00000000e+00 3.00000000e+00 3.00000000e+00]
[3.00000000e+00 5.00000000e+00 3.00000000e+00]
[3.00000000e+00 3.00000000e+00 4.20000000e+00]]]
The MatrixArray interacts with domain objects to tranform from Real- to Fourier- space and vice versa. Note that these operations occur in-place and this method does not return a new MatrixArray.
[9]:
print('Real-Space')
print(MA1.data,'\n')
domain.MatrixArray_to_fourier(MA1)
print('Fourier-Space')
print(MA1.data,'\n')
Real-Space
[[[991.68567821 0. 0. ]
[ 0. 9.99917183 0. ]
[ 0. 0. 1. ]]
[[967.35716814 0. 0. ]
[ 0. 9.99668766 0. ]
[ 0. 0. 1. ]]
[[928.76831425 0. 0. ]
[ 0. 9.9925496 0. ]
[ 0. 0. 1. ]]
...
[[ 1. 0. 0. ]
[ 0. 1. 0. ]
[ 0. 0. 1. ]]
[[ 1. 0. 0. ]
[ 0. 1. 0. ]
[ 0. 0. 1. ]]
[[ 1. 0. 0. ]
[ 0. 1. 0. ]
[ 0. 0. 1. ]]]
Fourier-Space
[[[ 2.46329939e+07 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 2.33022049e+07 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 2.13822840e+07]]
[[-2.51675413e+06 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 -3.77481824e+06 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 -5.34036209e+06]]
[[ 4.67049835e+06 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 3.51128971e+06 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 2.37581679e+06]]
...
[[-1.63564440e+01 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 -3.19218659e+01 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 -3.20627732e+01]]
[[ 4.77535180e+01 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 3.22033679e+01 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 3.20625988e+01]]
[[-1.63244229e+01 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 -3.18593685e+01 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 -3.20000000e+01]]]
MatrixArrays also have safety checks build in to avoid operating across arrays of different spaces. The cell below should throw an error as MA1 is now in Fourier-space and MA2 is still in Real-space.
[10]:
MA1.dot(MA2)
---------------------------------------------------------------------------
AssertionError Traceback (most recent call last)
<ipython-input-10-dcdd31d919dc> in <module>()
----> 1 MA1.dot(MA2)
~/software/pyPRISM/dev/src/pyPRISM/core/MatrixArray.py in dot(self, other, inplace)
315 '''
316 if isinstance(other,MatrixArray):
--> 317 assert (self.space == other.space) or (Space.NonSpatial in (self.space,other.space)),MatrixArray.SpaceError
318 if inplace:
319 self.data = np.einsum('lij,ljk->lik', self.data, other.data)
AssertionError: Attempting MatrixArray math in non-matching spaces
ValueTables and PairTables¶
While the MatrixArray is the workhorse of the PRISM mathematics, there are many times while setting up a calculation that a more flexible data structure is needed. In particular, it is often the case that data needs to be stored without mathematical features or the data is non-numerical in nature. pyPRISM provides two data structures towards this end: ValueTable and PairTable
To start we look at the ValueTable. One of the features of PRISM data and PRISM parameters is that they are almost always associated with sites. The ValueTable provides an efficient way to access and manipulate data via site-types. See the examples below for an overview:
[11]:
VT = pyPRISM.ValueTable(types=['A','B','C','D','E'],name='density')
# set the value for type A to be 0.25
VT['A'] = 0.25
# set the value for types B & C to be 0.35
VT[ ['B','C'] ] = 0.35
# set all other values to be 0.1
VT.setUnset(0.1)
for i,t,v in VT:
print('{}) {} for type {} is {}'.format(i,VT.name,t,v))
0) density for type A is 0.25
1) density for type B is 0.35
2) density for type C is 0.35
3) density for type D is 0.1
4) density for type E is 0.1
Alternatively, there are many data in PRISM that are identified by a pair of site types. This is what the PairTable was designed to handle. Note that in this example, we are storing strings rather than numbers.
[12]:
PT = pyPRISM.PairTable(['A','B','C'],name='potential')
# Set the 'A-A' pair
PT['A','A'] = 'Lennard-Jones'
# Set the 'B-A', 'A-B', 'B-B', 'B-C', and 'C-B' pairs
PT['B',['A','B','C']] = 'Weeks-Chandler-Andersen'
# Set the 'C-A', 'A-C', 'C-C' pairs
PT['C',['A','C'] ] = 'Exponential'
#print only independent pairs
print('Independent Pairs')
for i,t,v in PT.iterpairs():
print('{}) {} for pair {}-{} is {}'.format(i,VT.name,t[0],t[1],v))
#print all pairs
print('\n','All Pairs')
for i,t,v in PT.iterpairs(full=True):
print('{}) {} for pair {}-{} is {}'.format(i,VT.name,t[0],t[1],v))
Independent Pairs
(0, 0)) density for pair A-A is Lennard-Jones
(0, 1)) density for pair A-B is Weeks-Chandler-Andersen
(0, 2)) density for pair A-C is Exponential
(1, 1)) density for pair B-B is Weeks-Chandler-Andersen
(1, 2)) density for pair B-C is Weeks-Chandler-Andersen
(2, 2)) density for pair C-C is Exponential
All Pairs
(0, 0)) density for pair A-A is Lennard-Jones
(0, 1)) density for pair A-B is Weeks-Chandler-Andersen
(0, 2)) density for pair A-C is Exponential
(1, 0)) density for pair B-A is Weeks-Chandler-Andersen
(1, 1)) density for pair B-B is Weeks-Chandler-Andersen
(1, 2)) density for pair B-C is Weeks-Chandler-Andersen
(2, 0)) density for pair C-A is Exponential
(2, 1)) density for pair C-B is Weeks-Chandler-Andersen
(2, 2)) density for pair C-C is Exponential
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
[ ]: