def __init__(self, atoms_n_occu, coords, properties=None):
        """
        Create a *non-periodic* site.

        Args:
            atoms_n_occu: Species on the site. Can be:

                i.  A sequence of element / specie specified either as string
                    symbols, e.g. ["Li", "Fe2+", "P", ...] or atomic numbers,
                    e.g., (3, 56, ...) or actual Element or Specie objects.
                ii. List of dict of elements/species and occupancies, e.g.,
                    [{"Fe" : 0.5, "Mn":0.5}, ...]. This allows the setup of
                    disordered structures.
            coords: Cartesian coordinates of site.
            properties: Properties associated with the site as a dict, e.g.
                {"magmom": 5}. Defaults to None.
        """
        if issubclass(atoms_n_occu.__class__, collections.Mapping):
            self._species = Composition({get_el_sp(k): v
                                         for k, v in atoms_n_occu.items()})
            totaloccu = self._species.num_atoms
            if totaloccu > 1 + Composition.amount_tolerance:
                raise ValueError("Species occupancies sum to more than 1!")
            self._is_ordered = (totaloccu == 1 and len(self._species) == 1)
        else:
            self._species = Composition(
                {get_el_sp(atoms_n_occu): 1})
            self._is_ordered = True

        self._coords = coords
        self._properties = properties if properties else {}

 def _get_c_ranges(self, bonds):
     c_ranges = set()
     bonds = {(get_el_sp(s1), get_el_sp(s2)): dist for (s1, s2), dist in
              bonds.items()}
     for (sp1, sp2), bond_dist in bonds.items():
         for site in self.oriented_unit_cell:
             if sp1 in site.species_and_occu:
                 for nn, d in self.oriented_unit_cell.get_neighbors(
                         site, bond_dist):
                     if sp2 in nn.species_and_occu:
                         c_range = tuple(sorted([site.frac_coords[2],
                                                 nn.frac_coords[2]]))
                         if c_range[1] > 1:
                             # Takes care of PBC when c coordinate of site
                             # goes beyond the upper boundary of the cell
                             c_ranges.add((c_range[0], 1))
                             c_ranges.add((0, c_range[1] - 1))
                         elif c_range[0] < 0:
                             # Takes care of PBC when c coordinate of site
                             # is below the lower boundary of the unit cell
                             c_ranges.add((0, c_range[1]))
                             c_ranges.add((c_range[0] + 1, 1))
                         elif c_range[0] != c_range[1]:
                             c_ranges.add(c_range)
     return c_ranges

def get_bond_length(sp1, sp2, bond_order=1):
    """
    Get the bond length between two species.

    Args:
        sp1 (Specie): First specie.
        sp2 (Specie): Second specie.
        bond_order: For species with different possible bond orders,
            this allows one to obtain the bond length for a particular bond
            order. For example, to get the C=C bond length instead of the
            C-C bond length, this should be set to 2. Defaults to 1.

    Returns:
        Bond length in Angstrom. If no data is available, the sum of the atomic
        radii is used.
    """
    sp1 = get_el_sp(sp1)
    sp2 = get_el_sp(sp2)
    syms = tuple(sorted([sp1.symbol, sp2.symbol]))
    if syms in bond_lengths:
        all_lengths = bond_lengths[syms]
        if bond_order:
            return all_lengths.get(bond_order)
        else:
            return all_lengths.get(1)
    warnings.warn("No bond lengths for %s-%s found in database. Returning sum"
                  "of atomic radius." % (sp1, sp2))
    return sp1.atomic_radius + sp2.atomic_radius

def get_bond_length(sp1, sp2, bond_order=1):
    """
    Get the bond length between two species.

    Args:
        sp1 (Specie): First specie.
        sp2 (Specie): Second specie.
        bond_order: For species with different possible bond orders,
            this allows one to obtain the bond length for a particular bond
            order. For example, to get the C=C bond length instead of the
            C-C bond length, this should be set to 2. Defaults to 1.

    Returns:
        Bond length in Angstrom. None if no data is available.
    """

    syms = tuple(sorted([get_el_sp(sp1).symbol,
                         get_el_sp(sp2).symbol]))
    if syms in bond_lengths:
        all_lengths = bond_lengths[syms]
        if bond_order:
            return all_lengths.get(bond_order)
        else:
            return all_lengths.get(1)
    return None

    def __init__(
        self,
        sp,
        num=None,
        bond_lengths=None,
        interior_only=False,
        midpoints=True,
        min_solid_angle=0,
        site_distance_tol=0.1,
    ):
        """
        Args:
            sp:
                specie to add
            num:
                number to add (will end up with a disordered structure)
            bond_lengths:
                dictionary of bond lengths. For example, if {'O' : 1} is
                used, all points will be scaled from their voronoi locations
                to be exactly 1 angstrom from the oxygen (if the oxygen is
                at the center of the voronoi region). In the case of a point
                neighboring 2 oxygen atoms, 2 points will be added - one for
                each oxygen
            interior_only:
                only uses voronoi points that are inside the tetrahedra
                formed by their 4 nearest neighbors. This is useful because
                when this is not true, there is another point accessible
                without going through a bottleneck that is strictly larger.
                Currently this doesn't work well with bond lengths because 
                finding the correct association between sites and voronoi 
                points is difficult. Doesn't affect the delaunay midpoints
            midpoints:
                whether to add the points at the center of the voronoi
                faces (midpoints of the delaunay triangulation)
            min_solid_angle:
                Threshold for creating points at the delaunay midpoints. The
                solid angle is a metric for how much two sites are neighbors of
                each other (the solid angle swept out by the voronoi face 
                corresponding to the two sites). Typical values are 0 to Pi, with
                Pi corresponding to a tetrahedrally coordinated site.
            site_distance_tol:
                Tolerance passed to to_unit_cell(). When there are multiple
                sites within this tolerance, only one is used
        """
        self._sp = get_el_sp(sp)
        self._n = num
        self._interior = interior_only
        self._midpoints = midpoints
        self._msa = min_solid_angle
        self._sdt = site_distance_tol

        # transform keys to species
        if bond_lengths:
            self._bl = dict(zip(map(lambda x: get_el_sp(x), bond_lengths.keys()), bond_lengths.values()))
        else:
            self._bl = {}

 def apply_transformation(self, structure):
     species_map = {}
     for k, v in self._species_map.items():
         if isinstance(v, dict):
             value = {get_el_sp(x): y for x, y in v.items()}
         else:
             value = get_el_sp(v)
         species_map[get_el_sp(k)] = value
     s = Structure.from_sites(structure.sites)
     s.replace_species(species_map)
     return s

def reduce_formula(sym_amt, iupac_ordering=False):
    """
    Helper method to reduce a sym_amt dict to a reduced formula and factor.

    Args:
        sym_amt (dict): {symbol: amount}.
        iupac_ordering (bool, optional): Whether to order the
            formula by the iupac "electronegativity" series, defined in
            Table VI of "Nomenclature of Inorganic Chemistry (IUPAC
            Recommendations 2005)". This ordering effectively follows
            the groups and rows of the periodic table, except the
            Lanthanides, Actanides and hydrogen. Note that polyanions
            will still be determined based on the true electronegativity of
            the elements.

    Returns:
        (reduced_formula, factor).
    """
    syms = sorted(sym_amt.keys(), key=lambda x: [get_el_sp(x).X, x])

    syms = list(filter(
        lambda x: abs(sym_amt[x]) > Composition.amount_tolerance, syms))

    factor = 1
    # Enforce integers for doing gcd.
    if all((int(i) == i for i in sym_amt.values())):
        factor = abs(gcd(*(int(i) for i in sym_amt.values())))

    polyanion = []
    # if the composition contains a poly anion
    if len(syms) >= 3 and get_el_sp(syms[-1]).X - get_el_sp(syms[-2]).X < 1.65:
        poly_sym_amt = {syms[i]: sym_amt[syms[i]] / factor
                        for i in [-2, -1]}
        (poly_form, poly_factor) = reduce_formula(
            poly_sym_amt, iupac_ordering=iupac_ordering)

        if poly_factor != 1:
            polyanion.append("({}){}".format(poly_form, int(poly_factor)))

    syms = syms[:len(syms) - 2 if polyanion else len(syms)]

    if iupac_ordering:
        syms = sorted(syms,
                      key=lambda x: [get_el_sp(x).iupac_ordering, x])

    reduced_form = []
    for s in syms:
        normamt = sym_amt[s] * 1.0 / factor
        reduced_form.append(s)
        reduced_form.append(formula_double_format(normamt))

    reduced_form = "".join(reduced_form + polyanion)
    return reduced_form, factor

 def __getitem__(self, item):
     try:
         sp = get_el_sp(item)
         return self._data.get(sp, 0)
     except ValueError as ex:
         raise TypeError("Invalid key {}, {} for Composition\n"
                         "ValueError exception:\n{}".format(item, type(item), ex))

 def __contains__(self, item):
     try:
         sp = get_el_sp(item)
         return sp in self._data
     except ValueError:
         raise TypeError("Invalid key {}, {} for Composition\n"
                         "ValueError exception:\n{}".format(item, type(item), ex))

def get_conversion_factor(structure, species, temperature):
    """
    Conversion factor to convert between cm^2/s diffusivity measurements and
    mS/cm conductivity measurements based on number of atoms of diffusing
    species. Note that the charge is based on the oxidation state of the
    species (where available), or else the number of valence electrons
    (usually a good guess, esp for main group ions).

    Args:
        structure (Structure): Input structure.
        species (Element/Specie): Diffusing species.
        temperature (float): Temperature of the diffusion run in Kelvin.

    Returns:
        Conversion factor.
        Conductivity (in mS/cm) = Conversion Factor * Diffusivity (in cm^2/s)
    """
    df_sp = get_el_sp(species)
    if hasattr(df_sp, "oxi_state"):
        z = df_sp.oxi_state
    else:
        z = df_sp.full_electronic_structure[-1][2]

    n = structure.composition[species]

    vol = structure.volume * 1e-24  # units cm^3
    return 1000 * n / (vol * const.N_A) * z ** 2 * (const.N_A * const.e) ** 2 \
        / (const.R * temperature)

    def get_element_spd_dos(self, el):
        """
        Get element and spd projected Dos


        Args:
            el: Element in Structure.composition associated with LobsterCompleteDos

        Returns:
            dict of {Element: {"S": densities, "P": densities, "D": densities}}
        """
        el = get_el_sp(el)
        el_dos = {}
        for site, atom_dos in self.pdos.items():
            if site.specie == el:
                for orb, pdos in atom_dos.items():
                    orbital_type = _get_orb_type_lobster(orb)
                    if orbital_type not in el_dos:
                        el_dos[orbital_type] = pdos
                    else:
                        el_dos[orbital_type] = \
                            add_densities(el_dos[orbital_type], pdos)

        return {orb: Dos(self.efermi, self.energies, densities)
                for orb, densities in el_dos.items()}

    def __init__(self, atoms_n_occu, coords, properties=None):
        """
        Create a *non-periodic* site.

        Args:
            atoms_n_occu: Species on the site. Can be:
                i.  A Composition object (preferred)
                ii. An  element / specie specified either as a string
                    symbols, e.g. "Li", "Fe2+", "P" or atomic numbers,
                    e.g., 3, 56, or actual Element or Specie objects.
                iii.Dict of elements/species and occupancies, e.g.,
                    {"Fe" : 0.5, "Mn":0.5}. This allows the setup of
                    disordered structures.
            coords: Cartesian coordinates of site.
            properties: Properties associated with the site as a dict, e.g.
                {"magmom": 5}. Defaults to None.
        """
        if isinstance(atoms_n_occu, Composition):
            # Compositions are immutable, so don't need to copy (much faster)
            species = atoms_n_occu
        else:
            try:
                species = Composition({get_el_sp(atoms_n_occu): 1})
            except TypeError:
                species = Composition(atoms_n_occu)
        totaloccu = species.num_atoms
        if totaloccu > 1 + Composition.amount_tolerance:
            raise ValueError("Species occupancies sum to more than 1!")
        self._species = species
        self.coords = np.array(coords)
        self.properties = properties or {}

    def get_magnitude_of_effect_from_species(self, species, spin_state, motif):
        """
        Get magnitude of Jahn-Teller effect from provided species, spin state and motife.

        :param species: e.g. Fe2+
        :param spin_state (str): "high" or "low"
        :param motif (str): "oct" or "tet"

        :return (str):
        """

        magnitude = "none"

        sp = get_el_sp(species)

        # has to be Specie; we need to know the oxidation state
        if isinstance(sp, Specie) and sp.element.is_transition_metal:

            d_electrons = self._get_number_of_d_electrons(sp)

            if motif in self.spin_configs:
                if spin_state not in self.spin_configs[motif][d_electrons]:
                    spin_state = self.spin_configs[motif][d_electrons]['default']
                spin_config = self.spin_configs[motif][d_electrons][spin_state]
                magnitude = JahnTellerAnalyzer.get_magnitude_of_effect_from_spin_config(motif,
                                                                                        spin_config)
        else:
            warnings.warn("No data for this species.")

        return magnitude

    def __init__(self, *args, **kwargs):
        """
        Very flexible Composition construction, similar to the built-in Python
        dict(). Also extended to allow simple string init.

        Args:
            Any form supported by the Python built-in dict() function.

            1. A dict of either {Element/Specie: amount},

               {string symbol:amount}, or {atomic number:amount} or any mixture
               of these. E.g., {Element("Li"):2 ,Element("O"):1},
               {"Li":2, "O":1}, {3:2, 8:1} all result in a Li2O composition.
            2. Keyword arg initialization, similar to a dict, e.g.,

               Compostion(Li = 2, O = 1)

            In addition, the Composition constructor also allows a single
            string as an input formula. E.g., Composition("Li2O").
        """
        if len(args) == 1 and isinstance(args[0], basestring):
            elmap = self._parse_formula(args[0])
        else:
            elmap = dict(*args, **kwargs)
        for k, v in elmap.items():
            if v < -Composition.amount_tolerance:
                raise CompositionError("Amounts in Composition cannot be "
                                       "negative!")
            elif v < 0:
                del elmap[k]
        self._elmap = {get_el_sp(k): v for k, v in elmap.items()}
        self._natoms = sum(self._elmap.values())

    def __init__(self, entries, chempots, elements=None):
        """
        Standard constructor for grand potential phase diagram.

        Args:
            entries ([PDEntry]): A list of PDEntry-like objects having an
                energy, energy_per_atom and composition.
            chempots {Element: float}: Specify the chemical potentials
                of the open elements.
            elements ([Element]): Optional list of elements in the phase
                diagram. If set to None, the elements are determined from
                the the entries themselves.
        """
        if elements is None:
            elements = set()
            map(elements.update, [entry.composition.elements
                                  for entry in entries])
        self.chempots = {get_el_sp(el): u for el, u in chempots.items()}
        elements = set(elements).difference(self.chempots.keys())
        all_entries = [GrandPotPDEntry(e, self.chempots)
                       for e in entries
                       if (not e.is_element) or
                       e.composition.elements[0] in elements]

        super(GrandPotentialPhaseDiagram, self).__init__(all_entries, elements)

def reduce_formula(sym_amt):
    """
    Helper method to reduce a sym_amt dict to a reduced formula and factor.

    Args:
        sym_amt (dict): {symbol: amount}.

    Returns:
        (reduced_formula, factor).
    """
    syms = sorted(sym_amt.keys(),
                  key=lambda s: get_el_sp(s).X)

    syms = list(filter(lambda s: abs(sym_amt[s]) >
                       Composition.amount_tolerance, syms))
    num_el = len(syms)
    contains_polyanion = (num_el >= 3 and
                          get_el_sp(syms[num_el - 1]).X
                          - get_el_sp(syms[num_el - 2]).X < 1.65)

    factor = 1
    # Enforce integers for doing gcd.
    if all((int(i) == i for i in sym_amt.values())):
        factor = abs(gcd(*(int(i) for i in sym_amt.values())))

    reduced_form = []
    n = num_el - 2 if contains_polyanion else num_el
    for i in range(0, n):
        s = syms[i]
        normamt = sym_amt[s] * 1.0 / factor
        reduced_form.append(s)
        reduced_form.append(formula_double_format(normamt))

    if contains_polyanion:
        poly_sym_amt = {syms[i]: sym_amt[syms[i]] / factor
                        for i in range(n, num_el)}
        (poly_form, poly_factor) = reduce_formula(poly_sym_amt)

        if poly_factor != 1:
            reduced_form.append("({}){}".format(poly_form, int(poly_factor)))
        else:
            reduced_form.append(poly_form)

    reduced_form = "".join(reduced_form)

    return reduced_form, factor

    def apply_transformation(self, structure, return_ranked_list=False):
        #Make a mutable structure first
        mods = Structure.from_sites(structure)
        for sp, spin in self.mag_species_spin.items():
            sp = get_el_sp(sp)
            oxi_state = getattr(sp, "oxi_state", 0)
            if spin:
                up = Specie(sp.symbol, oxi_state, {"spin": abs(spin)})
                down = Specie(sp.symbol, oxi_state, {"spin": -abs(spin)})
                mods.replace_species(
                    {sp: Composition({up: self.order_parameter,
                                      down: 1 - self.order_parameter})})
            else:
                mods.replace_species(
                    {sp: Specie(sp.symbol, oxi_state, {"spin": spin})})

        if mods.is_ordered:
            return [mods] if return_ranked_list > 1 else mods

        enum_args = self.enum_kwargs

        enum_args["min_cell_size"] = max(int(
            MagOrderingTransformation.determine_min_cell(
                structure, self.mag_species_spin,
                self.order_parameter)),
            enum_args.get("min_cell_size"))

        max_cell = self.enum_kwargs.get('max_cell_size')
        if max_cell:
            if enum_args["min_cell_size"] > max_cell:
                raise ValueError('Specified max cell size is smaller'
                                 ' than the minimum enumerable cell size')
        else:
            enum_args["max_cell_size"] = enum_args["min_cell_size"]

        t = EnumerateStructureTransformation(**enum_args)

        alls = t.apply_transformation(mods,
                                      return_ranked_list=return_ranked_list)

        try:
            num_to_return = int(return_ranked_list)
        except ValueError:
            num_to_return = 1

        if num_to_return == 1 or not return_ranked_list:
            return alls[0]["structure"] if num_to_return else alls

        m = StructureMatcher(comparator=SpinComparator())

        grouped = m.group_structures([d["structure"] for d in alls])

        alls = [{"structure": g[0], "energy": self.emodel.get_energy(g[0])}
                for g in grouped]

        self._all_structures = sorted(alls, key=lambda d: d["energy"])

        return self._all_structures[0:num_to_return]

 def formula(self):
     """
     Returns a formula string, with elements sorted by electronegativity,
     e.g., Li4 Fe4 P4 O16.
     """
     sym_amt = self.get_el_amt_dict()
     syms = sorted(sym_amt.keys(), key=lambda sym: get_el_sp(sym).X)
     formula = [s + formula_double_format(sym_amt[s], False) for s in syms]
     return " ".join(formula)

    def list_prediction(self, species, to_this_composition=True):
        """
        Args:
            species:
                list of species
            to_this_composition:
                If true, substitutions with this as a final composition
                will be found. If false, substitutions with this as a
                starting composition will be found (these are slightly
                different)
        Returns:
            List of predictions in the form of dictionaries.
            If to_this_composition is true, the values of the dictionary
            will be from the list species. If false, the keys will be
            from that list.
        """
        for sp in species:
            if get_el_sp(sp) not in self.p.species:
                raise ValueError("the species {} is not allowed for the"
                                 "probability model you are using".format(sp))
        max_probabilities = []
        for s1 in species:
            if to_this_composition:
                max_p = max([self.p.cond_prob(s2, s1) for s2 in self.p.species])
            else:
                max_p = max([self.p.cond_prob(s1, s2) for s2 in self.p.species])
            max_probabilities.append(max_p)

        output = []

        def _recurse(output_prob, output_species):
            best_case_prob = list(max_probabilities)
            best_case_prob[:len(output_prob)] = output_prob
            if six.moves.reduce(mul, best_case_prob) > self.threshold:
                if len(output_species) == len(species):
                    odict = {
                        'probability': six.moves.reduce(mul, best_case_prob)}
                    if to_this_composition:
                        odict['substitutions'] = dict(
                            zip(output_species, species))
                    else:
                        odict['substitutions'] = dict(
                            zip(species, output_species))
                    if len(output_species) == len(set(output_species)):
                        output.append(odict)
                    return
                for sp in self.p.species:
                    i = len(output_prob)
                    if to_this_composition:
                        prob = self.p.cond_prob(sp, species[i])
                    else:
                        prob = self.p.cond_prob(species[i], sp)
                    _recurse(output_prob + [prob], output_species + [sp])

        _recurse([], [])
        logging.info('{} substitutions found'.format(len(output)))
        return output

 def __init__(self, *args, **kwargs):
     """
     Args:
         *args, **kwargs: any valid dict init arguments
     """
     d = dict(*args, **kwargs)
     super(ChemicalPotential, self).__init__((get_el_sp(k), v) for k, v in d.items())
     if len(d) != len(self):
         raise ValueError("Duplicate potential specified")