def __init__(self, lambda_table=None, alpha=-5):
        if lambda_table is not None:
            self._lambda_table = lambda_table
        else:
            module_dir = os.path.dirname(__file__)
            json_file = os.path.join(module_dir, 'data', 'lambda.json')
            with open(json_file) as f:
                self._lambda_table = json.load(f)

        # build map of specie pairs to lambdas
        self.alpha = alpha
        self._l = {}
        self.species = set()
        for row in self._lambda_table:
            if 'D1+' not in row:
                s1 = Specie.from_string(row[0])
                s2 = Specie.from_string(row[1])
                self.species.add(s1)
                self.species.add(s2)
                self._l[frozenset([s1, s2])] = float(row[2])

        # create Z and px
        self.Z = 0
        self._px = defaultdict(float)
        for s1, s2 in itertools.product(self.species, repeat=2):
            value = math.exp(self.get_lambda(s1, s2))
            self._px[s1] += value / 2
            self._px[s2] += value / 2
            self.Z += value

 def as_dict(self):
     """
     JSON serialization of a SubStructureSpecie
     """
     sp = Specie(self.specie.symbol, self.specie.oxi_state, self.specie._properties)
     return {"@module": self.__class__.__module__,
             "@class": self.__class__.__name__,
             "specie": sp.as_dict(),
             "weight": self.weight}

 def test_to_from_string(self):
     fe3 = Specie("Fe", 3, {"spin": 5})
     self.assertEqual(str(fe3), "Fe3+spin=5")
     fe = Specie.from_string("Fe3+spin=5")
     self.assertEqual(fe.spin, 5)
     mo0 = Specie("Mo", 0, {"spin": 5})
     self.assertEqual(str(mo0), "Mo0+spin=5")
     mo = Specie.from_string("Mo0+spin=4")
     self.assertEqual(mo.spin, 4)

    def from_dict(cls, d, lattice=None):
        """
        Create PeriodicSite from dict representation.

        Args:
            d (dict): dict representation of PeriodicSite
            lattice: Optional lattice to override lattice specified in d.
                Useful for ensuring all sites in a structure share the same
                lattice.

        Returns:
            PeriodicSite
        """
        atoms_n_occu = {}
        for sp_occu in d["species"]:
            if "oxidation_state" in sp_occu and Element.is_valid_symbol(
                    sp_occu["element"]):
                sp = Specie.from_dict(sp_occu)
            elif "oxidation_state" in sp_occu:
                sp = DummySpecie.from_dict(sp_occu)
            else:
                sp = Element(sp_occu["element"])
            atoms_n_occu[sp] = sp_occu["occu"]
        props = d.get("properties", None)
        lattice = lattice if lattice else Lattice.from_dict(d["lattice"])
        return cls(atoms_n_occu, d["abc"], lattice, properties=props)

 def from_dict(d):
     """
     Create Site from dict representation
     """
     atoms_n_occu = {}
     for sp_occu in d["species"]:
         sp = Specie.from_dict(sp_occu) if "oxidation_state" in sp_occu \
             else Element(sp_occu["element"])
         atoms_n_occu[sp] = sp_occu["occu"]
     props = d.get("properties", None)
     return Site(atoms_n_occu, d["xyz"], properties=props)

    def __init__(self, lambda_table=None, alpha=-5):
        #store the input table for the to_dict method
        self._lambda_table = lambda_table

        if not lambda_table:
            module_dir = os.path.dirname(__file__)
            json_file = os.path.join(module_dir, 'data', 'lambda.json')
            with open(json_file) as f:
                lambda_table = json.load(f)

        #build map of specie pairs to lambdas
        l = {}
        for row in lambda_table:
            if not row[0] == 'D1+' and not row[1] == 'D1+':
                s1 = Specie.from_string(row[0])
                s2 = Specie.from_string(row[1])
                l[frozenset([s1, s2])] = float(row[2])

        self._lambda = l
        self._alpha = alpha

        #create the partition functions Z and px
        sp_set = set()
        for key in self._lambda.keys():
            sp_set.update(key)
        px = dict.fromkeys(sp_set, 0.)
        Z = 0
        for s1, s2 in itertools.product(sp_set, repeat=2):
            value = math.exp(self._lambda.get(frozenset([s1, s2]),
                                              self._alpha))
            #not sure why the factor of 2 is here but it matches up
            #with BURP. BURP may actually be missing a factor of 2,
            #but it doesn't have a huge effect
            px[s1] += value / 2
            px[s2] += value / 2
            Z += value

        self._Z = Z
        self._px = px
        self.species_list = list(sp_set)

def find_connected_atoms(struct, tolerance=0.45, ldict=JmolNN().el_radius):
    """
    Finds the list of bonded atoms.

    Author: "Gowoon Cheon"
    Email: "[email protected]"

    Args:
        struct (Structure): Input structure
        tolerance: length in angstroms used in finding bonded atoms. Two atoms
            are considered bonded if (radius of atom 1) + (radius of atom 2) +
            (tolerance) < (distance between atoms 1 and 2). Default
            value = 0.45, the value used by JMol and Cheon et al.
        ldict: dictionary of bond lengths used in finding bonded atoms. Values
            from JMol are used as default

    Returns:
        (np.ndarray): A numpy array of shape (number of bonded pairs, 2); each
        row of is of the form [atomi, atomj]. atomi and atomj are the indices of
        the atoms in the input structure. If any image of atomj is bonded to
        atomi with periodic boundary conditions, [atomi, atomj] is included in
        the list. If atomi is bonded to multiple images of atomj, it is only
        counted once.
    """
    n_atoms = len(struct.species)
    fc = np.array(struct.frac_coords)
    species = list(map(str, struct.species))

    # in case of charged species
    for i, item in enumerate(species):
        if item not in ldict.keys():
            species[i] = str(Specie.from_string(item).element)
    latmat = struct.lattice.matrix
    connected_list = []

    for i in range(n_atoms):
        for j in range(i + 1, n_atoms):
            max_bond_length = ldict[species[i]] + ldict[species[j]] + tolerance
            add_ij = False
            for move_cell in itertools.product(
                    [0, 1, -1], [0, 1, -1], [0, 1, -1]):
                if not add_ij:
                    frac_diff = fc[j] + move_cell - fc[i]
                    distance_ij = np.dot(latmat.T, frac_diff)
                    if np.linalg.norm(distance_ij) < max_bond_length:
                        add_ij = True
            if add_ij:
                connected_list.append([i, j])
    return np.array(connected_list)

    def test_get_shannon_radius(self):
        self.assertEqual(Specie("Li", 1).get_shannon_radius("IV"), 0.59)
        mn2 = Specie("Mn", 2)
        self.assertEqual(mn2.get_shannon_radius("IV", "High Spin"), 0.66)
        self.assertEqual(mn2.get_shannon_radius("V", "High Spin"), 0.75)

        with warnings.catch_warnings(record=True) as w:
            warnings.simplefilter("always")
            # Trigger a warning.
            r = mn2.get_shannon_radius("V")
            # Verify some things
            self.assertEqual(len(w), 1)
            self.assertIs(w[-1].category, UserWarning)
            self.assertEqual(r, 0.75)

        self.assertEqual(mn2.get_shannon_radius("VI", "Low Spin"), 0.67)
        self.assertEqual(mn2.get_shannon_radius("VI", "High Spin"), 0.83)
        self.assertEqual(mn2.get_shannon_radius("VII", "High Spin"), 0.9)
        self.assertEqual(mn2.get_shannon_radius("VIII"), 0.96)

def find_connected_atoms(struct, tolerance=0.45, ldict=JmolNN().el_radius):
    """
    Finds bonded atoms and returns a adjacency matrix of bonded atoms.

    Author: "Gowoon Cheon"
    Email: "[email protected]"

    Args:
        struct (Structure): Input structure
        tolerance: length in angstroms used in finding bonded atoms. Two atoms
            are considered bonded if (radius of atom 1) + (radius of atom 2) +
            (tolerance) < (distance between atoms 1 and 2). Default
            value = 0.45, the value used by JMol and Cheon et al.
        ldict: dictionary of bond lengths used in finding bonded atoms. Values
            from JMol are used as default

    Returns:
        (np.ndarray): A numpy array of shape (number of atoms, number of atoms);
        If any image of atom j is bonded to atom i with periodic boundary
        conditions, the matrix element [atom i, atom j] is 1.
    """
    n_atoms = len(struct.species)
    fc = np.array(struct.frac_coords)
    fc_copy = np.repeat(fc[:, :, np.newaxis], 27, axis=2)
    neighbors = np.array(list(itertools.product([0, 1, -1], [0, 1, -1], [0, 1, -1]))).T
    neighbors = np.repeat(neighbors[np.newaxis, :, :], 1, axis=0)
    fc_diff = fc_copy - neighbors
    species = list(map(str, struct.species))
    # in case of charged species
    for i, item in enumerate(species):
        if not item in ldict.keys():
            species[i] = str(Specie.from_string(item).element)
    latmat = struct.lattice.matrix
    connected_matrix = np.zeros((n_atoms,n_atoms))

    for i in range(n_atoms):
        for j in range(i + 1, n_atoms):
            max_bond_length = ldict[species[i]] + ldict[species[j]] + tolerance
            frac_diff = fc_diff[j] - fc_copy[i]
            distance_ij = np.dot(latmat.T, frac_diff)
            # print(np.linalg.norm(distance_ij,axis=0))
            if sum(np.linalg.norm(distance_ij, axis=0) < max_bond_length) > 0:
                connected_matrix[i, j] = 1
                connected_matrix[j, i] = 1
    return connected_matrix

    def coupling_constant(self, specie):
        """
        Computes the couplling constant C_q as defined in:
            Wasylishen R E, Ashbrook S E, Wimperis S. NMR of quadrupolar nuclei
            in solid materials[M]. John Wiley & Sons, 2012. (Chapter 3.2)

        C_q for a specific atom type for this electric field tensor:
                C_q=e*Q*V_zz/h
            h: planck's constant
            Q: nuclear electric quadrupole moment in mb (millibarn
            e: elementary proton charge

        Args:
            specie: flexible input to specify the species at this site.
                    Can take a isotope or element string, Specie object,
                    or Site object

        Return:

            the coupling constant as a FloatWithUnit in MHz
        """
        planks_constant=FloatWithUnit(6.62607004E-34, "m^2 kg s^-1")
        Vzz=FloatWithUnit(self.V_zz, "V ang^-2")
        e=FloatWithUnit(-1.60217662E-19, "C")

        # Convert from string to Specie object
        if isinstance(specie, str):
            # isotope was provided in string format
            if len(specie.split("-")) > 1:
                isotope=str(specie)
                specie=Specie(specie.split("-")[0])
                Q=specie.get_nmr_quadrupole_moment(isotope)
            else:
                specie=Specie(specie)
                Q=specie.get_nmr_quadrupole_moment()
        elif isinstance(specie, Site):
            specie=specie.specie
            Q=specie.get_nmr_quadrupole_moment()
        elif isinstance(specie, Specie):
            Q=specie.get_nmr_quadrupole_moment()
        else:
            raise ValueError("Invalid speciie provided for quadrupolar coupling constant calcuations")

        return (e * Q * Vzz / planks_constant).to("MHz")

 def setUp(self):
     self.specie1 = Specie.from_string("Fe2+")
     self.specie2 = Specie("Fe", 3)
     self.specie3 = Specie("Fe", 2)
     self.specie4 = Specie("Fe", 2, {"spin": 5})

    def _get_oxid_state_guesses(self, all_oxi_states, max_sites, oxi_states_override, target_charge):
        """
        Utility operation for guessing oxidation states.

        See `oxi_state_guesses` for full details. This operation does the calculation
        of the most likely oxidation states

        Args:
            oxi_states_override (dict): dict of str->list to override an
                element's common oxidation states, e.g. {"V": [2,3,4,5]}
            target_charge (int): the desired total charge on the structure.
                Default is 0 signifying charge balance.
            all_oxi_states (bool): if True, an element defaults to
                all oxidation states in pymatgen Element.icsd_oxidation_states.
                Otherwise, default is Element.common_oxidation_states. Note
                that the full oxidation state list is *very* inclusive and
                can produce nonsensical results.
            max_sites (int): if possible, will reduce Compositions to at most
                this many many sites to speed up oxidation state guesses. Set
                to -1 to just reduce fully.
        Returns:
            A list of dicts - each dict reports an element symbol and average
                oxidation state across all sites in that composition. If the
                composition is not charge balanced, an empty list is returned.
            A list of dicts - each dict maps the element symbol to a list of
                oxidation states for each site of that element. For example, Fe3O4 could
                return a list of [2,2,2,3,3,3] for the oxidation states of If the composition
                is

            """
        comp = self.copy()
        # reduce Composition if necessary
        if max_sites == -1:
            comp = self.reduced_composition

        elif max_sites and comp.num_atoms > max_sites:
            reduced_comp, reduced_factor = self. \
                get_reduced_composition_and_factor()
            if reduced_factor > 1:
                reduced_comp *= max(1, int(max_sites / reduced_comp.num_atoms))
                comp = reduced_comp  # as close to max_sites as possible
            if comp.num_atoms > max_sites:
                raise ValueError("Composition {} cannot accommodate max_sites "
                                 "setting!".format(comp))

        # Load prior probabilities of oxidation states, used to rank solutions
        if not Composition.oxi_prob:
            module_dir = os.path.join(os.path.
                                      dirname(os.path.abspath(__file__)))
            all_data = loadfn(os.path.join(module_dir, "..",
                                           "analysis", "icsd_bv.yaml"))
            Composition.oxi_prob = {Specie.from_string(sp): data
                                    for sp, data in
                                    all_data["occurrence"].items()}
        oxi_states_override = oxi_states_override or {}
        # assert: Composition only has integer amounts
        if not all(amt == int(amt) for amt in comp.values()):
            raise ValueError("Charge balance analysis requires integer "
                             "values in Composition!")

        # for each element, determine all possible sum of oxidations
        # (taking into account nsites for that particular element)
        el_amt = comp.get_el_amt_dict()
        els = el_amt.keys()
        el_sums = []  # matrix: dim1= el_idx, dim2=possible sums
        el_sum_scores = defaultdict(set)  # dict of el_idx, sum -> score
        el_best_oxid_combo = {}  # dict of el_idx, sum -> oxid combo with best score
        for idx, el in enumerate(els):
            el_sum_scores[idx] = {}
            el_best_oxid_combo[idx] = {}
            el_sums.append([])
            if oxi_states_override.get(el):
                oxids = oxi_states_override[el]
            elif all_oxi_states:
                oxids = Element(el).oxidation_states
            else:
                oxids = Element(el).icsd_oxidation_states or \
                        Element(el).oxidation_states

            # get all possible combinations of oxidation states
            # and sum each combination
            for oxid_combo in combinations_with_replacement(oxids,
                                                            int(el_amt[el])):

                # List this sum as a possible option
                oxid_sum = sum(oxid_combo)
                if oxid_sum not in el_sums[idx]:
                    el_sums[idx].append(oxid_sum)

                # Determine how probable is this combo?
                score = sum([Composition.oxi_prob.get(Specie(el, o), 0) for
                             o in oxid_combo])

                # If it is the most probable combo for a certain sum,
                #   store the combination
                if oxid_sum not in el_sum_scores[idx] or score > el_sum_scores[idx].get(oxid_sum, 0):
                    el_sum_scores[idx][oxid_sum] = score
                    el_best_oxid_combo[idx][oxid_sum] = oxid_combo

        # Determine which combination of oxidation states for each element
#.........这里部分代码省略.........

 def from_dict(cls, d):
     return SubStructureSpecie.from_specie(Specie.from_dict(d['specie']), d.get('weight', 0.0))

    def oxi_state_guesses(self, oxi_states_override=None, target_charge=0,
                          all_oxi_states=False, max_sites=None):
        """
        Checks if the composition is charge-balanced and returns back all
        charge-balanced oxidation state combinations. Composition must have
        integer values. Note that more num_atoms in the composition gives
        more degrees of freedom. e.g., if possible oxidation states of
        element X are [2,4] and Y are [-3], then XY is not charge balanced
        but X2Y2 is. Results are returned from most to least probable based
        on ICSD statistics. Use max_sites to improve performance if needed.

        Args:
            oxi_states_override (dict): dict of str->list to override an
                element's common oxidation states, e.g. {"V": [2,3,4,5]}
            target_charge (int): the desired total charge on the structure.
                Default is 0 signifying charge balance.
            all_oxi_states (bool): if True, an element defaults to
                all oxidation states in pymatgen Element.icsd_oxidation_states.
                Otherwise, default is Element.common_oxidation_states. Note
                that the full oxidation state list is *very* inclusive and
                can produce nonsensical results.
            max_sites (int): if possible, will reduce Compositions to at most
                this many many sites to speed up oxidation state guesses. Set
                to -1 to just reduce fully.

        Returns:
            A list of dicts - each dict reports an element symbol and average
                oxidation state across all sites in that composition. If the
                composition is not charge balanced, an empty list is returned.
        """

        comp = self.copy()

        # reduce Composition if necessary
        if max_sites == -1:
            comp = self.reduced_composition

        elif max_sites and comp.num_atoms > max_sites:
            reduced_comp, reduced_factor = self.\
                get_reduced_composition_and_factor()
            if reduced_factor > 1:
                reduced_comp *= max(1, int(max_sites / reduced_comp.num_atoms))
                comp = reduced_comp  # as close to max_sites as possible
            if comp.num_atoms > max_sites:
                raise ValueError("Composition {} cannot accommodate max_sites "
                                 "setting!".format(comp))

        # Load prior probabilities of oxidation states, used to rank solutions
        if not Composition.oxi_prob:
            module_dir = os.path.join(os.path.
                                      dirname(os.path.abspath(__file__)))
            all_data = loadfn(os.path.join(module_dir, "..",
                                           "analysis", "icsd_bv.yaml"))
            Composition.oxi_prob = {Specie.from_string(sp): data
                                    for sp, data in
                                    all_data["occurrence"].items()}

        oxi_states_override = oxi_states_override or {}

        # assert: Composition only has integer amounts
        if not all(amt == int(amt) for amt in comp.values()):
            raise ValueError("Charge balance analysis requires integer "
                             "values in Composition!")

        # for each element, determine all possible sum of oxidations
        # (taking into account nsites for that particular element)
        el_amt = comp.get_el_amt_dict()
        els = el_amt.keys()
        el_sums = []  # matrix: dim1= el_idx, dim2=possible sums
        el_sum_scores = defaultdict(set)  # dict of el_idx, sum -> score
        for idx, el in enumerate(els):
            el_sum_scores[idx] = {}
            el_sums.append([])
            if oxi_states_override.get(el):
                oxids = oxi_states_override[el]
            elif all_oxi_states:
                oxids = Element(el).oxidation_states
            else:
                oxids = Element(el).icsd_oxidation_states or \
                        Element(el).oxidation_states

            # get all possible combinations of oxidation states
            # and sum each combination
            for oxid_combo in combinations_with_replacement(oxids,
                                                            int(el_amt[el])):
                if sum(oxid_combo) not in el_sums[idx]:
                    el_sums[idx].append(sum(oxid_combo))
                    score = sum([Composition.oxi_prob.get(Specie(el, o), 0) for
                                 o in oxid_combo])  # how probable is this combo?
                    el_sum_scores[idx][sum(oxid_combo)] = max(
                        el_sum_scores[idx].get(sum(oxid_combo), 0), score)


        all_sols = []  # will contain all solutions
        all_scores = []  # will contain a score for each solution
        for x in product(*el_sums):
            # each x is a trial of one possible oxidation sum for each element
            if sum(x) == target_charge:  # charge balance condition
                el_sum_sol = dict(zip(els, x))  # element->oxid_sum
                # normalize oxid_sum by amount to get avg oxid state
#.........这里部分代码省略.........

                           "N", "P", "As", "Sb",
                           "O", "S", "Se", "Te",
                           "F", "Cl", "Br", "I"])

module_dir = os.path.dirname(os.path.abspath(__file__))

#Read in BV parameters.
BV_PARAMS = {}
with open(os.path.join(module_dir, "bvparam_1991.json"), "r") as f:
    for k, v in json.load(f).items():
        BV_PARAMS[Element(k)] = v

#Read in json containing data-mined ICSD BV data.
with open(os.path.join(module_dir, "icsd_bv.json"), "r") as f:
    all_data = json.load(f)
    ICSD_BV_DATA = {Specie.from_string(sp): data
                    for sp, data in all_data["bvsum"].items()}
    PRIOR_PROB = {Specie.from_string(sp): data
                  for sp, data in all_data["occurrence"].items()}


def calculate_bv_sum(site, nn_list, scale_factor=1.0):
    """
    Calculates the BV sum of a site.

    Args:
        site:
            The site
        nn_list:
            List of nearest neighbors in the format [(nn_site, dist), ...].
        anion_el: