4075885785

Build Type

push

github

Committed by Shyue Ping Ong

Commit Message

Merge branch 'master' of github.com:materialsproject/pymatgen

Run Details

96 of 96 new or added lines in 27 files covered. (100.0%)

81013 of 102710 relevant lines covered (78.88%)

0.79 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

55.72

/pymatgen/analysis/bond_valence.py

# Copyright (c) Pymatgen Development Team.
# Distributed under the terms of the MIT License.

"""
This module implements classes to perform bond valence analyses.
"""

from __future__ import annotations

import collections
import functools
import operator
import os
from math import exp, sqrt

import numpy as np
from monty.serialization import loadfn

from pymatgen.core.periodic_table import Element, Species, get_el_sp
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer

# Let's initialize some module level properties.

# List of electronegative elements specified in M. O'Keefe, & N. Brese,
# JACS, 1991, 113(9), 3226-3229. doi:10.1021/ja00009a002.
ELECTRONEG = [
    Element(sym)
    for sym in [
        "H",
        "B",
        "C",
        "Si",
        "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 = {}
for k, v in loadfn(os.path.join(module_dir, "bvparam_1991.yaml")).items():
    BV_PARAMS[Element(k)] = v

# Read in yaml containing data-mined ICSD BV data.
all_data = loadfn(os.path.join(module_dir, "icsd_bv.yaml"))
ICSD_BV_DATA = {Species.from_string(sp): data for sp, data in all_data["bvsum"].items()}
PRIOR_PROB = {Species.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 (PeriodicSite): The central site to calculate the bond valence
        nn_list ([Neighbor]): A list of namedtuple Neighbors having "distance"
            and "site" attributes
        scale_factor (float): A scale factor to be applied. This is useful for
            scaling distance, esp in the case of calculation-relaxed structures
            which may tend to under (GGA) or over bind (LDA).
    """
    el1 = Element(site.specie.symbol)
    bvsum = 0
    for nn in nn_list:
        el2 = Element(nn.specie.symbol)
        if (el1 in ELECTRONEG or el2 in ELECTRONEG) and el1 != el2:
            r1 = BV_PARAMS[el1]["r"]
            r2 = BV_PARAMS[el2]["r"]
            c1 = BV_PARAMS[el1]["c"]
            c2 = BV_PARAMS[el2]["c"]
            R = r1 + r2 - r1 * r2 * (sqrt(c1) - sqrt(c2)) ** 2 / (c1 * r1 + c2 * r2)
            vij = exp((R - nn.nn_distance * scale_factor) / 0.31)
            bvsum += vij * (1 if el1.X < el2.X else -1)
    return bvsum


def calculate_bv_sum_unordered(site, nn_list, scale_factor=1):
    """
    Calculates the BV sum of a site for unordered structures.

    Args:
        site (PeriodicSite): The central site to calculate the bond valence
        nn_list ([Neighbor]): A list of namedtuple Neighbors having "distance"
            and "site" attributes
        scale_factor (float): A scale factor to be applied. This is useful for
            scaling distance, esp in the case of calculation-relaxed structures
            which may tend to under (GGA) or over bind (LDA).
    """
    # If the site "site" has N partial occupations as : f_{site}_0,
    # f_{site}_1, ... f_{site}_N of elements
    # X_{site}_0, X_{site}_1, ... X_{site}_N, and each neighbors nn_i in nn
    # has N_{nn_i} partial occupations as :
    # f_{nn_i}_0, f_{nn_i}_1, ..., f_{nn_i}_{N_{nn_i}}, then the bv sum of
    # site "site" is obtained as :
    # \sum_{nn} \sum_j^N \sum_k^{N_{nn}} f_{site}_j f_{nn_i}_k vij_full
    # where vij_full is the valence bond of the fully occupied bond
    bvsum = 0
    for specie1, occu1 in site.species.items():
        el1 = Element(specie1.symbol)
        for nn in nn_list:
            for specie2, occu2 in nn.species.items():
                el2 = Element(specie2.symbol)
                if (el1 in ELECTRONEG or el2 in ELECTRONEG) and el1 != el2:
                    r1 = BV_PARAMS[el1]["r"]
                    r2 = BV_PARAMS[el2]["r"]
                    c1 = BV_PARAMS[el1]["c"]
                    c2 = BV_PARAMS[el2]["c"]
                    R = r1 + r2 - r1 * r2 * (sqrt(c1) - sqrt(c2)) ** 2 / (c1 * r1 + c2 * r2)
                    vij = exp((R - nn.nn_distance * scale_factor) / 0.31)
                    bvsum += occu1 * occu2 * vij * (1 if el1.X < el2.X else -1)
    return bvsum


class BVAnalyzer:
    """
    This class implements a maximum a posteriori (MAP) estimation method to
    determine oxidation states in a structure. The algorithm is as follows:
    1) The bond valence sum of all symmetrically distinct sites in a structure
    is calculated using the element-based parameters in M. O'Keefe, & N. Brese,
    JACS, 1991, 113(9), 3226-3229. doi:10.1021/ja00009a002.
    2) The posterior probabilities of all oxidation states is then calculated
    using: P(oxi_state/BV) = K * P(BV/oxi_state) * P(oxi_state), where K is
    a constant factor for each element. P(BV/oxi_state) is calculated as a
    Gaussian with mean and std deviation determined from an analysis of
    the ICSD. The posterior P(oxi_state) is determined from a frequency
    analysis of the ICSD.
    3) The oxidation states are then ranked in order of decreasing probability
    and the oxidation state combination that result in a charge neutral cell
    is selected.
    """

    CHARGE_NEUTRALITY_TOLERANCE = 0.00001

    def __init__(
        self,
        symm_tol=0.1,
        max_radius=4,
        max_permutations=100000,
        distance_scale_factor=1.015,
        charge_neutrality_tolerance=CHARGE_NEUTRALITY_TOLERANCE,
        forbidden_species=None,
    ):
        """
        Initializes the BV analyzer, with useful defaults.

        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.
            charge_neutrality_tolerance:
                Tolerance on the charge neutrality when unordered structures
                are at stake.
            forbidden_species:
                List of species that are forbidden (example : ["O-"] cannot be
                used) It is used when e.g. someone knows that some oxidation
                state cannot occur for some atom in a structure or list of
                structures.
        """
        self.symm_tol = symm_tol
        self.max_radius = max_radius
        self.max_permutations = max_permutations
        self.dist_scale_factor = distance_scale_factor
        self.charge_neutrality_tolerance = charge_neutrality_tolerance
        forbidden_species = [get_el_sp(sp) for sp in forbidden_species] if forbidden_species else []
        self.icsd_bv_data = (
            {get_el_sp(specie): data for specie, data in ICSD_BV_DATA.items() if specie not in forbidden_species}
            if len(forbidden_species) > 0
            else ICSD_BV_DATA
        )

    def _calc_site_probabilities(self, site, nn):
        el = site.specie.symbol
        bv_sum = calculate_bv_sum(site, nn, scale_factor=self.dist_scale_factor)
        prob = {}
        for sp, data in self.icsd_bv_data.items():
            if sp.symbol == el and sp.oxi_state != 0 and data["std"] > 0:
                u = data["mean"]
                sigma = data["std"]
                # Calculate posterior probability. Note that constant
                # factors are ignored. They have no effect on the results.
                prob[sp.oxi_state] = exp(-((bv_sum - u) ** 2) / 2 / (sigma**2)) / sigma * PRIOR_PROB[sp]
        # Normalize the probabilities
        try:
            prob = {k: v / sum(prob.values()) for k, v in prob.items()}
        except ZeroDivisionError:
            prob = {k: 0.0 for k in prob}
        return prob

    def _calc_site_probabilities_unordered(self, site, nn):
        bv_sum = calculate_bv_sum_unordered(site, nn, scale_factor=self.dist_scale_factor)
        prob = {}
        for specie in site.species:
            el = specie.symbol

            prob[el] = {}
            for sp, data in self.icsd_bv_data.items():
                if sp.symbol == el and sp.oxi_state != 0 and data["std"] > 0:
                    u = data["mean"]
                    sigma = data["std"]
                    # Calculate posterior probability. Note that constant
                    # factors are ignored. They have no effect on the results.
                    prob[el][sp.oxi_state] = exp(-((bv_sum - u) ** 2) / 2 / (sigma**2)) / sigma * PRIOR_PROB[sp]
            # Normalize the probabilities
            try:
                prob[el] = {k: v / sum(prob[el].values()) for k, v in prob[el].items()}
            except ZeroDivisionError:
                prob[el] = {k: 0.0 for k in prob[el]}
        return prob

    def get_valences(self, structure):
        """
        Returns a list of valences for the structure. This currently works only
        for ordered structures only.

        Args:
            structure: Structure to analyze

        Returns:
            A list of valences for each site in the structure (for an ordered
            structure), e.g., [1, 1, -2] or a list of lists with the
            valences for each fractional element of each site in the
            structure (for an unordered structure),
            e.g., [[2, 4], [3], [-2], [-2], [-2]]

        Raises:
            A ValueError if the valences cannot be determined.
        """
        els = [Element(el.symbol) for el in structure.composition.elements]

        if not set(els).issubset(set(BV_PARAMS)):
            raise ValueError("Structure contains elements not in set of BV parameters!")

        # Perform symmetry determination and get sites grouped by symmetry.
        if self.symm_tol:
            finder = SpacegroupAnalyzer(structure, self.symm_tol)
            symm_structure = finder.get_symmetrized_structure()
            equi_sites = symm_structure.equivalent_sites
        else:
            equi_sites = [[site] for site in structure]

        # Sort the equivalent sites by decreasing electronegativity.
        equi_sites = sorted(equi_sites, key=lambda sites: -sites[0].species.average_electroneg)

        # Get a list of valences and probabilities for each symmetrically
        # distinct site.
        valences = []
        all_prob = []
        if structure.is_ordered:
            for sites in equi_sites:
                test_site = sites[0]
                nn = structure.get_neighbors(test_site, self.max_radius)
                prob = self._calc_site_probabilities(test_site, nn)
                all_prob.append(prob)
                val = list(prob)
                # Sort valences in order of decreasing probability.
                val = sorted(val, key=lambda v: -prob[v])
                # Retain probabilities that are at least 1/100 of highest prob.
                valences.append(list(filter(lambda v: prob[v] > 0.01 * prob[val[0]], val)))
        else:
            full_all_prob = []
            for sites in equi_sites:
                test_site = sites[0]
                nn = structure.get_neighbors(test_site, self.max_radius)
                prob = self._calc_site_probabilities_unordered(test_site, nn)
                all_prob.append(prob)
                full_all_prob.extend(prob.values())
                vals = []
                for elsp, _ in get_z_ordered_elmap(test_site.species):
                    val = list(prob[elsp.symbol])
                    # Sort valences in order of decreasing probability.
                    val = sorted(val, key=lambda v: -prob[elsp.symbol][v])
                    # Retain probabilities that are at least 1/100 of highest
                    # prob.
                    vals.append(
                        list(
                            filter(
                                lambda v: prob[elsp.symbol][v] > 0.001 * prob[elsp.symbol][val[0]],
                                val,
                            )
                        )
                    )
                valences.append(vals)

        # make variables needed for recursion
        if structure.is_ordered:
            nsites = np.array([len(i) for i in equi_sites])
            vmin = np.array([min(i) for i in valences])
            vmax = np.array([max(i) for i in valences])

            self._n = 0
            self._best_score = 0
            self._best_vset = None

            def evaluate_assignment(v_set):
                el_oxi = collections.defaultdict(list)
                for i, sites in enumerate(equi_sites):
                    el_oxi[sites[0].specie.symbol].append(v_set[i])
                max_diff = max(max(v) - min(v) for v in el_oxi.values())
                if max_diff > 1:
                    return
                score = functools.reduce(operator.mul, [all_prob[i][v] for i, v in enumerate(v_set)])
                if score > self._best_score:
                    self._best_vset = v_set
                    self._best_score = score

            def _recurse(assigned=None):
                # recurses to find permutations of valences based on whether a
                # charge balanced assignment can still be found
                if self._n > self.max_permutations:
                    return None
                if assigned is None:
                    assigned = []

                i = len(assigned)
                highest = vmax.copy()
                highest[:i] = assigned
                highest *= nsites
                highest = np.sum(highest)

                lowest = vmin.copy()
                lowest[:i] = assigned
                lowest *= nsites
                lowest = np.sum(lowest)

                if highest < 0 or lowest > 0:
                    self._n += 1
                    return None

                if i == len(valences):
                    evaluate_assignment(assigned)
                    self._n += 1
                    return None
                for v in valences[i]:
                    new_assigned = list(assigned)
                    _recurse(new_assigned + [v])
                return None

        else:
            nsites = np.array([len(i) for i in equi_sites])
            tmp = []
            attrib = []
            for insite, nsite in enumerate(nsites):
                for _ in valences[insite]:
                    tmp.append(nsite)
                    attrib.append(insite)
            new_nsites = np.array(tmp)
            fractions = []
            elements = []
            for sites in equi_sites:
                for sp, occu in get_z_ordered_elmap(sites[0].species):
                    elements.append(sp.symbol)
                    fractions.append(occu)
            fractions = np.array(fractions, np.float_)
            new_valences = []
            for vals in valences:
                for val in vals:
                    new_valences.append(val)
            vmin = np.array([min(i) for i in new_valences], np.float_)
            vmax = np.array([max(i) for i in new_valences], np.float_)

            self._n = 0
            self._best_score = 0
            self._best_vset = None

            def evaluate_assignment(v_set):
                el_oxi = collections.defaultdict(list)
                jj = 0
                for sites in equi_sites:
                    for specie, _ in get_z_ordered_elmap(sites[0].species):
                        el_oxi[specie.symbol].append(v_set[jj])
                        jj += 1
                max_diff = max(max(v) - min(v) for v in el_oxi.values())
                if max_diff > 2:
                    return

                score = functools.reduce(
                    operator.mul,
                    [all_prob[attrib[iv]][elements[iv]][vv] for iv, vv in enumerate(v_set)],
                )
                if score > self._best_score:
                    self._best_vset = v_set
                    self._best_score = score

            def _recurse(assigned=None):
                # recurses to find permutations of valences based on whether a
                # charge balanced assignment can still be found
                if self._n > self.max_permutations:
                    return None
                if assigned is None:
                    assigned = []

                i = len(assigned)
                highest = vmax.copy()
                highest[:i] = assigned
                highest *= new_nsites
                highest *= fractions
                highest = np.sum(highest)

                lowest = vmin.copy()
                lowest[:i] = assigned
                lowest *= new_nsites
                lowest *= fractions
                lowest = np.sum(lowest)

                if highest < -self.charge_neutrality_tolerance or lowest > self.charge_neutrality_tolerance:
                    self._n += 1
                    return None

                if i == len(new_valences):
                    evaluate_assignment(assigned)
                    self._n += 1
                    return None

                for v in new_valences[i]:
                    new_assigned = list(assigned)
                    _recurse(new_assigned + [v])

                return None

        _recurse()

        if self._best_vset:
            if structure.is_ordered:
                assigned = {}
                for val, sites in zip(self._best_vset, equi_sites):
                    for site in sites:
                        assigned[site] = val

                return [int(assigned[site]) for site in structure]
            assigned = {}
            new_best_vset = []
            for _ in equi_sites:
                new_best_vset.append([])
            for ival, val in enumerate(self._best_vset):
                new_best_vset[attrib[ival]].append(val)
            for val, sites in zip(new_best_vset, equi_sites):
                for site in sites:
                    assigned[site] = val

            return [[int(frac_site) for frac_site in assigned[site]] for site in structure]
        raise ValueError("Valences cannot be assigned!")

    def get_oxi_state_decorated_structure(self, structure):
        """
        Get an oxidation state decorated structure. This currently works only
        for ordered structures only.

        Args:
            structure: Structure to analyze

        Returns:
            A modified structure that is oxidation state decorated.

        Raises:
            ValueError if the valences cannot be determined.
        """
        s = structure.copy()
        if s.is_ordered:
            valences = self.get_valences(s)
            s.add_oxidation_state_by_site(valences)
        else:
            valences = self.get_valences(s)
            s = add_oxidation_state_by_site_fraction(s, valences)
        return s


def get_z_ordered_elmap(comp):
    """
    Arbitrary ordered elmap on the elements/species of a composition of a
    given site in an unordered structure. Returns a list of tuples (
    element_or_specie: occupation) in the arbitrary order.

    The arbitrary order is based on the Z of the element and the smallest
    fractional occupations first.
    Example : {"Ni3+": 0.2, "Ni4+": 0.2, "Cr3+": 0.15, "Zn2+": 0.34,
    "Cr4+": 0.11} will yield the species in the following order :
    Cr4+, Cr3+, Ni3+, Ni4+, Zn2+ ... or
    Cr4+, Cr3+, Ni4+, Ni3+, Zn2+
    """
    return sorted((elsp, comp[elsp]) for elsp in comp)


def add_oxidation_state_by_site_fraction(structure, oxidation_states):
    """
    Add oxidation states to a structure by fractional site.

    Args:
        oxidation_states (list): List of list of oxidation states for each
            site fraction for each site.
            E.g., [[2, 4], [3], [-2], [-2], [-2]]
    """
    try:
        for i, site in enumerate(structure):
            new_sp = collections.defaultdict(float)
            for j, (el, occu) in enumerate(get_z_ordered_elmap(site.species)):
                specie = Species(el.symbol, oxidation_states[i][j])
                new_sp[specie] += occu
            structure[i] = new_sp
        return structure
    except IndexError:
        raise ValueError("Oxidation state of all sites must be specified in the list.")

1	# Copyright (c) Pymatgen Development Team.
2	# Distributed under the terms of the MIT License.
3
4	"""	1✔
5	This module implements classes to perform bond valence analyses.
6	"""
7
8	from __future__ import annotations	1✔
9
10	import collections	1✔
11	import functools	1✔
12	import operator	1✔
13	import os	1✔
14	from math import exp, sqrt	1✔
15
16	import numpy as np	1✔
17	from monty.serialization import loadfn	1✔
18
19	from pymatgen.core.periodic_table import Element, Species, get_el_sp	1✔
20	from pymatgen.symmetry.analyzer import SpacegroupAnalyzer	1✔
21
22	# Let's initialize some module level properties.
23
24	# List of electronegative elements specified in M. O'Keefe, & N. Brese,
25	# JACS, 1991, 113(9), 3226-3229. doi:10.1021/ja00009a002.
26	ELECTRONEG = [	1✔
27	Element(sym)
28	for sym in [
29	"H",
30	"B",
31	"C",
32	"Si",
33	"N",
34	"P",
35	"As",
36	"Sb",
37	"O",
38	"S",
39	"Se",
40	"Te",
41	"F",
42	"Cl",
43	"Br",
44	"I",
45	]
46	]
47
48	module_dir = os.path.dirname(os.path.abspath(__file__))	1✔
49
50	# Read in BV parameters.
51	BV_PARAMS = {}	1✔
52	for k, v in loadfn(os.path.join(module_dir, "bvparam_1991.yaml")).items():	1✔
53	BV_PARAMS[Element(k)] = v	1✔
54
55	# Read in yaml containing data-mined ICSD BV data.
56	all_data = loadfn(os.path.join(module_dir, "icsd_bv.yaml"))	1✔
57	ICSD_BV_DATA = {Species.from_string(sp): data for sp, data in all_data["bvsum"].items()}	1✔
58	PRIOR_PROB = {Species.from_string(sp): data for sp, data in all_data["occurrence"].items()}	1✔
59
60
61	def calculate_bv_sum(site, nn_list, scale_factor=1.0):	1✔
62	"""
63	Calculates the BV sum of a site.
64
65	Args:
66	site (PeriodicSite): The central site to calculate the bond valence
67	nn_list ([Neighbor]): A list of namedtuple Neighbors having "distance"
68	and "site" attributes
69	scale_factor (float): A scale factor to be applied. This is useful for
70	scaling distance, esp in the case of calculation-relaxed structures
71	which may tend to under (GGA) or over bind (LDA).
72	"""
73	el1 = Element(site.specie.symbol)	1✔
74	bvsum = 0	1✔
75	for nn in nn_list:	1✔
76	el2 = Element(nn.specie.symbol)	1✔
77	if (el1 in ELECTRONEG or el2 in ELECTRONEG) and el1 != el2:	1✔
78	r1 = BV_PARAMS[el1]["r"]	1✔
79	r2 = BV_PARAMS[el2]["r"]	1✔
80	c1 = BV_PARAMS[el1]["c"]	1✔
81	c2 = BV_PARAMS[el2]["c"]	1✔
82	R = r1 + r2 - r1 * r2 * (sqrt(c1) - sqrt(c2)) ** 2 / (c1 * r1 + c2 * r2)	1✔
83	vij = exp((R - nn.nn_distance * scale_factor) / 0.31)	1✔
84	bvsum += vij * (1 if el1.X < el2.X else -1)	1✔
85	return bvsum	1✔
86
87
88	def calculate_bv_sum_unordered(site, nn_list, scale_factor=1):	1✔
89	"""
90	Calculates the BV sum of a site for unordered structures.
91
92	Args:
93	site (PeriodicSite): The central site to calculate the bond valence
94	nn_list ([Neighbor]): A list of namedtuple Neighbors having "distance"
95	and "site" attributes
96	scale_factor (float): A scale factor to be applied. This is useful for
97	scaling distance, esp in the case of calculation-relaxed structures
98	which may tend to under (GGA) or over bind (LDA).
99	"""
100	# If the site "site" has N partial occupations as : f_{site}_0,
101	# f_{site}_1, ... f_{site}_N of elements
102	# X_{site}_0, X_{site}_1, ... X_{site}_N, and each neighbors nn_i in nn
103	# has N_{nn_i} partial occupations as :
104	# f_{nn_i}_0, f_{nn_i}_1, ..., f_{nn_i}_{N_{nn_i}}, then the bv sum of
105	# site "site" is obtained as :
106	# \sum_{nn} \sum_j^N \sum_k^{N_{nn}} f_{site}_j f_{nn_i}_k vij_full
107	# where vij_full is the valence bond of the fully occupied bond
108	bvsum = 0	1✔
109	for specie1, occu1 in site.species.items():	1✔
110	el1 = Element(specie1.symbol)	1✔
111	for nn in nn_list:	1✔
112	for specie2, occu2 in nn.species.items():	1✔
113	el2 = Element(specie2.symbol)	1✔
114	if (el1 in ELECTRONEG or el2 in ELECTRONEG) and el1 != el2:	1✔
115	r1 = BV_PARAMS[el1]["r"]	1✔
116	r2 = BV_PARAMS[el2]["r"]	1✔
117	c1 = BV_PARAMS[el1]["c"]	1✔
118	c2 = BV_PARAMS[el2]["c"]	1✔
119	R = r1 + r2 - r1 * r2 * (sqrt(c1) - sqrt(c2)) ** 2 / (c1 * r1 + c2 * r2)	1✔
120	vij = exp((R - nn.nn_distance * scale_factor) / 0.31)	1✔
121	bvsum += occu1 * occu2 * vij * (1 if el1.X < el2.X else -1)	1✔
122	return bvsum	1✔
123
124
125	class BVAnalyzer:	1✔
126	"""
127	This class implements a maximum a posteriori (MAP) estimation method to
128	determine oxidation states in a structure. The algorithm is as follows:
129	1) The bond valence sum of all symmetrically distinct sites in a structure
130	is calculated using the element-based parameters in M. O'Keefe, & N. Brese,
131	JACS, 1991, 113(9), 3226-3229. doi:10.1021/ja00009a002.
132	2) The posterior probabilities of all oxidation states is then calculated
133	using: P(oxi_state/BV) = K * P(BV/oxi_state) * P(oxi_state), where K is
134	a constant factor for each element. P(BV/oxi_state) is calculated as a
135	Gaussian with mean and std deviation determined from an analysis of
136	the ICSD. The posterior P(oxi_state) is determined from a frequency
137	analysis of the ICSD.
138	3) The oxidation states are then ranked in order of decreasing probability
139	and the oxidation state combination that result in a charge neutral cell
140	is selected.
141	"""
142
143	CHARGE_NEUTRALITY_TOLERANCE = 0.00001	1✔
144
145	def __init__(	1✔
146	self,
147	symm_tol=0.1,
148	max_radius=4,
149	max_permutations=100000,
150	distance_scale_factor=1.015,
151	charge_neutrality_tolerance=CHARGE_NEUTRALITY_TOLERANCE,
152	forbidden_species=None,
153	):
154	"""
155	Initializes the BV analyzer, with useful defaults.
156
157	Args:
158	symm_tol:
159	Symmetry tolerance used to determine which sites are
160	symmetrically equivalent. Set to 0 to turn off symmetry.
161	max_radius:
162	Maximum radius in Angstrom used to find nearest neighbors.
163	max_permutations:
164	The maximum number of permutations of oxidation states to test.
165	distance_scale_factor:
166	A scale factor to be applied. This is useful for scaling
167	distances, esp in the case of calculation-relaxed structures
168	which may tend to under (GGA) or over bind (LDA). The default
169	of 1.015 works for GGA. For experimental structure, set this to
170	1.
171	charge_neutrality_tolerance:
172	Tolerance on the charge neutrality when unordered structures
173	are at stake.
174	forbidden_species:
175	List of species that are forbidden (example : ["O-"] cannot be
176	used) It is used when e.g. someone knows that some oxidation
177	state cannot occur for some atom in a structure or list of
178	structures.
179	"""
180	self.symm_tol = symm_tol	1✔
181	self.max_radius = max_radius	1✔
182	self.max_permutations = max_permutations	1✔
183	self.dist_scale_factor = distance_scale_factor	1✔
184	self.charge_neutrality_tolerance = charge_neutrality_tolerance	1✔
185	forbidden_species = [get_el_sp(sp) for sp in forbidden_species] if forbidden_species else []	1✔
186	self.icsd_bv_data = (	1✔
187	{get_el_sp(specie): data for specie, data in ICSD_BV_DATA.items() if specie not in forbidden_species}
188	if len(forbidden_species) > 0
189	else ICSD_BV_DATA
190	)
191
192	def _calc_site_probabilities(self, site, nn):	1✔
193	el = site.specie.symbol	1✔
194	bv_sum = calculate_bv_sum(site, nn, scale_factor=self.dist_scale_factor)	1✔
195	prob = {}	1✔
196	for sp, data in self.icsd_bv_data.items():	1✔
197	if sp.symbol == el and sp.oxi_state != 0 and data["std"] > 0:	1✔
198	u = data["mean"]	1✔
199	sigma = data["std"]	1✔
200	# Calculate posterior probability. Note that constant
201	# factors are ignored. They have no effect on the results.
202	prob[sp.oxi_state] = exp(-((bv_sum - u) 2) / 2 / (sigma2)) / sigma * PRIOR_PROB[sp]	1✔
203	# Normalize the probabilities
204	try:	1✔
205	prob = {k: v / sum(prob.values()) for k, v in prob.items()}	1✔
206	except ZeroDivisionError:	×
207	prob = {k: 0.0 for k in prob}	×
208	return prob	1✔
209
210	def _calc_site_probabilities_unordered(self, site, nn):	1✔
211	bv_sum = calculate_bv_sum_unordered(site, nn, scale_factor=self.dist_scale_factor)	×
212	prob = {}	×
213	for specie in site.species:	×
214	el = specie.symbol	×
215
216	prob[el] = {}	×
217	for sp, data in self.icsd_bv_data.items():	×
218	if sp.symbol == el and sp.oxi_state != 0 and data["std"] > 0:	×
219	u = data["mean"]	×
220	sigma = data["std"]	×
221	# Calculate posterior probability. Note that constant
222	# factors are ignored. They have no effect on the results.
223	prob[el][sp.oxi_state] = exp(-((bv_sum - u) 2) / 2 / (sigma2)) / sigma * PRIOR_PROB[sp]	×
224	# Normalize the probabilities
225	try:	×
226	prob[el] = {k: v / sum(prob[el].values()) for k, v in prob[el].items()}	×
227	except ZeroDivisionError:	×
228	prob[el] = {k: 0.0 for k in prob[el]}	×
229	return prob	×
230
231	def get_valences(self, structure):	1✔
232	"""
233	Returns a list of valences for the structure. This currently works only
234	for ordered structures only.
235
236	Args:
237	structure: Structure to analyze
238
239	Returns:
240	A list of valences for each site in the structure (for an ordered
241	structure), e.g., [1, 1, -2] or a list of lists with the
242	valences for each fractional element of each site in the
243	structure (for an unordered structure),
244	e.g., [[2, 4], [3], [-2], [-2], [-2]]
245
246	Raises:
247	A ValueError if the valences cannot be determined.
248	"""
249	els = [Element(el.symbol) for el in structure.composition.elements]	1✔
250
251	if not set(els).issubset(set(BV_PARAMS)):	1✔
252	raise ValueError("Structure contains elements not in set of BV parameters!")	×
253
254	# Perform symmetry determination and get sites grouped by symmetry.
255	if self.symm_tol:	1✔
256	finder = SpacegroupAnalyzer(structure, self.symm_tol)	1✔
257	symm_structure = finder.get_symmetrized_structure()	1✔
258	equi_sites = symm_structure.equivalent_sites	1✔
259	else:
260	equi_sites = [[site] for site in structure]	1✔
261
262	# Sort the equivalent sites by decreasing electronegativity.
263	equi_sites = sorted(equi_sites, key=lambda sites: -sites[0].species.average_electroneg)	1✔
264
265	# Get a list of valences and probabilities for each symmetrically
266	# distinct site.
267	valences = []	1✔
268	all_prob = []	1✔
269	if structure.is_ordered:	1✔
270	for sites in equi_sites:	1✔
271	test_site = sites[0]	1✔
272	nn = structure.get_neighbors(test_site, self.max_radius)	1✔
273	prob = self._calc_site_probabilities(test_site, nn)	1✔
274	all_prob.append(prob)	1✔
275	val = list(prob)	1✔
276	# Sort valences in order of decreasing probability.
277	val = sorted(val, key=lambda v: -prob[v])	1✔
278	# Retain probabilities that are at least 1/100 of highest prob.
279	valences.append(list(filter(lambda v: prob[v] > 0.01 * prob[val[0]], val)))	1✔
280	else:
281	full_all_prob = []	×
282	for sites in equi_sites:	×
283	test_site = sites[0]	×
284	nn = structure.get_neighbors(test_site, self.max_radius)	×
285	prob = self._calc_site_probabilities_unordered(test_site, nn)	×
286	all_prob.append(prob)	×
287	full_all_prob.extend(prob.values())	×
288	vals = []	×
289	for elsp, _ in get_z_ordered_elmap(test_site.species):	×
290	val = list(prob[elsp.symbol])	×
291	# Sort valences in order of decreasing probability.
292	val = sorted(val, key=lambda v: -prob[elsp.symbol][v])	×
293	# Retain probabilities that are at least 1/100 of highest
294	# prob.
295	vals.append(	×
296	list(
297	filter(
298	lambda v: prob[elsp.symbol][v] > 0.001 * prob[elsp.symbol][val[0]],
299	val,
300	)
301	)
302	)
303	valences.append(vals)	×
304
305	# make variables needed for recursion
306	if structure.is_ordered:	1✔
307	nsites = np.array([len(i) for i in equi_sites])	1✔
308	vmin = np.array([min(i) for i in valences])	1✔
309	vmax = np.array([max(i) for i in valences])	1✔
310
311	self._n = 0	1✔
312	self._best_score = 0	1✔
313	self._best_vset = None	1✔
314
315	def evaluate_assignment(v_set):	1✔
316	el_oxi = collections.defaultdict(list)	1✔
317	for i, sites in enumerate(equi_sites):	1✔
318	el_oxi[sites[0].specie.symbol].append(v_set[i])	1✔
319	max_diff = max(max(v) - min(v) for v in el_oxi.values())	1✔
320	if max_diff > 1:	1✔
321	return	1✔
322	score = functools.reduce(operator.mul, [all_prob[i][v] for i, v in enumerate(v_set)])	1✔
323	if score > self._best_score:	1✔
324	self._best_vset = v_set	1✔
325	self._best_score = score	1✔
326
327	def _recurse(assigned=None):	1✔
328	# recurses to find permutations of valences based on whether a
329	# charge balanced assignment can still be found
330	if self._n > self.max_permutations:	1✔
331	return None	×
332	if assigned is None:	1✔
333	assigned = []	1✔
334
335	i = len(assigned)	1✔
336	highest = vmax.copy()	1✔
337	highest[:i] = assigned	1✔
338	highest *= nsites	1✔
339	highest = np.sum(highest)	1✔
340
341	lowest = vmin.copy()	1✔
342	lowest[:i] = assigned	1✔
343	lowest *= nsites	1✔
344	lowest = np.sum(lowest)	1✔
345
346	if highest < 0 or lowest > 0:	1✔
347	self._n += 1	1✔
348	return None	1✔
349
350	if i == len(valences):	1✔
351	evaluate_assignment(assigned)	1✔
352	self._n += 1	1✔
353	return None	1✔
354	for v in valences[i]:	1✔
355	new_assigned = list(assigned)	1✔
356	_recurse(new_assigned + [v])	1✔
357	return None	1✔
358
359	else:
360	nsites = np.array([len(i) for i in equi_sites])	×
361	tmp = []	×
362	attrib = []	×
363	for insite, nsite in enumerate(nsites):	×
364	for _ in valences[insite]:	×
365	tmp.append(nsite)	×
366	attrib.append(insite)	×
367	new_nsites = np.array(tmp)	×
368	fractions = []	×
369	elements = []	×
370	for sites in equi_sites:	×
371	for sp, occu in get_z_ordered_elmap(sites[0].species):	×
372	elements.append(sp.symbol)	×
373	fractions.append(occu)	×
374	fractions = np.array(fractions, np.float_)	×
375	new_valences = []	×
376	for vals in valences:	×
377	for val in vals:	×
378	new_valences.append(val)	×
379	vmin = np.array([min(i) for i in new_valences], np.float_)	×
380	vmax = np.array([max(i) for i in new_valences], np.float_)	×
381
382	self._n = 0	×
383	self._best_score = 0	×
384	self._best_vset = None	×
385
386	def evaluate_assignment(v_set):	×
387	el_oxi = collections.defaultdict(list)	×
388	jj = 0	×
389	for sites in equi_sites:	×
390	for specie, _ in get_z_ordered_elmap(sites[0].species):	×
391	el_oxi[specie.symbol].append(v_set[jj])	×
392	jj += 1	×
393	max_diff = max(max(v) - min(v) for v in el_oxi.values())	×
394	if max_diff > 2:	×
395	return	×
396
397	score = functools.reduce(	×
398	operator.mul,
399	[all_prob[attrib[iv]][elements[iv]][vv] for iv, vv in enumerate(v_set)],
400	)
401	if score > self._best_score:	×
402	self._best_vset = v_set	×
403	self._best_score = score	×
404
405	def _recurse(assigned=None):	×
406	# recurses to find permutations of valences based on whether a
407	# charge balanced assignment can still be found
408	if self._n > self.max_permutations:	×
409	return None	×
410	if assigned is None:	×
411	assigned = []	×
412
413	i = len(assigned)	×
414	highest = vmax.copy()	×
415	highest[:i] = assigned	×
416	highest *= new_nsites	×
417	highest *= fractions	×
418	highest = np.sum(highest)	×
419
420	lowest = vmin.copy()	×
421	lowest[:i] = assigned	×
422	lowest *= new_nsites	×
423	lowest *= fractions	×
424	lowest = np.sum(lowest)	×
425
426	if highest < -self.charge_neutrality_tolerance or lowest > self.charge_neutrality_tolerance:	×
427	self._n += 1	×
428	return None	×
429
430	if i == len(new_valences):	×
431	evaluate_assignment(assigned)	×
432	self._n += 1	×
433	return None	×
434
435	for v in new_valences[i]:	×
436	new_assigned = list(assigned)	×
437	_recurse(new_assigned + [v])	×
438
439	return None	×
440
441	_recurse()	1✔
442
443	if self._best_vset:	1✔
444	if structure.is_ordered:	1✔
445	assigned = {}	1✔
446	for val, sites in zip(self._best_vset, equi_sites):	1✔
447	for site in sites:	1✔
448	assigned[site] = val	1✔
449
450	return [int(assigned[site]) for site in structure]	1✔
451	assigned = {}	×
452	new_best_vset = []	×
453	for _ in equi_sites:	×
454	new_best_vset.append([])	×
455	for ival, val in enumerate(self._best_vset):	×
456	new_best_vset[attrib[ival]].append(val)	×
457	for val, sites in zip(new_best_vset, equi_sites):	×
458	for site in sites:	×
459	assigned[site] = val	×
460
461	return [[int(frac_site) for frac_site in assigned[site]] for site in structure]	×
462	raise ValueError("Valences cannot be assigned!")	1✔
463
464	def get_oxi_state_decorated_structure(self, structure):	1✔
465	"""
466	Get an oxidation state decorated structure. This currently works only
467	for ordered structures only.
468
469	Args:
470	structure: Structure to analyze
471
472	Returns:
473	A modified structure that is oxidation state decorated.
474
475	Raises:
476	ValueError if the valences cannot be determined.
477	"""
478	s = structure.copy()	1✔
479	if s.is_ordered:	1✔
480	valences = self.get_valences(s)	1✔
481	s.add_oxidation_state_by_site(valences)	1✔
482	else:
483	valences = self.get_valences(s)	×
484	s = add_oxidation_state_by_site_fraction(s, valences)	×
485	return s	1✔
486
487
488	def get_z_ordered_elmap(comp):	1✔
489	"""
490	Arbitrary ordered elmap on the elements/species of a composition of a
491	given site in an unordered structure. Returns a list of tuples (
492	element_or_specie: occupation) in the arbitrary order.
493
494	The arbitrary order is based on the Z of the element and the smallest
495	fractional occupations first.
496	Example : {"Ni3+": 0.2, "Ni4+": 0.2, "Cr3+": 0.15, "Zn2+": 0.34,
497	"Cr4+": 0.11} will yield the species in the following order :
498	Cr4+, Cr3+, Ni3+, Ni4+, Zn2+ ... or
499	Cr4+, Cr3+, Ni4+, Ni3+, Zn2+
500	"""
501	return sorted((elsp, comp[elsp]) for elsp in comp)	×
502
503
504	def add_oxidation_state_by_site_fraction(structure, oxidation_states):	1✔
505	"""
506	Add oxidation states to a structure by fractional site.
507
508	Args:
509	oxidation_states (list): List of list of oxidation states for each
510	site fraction for each site.
511	E.g., [[2, 4], [3], [-2], [-2], [-2]]
512	"""
513	try:	×
514	for i, site in enumerate(structure):	×
515	new_sp = collections.defaultdict(float)	×
516	for j, (el, occu) in enumerate(get_z_ordered_elmap(site.species)):	×
517	specie = Species(el.symbol, oxidation_states[i][j])	×
518	new_sp[specie] += occu	×
519	structure[i] = new_sp	×
520	return structure	×
521	except IndexError:	×
522	raise ValueError("Oxidation state of all sites must be specified in the list.")	×

materialsproject / pymatgen / 4075885785

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous