Merge pull request #14 from bjmorgan/master

Fix for RDF coordination number calculation

pymatgen_diffusion/aimd/tests/test_van_hove.py

@@ -3,6 +3,7 @@
 # Distributed under the terms of the BSD License.
 
 from __future__ import division, unicode_literals, print_function
+from pymatgen import Structure, Lattice
 
 __author__ = "Iek-Heng Chu"
 __version__ = 1.0
@@ -56,6 +57,14 @@ def test_rdf(self):
         self.assertTrue(check)
         self.assertAlmostEqual(obj.rdf.max(), 1.6831, 4)
 
+    def test_rdf_coordination_number(self):
+        # create a simple cubic lattice
+        coords = np.array( [ [ 0.5, 0.5, 0.5 ] ] )
+        atom_list = [ 'S' ]
+        lattice = Lattice.from_parameters( a=1.0, b=1.0, c=1.0, alpha=90, beta=90, gamma=90 )
+        structure = Structure( lattice, atom_list, coords )
+        rdf = RadialDistributionFunction( structures=[ structure ], species=['S'], rmax=5.0, sigma=0.1, ngrid=500 )
+        self.assertEqual( rdf.coordination_number[100], 6.0 )
 
 class EvolutionAnalyzerTest(unittest.TestCase):
     def test_get_df(self):

pymatgen_diffusion/aimd/van_hove.py

@@ -326,6 +326,7 @@ def __init__(self, structures, ngrid=101, rmax=10.0, cellrange=1, sigma=0.1,
         dr = rmax / (ngrid - 1)
         interval = np.linspace(0.0, rmax, ngrid)
         rdf = np.zeros((ngrid), dtype=np.double)
+        raw_rdf = np.zeros((ngrid), dtype=np.double)
         dns = Counter()
 
         # generate the translational vectors
@@ -364,9 +365,11 @@ def __init__(self, structures, ngrid=101, rmax=10.0, cellrange=1, sigma=0.1,
             rdf[:] += stats.norm.pdf(interval, interval[indx], sigma) * dn \
                       / float(len(ref_indices)) / ff / self.rho / len(
                 fcoords_list)
+            raw_rdf[indx] += dn / float(len(ref_indices)) / ff / self.rho / len(fcoords_list)
 
         self.structures = structures
         self.rdf = rdf
+        self.raw_rdf = raw_rdf
         self.interval = interval
         self.cellrange = cellrange
         self.rmax = rmax
@@ -382,7 +385,7 @@ def coordination_number(self):
         Returns:
             numpy array
         """
-        return np.cumsum(self.rdf * self.rho * 4.0 * np.pi * self.interval ** 2)
+        return np.cumsum(self.raw_rdf * self.rho * 4.0 * np.pi * self.interval ** 2)
 
     def get_rdf_plot(self, label=None, xlim=(0.0, 8.0), ylim=(-0.005, 3.0)):
         """