def tersoff_potential(self, structure):
        """
        Generate the species, tersoff potential lines for an oxide structure

        Args:
            structure:
                pymatgen.core.structure.Structure
        """
        bv = BVAnalyzer()
        el = [site.species_string for site in structure.sites]
        valences = bv.get_valences(structure)
        el_val_dict = dict(zip(el, valences))

        gin = "species \n"
        qerfstring = "qerfc\n"

        for key in el_val_dict.keys():
            if key != "O" and el_val_dict[key] % 1 != 0:
                raise SystemError("Oxide has mixed valence on metal")
            specie_string = key + " core " + str(el_val_dict[key]) + "\n"
            gin += specie_string
            qerfstring += key + " " + key + " 0.6000 10.0000 \n"

        gin += "# noelectrostatics \n Morse \n"
        met_oxi_ters = Tersoff_pot().data
        for key in el_val_dict.keys():
            if key != "O":
                metal = key + "(" + str(int(el_val_dict[key])) + ")"
                ters_pot_str = met_oxi_ters[metal]
                gin += ters_pot_str

        gin += qerfstring
        return gin

    def __init__(self, defect):
        """
        Args:
            defect(Defect): pymatgen Defect object
        """
        self.defect = defect

        try:
            bv = BVAnalyzer()
            struct_valences = bv.get_valences(self.defect.bulk_structure)
            site_index = self.defect.bulk_structure.get_sites_in_sphere(
                self.defect.site.coords, 0.1, include_index=True)[0][2]
            def_site_valence = struct_valences[site_index]
        except Exception:  # sometimes valences cant be assigned
            def_site_valence = 0

        if isinstance(defect, Vacancy):
            self.charges = [-1 * def_site_valence]
        elif isinstance(defect, Substitution):
            #(minimize difference with host site specie)
            probable_chgs = [ox - def_site_valence for ox in self.defect.site.specie.oxidation_states]
            self.charges = [min(probable_chgs, key=abs)]
        elif isinstance(defect, Interstitial):
            self.charges = [0]
        else:
            raise ValueError("Defect Type not recognized.")

def get_basic_analysis_and_error_checks(d):
    initial_vol = d["input"]["crystal"]["lattice"]["volume"]
    final_vol = d["output"]["crystal"]["lattice"]["volume"]
    delta_vol = final_vol - initial_vol
    percent_delta_vol = delta_vol / initial_vol
    coord_num = get_coordination_numbers(d)
    calc = d["calculations"][-1]
    gap = calc["output"]["bandgap"]
    cbm = calc["output"]["cbm"]
    vbm = calc["output"]["vbm"]
    is_direct = calc["output"]["is_gap_direct"]

    if abs(percent_delta_vol) > 0.20:
        warning_msgs = ["Volume change > 20%"]
    else:
        warning_msgs = []

    bv_struct = Structure.from_dict(d["output"]["crystal"])
    try:
        bva = BVAnalyzer()
        bv_struct = bva.get_oxi_state_decorated_structure(bv_struct)
    except ValueError as e:
        logger.error("Valence cannot be determined due to {e}."
                     .format(e=e))
    except Exception as ex:
        logger.error("BVAnalyzer error {e}.".format(e=str(ex)))

    return {"delta_volume": delta_vol,
            "percent_delta_volume": percent_delta_vol,
            "warnings": warning_msgs, "coordination_numbers": coord_num,
            "bandgap": gap, "cbm": cbm, "vbm": vbm,
            "is_gap_direct": is_direct,
            "bv_structure": bv_struct.to_dict}

def get_basic_analysis_and_error_checks(d, max_force_threshold=0.5,
                                        volume_change_threshold=0.2):

    initial_vol = d["input"]["crystal"]["lattice"]["volume"]
    final_vol = d["output"]["crystal"]["lattice"]["volume"]
    delta_vol = final_vol - initial_vol
    percent_delta_vol = delta_vol / initial_vol
    coord_num = get_coordination_numbers(d)
    calc = d["calculations"][-1]
    gap = calc["output"]["bandgap"]
    cbm = calc["output"]["cbm"]
    vbm = calc["output"]["vbm"]
    is_direct = calc["output"]["is_gap_direct"]

    warning_msgs = []
    error_msgs = []

    if abs(percent_delta_vol) > volume_change_threshold:
        warning_msgs.append("Volume change > {}%"
                            .format(volume_change_threshold * 100))

    bv_struct = Structure.from_dict(d["output"]["crystal"])
    try:
        bva = BVAnalyzer()
        bv_struct = bva.get_oxi_state_decorated_structure(bv_struct)
    except ValueError as e:
        logger.error("Valence cannot be determined due to {e}."
                     .format(e=e))
    except Exception as ex:
        logger.error("BVAnalyzer error {e}.".format(e=str(ex)))

    max_force = None
    if d["state"] == "successful" and \
            d["calculations"][0]["input"]["parameters"].get("NSW", 0) > 0:
        # handle the max force and max force error
        max_force = max([np.linalg.norm(a)
                        for a in d["calculations"][-1]["output"]
                        ["ionic_steps"][-1]["forces"]])

        if max_force > max_force_threshold:
            error_msgs.append("Final max force exceeds {} eV"
                              .format(max_force_threshold))
            d["state"] = "error"

        s = Structure.from_dict(d["output"]["crystal"])
        if not s.is_valid():
            error_msgs.append("Bad structure (atoms are too close!)")
            d["state"] = "error"

    return {"delta_volume": delta_vol,
            "max_force": max_force,
            "percent_delta_volume": percent_delta_vol,
            "warnings": warning_msgs,
            "errors": error_msgs,
            "coordination_numbers": coord_num,
            "bandgap": gap, "cbm": cbm, "vbm": vbm,
            "is_gap_direct": is_direct,
            "bv_structure": bv_struct.as_dict()}

 def _get_valences(self):
     """
     Computes ionic valences of elements for all sites in the structure.
     """
     bv = BVAnalyzer()
     valences = bv.get_valences(self._structure)
     el = [site.species_string for site in self.structure.sites]
     valence_dict = dict(zip(el, valences))
     return valence_dict

 def setUp(self):
     filepath1 = os.path.join(test_dir, 'Li2O.cif')
     p = CifParser(filepath1).get_structures(False)[0]
     bv = BVAnalyzer()
     valences = bv.get_valences(p)
     el = [site.species_string for site in p.sites]
     val_dict = dict(zip(el, valences))
     self._radii = {}
     for k,v in val_dict.items():
         k1 = re.sub('[1-9,+,\-]', '', k)
         self._radii[k] = Specie(k1, v).ionic_radius
     p.remove(0)
     self._vac_struct = p

 def setUp(self):
     filepath = os.path.join(test_dir, 'POSCAR')
     p = Poscar.from_file(filepath)
     self.structure = p.structure
     bv = BVAnalyzer()
     self.structure = bv.get_oxi_state_decorated_structure(self.structure)
     valences = bv.get_valences(self.structure)
     radii = []
     for i in range(len(valences)):
         el = self.structure.sites[i].specie.symbol
         radius = Specie(el, valences[i]).ionic_radius
         radii.append(radius)
     el = [site.species_string for site in self.structure.sites]
     self.rad_dict = dict(zip(el, radii))
     for el in self.rad_dict.keys():
         print((el, self.rad_dict[el].real))

    def _get_valences(self):
        """
        Computes ionic valences of elements for all sites in the structure.
        """
        bv = BVAnalyzer()
        try:
            valences = bv.get_valences(self._structure)
        except:
            err_str = "BVAnalyzer failed. The defect effective charge, and"
            err_str += " volume and surface area may not work"
            print err_str
            raise LookupError()

        el = [site.species_string for site in self.structure.sites]
        valence_dict = dict(zip(el, valences))
        return valence_dict

    def setUp(self):
        """
        Setup MgO rocksalt structure for testing Vacancy
        """
        mgo_latt = [[4.212, 0, 0], [0, 4.212, 0], [0, 0, 4.212]]
        mgo_specie = ["Mg"] * 4 + ["O"] * 4
        mgo_frac_cord = [[0, 0, 0], [0.5, 0.5, 0], [0.5, 0, 0.5], [0, 0.5, 0.5],
                         [0.5, 0, 0], [0, 0.5, 0], [0, 0, 0.5], [0.5, 0.5, 0.5]]
        self._mgo_uc = Structure(mgo_latt, mgo_specie, mgo_frac_cord, True,
                                 True)

        bv = BVAnalyzer()
        self._mgo_uc = bv.get_oxi_state_decorated_structure(self._mgo_uc)
        self._mgo_val_rad_eval = ValenceIonicRadiusEvaluator(self._mgo_uc)
        self._mgo_val = self._mgo_val_rad_eval.valences
        self._mgo_rad = self._mgo_val_rad_eval.radii
        self._mgo_vac = Vacancy(self._mgo_uc, self._mgo_val, self._mgo_rad)

 def __init__(self, symm_tol=0.1, max_radius=4, max_permutations=100000,
              distance_scale_factor=1.015):
     self.symm_tol = symm_tol
     self.max_radius = max_radius
     self.max_permutations = max_permutations
     self.distance_scale_factor = distance_scale_factor
     self.analyzer = BVAnalyzer(symm_tol, max_radius, max_permutations,
                                distance_scale_factor)

    def predict(self, structure, ref_structure, test_isostructural=True):
        """
        Given a structure, returns back the predicted volume.

        Args:
            structure (Structure): structure w/unknown volume
            ref_structure (Structure): A reference structure with a similar
                structure but different species.
            test_isostructural (bool): Whether to test that the two
                structures are isostructural. This algo works best for
                isostructural compounds. Defaults to True.

        Returns:
            a float value of the predicted volume
        """
        if not is_ox(structure):
            a = BVAnalyzer()
            structure = a.get_oxi_state_decorated_structure(structure)
        if not is_ox(ref_structure):
            a = BVAnalyzer()
            ref_structure = a.get_oxi_state_decorated_structure(ref_structure)

        if test_isostructural:
            m = StructureMatcher()
            mapping = m.get_best_electronegativity_anonymous_mapping(structure, ref_structure)
            if mapping is None:
                raise ValueError("Input structures do not match!")

        comp = structure.composition
        ref_comp = ref_structure.composition

        numerator = 0
        denominator = 0

        # Here, the 1/3 factor on the composition accounts for atomic
        # packing. We want the number per unit length.

        # TODO: AJ doesn't understand the (1/3). It would make sense to him
        # if you were doing atomic volume and not atomic radius
        for k, v in comp.items():
            numerator += k.ionic_radius * v ** (1 / 3)
        for k, v in ref_comp.items():
            denominator += k.ionic_radius * v ** (1 / 3)

        # The scaling factor is based on lengths. We apply a power of 3.
        return ref_structure.volume * (numerator / denominator) ** 3

    def __init__(self, structure, valences, radii):
        """
        Given a structure, generate symmetrically distinct interstitial sites.

        Args:
            structure: pymatgen.core.structure.Structure
            valences: Dictionary of oxidation states of elements in {
                El:valence} form
            radii: Radii of elemnts in the structure
        """

        bv = BVAnalyzer()
        self._structure = bv.get_oxi_state_decorated_structure(structure)
        #self._structure = structure
        self._valence_dict = valences
        self._rad_dict = radii

        #Use Zeo++ to obtain the voronoi nodes. Apply symmetry reduction and
        #the symmetry reduced voronoi nodes
        #are possible candidates for interstitial sites
        #try:
        possible_interstitial_sites = symmetry_reduced_voronoi_nodes(
                self._structure, self._rad_dict)
        #except:
        #    raise ValueError("Symmetry_reduced_voronoi_nodes failed")

        #Do futher processing on possibleInterstitialSites to obtain
        #interstitial sites
        self._defect_sites = possible_interstitial_sites
        self._defectsite_coord_no = []
        self._defect_coord_sites = []
        self._defect_coord_charge = []
        self._radii = []

        for site in self._defect_sites:
            coord_no, coord_sites, chrg = self._get_coord_no_sites_chrg(site)
            self._defectsite_coord_no.append(coord_no)
            self._defect_coord_sites.append(coord_sites)
            self._defect_coord_charge.append(chrg)

        for site in self._defect_sites:
            self._radii.append(float(site.properties['voronoi_radius']))

    def _get_valences(self):
        """
        Computes ionic valences of elements for all sites in the structure.
        """
        bv = BVAnalyzer()
        self._structure = bv.get_oxi_state_decorated_structure(self._structure)
        try:
            valences = bv.get_valences(self._structure)
        except:
            try:
                valences = bv.get_valences(self._structure, symm_tol=0.0)
            except:
                raise

        #print valences
        #el = [site.specie.symbol for site in self._structure.sites]
        #el = [site.species_string for site in self._structure.sites]
        #el = [site.specie for site in self._structure.sites]
        #valence_dict = dict(zip(el, valences))
        #print valence_dict
        return valences

class BVAnalyzerTest(PymatgenTest):

    def setUp(self):
        self.analyzer = BVAnalyzer()

    def test_get_valence(self):
        s = Structure.from_file(os.path.join(test_dir, "LiMn2O4.json"))
        ans = [1, 1, 3, 3, 4, 4, -2, -2, -2, -2, -2, -2, -2, -2]
        self.assertEqual(self.analyzer.get_valences(s), ans)
        s = self.get_structure("LiFePO4")
        ans = [1, 1, 1, 1, 2, 2, 2, 2, 5, 5, 5, 5, -2, -2, -2, -2, -2, -2, -2,
               - 2, -2, -2, -2, -2, -2, -2, -2, -2]
        self.assertEqual(self.analyzer.get_valences(s), ans)
        s = self.get_structure("Li3V2(PO4)3")
        ans = [1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 5, 5, 5, 5, 5, 5, -2, -2, -2, -2,
               - 2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2,
               - 2, -2, -2, -2]
        self.assertEqual(self.analyzer.get_valences(s), ans)
        s = Structure.from_file(os.path.join(test_dir, "Li4Fe3Mn1(PO4)4.json"))
        ans = [1, 1, 1, 1, 2, 2, 2, 2, 5, 5, 5, 5, -2, -2, -2, -2, -2, -2, -2,
               - 2, -2, -2, -2, -2, -2, -2, -2, -2]
        self.assertEqual(self.analyzer.get_valences(s), ans)
        s = self.get_structure("NaFePO4")
        ans = [1, 1, 1, 1, 2, 2, 2, 2, 5, 5, 5, 5, -2, -2, -2, -2, -2, -2, -2,
               - 2, -2, -2, -2, -2, -2, -2, -2, -2]
        self.assertEqual(self.analyzer.get_valences(s), ans)

    def test_get_oxi_state_structure(self):
        s = Structure.from_file(os.path.join(test_dir, "LiMn2O4.json"))
        news = self.analyzer.get_oxi_state_decorated_structure(s)
        self.assertIn(Specie("Mn", 3), news.composition.elements)
        self.assertIn(Specie("Mn", 4), news.composition.elements)

    def get_coordsites_min_max_charge(self, n):
        """
        Minimum and maximum charge of sites surrounding the vacancy site.

        Args:
            n: Index of vacancy list
        """
        bv = BVAnalyzer()
        struct_valences = bv.get_valences(self._structure)
        coordinated_site_valences = []

        def _get_index(site):
            for i in range(len(self._structure.sites)):
                if site.is_periodic_image(self._structure.sites[i]):
                    return i
            raise ValueError("Site not found")

        for site in self._defect_coord_sites[n]:
            ind = _get_index(site)
            coordinated_site_valences.append(struct_valences[ind])
        coordinated_site_valences.sort()
        return coordinated_site_valences[0], coordinated_site_valences[-1]

class AutoOxiStateDecorationTransformation(AbstractTransformation):
    """
    This transformation automatically decorates a structure with oxidation
    states using a bond valence approach.
    """

    def __init__(self, symm_tol=0.1, max_radius=4, max_permutations=100000, distance_scale_factor=1.015):
        """
        Args:
            symm_tol:
                Symmetry tolerance used to determine which sites are
                symmetrically equivalent. Set to 0 to turn off symmetry.
            max_radius:
                Maximum radius in Angstrom used to find nearest neighbors.
            max_permutations:
                The maximum number of permutations of oxidation states to test.
            distance_scale_factor:
                A scale factor to be applied. This is useful for scaling
                distances, esp in the case of calculation-relaxed structures
                which may tend to under (GGA) or over bind (LDA). The default
                of 1.015 works for GGA. For experimental structure, set this to
                1.
        """
        self.analyzer = BVAnalyzer(symm_tol, max_radius, max_permutations, distance_scale_factor)

    def apply_transformation(self, structure):
        return self.analyzer.get_oxi_state_decorated_structure(structure)

    @property
    def inverse(self):
        return None

    @property
    def is_one_to_many(self):
        return False

    @property
    def to_dict(self):
        return {
            "name": self.__class__.__name__,
            "version": __version__,
            "init_args": {
                "symm_tol": self.analyzer.symm_tol,
                "max_radius": self.analyzer.max_radius,
                "max_permutations": self.analyzer.max_permutations,
                "distance_scale_factor": self.analyzer.dist_scale_factor,
            },
            "@module": self.__class__.__module__,
            "@class": self.__class__.__name__,
        }

 def __init__(self, symm_tol=0.1, max_radius=4, max_permutations=100000, distance_scale_factor=1.015):
     """
     Args:
         symm_tol:
             Symmetry tolerance used to determine which sites are
             symmetrically equivalent. Set to 0 to turn off symmetry.
         max_radius:
             Maximum radius in Angstrom used to find nearest neighbors.
         max_permutations:
             The maximum number of permutations of oxidation states to test.
         distance_scale_factor:
             A scale factor to be applied. This is useful for scaling
             distances, esp in the case of calculation-relaxed structures
             which may tend to under (GGA) or over bind (LDA). The default
             of 1.015 works for GGA. For experimental structure, set this to
             1.
     """
     self.analyzer = BVAnalyzer(symm_tol, max_radius, max_permutations, distance_scale_factor)

class AutoOxiStateDecorationTransformation(AbstractTransformation):
    """
    This transformation automatically decorates a structure with oxidation
    states using a bond valence approach.

    Args:
        symm_tol (float): Symmetry tolerance used to determine which sites are
            symmetrically equivalent. Set to 0 to turn off symmetry.
        max_radius (float): Maximum radius in Angstrom used to find nearest
            neighbors.
        max_permutations (int): Maximum number of permutations of oxidation
            states to test.
        distance_scale_factor (float): A scale factor to be applied. This is
            useful for scaling distances, esp in the case of
            calculation-relaxed structures, which may tend to under (GGA) or
            over bind (LDA). The default of 1.015 works for GGA. For
            experimental structure, set this to 1.
    """

    def __init__(self, symm_tol=0.1, max_radius=4, max_permutations=100000,
                 distance_scale_factor=1.015):
        self.symm_tol = symm_tol
        self.max_radius = max_radius
        self.max_permutations = max_permutations
        self.distance_scale_factor = distance_scale_factor
        self.analyzer = BVAnalyzer(symm_tol, max_radius, max_permutations,
                                   distance_scale_factor)

    def apply_transformation(self, structure):
        return self.analyzer.get_oxi_state_decorated_structure(structure)

    @property
    def inverse(self):
        return None

    @property
    def is_one_to_many(self):
        return False

class BVAnalyzerTest(unittest.TestCase):

    def setUp(self):
        self.analyzer = BVAnalyzer()

    def test_get_valence(self):
        parser = CifParser(os.path.join(test_dir, "LiMn2O4.cif"))
        s = parser.get_structures()[0]
        ans = [1, 1, 3, 3, 4, 4, -2, -2, -2, -2, -2, -2, -2, -2]
        self.assertEqual(self.analyzer.get_valences(s), ans)
        parser = CifParser(os.path.join(test_dir, "LiFePO4.cif"))
        s = parser.get_structures()[0]
        ans = [1, 1, 1, 1, 2, 2, 2, 2, 5, 5, 5, 5, -2, -2, -2, -2, -2, -2, -2,
               - 2, -2, -2, -2, -2, -2, -2, -2, -2]
        self.assertEqual(self.analyzer.get_valences(s), ans)
        parser = CifParser(os.path.join(test_dir, "Li3V2(PO4)3.cif"))
        s = parser.get_structures()[0]
        ans = [1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 5, 5, 5, 5, 5, 5, -2, -2, -2, -2,
               - 2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2,
               - 2, -2, -2, -2]
        self.assertEqual(self.analyzer.get_valences(s), ans)
        parser = CifParser(os.path.join(test_dir, "Li4Fe3Mn1(PO4)4.cif"))
        s = parser.get_structures()[0]
        ans = [1, 1, 1, 1, 2, 2, 2, 2, 5, 5, 5, 5, -2, -2, -2, -2, -2, -2, -2,
               - 2, -2, -2, -2, -2, -2, -2, -2, -2]
        self.assertEqual(self.analyzer.get_valences(s), ans)
        parser = CifParser(os.path.join(test_dir, "NaFePO4.cif"))
        s = parser.get_structures()[0]
        ans = [1, 1, 1, 1, 2, 2, 2, 2, 5, 5, 5, 5, -2, -2, -2, -2, -2, -2, -2,
               - 2, -2, -2, -2, -2, -2, -2, -2, -2]
        self.assertEqual(self.analyzer.get_valences(s), ans)

    def test_get_oxi_state_structure(self):
        parser = CifParser(os.path.join(test_dir, "LiMn2O4.cif"))
        s = parser.get_structures()[0]
        news = self.analyzer.get_oxi_state_decorated_structure(s)
        self.assertIn(Specie("Mn", 3), news.composition.elements)
        self.assertIn(Specie("Mn", 4), news.composition.elements)

    def __init__(self, structure, valences, radii, site_type='voronoi_vertex',
                 accuracy='Normal', symmetry_flag=True):
        """
        Given a structure, generate symmetrically distinct interstitial sites.

        Args:
            structure: pymatgen.core.structure.Structure
            valences: Dictionary of oxidation states of elements in 
                {el:valence} form
            radii: Radii of elemnts in the structure
            site_type: "voronoi_vertex" uses voronoi nodes
                "voronoi_facecenter" uses voronoi polyhedra face centers
                Default is "voronoi_vertex"
            accuracy: Flag denoting whether to use high accuracy version 
                of Zeo++. Options are "Normal" and "High". Default is normal.
        """
        try:
            bv = BVAnalyzer()
            self._structure = bv.get_oxi_state_decorated_structure(structure)
        except:
            try:
                bv = BVAnalyzer(symm_tol=0.0)
                self._structure = bv.get_oxi_state_decorated_structure(
                        structure
                        )
            except:
                raise
        self._valence_dict = valences
        self._rad_dict = radii

        """
        Use Zeo++ to obtain the voronoi nodes. Apply symmetry reduction
        and the symmetry reduced voronoi nodes are possible candidates
        for interstitial sites.
        """
        if accuracy == "Normal":
            high_accuracy_flag = False
        elif accuracy == "High":
            high_accuracy_flag = True
        else:
            raise ValueError("Accuracy setting not understood.")

        vor_node_sites, vor_facecenter_sites = symmetry_reduced_voronoi_nodes(
                self._structure, self._rad_dict, high_accuracy_flag, symmetry_flag
                )
        
        if site_type == 'voronoi_vertex':
            possible_interstitial_sites = vor_node_sites
        elif site_type == 'voronoi_facecenter':
            possible_interstitial_sites = vor_facecenter_sites
        else:
            raise ValueError("Input site type not implemented")

        #Do futher processing on possibleInterstitialSites to obtain
        #interstitial sites
        self._defect_sites = possible_interstitial_sites
        self._defectsite_coord_no = []
        self._defect_coord_sites = []
        self._defect_coord_charge = []
        self._radii = []

        for site in self._defect_sites:
            coord_no, coord_sites, chrg = self._get_coord_no_sites_chrg(site)
            self._defectsite_coord_no.append(coord_no)
            self._defect_coord_sites.append(coord_sites)
            self._defect_coord_charge.append(chrg)

        for site in self._defect_sites:
            self._radii.append(float(site.properties['voronoi_radius']))