Generation of low-entropy structures

[yutaoki]

Is it possible to generate top N low-entropy structures before generating all the possible structures?
It is related to #13 but this issue would cover wider concept.

Low-entropy structures can be categorized into the following two types:

1. High-symmetry structures. It is already solved in Generate only higher-symmetry substructures #13.
2. Impurity-clustered structures. It is important to model segregated structures in alloy systems.

Degree of entropy can be defined by the multiplicity (the number of the symmetrically equivalent configurations) of the structure. I guess it is mentioned as "weight" in shry.

When we model the segregation in alloy systems, we usually employ a relatively huge supercell so if we generate all the symmetrically unique structures, we often encounter the out-of-memory error. So this function would be very helpful.

[kousuke-nakano]

@yutaoki

Thank you for the issue!

Degree of entropy can be defined by the multiplicity (the number of the symmetrically equivalent configurations) of the structure. I guess it is mentioned as "weight" in shry.

Yes. The weight refers to the multiplicity in SHRY.

When we model the segregation in alloy systems, we usually employ a relatively huge supercell so if we generate all the symmetrically unique structures, we often encounter the out-of-memory error. So this function would be very helpful.

Umm, I do not think the current implementation of SHRY is well suited for this purpose. The SHRY's enumeration relies on a tree search that is pruned using global crystal symmetry and local site symmetries, which cannot be directly applied for truncating cluster searches. For your purpose, a sampling-based approach (e.g., randomized search) is usually more practical than exhaustive enumeration. If there is no suitable existing software for this use case, I’d be happy to consider implementing such a feature in SHRY. Do you know of any alternative software packages that already provide this functionality?

[yutaoki]

@kousuke-nakano
Thank you for your comment.

Do you know of any alternative software packages that already provide this functionality?

I don't know actually. If we pursue an approach which doesn't employ exhaustive enumeration, aevo might be a solution.

https://github.com/atomisticnet/aevo

It searches the global optimized substitution/vacancy configurations using an evolutionary algorithm. If we set a fitness function related to the multiplicity, it may work. However, evaluating the multiplicity directly for a given structure doesn't seem to be so easy so I think utilizing related functions in enumerating software such as SHRY would be a better way. Also, it expect that it would be possible to rank the structure based on the multiplicity and it would be more convenient.

[kousuke-nakano]

@yutaoki Sure. Thank you. I am considering it. Is SQS also useful? To implement what you requested, we need not only exhaustive enumeration but also some heuristic search algorithms, such as SA. SQS is along the context.

[yutaoki]

@kousuke-nakano Thank you. In my context, I'm not sure whether SQS is useful because the desired structures are low-entropy structures. In my understanding, SQS is a representative random structure so it sounds opposite feature. But I agree that we need some heuristic algorithm when the search space is tremendously huge.

[kousuke-nakano]

Dear @yutaoki

I think the following script can write top-N structures for 2)

"""Select impurity-clustered structures using a distance-matrix score.

This example streams structures from SHRY without writing intermediate CIFs,
scores each structure by how strongly a target species is clustered, and keeps
only TOP_N best structures using a heap.

Primary score (smaller is better):
    mean nearest-neighbor distance among target sites.

Tie-breakers:
    1) mean pairwise distance among target sites (smaller is better),
    2) SHRY weight (smaller is better / higher symmetry),
    3) generation index.
"""

import math
import os
from heapq import heappush, heapreplace

import numpy as np
from pymatgen.core import Structure
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer

from shry import _ScriptHelper


# ---------------------------------------------------------------------------
# User settings
# ---------------------------------------------------------------------------
cif_file = "PbSnTe.cif"  # Replace with your CIF if desired

helper = _ScriptHelper(
    structure_file=cif_file,
    from_species=("Sn",),
    to_species=("Sn0.25Sb0.75",),
    scaling_matrix=(2, 2, 2),
)

TARGET_SYMBOL = "Sn"  # Species to evaluate clustering for
TOP_N = 5

output_dir = "output_cluster_top"
os.makedirs(output_dir, exist_ok=True)


# Heap of fixed size TOP_N, keeping "worst among current best" at heap[0]
# by storing negated score keys.
best = []

for i, packet in enumerate(helper.substitutor.quantities(("weight", "structure"))):
    weight = float(packet["weight"])
    structure = packet["structure"]

    # Collect target-species site indices.
    target_indices = []
    for site_index, site in enumerate(structure.sites):
        try:
            symbol = site.specie.symbol
        except Exception:
            symbol = site.species_string
        if symbol == TARGET_SYMBOL:
            target_indices.append(site_index)

    if len(target_indices) < 2:
        nn_mean = math.inf
        pair_mean = math.inf
    else:
        dm = np.array(structure.distance_matrix)[np.ix_(target_indices, target_indices)]

        # Primary score: mean nearest-neighbor distance among target atoms.
        dm_nn = dm.copy()
        np.fill_diagonal(dm_nn, np.inf)
        nn_mean = float(np.min(dm_nn, axis=1).mean())

        # Tie-break score: mean pairwise distance among target atoms.
        iu = np.triu_indices_from(dm, k=1)
        pair_mean = float(dm[iu].mean())

    key = (nn_mean, pair_mean, weight, i)
    neg_key = (-nn_mean, -pair_mean, -weight, -i)

    record = (neg_key, key, i, structure)
    if len(best) < TOP_N:
        heappush(best, record)
    elif neg_key > best[0][0]:
        heapreplace(best, record)


for rank, (_neg_key, key, idx, structure) in enumerate(sorted(best, key=lambda x: x[1])):
    nn_mean, pair_mean, weight, _ = key

    # Drop site properties so SpacegroupAnalyzer can hash/cache safely.
    plain_structure = Structure(
        lattice=structure.lattice,
        species=structure.species,
        coords=structure.frac_coords,
        coords_are_cartesian=False,
    )
    sga = SpacegroupAnalyzer(plain_structure)
    sg_number = sga.get_space_group_number()
    sg_symbol = sga.get_space_group_symbol().replace("/", "_")

    filename = os.path.join(
        output_dir,
        (
            f"cluster_rank{rank:03d}"
            f"_nn{nn_mean:.4f}"
            f"_pair{pair_mean:.4f}"
            f"_w{int(weight)}"
            f"_{sg_number}_{sg_symbol}"
            f"_idx{idx}.cif"
        ),
    )
    structure.to(fmt="cif", filename=filename)
    print(f"Saved {filename}")

[yutaoki]

@kousuke-nakano
Thank you! Let met try.

[kousuke-nakano]

@yutaoki. Thank you. Let me know if it works. I will merge it into the script helper if it works as you thought.

[yutaoki]

@kousuke-nakano
Hi, I have tested it for Al3Cu alloy based on the fcc Al 2x2x2 supercell. I expected that the structure with the planar Cu defect would be ranked as the top but the result was a bit different as attached (but still very close to what I expected!)

Maybe the primary scores for at least top 20 results are totally same. The definition of the tie-break score leads to slightly different behavior from what I originally expected.

Anyway, I think that the current definition for the tie-break score is very suitable for the three-dimensional clustering behavior (ball-like clustering). Maybe we can consider several different tie-break scores such as for planar defects and ball-like clusters and make users can switch them. Deviation of the pair-wise distance among the target elements may be a good tie-break score for planar defects.

I would test the following metrics as the tie-break score and let you know the result.
pair_std = float(dm[iu].std())

[kousuke-nakano]

@yutaoki

I would test the following metrics as the tie-break score and let you know the result.

Thank you! Please let me know if it works.