# Imports of required packages
from rdftool import InterRDF
import MDAnalysis as mda
%matplotlib inline
import matplotlib.pyplot as plt

"""
Explanation: Making RDFs
This notebook demonstrates how to build an intermolecular RDF using tools built upon MDAnalysis.
Requirements:
 - MDAnalysis version 0.10.0
 - numpy
End of explanation
"""


u = mda.Universe('conf.gro')
print u
print len(u.trajectory)
print set(u.atoms.names())
print set(u.atoms.types())

"""
Explanation: Firstly, let's load up a Universe and inspect what we have in our system.  The test system here is a single frame of octanol.
End of explanation
"""


oxygen = u.selectAtoms('name HO')
hydrogen = u.selectAtoms('name HO')

print len(oxygen)
print len(hydrogen)

"""
Explanation: Making selections
To make the RDF, we want to select two different species to calculate the distribution of pairwise distances between.  We will do this based upon the atom's name.
This gives us 2 AtomGroups, both with 512 atoms in.
End of explanation
"""


rdf = InterRDF(oxygen, hydrogen, bins=50, range=(0.0, 12.0))

rdf.run()

plt.plot(rdf.bins, rdf.rdf)

"""
Explanation: Next we create an RDF making object, and pass it our two AtomGroup selections
End of explanation
"""


rdf = InterRDF(oxygen, hydrogen, bins=50, range=(0.0, 12.0), exclusion_block=(1,1))
rdf.run()

plt.plot(rdf.bins, rdf.rdf)

"""
Explanation: Oops, that didn't work.  We've got a single huge peak at an incredibly small distance which obscures everything.  This is caused by the OH and HO bonded pairs which are being included.
These can be excluded using the "exclusion_block" keyword.
End of explanation
"""


