

In [ ]:

    
# One last time
import numpy as np
import scipy as sp
import pandas as pd

a = 12.55               # Cell size

Load the twobody data



In [ ]:

    
xyz = pd.read_hdf('xyz.hdf5', 'xyz')
twobody = pd.read_hdf('twobody.hdf5', 'twobody')

The Pair Correlation Function (or Radial Distribution Function)

Nice picture/description here

Basically we want to compute the following for a given pair of atom symbols ($A, B$):

\begin{equation} g_{AB}\left(r\right) = \frac{V}{4\pi r^{2}\Delta r MN_{A}N_{B}}\sum_{m=1}^{M}\sum_{a=1}^{N_{A}}\sum_{b=1}^{N_{B}}Q_{m}\left(r_{a}, r_{b}; r, \Delta r\right) \end{equation}\begin{equation} Q_{m}\left(r_{a}, r_{b}; r, \Delta r\right) = \begin{cases} 1\ \ if\ r - \frac{\Delta r}{2} \le \left|r_{a} - r_{b}\right|\lt r + \frac{\Delta r}{2}\\ 0\ \ otherwise \end{cases} \end{equation}

Note that that is the analytical form of the equation (meaning continuous values for r). As a consequence the denominator is simplified using an approximation for a volume of a spherical shell when $\Delta r$ is small.

Note:

\begin{equation} \frac{4}{3}\pi\left(r_{i+1}^{3} - r_{i}^{3}\right) \approx 4\pi r_{i}^{2}\Delta r \end{equation}

Computationally things will be a bit simpler...the summations are simply a histogram and there is no need to make the approximation above.

Algorithm:

Select the distances of interest
Compute the distance histogram
Multiply by the normalization constant

\begin{equation} \frac{V}{4\pi r^{2}\Delta r MN_{A}N_{B}} \equiv \frac{volume}{\left(distance\ count\right)\left(4 / 3 \pi\right)\left(r_{i+1}^{3} - r_{i}^{3}\right)} \end{equation}

Let's also compute the normalized integration of $g_{AB}(r)$ which returns the pairwise count with respect to distance:

\begin{equation} n(r) = \rho 4\pi r^2 g_{AB}(r) \end{equation}



In [ ]:

    
from scipy.integrate import cumtrapz



In [ ]:

    
def pcf(A, B, a, twobody, dr=0.05, start=0.5, end=7.5):
    '''
    Pair correlation function between two atom types.
    '''
    distances = twobody.loc[(twobody['symbols'] == A + B) |
                            (twobody['symbols'] == B + A), 'distance'].values
    bins = np.arange(start, end, dr)
    bins = np.append(bins, bins[-1] + dr)
    hist, bins = np.histogram(distances, bins)
    #...
    #...
    # r = ?
    # g = ?
    # n = ?
    return pd.DataFrame.from_dict({'r': None, 'g': None, 'n': None})



In [ ]:

    
%load -s pcf, snippets/pcf.py

Compute!



In [ ]:

    
A = 'O'
B = 'O'
df = pcf(A, B, a, twobody)

Plot!



In [ ]:

    
import seaborn as sns
sns.set_context('poster', font_scale=1.3)
sns.set_style('white')
sns.set_palette('colorblind')



In [ ]:

    
# Lets modify a copy of the data for plotting
plotdf = df.set_index('r')
plotdf.columns = ['PCF', 'Pair Count']



In [ ]:

    
# Generate the plot
ax = plotdf.plot(secondary_y='Pair Count')
ax.set_ylabel('Pair Correlation Function ({0}, {1})'.format(A, B))
ax.right_ax.set_ylabel('Pair Count ({0}, {1})'.format(A, B))
ax.set_xlabel('Distance ($\AA$)')
patches, labels = ax.get_legend_handles_labels()
patches2, labels2 = ax.right_ax.get_legend_handles_labels()
legend = ax.legend(patches+patches2, labels+labels2, loc='upper center', frameon=True)
frame = legend.get_frame()
frame.set_facecolor('white')
frame.set_edgecolor('black')

Save the everything for later

First generate a beautiful graph...



In [ ]:

    
df1 = pcf('O', 'O', a, twobody)
df2 = pcf('O', 'H', a, twobody)
df3 = pcf('H', 'H', a, twobody)
df = pd.concat((df1, df2, df3), axis=1)
df.columns = ['$g_{OO}$', '$n_{OO}$', '$r$', '$g_{OH}$', '$n_{OH}$', 'del1', '$g_{HH}$', '$n_{HH}$', 'del2']
del df['del1']
del df['del2']
df.set_index('$r$', inplace=True)
ax = df.plot(secondary_y=['$n_{OO}$', '$n_{OH}$', '$n_{HH}$'])
ax.set_ylabel('Pair Correlation Function ($g_{AB}$)')
ax.right_ax.set_ylabel('Pairwise Count ($n_{AB}$)')
ax.set_xlabel('Distance ($\AA$)')
ax.set_ylim(0, 5)
ax.right_ax.set_ylim(0, 20)
patches, labels = ax.get_legend_handles_labels()
patches2, labels2 = ax.right_ax.get_legend_handles_labels()
legend = ax.legend(patches+patches2, labels+labels2, loc='upper right', frameon=True)
frame = legend.get_frame()
frame.set_facecolor('white')
frame.set_edgecolor('black')



In [ ]:

    
# Save the figure
fig = ax.get_figure()
fig.savefig('pcf.pdf')



In [ ]:

    
# Save the pcf data
store = pd.HDFStore('pcf.hdf5', mode='w')
store.put('pcf', df)
store.close()