Metrics and analysis
The evaluation layer defines structure validity, uniqueness, novelty, XRD agreement, property targets, and stability metrics for generated CIFs.
The validity predicates define the structural checks used by the evaluation workflows. The scoring CLIs operate on post-processed generation results stored in parquet files.
Validity predicates
These predicates return booleans and can be applied to model output that may fail to parse.
is_valid(cif_str, bond_length_acceptability_cutoff=1.0, allow_stated_p1_mismatch=False, debug=False)
Check whether a generated CIF passes the structural validity checks.
The checks are applied in order and stop at the first failure. The formula must match the atom
sites, atom-site multiplicities must be self-consistent, the bond-length score must reach
bond_length_acceptability_cutoff, and the detected symmetry must match the declared space
group. allow_stated_p1_mismatch permits a declared P1 when the structure has higher symmetry.
is_sensible is not part of this composite check.
Source code in _utils/validity.py
| def is_valid(cif_str: str, bond_length_acceptability_cutoff: float=1.0, allow_stated_p1_mismatch: bool=False, debug: bool=False) -> bool:
"""Check whether a generated CIF passes the structural validity checks.
The checks are applied in order and stop at the first failure. The formula must match the atom
sites, atom-site multiplicities must be self-consistent, the bond-length score must reach
`bond_length_acceptability_cutoff`, and the detected symmetry must match the declared space
group. `allow_stated_p1_mismatch` permits a declared P1 when the structure has higher symmetry.
`is_sensible` is not part of this composite check.
"""
if not is_formula_consistent(cif_str):
if debug:
print(f"Formula is inconsistent for {cif_str}")
return False
if not is_atom_site_multiplicity_consistent(cif_str):
if debug:
print(f"Atom site multiplicity is inconsistent for {cif_str}")
return False
bond_length_score = bond_length_reasonableness_score(cif_str)
if bond_length_score is not None and bond_length_score < bond_length_acceptability_cutoff:
if debug:
print(f"Bond length is unreasonable for {cif_str}")
return False
if not is_space_group_consistent(cif_str, allow_stated_p1_mismatch=allow_stated_p1_mismatch):
if debug:
print(f"Space group is inconsistent for {cif_str}")
return False
return True
|
is_sensible(cif_str, length_lo=0.5, length_hi=1000.0, angle_lo=10.0, angle_hi=170.0)
Check whether the unit-cell dimensions fall within the supplied physical bounds.
The check uses only cell parameters parsed from the CIF text. All lengths must lie between
length_lo and length_hi in Å, and all angles must lie between angle_lo and angle_hi in
degrees. This check is separate from is_valid.
Source code in _utils/validity.py
| def is_sensible(cif_str: str, length_lo: float=0.5, length_hi: float=1000., angle_lo: float=10., angle_hi: float=170.) -> bool:
"""Check whether the unit-cell dimensions fall within the supplied physical bounds.
The check uses only cell parameters parsed from the CIF text. All lengths must lie between
`length_lo` and `length_hi` in Å, and all angles must lie between `angle_lo` and `angle_hi` in
degrees. This check is separate from `is_valid`.
"""
cell_length_pattern = re.compile(r"_cell_length_[abc]\s+([\d\.]+)")
cell_angle_pattern = re.compile(r"_cell_angle_(alpha|beta|gamma)\s+([\d\.]+)")
cell_lengths = cell_length_pattern.findall(cif_str)
for length_str in cell_lengths:
length = float(length_str)
if length < length_lo or length > length_hi:
return False
cell_angles = cell_angle_pattern.findall(cif_str)
for _, angle_str in cell_angles:
angle = float(angle_str)
if angle < angle_lo or angle > angle_hi:
return False
return True
|
is_formula_consistent(cif_str)
Check whether the chemical formula declared by the CIF matches its atom-site composition.
Compared up to cell scale, so a supercell matches the formula it is a supercell of.
Source code in _utils/validity.py
| def is_formula_consistent(cif_str: str) -> bool:
"""Check whether the chemical formula declared by the CIF matches its atom-site composition.
Compared up to cell scale, so a supercell matches the formula it is a supercell of.
"""
try:
parser = cif_parser(cif_str)
cif_data = parser.as_dict()
key = list(cif_data.keys())[0]
formula_data = Composition(extract_data_formula(cif_str))
formula_sum = Composition(cif_data[key].get("_chemical_formula_sum", ""))
formula_structural = Composition(cif_data[key].get("_chemical_formula_structural", ""))
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=UserWarning)
try:
structure = parser.parse_structures(primitive=False)[0]
except Exception:
# Some valid disordered CIFs need occupancy rescaling to parse.
parser = cif_parser(cif_str, occupancy_tolerance=2.0)
structure = parser.parse_structures(primitive=False)[0]
formula_geometry = structure.composition
names_match = (
formula_data.reduced_formula == formula_sum.reduced_formula ==
formula_structural.reduced_formula
)
return names_match and _compositions_match(formula_sum, formula_geometry)
except Exception:
return False
|
is_space_group_consistent(cif_str, allow_stated_p1_mismatch=False)
Check whether a CIF's structure matches its declared space group.
allow_stated_p1_mismatch permits a CIF declaring P1 when the detected structure has higher
symmetry.
Source code in _utils/validity.py
| def is_space_group_consistent(cif_str: str, allow_stated_p1_mismatch: bool=False) -> bool:
"""Check whether a CIF's structure matches its declared space group.
`allow_stated_p1_mismatch` permits a CIF declaring P1 when the detected structure has higher
symmetry.
"""
structure = cif_parser(cif_str).parse_structures(primitive=False)[0]
parser = cif_parser(cif_str)
cif_data = parser.as_dict()
# Extract the stated space group from the CIF file
stated_space_group = cif_data[list(cif_data.keys())[0]]['_symmetry_space_group_name_H-M']
# Analyze the symmetry of the structure
spacegroup_analyzer = SpacegroupAnalyzer(structure, symprec=0.1)
# Get the detected space group
detected_space_group = spacegroup_analyzer.get_space_group_symbol()
# Check if the detected space group matches the stated space group
is_match = stated_space_group.strip() == detected_space_group.strip()
if not is_match and allow_stated_p1_mismatch:
stated_normalized = stated_space_group.replace(" ", "").upper()
if stated_normalized == "P1":
return True
return is_match
|
is_atom_site_multiplicity_consistent(cif_str)
Check whether each atom site's stated multiplicity matches its symmetry orbit.
Source code in _utils/validity.py
| def is_atom_site_multiplicity_consistent(cif_str: str) -> bool:
# Parse the CIF string
"""Check whether each atom site's stated multiplicity matches its symmetry orbit."""
parser = cif_parser(cif_str)
cif_data = parser.as_dict()
# Extract the chemical formula sum from the CIF data
formula_sum = cif_data[list(cif_data.keys())[0]]["_chemical_formula_sum"]
# Convert the formula sum into a dictionary
expected_atoms = Composition(formula_sum).as_dict()
# Count the atoms provided in the _atom_site_type_symbol section
actual_atoms = {}
for key in cif_data:
if "_atom_site_type_symbol" in cif_data[key] and "_atom_site_symmetry_multiplicity" in cif_data[key]:
for atom_type, multiplicity in zip(cif_data[key]["_atom_site_type_symbol"],
cif_data[key]["_atom_site_symmetry_multiplicity"]):
if atom_type in actual_atoms:
actual_atoms[atom_type] += int(multiplicity)
else:
actual_atoms[atom_type] = int(multiplicity)
# Validate if the expected and actual atom counts match
return expected_atoms == actual_atoms
|
bond_length_reasonableness_score(cif_str, tolerance=0.32, h_factor=2.5)
Score bond lengths against the sum of covalent radii.
The score is the fraction of bonds whose length lies within tolerance of the expected
covalent-radii sum. A score of 1.0 means every counted bond passes. Neighbours are obtained
with pymatgen's CrystalNN. Bonds involving hydrogen use a widened tolerance controlled by
h_factor.
Returns None for a disordered structure, meaning not checked rather than passed.
Source code in _utils/validity.py
| def bond_length_reasonableness_score(cif_str: str, tolerance: float=0.32, h_factor: float=2.5) -> float | None:
"""Score bond lengths against the sum of covalent radii.
The score is the fraction of bonds whose length lies within `tolerance` of the expected
covalent-radii sum. A score of `1.0` means every counted bond passes. Neighbours are obtained
with pymatgen's `CrystalNN`. Bonds involving hydrogen use a widened tolerance controlled by
`h_factor`.
Returns `None` for a disordered structure, meaning not checked rather than passed.
"""
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=UserWarning)
structure = cif_parser(cif_str).parse_structures(primitive=False)[0]
# TODO: make function work for disordered materials. CrystalNN reads site.specie,
# which raises on partial occupancies, as do the radii lookups below.
if not structure.is_ordered:
return None
crystal_nn = CrystalNN()
min_ratio = 1 - tolerance
max_ratio = 1 + tolerance
# calculate the score based on bond lengths and covalent radii
score = 0
bond_count = 0
for i, site in enumerate(structure):
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=UserWarning)
bonded_sites = crystal_nn.get_nn_info(structure, i)
for connected_site_info in bonded_sites:
j = connected_site_info['site_index']
if i == j: # skip if they're the same site
continue
connected_site = connected_site_info['site']
bond_length = site.distance(connected_site)
is_hydrogen_bond = "H" in [site.specie.symbol, connected_site.specie.symbol]
electronegativity_diff = abs(site.specie.X - connected_site.specie.X)
"""
According to the Pauling scale, when the electronegativity difference
between two bonded atoms is less than 1.7, the bond can be considered
to have predominantly covalent character, while a difference greater
than or equal to 1.7 indicates that the bond has significant ionic
character.
"""
if electronegativity_diff >= 1.7:
# use ionic radii
if site.specie.X < connected_site.specie.X:
expected_length = site.specie.average_cationic_radius + connected_site.specie.average_anionic_radius
else:
expected_length = site.specie.average_anionic_radius + connected_site.specie.average_cationic_radius
else:
expected_length = site.specie.atomic_radius + connected_site.specie.atomic_radius
bond_ratio = bond_length / expected_length
# Penalise bond lengths that are too short or too long.
# Hydrogen bonds use a looser tolerance.
if is_hydrogen_bond:
if bond_ratio < h_factor:
score += 1
else:
if min_ratio < bond_ratio < max_ratio:
score += 1
bond_count += 1
normalized_score = score / bond_count if bond_count > 0 else 0
return normalized_score
|
Density is calculated from the parsed structure.
Compute the crystallographic density of a CIF in g/cm³.
Returns NaN when the CIF cannot be parsed or pymatgen reports incorrect stoichiometry, so
batch scoring keeps invalid rows instead of raising.
Source code in _utils/validity.py
| def get_density(cif: str) -> float:
"""Compute the crystallographic density of a CIF in g/cm³.
Returns `NaN` when the CIF cannot be parsed or pymatgen reports incorrect stoichiometry, so
batch scoring keeps invalid rows instead of raising.
"""
try:
with warnings.catch_warnings(record=True) as w:
warnings.simplefilter("always")
structure = cif_parser(cif).parse_structures(primitive=False)[0]
for warn in w:
if ("Incorrect stoichiometry" in str(warn.message)):
return np.nan
return structure.density
except Exception:
return np.nan
|
Evaluation CLIs
Validity evaluation
Evaluate generated CIFs for validity, consistency, and crystallographic sanity.
Generated CIFs can be evaluated directly without postprocessing. The valid subset can optionally be
written to a new parquet.
From original repo: https://github.com/lantunes/CrystaLLM/blob/main/bin/evaluate_cifs.py
Usage
python _utils/_generating/evaluate_cifs.py --input_parquet gen.parquet \
--output_parquet valid.parquet
Validity, uniqueness, and novelty
Score generated CIFs for validity, uniqueness, and novelty.
The VUN metrics are computed over generated structures. Compositional novelty requires the training
data, supplied through --huggingface_dataset or a processed cache loaded with
--load_processed_data.
Usage
python _utils/_scoring/vun_metrics.py --input_parquet gen.parquet \
--huggingface_dataset c-bone/mpdb-2prop_clean --output_csv vun.csv
Count structures that pass each VUN criterion.
Source code in _utils/_scoring/vun_metrics.py
| def compute_vun_metrics(df: pd.DataFrame) -> dict[str, int]:
"""Count structures that pass each VUN criterion."""
metrics = {
"total": len(df),
"valid": int(df['is_valid'].sum()),
"unique": int(df['is_unique'].sum()),
"novel": int(df['is_novel'].sum()),
}
if 'is_comp_novel' in df.columns:
metrics['comp_novel'] = int(df['is_comp_novel'].sum())
return metrics
|
XRD structure match
Compute structure-matching metrics for XRD-conditioned generation.
The workflow reports match rate and mean RMS distance over matched pairs using DiffCSP-compliant
validity checks. It also reports lattice-parameter differences and matched counts over all generated
structures.
Benchmark metrics are adapted from: https://github.com/jiaor17/DiffCSP/tree/main as well as the
original https://github.com/lantunes/CrystaLLM/tree/main repo
Usage
python _utils/_scoring/xrd_metrics.py --input_parquet gen.parquet \
--ref_parquet reference.parquet --num_gens 20 --output_parquet xrd_metrics.parquet
get_match_rate_and_rms(gen_structs, true_structs, matcher, args, score_data=None, num_workers=None, atom_counts=None, novel_data=None, material_ids_order=None)
Compute XRD match rate, RMS distance, and related metrics.
Only structures passing SMACT and structure-validity checks are considered. The minimum RMS
distance among valid matches is used for each material. Match rate is the fraction of materials
with at least one valid StructureMatcher match. Lattice-parameter matching is used only as an
analysis fallback and does not contribute to match rate.
Adapted from:
https://github.com/jiaor17/DiffCSP/blob/ee131b03a1c6211828e8054d837caa8f1a980c3e/scripts/compute_metrics.py
Parameters:
| Name |
Type |
Description |
Default |
score_data
|
dict | None
|
Mapping from material id to scores corresponding to generated structures.
|
None
|
num_workers
|
int | None
|
Number of parallel workers. Defaults to half the available CPU count.
|
None
|
atom_counts
|
list | None
|
True structure atom counts, one per material in the same order as the generated
and reference structures.
|
None
|
novel_data
|
dict | None
|
Additional columns mapped from material id to values.
|
None
|
material_ids_order
|
list | None
|
Material ids in the same order as the generated and reference
structures.
|
None
|
Source code in _utils/_scoring/xrd_metrics.py
| def get_match_rate_and_rms(gen_structs, true_structs, matcher, args, score_data: dict | None = None,
num_workers: int | None = None, atom_counts: list | None = None,
novel_data: dict | None = None, material_ids_order: list | None = None):
"""Compute XRD match rate, RMS distance, and related metrics.
Only structures passing SMACT and structure-validity checks are considered. The minimum RMS
distance among valid matches is used for each material. Match rate is the fraction of materials
with at least one valid `StructureMatcher` match. Lattice-parameter matching is used only as an
analysis fallback and does not contribute to match rate.
Adapted from:
https://github.com/jiaor17/DiffCSP/blob/ee131b03a1c6211828e8054d837caa8f1a980c3e/scripts/compute_metrics.py
Args:
score_data: Mapping from material id to scores corresponding to generated structures.
num_workers: Number of parallel workers. Defaults to half the available CPU count.
atom_counts: True structure atom counts, one per material in the same order as the generated
and reference structures.
novel_data: Additional columns mapped from material id to values.
material_ids_order: Material ids in the same order as the generated and reference
structures.
"""
if num_workers is None:
num_workers = max(1, multiprocessing.cpu_count() // 2)
# Prepare data for parallel processing
material_ids = list(score_data.keys()) if score_data else []
matcher_params = {
'stol': matcher.stol,
'angle_tol': matcher.angle_tol,
'ltol': matcher.ltol
}
# Build work items
work_items = []
for i in range(len(gen_structs)):
current_material_id = material_ids[i] if i < len(material_ids) else None
current_scores = score_data.get(current_material_id, []) if score_data and current_material_id else []
work_items.append((
i, gen_structs[i], true_structs[i], current_material_id, current_scores,
args.validity_check, matcher_params
))
# Process in parallel
results = []
if num_workers == 1:
# Serial processing for debugging or single-core systems
for work_item in tqdm(work_items, desc="Comparing structures"):
results.append(_process_material_comparison(work_item))
else:
# Parallel processing
with ProcessPoolExecutor(max_workers=num_workers) as executor:
results = list(tqdm(
executor.map(_process_material_comparison, work_items),
total=len(work_items),
desc="Comparing structures"
))
# Collect results in original order
rms_dists = [None] * len(gen_structs)
a_diffs = [None] * len(gen_structs)
b_diffs = [None] * len(gen_structs)
c_diffs = [None] * len(gen_structs)
rows = [None] * len(gen_structs)
for result in results:
if result and len(result) == 7:
i, rms_dist, a_diff, b_diff, c_diff, row, n_matched_struct = result
if i is not None and 0 <= i < len(gen_structs):
rms_dists[i] = rms_dist
a_diffs[i] = a_diff
b_diffs[i] = b_diff
c_diffs[i] = c_diff
rows[i] = row
# Handle any failed results by creating empty rows
for i in range(len(gen_structs)):
if rows[i] is None:
current_material_id = material_ids[i] if i < len(material_ids) else None
current_scores = score_data.get(current_material_id, []) if score_data and current_material_id else []
best_score = current_scores[0] if current_scores else None
rows[i] = _create_empty_row(true_structs[i], len(gen_structs[i]), best_score)
stats_df = pd.DataFrame(rows)
if atom_counts is not None:
stats_df["atom_counts"] = atom_counts
if novel_data and material_ids_order:
for col_name, mid_to_val in novel_data.items():
stats_df[col_name] = [mid_to_val.get(mid) for mid in material_ids_order]
return _calculate_metrics(rms_dists, a_diffs, b_diffs, c_diffs, gen_structs), stats_df
|
Check benchmark validity using both composition and structure validity.
Adapted from:
https://github.com/jiaor17/DiffCSP/blob/ee131b03a1c6211828e8054d837caa8f1a980c3e/scripts/compute_metrics.py
Source code in _utils/_scoring/xrd_metrics.py
| def is_valid_bench(struct):
"""Check benchmark validity using both composition and structure validity.
Adapted from:
https://github.com/jiaor17/DiffCSP/blob/ee131b03a1c6211828e8054d837caa8f1a980c3e/scripts/compute_metrics.py
"""
elem_counter = Counter([specie.Z for specie in struct.species])
elems = [(elem, elem_counter[elem]) for elem in sorted(elem_counter.keys())]
comp, elem_counts = list(zip(*elems))
elem_counts = np.array(elem_counts)
elem_counts = elem_counts / np.gcd.reduce(elem_counts)
count = tuple(elem_counts.astype("int").tolist())
comp_valid = smact_validity(comp, count)
struct_valid = structure_validity(struct)
return comp_valid and struct_valid
|
Property targets
Score generated structures against their conditioned property targets.
Targets are normalised using the training convention, properties are predicted from the generated
structures, and errors are grouped by condition vector so each target value is scored separately.
Usage
python _utils/_scoring/property_metrics.py --post_parquet gen_post.parquet \
--output_metrics density_metrics.csv --property_targets density
Stability
Compute energy above hull using the MACE potential and MP2020 corrections.
Each generated structure is relaxed with MACE and placed on the Materials Project convex hull.
MP2020 compatibility corrections make the resulting energies comparable with Materials Project
entries.
Inspired by: https://github.com/facebookresearch/crystal-text-llm/blob/main/e_above_hull.py
Usage
python _utils/_scoring/mace_ehull.py --post_parquet gen_post.parquet \
--output_parquet ehull.parquet --mp_data mp_data.pkl --num_workers 4
Compute energy above hull and formation energy from DFT energies.
Input energies are rescaled against Materials Project reference data so hull distances are
comparable with Materials Project entries. Novelty and density can also be reported.
Usage
python _utils/_scoring/dft_ehull.py --input_csv dft_energies.csv \
--output_parquet ehull.parquet --num_workers 8
Materials Project reference data
Materials Project entries are downloaded and cached for the stability metrics.
Provide Materials Project data with on-demand phase-diagram construction.
Source code in _utils/mp_data.py
| class MPDataProvider:
"""Provide Materials Project data with on-demand phase-diagram construction."""
def __init__(self, mp_data_path: str) -> None:
print("Loading MP entries...")
df_mp = pd.read_json(mp_data_path)
# Filter to only GGA entries (like reference script)
if 'index' in df_mp.columns:
df_mp = df_mp[df_mp['index'].str.contains("GGA")]
print(f"Found {len(df_mp)} MP entries")
# Convert to ComputedEntry objects and organize by chemical system
print("Organizing MP entries by chemical system...")
self.entries_by_chemsys = {}
for entry_dict in tqdm(df_mp.entry, desc="Processing MP entries"):
if 'GGA' in entry_dict['parameters']['run_type']:
entry = ComputedEntry.from_dict(entry_dict)
# Filter out R2SCAN entries (exactly like mace_ehull_copy.py)
if not np.any(['R2SCAN' in a.name for a in entry.energy_adjustments]):
# Get chemical system (sorted elements)
elements = sorted(entry.composition.elements)
chemsys = tuple(str(el) for el in elements)
if chemsys not in self.entries_by_chemsys:
self.entries_by_chemsys[chemsys] = []
self.entries_by_chemsys[chemsys].append(entry)
# Cache for built phase diagrams
self.pd_cache = {}
print(f"Organized {sum(len(v) for v in self.entries_by_chemsys.values())} entries across {len(self.entries_by_chemsys)} chemical systems")
def get_phase_diagram(self, elements: frozenset) -> object:
"""Return or construct the phase diagram for a chemical system."""
chemsys = tuple(sorted(str(el) for el in elements))
if chemsys in self.pd_cache:
return self.pd_cache[chemsys]
# Get ALL entries that contain ONLY these elements (like MP API's get_entries_in_chemsys)
# This includes pure element entries and all combinations
relevant_entries = []
element_set = set(str(el) for el in elements)
for system, entries in self.entries_by_chemsys.items():
# Check if this chemical system is a subset of our target elements
system_elements = set(system)
if system_elements.issubset(element_set):
relevant_entries.extend(entries)
if len(relevant_entries) < 2:
# Not enough entries to build a meaningful phase diagram
self.pd_cache[chemsys] = None
return None
# Build phase diagram with all relevant entries
pd_sys = PhaseDiagram(relevant_entries)
self.pd_cache[chemsys] = pd_sys
return pd_sys
def compute_ehull_and_eform(self, structure: Structure, energy_eV: float) -> tuple:
"""Compute energy above hull and formation energy."""
elements = sorted({el.symbol for el in structure.composition.elements})
pd_sys = self.get_phase_diagram(elements)
if pd_sys is None:
return np.nan, np.nan
entry = ComputedEntry(composition=structure.composition, energy=energy_eV)
# Apply MP2020 corrections only to user entry (same as mace_ehull_copy.py)
compat = MaterialsProject2020Compatibility(check_potcar=False)
try:
entry.parameters["software"] = "non-vasp"
entry.parameters["run_type"] = "GGA"
entry = compat.process_entry(entry, clean=True)
if entry is None:
return np.nan, np.nan
except Exception:
return np.nan, np.nan
eh = pd_sys.get_e_above_hull(entry, allow_negative=True)
eform_pa = pd_sys.get_form_energy_per_atom(entry)
return eh, eform_pa
|
get_phase_diagram(elements)
Return or construct the phase diagram for a chemical system.
Source code in _utils/mp_data.py
| def get_phase_diagram(self, elements: frozenset) -> object:
"""Return or construct the phase diagram for a chemical system."""
chemsys = tuple(sorted(str(el) for el in elements))
if chemsys in self.pd_cache:
return self.pd_cache[chemsys]
# Get ALL entries that contain ONLY these elements (like MP API's get_entries_in_chemsys)
# This includes pure element entries and all combinations
relevant_entries = []
element_set = set(str(el) for el in elements)
for system, entries in self.entries_by_chemsys.items():
# Check if this chemical system is a subset of our target elements
system_elements = set(system)
if system_elements.issubset(element_set):
relevant_entries.extend(entries)
if len(relevant_entries) < 2:
# Not enough entries to build a meaningful phase diagram
self.pd_cache[chemsys] = None
return None
# Build phase diagram with all relevant entries
pd_sys = PhaseDiagram(relevant_entries)
self.pd_cache[chemsys] = pd_sys
return pd_sys
|
compute_ehull_and_eform(structure, energy_eV)
Compute energy above hull and formation energy.
Source code in _utils/mp_data.py
| def compute_ehull_and_eform(self, structure: Structure, energy_eV: float) -> tuple:
"""Compute energy above hull and formation energy."""
elements = sorted({el.symbol for el in structure.composition.elements})
pd_sys = self.get_phase_diagram(elements)
if pd_sys is None:
return np.nan, np.nan
entry = ComputedEntry(composition=structure.composition, energy=energy_eV)
# Apply MP2020 corrections only to user entry (same as mace_ehull_copy.py)
compat = MaterialsProject2020Compatibility(check_potcar=False)
try:
entry.parameters["software"] = "non-vasp"
entry.parameters["run_type"] = "GGA"
entry = compat.process_entry(entry, clean=True)
if entry is None:
return np.nan, np.nan
except Exception:
return np.nan, np.nan
eh = pd_sys.get_e_above_hull(entry, allow_negative=True)
eform_pa = pd_sys.get_form_energy_per_atom(entry)
return eh, eform_pa
|
Virtual crystals
Convert ordered structures to disordered virtual crystals and promote them to higher-symmetry parent structures.
Virtualise selected element groups in an ordered crystal structure.
The tool replaces each selected group (two or more elements) on a shared sublattice with fractional
mixed occupancy, refines the resulting structure to higher symmetry, and writes the virtual crystal
as a CIF. Vacancies are not handled explicitly, and all sites containing members of a group are
assumed to belong to the same sublattice.
Contribution by Dr Ricardo Grau-Crespo
https://github.com/rgraucrespo
Usage
python _utils/_virtualiser/virtualiser.py --in Mg3ZnO4.cif --config config.yaml --out virtual.cif
virtualise_structure(struct, virtual_pairs)
Merge grouped elements onto shared sites with fractional occupancy.
Each site containing a member of a selected group is replaced by a mixed site weighted by the
group fractions. Absent groups are skipped and existing disordered sites are preserved.
Oxidation states are removed from the result before symmetry processing.
Source code in _utils/_virtualiser/virtualiser.py
| def virtualise_structure(struct: Structure, virtual_pairs: list[tuple[str, ...]]) -> Structure:
# Build a mapping from elements that are in any group to their group-fractions
"""Merge grouped elements onto shared sites with fractional occupancy.
Each site containing a member of a selected group is replaced by a mixed site weighted by the
group fractions. Absent groups are skipped and existing disordered sites are preserved.
Oxidation states are removed from the result before symmetry processing.
"""
replace_map: dict[str, dict[str, float]] = {}
for group in virtual_pairs:
fracs = compute_pair_fractions(struct, group)
if all(f == 0.0 for f in fracs.values()):
continue
for el in group:
replace_map[el] = fracs
new_species = []
new_coords = []
for site in struct.sites:
if len(site.species) == 1:
el = list(site.species.as_dict().keys())[0]
el = str(el)
if el in replace_map:
fracs = replace_map[el]
spec_map = {Element(k): float(v) for k, v in fracs.items() if v > 0.0}
new_species.append(spec_map)
else:
new_species.append(site.species)
else:
# already disordered, keep as-is
new_species.append(site.species)
new_coords.append(site.frac_coords)
virt = Structure(struct.lattice, new_species, new_coords, coords_are_cartesian=False,
site_properties=struct.site_properties if struct.site_properties else None)
virt.remove_oxidation_states() # spglib requires species without oxidation states
return virt
|
promote_symmetry(struct, symprec, angle_tol)
Refine a structure to the highest symmetry consistent with the supplied tolerances.
If refinement fails, the conventional standard cell is used instead.
Source code in _utils/_virtualiser/virtualiser.py
| def promote_symmetry(struct: Structure, symprec: float, angle_tol: float) -> Structure:
"""Refine a structure to the highest symmetry consistent with the supplied tolerances.
If refinement fails, the conventional standard cell is used instead.
"""
sga = SpacegroupAnalyzer(struct, symprec=symprec, angle_tolerance=angle_tol)
try:
refined = sga.get_refined_structure()
except Exception:
refined = sga.get_conventional_standard_structure()
return refined
|