def _get_coord_no_sites_chrg(self, site):
        """
        Compute the coordination number and coordination charge

        Args:
            site:
                pymatgen.core.sites.Site
        """
        struct = self._structure.copy()
        struct.append(site.specie.symbol, site.frac_coords)
        coord_finder = VoronoiCoordFinder(struct)
        coord_no = coord_finder.get_coordination_number(-1)
        coord_sites = coord_finder.get_coordinated_sites(-1)

        # In some cases coordination sites to interstitials include
        # interstitials also. Filtering them.
        def no_inter(site):
            return not site.specie.symbol == 'X'
        coord_sites = filter(no_inter, coord_sites)

        coord_chrg = 0
        for site, weight in coord_finder.get_voronoi_polyhedra(-1).items():
            if not site.specie.symbol == 'X':
                coord_chrg += weight * self._valence_dict[site.species_string]

        return coord_no, coord_sites, coord_chrg

    def _get_coord_no_sites_chrg(self, site):
        """
        Compute the coordination number and coordination charge

        Args:
            site:
                pymatgen.core.sites.Site
        """
        struct = self._structure.copy()
        struct.append(site.species_string, site.frac_coords)
        coord_finder = VoronoiCoordFinder(struct)
        coord_no = coord_finder.get_coordination_number(-1)
        coord_sites = coord_finder.get_coordinated_sites(-1)

        # In some cases coordination sites to interstitials include 
        # interstitials also. 
        sites_to_be_deleted = []
        for i in range(len(coord_sites)):
            if coord_sites[i].species_string == 'X':
                sites_to_be_deleted.append(i)
        sites_to_be_deleted.reverse()  # So index won't change in between
        for ind in sites_to_be_deleted:
            del coord_sites[ind]

        coord_chrg = 0
        for site, weight in coord_finder.get_voronoi_polyhedra(-1).items():
            if not site.species_string == 'X':
                coord_chrg += weight * self._valence_dict[site.species_string]

        return coord_no, coord_sites, coord_chrg

    def from_bulk_and_miller(cls, structure, miller_index, min_slab_size=5.0,
                             min_vacuum_size=10.0, max_normal_search=None, 
                             center_slab = True, selective_dynamics=False,
                             undercoord_threshold = 0.09):
        """
        This method constructs the adsorbate site finder from a bulk 
        structure and a miller index, which allows the surface sites 
        to be determined from the difference in bulk and slab coordination, 
        as opposed to the height threshold.
        
        Args:
            structure (Structure): structure from which slab
                input to the ASF is constructed
            miller_index (3-tuple or list): miller index to be used
            min_slab_size (float): min slab size for slab generation
            min_vacuum_size (float): min vacuum size for slab generation
            max_normal_search (int): max normal search for slab generation
            center_slab (bool): whether to center slab in slab generation
            selective dynamics (bool): whether to assign surface sites
                to selective dynamics
            undercoord_threshold (float): threshold of "undercoordation"
                to use for the assignment of surface sites.  Default is
                0.1, for which surface sites will be designated if they
                are 10% less coordinated than their bulk counterpart
        """
        # TODO: for some reason this works poorly with primitive cells
        vcf_bulk = VoronoiCoordFinder(structure)
        bulk_coords = [len(vcf_bulk.get_coordinated_sites(n))
                       for n in range(len(structure))]
        struct = structure.copy(site_properties = {'bulk_coordinations':bulk_coords})
        slabs = generate_all_slabs(struct, max_index=max(miller_index), 
                                   min_slab_size=min_slab_size,
                                   min_vacuum_size=min_vacuum_size,
                                   max_normal_search = max_normal_search,
                                   center_slab = center_slab)

        slab_dict = {slab.miller_index:slab for slab in slabs}
        
        if miller_index not in slab_dict:
            raise ValueError("Miller index not in slab dict")

        this_slab = slab_dict[miller_index]

        vcf_surface = VoronoiCoordFinder(this_slab, allow_pathological=True)
        surf_props = []
        this_mi_vec = get_mi_vec(this_slab.miller_index)
        mi_mags = [np.dot(this_mi_vec, site.coords) for site in this_slab]
        average_mi_mag = np.average(mi_mags)
        for n, site in enumerate(this_slab):
            bulk_coord = this_slab.site_properties['bulk_coordinations'][n]
            slab_coord = len(vcf_surface.get_coordinated_sites(n))
            mi_mag = np.dot(this_mi_vec, site.coords)
            undercoord = (bulk_coord - slab_coord)/bulk_coord
            if undercoord > undercoord_threshold and mi_mag > average_mi_mag:
                surf_props += ['surface']
            else:
                surf_props += ['subsurface']
        new_site_properties = {'surface_properties':surf_props}
        new_slab = this_slab.copy(site_properties=new_site_properties)
        return cls(new_slab, selective_dynamics)

    def _get_ionic_radii(self):
        """
        Computes ionic radii of elements for all sites in the structure.
        If valence is zero, atomic radius is used.
        """
        radii = []
        coord_finder = VoronoiCoordFinder(self._structure)

        def nearest_key(sorted_vals, key):
            i = bisect_left(sorted_vals, key)
            if i == len(sorted_vals):
                i = -1
            return sorted_vals[i]

        for i in range(len(self._structure.sites)):
            site = self._structure.sites[i]
            el = site.specie.symbol
            oxi_state = int(round(site.specie.oxi_state))
            coord_no = int(round(coord_finder.get_coordination_number(i)))
            try:
                tab_oxi_states = map(int, _ion_radii[el].keys())
                tab_oxi_states.sort()
                oxi_state = nearest_key(tab_oxi_states, oxi_state)
                radius = _ion_radii[el][str(oxi_state)][str(coord_no)]
            except KeyError:
                if coord_finder.get_coordination_number(i)-coord_no > 0:
                    new_coord_no = coord_no + 1
                else:
                    new_coord_no = coord_no - 1
                try:
                    radius = _ion_radii[el][str(oxi_state)][str(new_coord_no)]
                    coord_no = new_coord_no
                except:
                    tab_coords = map(int, _ion_radii[el][str(oxi_state)].keys())
                    tab_coords.sort()
                    new_coord_no = nearest_key(tab_coords, coord_no)
                    i = 0
                    for val in tab_coords:
                        if  val > coord_no:
                            break
                        i = i + 1
                    if i == len(tab_coords):
                        key = str(tab_coords[-1])
                        radius = _ion_radii[el][str(oxi_state)][key]
                    elif i == 0:
                        key = str(tab_coords[0])
                        radius = _ion_radii[el][str(oxi_state)][key]
                    else:
                        key = str(tab_coords[i-1])
                        radius1 = _ion_radii[el][str(oxi_state)][key]
                        key = str(tab_coords[i])
                        radius2 = _ion_radii[el][str(oxi_state)][key]
                        radius = (radius1+radius2)/2

            #implement complex checks later
            radii.append(radius)
        return radii

 def compute(self, structure, n):
     params = self._params
     vor = VoronoiCoordFinder(structure, **params)
     vorp = vor.get_voronoi_polyhedra(n)
     cdict = {}
     for i in vorp:
         if i.species_string not in cdict:
             cdict[i.species_string] = vorp[i]
         else:
             cdict[i.species_string] += vorp[i]
     return cdict

 def _coord_find(self):
     """
     calls VoronoiCoordFinder to compute the coordination number,
     coordination charge
     """
     for i in range(self.defect_count()):
         struct = self._structs[i].copy()
         coord_finder = VoronoiCoordFinder(struct)
         self._coord_no.append(coord_finder.get_coordination_number(-1))
         self._coord_sites.append(coord_finder.get_coordinated_sites(-1))
         coord_chrg = 0
         for site, weight in coord_finder.get_voronoi_polyhedra(-1).items():
             coord_chrg += weight * self._valence_dict[site.species_string]
         self._coord_charge_no.append(coord_chrg)

class VoronoiCoordFinderTest(PymatgenTest):
    def setUp(self):
        s = self.get_structure('LiFePO4')
        self.finder = VoronoiCoordFinder(s, [Element("O")])

    def test_get_voronoi_polyhedra(self):
        self.assertEqual(len(self.finder.get_voronoi_polyhedra(0).items()), 8)

    def test_get_coordination_number(self):
        self.assertAlmostEqual(self.finder.get_coordination_number(0),
                               5.809265748999465, 7)

    def test_get_coordinated_sites(self):
        self.assertEqual(len(self.finder.get_coordinated_sites(0)), 8)

    def __init__(self, structure, valences, radii):
        """
        Args:
            structure:
                pymatgen.core.structure.Structure
            valences:
                valences of elements as a dictionary 
            radii:
                Radii of elements as a dictionary
        """

        self._structure = structure
        self._valence_dict = valences
        self._rad_dict = radii
        # Store symmetrically distinct sites, their coordination numbers
        # coordinated_sites, effective charge
        symm_finder = SymmetryFinder(self._structure)
        symm_structure = symm_finder.get_symmetrized_structure()
        equiv_site_seq = symm_structure.equivalent_sites

        self._defect_sites = []
        for equiv_sites in equiv_site_seq:
            self._defect_sites.append(equiv_sites[0])

        self._vac_site_indices = []
        for site in self._defect_sites:
            for i in range(len(self._structure.sites)):
                if site == self._structure[i]:
                    self._vac_site_indices.append(i)

        coord_finder = VoronoiCoordFinder(self._structure)
        self._defectsite_coord_no = []
        self._defect_coord_sites = []
        for i in self._vac_site_indices:
            self._defectsite_coord_no.append(
                coord_finder.get_coordination_number(i)
            )
            self._defect_coord_sites.append(
                coord_finder.get_coordinated_sites(i)
            )

        # Store the ionic radii for the elements in the structure
        # (Used to  computing the surface are and volume)
        # Computed based on valence of each element

        self._vac_eff_charges = None
        self._vol = None
        self._sa = None

class VoronoiCoordFinderTest(unittest.TestCase):

    def setUp(self):
        filepath = os.path.join(test_dir, 'vasprun.xml')
        reader = Vasprun(filepath)
        s = reader.final_structure
        self.finder = VoronoiCoordFinder(s,[Element("O")])
    
    def test_get_voronoi_polyhedra(self):
        self.assertEqual(len(self.finder.get_voronoi_polyhedra(0).items()),10, "Incorrect number of results returned for get_voronoi_polyhedra")
        
    def test_get_coordination_number(self):
        print self.finder.get_coordination_number(0)
        self.assertAlmostEqual(self.finder.get_coordination_number(0), 5.60588600732, 7, "Incorrect coordination number returned!")                

    def test_get_coordinated_sites(self):
        self.assertEqual(len(self.finder.get_coordinated_sites(0)), 10)

def substructures_from_structure(structure, weight_cutoff=1e-2):
    """
    Helper method to calculate substructural components from a
    pymatgen structure object.

    Args:
        structure (Structure): Input structure
        weight_cutoff (float): minimum solid angle weight for
                               inclusion in a substructure

    Returns:
            A list of Substructure objects for the given structure
    """

    vcf = VoronoiCoordFinder(structure)
    substructures = []

    for i, site in enumerate(structure):

        charge = sum([getattr(specie, "oxi_state", 0) * amt
                      for specie, amt in site.species_and_occu.items()])

        central_subspecies = SubStructureSpecie(site.specie.symbol,
                                         oxidation_state=charge,
                                         properties=site.properties,
                                         weight=1.0)

        peripheral_subspecies = []
        for peripheral_site, weight in vcf.get_voronoi_polyhedra(i).iteritems():

            if weight > weight_cutoff:
                charge = sum([getattr(specie, "oxi_state", 0) * amt
                              for specie, amt in
                              peripheral_site.species_and_occu.items()])

                peripheral_subspecies.append(
                    SubStructureSpecie(peripheral_site.specie.symbol,
                                       oxidation_state=charge,
                                       properties=peripheral_site.properties,
                                       weight=weight))

        substructures.append(SubStructure(peripheral_subspecies, central_sub_species=central_subspecies))

    return substructures

class VoronoiCoordFinderTest(unittest.TestCase):

    def setUp(self):
        filepath = os.path.join(test_dir, 'LiFePO4.cif')
        parser = CifParser(filepath)
        s = parser.get_structures()[0]
        self.finder = VoronoiCoordFinder(s, [Element("O")])

    def test_get_voronoi_polyhedra(self):
        self.assertEqual(len(self.finder.get_voronoi_polyhedra(0).items()), 8,
                         "Incorrect number of results returned for " +
                         "get_voronoi_polyhedra")

    def test_get_coordination_number(self):
        self.assertAlmostEqual(self.finder.get_coordination_number(0),
                               5.809265748999465, 7)

    def test_get_coordinated_sites(self):
        self.assertEqual(len(self.finder.get_coordinated_sites(0)), 8)

def get_coordination_numbers(d):
    """
    Helper method to get the coordination number of all sites in the final
    structure from a run.

    Args:
        d:
            Run dict generated by VaspToDbTaskDrone.

    Returns:
        Coordination numbers as a list of dict of [{"site": site_dict,
        "coordination": number}, ...].
    """
    structure = Structure.from_dict(d["output"]["crystal"])
    f = VoronoiCoordFinder(structure)
    cn = []
    for i, s in enumerate(structure.sites):
        try:
            n = f.get_coordination_number(i)
            number = int(round(n))
            cn.append({"site": s.to_dict, "coordination": number})
        except Exception:
            logger.error("Unable to parse coordination errors")
    return cn

 def setUp(self):
     s = self.get_structure('LiFePO4')
     self.finder = VoronoiCoordFinder(s, [Element("O")])

 def setUp(self):
     filepath = os.path.join(test_dir, 'LiFePO4.cif')
     parser = CifParser(filepath)
     s = parser.get_structures()[0]
     self.finder = VoronoiCoordFinder(s, [Element("O")])

 def setUp(self):
     filepath = os.path.join(test_dir, 'vasprun.xml')
     reader = Vasprun(filepath)
     s = reader.final_structure
     self.finder = VoronoiCoordFinder(s,[Element("O")])