<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Load the twobody data
Step1:  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$)
Step2:  Compute!
Step3:  Plot!
Step4:  Save the everything for later
First generate a beautiful graph...
<ASSISTANT_TASK:>
Python Code:
xyz = pd.read_hdf('xyz.hdf5', 'xyz')
twobody = pd.read_hdf('twobody.hdf5', 'twobody')
Explanation: Load the twobody data
End of explanation
from scipy.integrate import cumtrapz
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})
%load -s pcf, snippets/pcf.py
Explanation: 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}
End of explanation
A = 'O'
B = 'O'
df = pcf(A, B, a, twobody)
Explanation: Compute!
End of explanation
import seaborn as sns
sns.set_context('poster', font_scale=1.3)
sns.set_style('white')
sns.set_palette('colorblind')
# Lets modify a copy of the data for plotting
plotdf = df.set_index('r')
plotdf.columns = ['PCF', 'Pair Count']
# 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')
Explanation: Plot!
End of explanation
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')
# Save the figure
fig = ax.get_figure()
fig.savefig('pcf.pdf')
# Save the pcf data
store = pd.HDFStore('pcf.hdf5', mode='w')
store.put('pcf', df)
store.close()
Explanation: Save the everything for later
First generate a beautiful graph...
End of explanation
<END_TASK>