Tools#

Mapping structures#

icet.tools.map_structure_to_reference(structure, reference, *, inert_species=None, tol_positions=0.0001, tol_cell=0.25, symprec=1e-05, find_translation=True, trigger_levels=None, suppress_warnings=False, assume_no_cell_relaxation=False)[source]#

Maps a structure onto a reference structure. This is often desirable when, for example, a structure has been relaxed using DFT, and one wants to use it as a training structure in a cluster expansion.

The function returns a tuple comprising the ideal supercell most closely matching the input structure and a StructureMapping carrying the supplementary information about the mapping. The latter reports the relaxation distances, the transformation matrix, any rigid offset that was removed, the sites whose assignment is ambiguous, the warnings that were triggered, the strain and its separation into a change of volume and a change of shape, and the rotation. Its attributes are documented there.

The rotation is reported as disorientation_angle and disorientation_axis, reduced over the symmetry of the reference structure, because a supercell related to the returned one by a symmetry operation describes the same crystal, so only the smallest rotation among those descriptions means anything on its own. Successive rotations of a structure are not reported separately, since any number of them compose into a single rotation about a single axis.

The returned Atoms object carries the following per-atom arrays, each with one entry per site of the returned supercell rather than per atom of the input structure.

Displacement, shape (n, 3), the displacement in Ångstrom pointing from the site to the position of the atom that occupies it,
Displacement_Magnitude, shape (n,), its norm,
Minimum_Distances, shape (n, 3), the distances in Ångstrom to the three closest sites, irrespective of whether those sites are occupied, padded if the reference structure has fewer than three sites,
IndexMapping, shape (n,), the index of the atom of the input structure that occupies the site.

At a vacant site the three arrays of distances are undefined and IndexMapping is -1.

Notes

Only structure and reference can be given positionally. Everything else has to be given by keyword.

The input structure may be rotated relative to the reference structure. The transformation matrix relating the two is then determined by a search, which is more expensive than the rounding that suffices when the two structures adhere to the same crystallographic setting. The returned structure is always expressed in the setting of the reference structure.

A rigid translation of the input structure relative to the reference structure is detected and removed, unless find_translation is set to False.

The reference structure has to be periodic along at least two directions. A slab is therefore mapped, with the amount and the direction of its vacuum having no bearing on the outcome, whereas a wire or a molecule is refused: the cell of a structure periodic along one direction is unchanged by a rotation about that direction, and a structure periodic along none has no cell vector that means anything, so in neither case does the cell determine the orientation.

Vacancies. The reference structure provides the sites and the input structure provides the atoms. If there are fewer atoms than sites, the sites that no atom is assigned to are vacant and come out carrying the species X, which is how a vacancy is denoted in icet. Which sites end up vacant is decided by the assignment and is not under the control of the caller: leaving a site vacant carries no direct cost for the underlying Hungarian algorithm, so the vacancies are placed wherever they leave the remaining atoms closest to their sites. A vacant site carries no displacement, so Displacement, Displacement_Magnitude and Minimum_Distances are undefined there and IndexMapping is -1; accordingly drmax and dravg are taken over the atoms that are present rather than over the sites. This is also the reason for inert_species: an input structure with vacancies has fewer atoms than the reference structure has sites, so rescaling the volume by the total number of atoms would be wrong and only the species that are never substituted for a vacancy may be counted.

Interstitials. Interstitials are handled by making the interstitial sites part of the reference structure. The sites that are not occupied in a particular input structure then come out as vacancies, which is how, for example, hydrogen in a metal is treated. The converse is not supported: an input structure cannot contain atoms for which the reference structure provides no site, and attempting it raises a ValueError. There is no mechanism for introducing a site that the reference structure does not define.

Parameters:

structure (Atoms) – Input structure, typically a relaxed structure.
reference (Atoms) – Reference structure, which can but need not be the primitive structure.
inert_species (list[str] | None) – List of chemical symbols (e.g., ['Au', 'Pd']) that are never substituted for a vacancy. The number of inert sites is used to rescale the volume of the input structure to match the reference structure.
tol_positions (float) – Tolerance factor applied when scanning for overlapping positions in Ångstrom (forwarded to ase.build.make_supercell()).
tol_cell (float) – Largest acceptable deviation of the cell metric of the input structure from an integer transformation of the reference cell metric, measured as the largest absolute eigenvalue of the residual Biot strain tensor. This measure is invariant under a rotation of either cell and independent of the size of the cells. For a structure that is not periodic along every direction it covers the subspace spanned by the periodic cell vectors alone, since the extent of the cell along the other directions is arbitrary. A strongly strained input structure may require a larger value.
symprec (float) – Tolerance imposed when analyzing the symmetry of the reference structure using spglib. The symmetry is used to discard candidate transformation matrices that describe the same mapping.
find_translation (bool) – If True a rigid offset of the input structure relative to the reference structure is removed before the atoms are assigned to the sites. The offset is only removed if doing so reduces the average relaxation distance by at least a factor of two, so that a structure that is already aligned is left alone. The offset that was removed is reported as translation.
trigger_levels (dict[str, float] | None) –
Levels at which the warnings are triggered, for those that are to deviate from the defaults. Levels that are not given keep their default value, and an unknown key raises a ValueError that lists the valid ones. The valid keys and their default values are
- volumetric_strain (0.25), the absolute volumetric strain,
- anisotropic_strain (0.1), the difference between the largest and the smallest eigenvalue of the strain tensor,
- maximum_displacement (1.0 Å), the largest distance between a relaxed position and its site,
- average_displacement (0.5 Å), the average distance between the relaxed positions and their sites,
- maximum_displacement_fraction (0.3), the same as a fraction of the distance between two neighbouring sites, since whether a displacement is large depends on the lattice,
- average_displacement_fraction (0.15), likewise for the average,
- ambiguity_tolerance (1e-6 Å), the slack allowed when deciding whether an atom was assigned to a site further away than the closest one,
- ambiguity_ratio (0.9), the fraction of the distance to the second closest site beyond which an atom counts as almost equally far from both.
Setting assume_no_cell_relaxation to True tightens the defaults for volumetric_strain and anisotropic_strain to 1e-3, which an entry given here then overrides in turn.
suppress_warnings (bool) – If True no warnings are printed. They are still reported in the supplementary information, so that they can be inspected programmatically.
assume_no_cell_relaxation (bool) – If True the volume and cell metric of the input structure are not rescaled to match the reference structure. Skipping the rescaling can be advantageous for some structures, e.g., with many vacancies. Note that the input structure must then be obtainable via an integer transformation matrix from the reference cell metric, i.e. it should not involve relaxations of the volume or of the cell metric.

Raises:

ValueError – If the boundary conditions of the two structures do not match.
ValueError – If either structure contains no atoms or has a cell of less than full rank.
ValueError – If a species given as inert_species occupies no site of one of them.
ValueError – If no integer transformation of the reference cell reproduces the input cell.
ValueError – If the input structure has more atoms than the reference structure has sites.
ValueError – If trigger_levels contains a key that is not a known trigger level.

Return type:

tuple[Atoms, StructureMapping]

Example

The following code snippet illustrates the general usage. It first creates a primitive FCC cell, which is latter used as reference structure. To emulate a relaxed structure obtained from, e.g., a density functional theory calculation, the code then creates a 4x4x4 conventional FCC supercell, which is populated with two different atom types, has distorted cell vectors, and random displacements to the atoms. Finally, the present function is used to map the structure back the ideal lattice:

>>> from ase.build import bulk
>>> reference = bulk('Au', a=4.09)
>>> structure = bulk('Au', cubic=True, a=4.09).repeat(4)
>>> structure.symbols = 10 * ['Ag'] + (len(structure) - 10) * ['Au']
>>> structure.set_cell(structure.cell * 1.02, scale_atoms=True)
>>> structure.rattle(0.1, seed=42)
>>> mapped_structure, info = map_structure_to_reference(structure, reference)
>>> print('{:.4f}'.format(info['dravg']))
0.1541

class icet.tools.StructureMapping(drmax, dravg, transformation_matrix, translation, ambiguous_sites, warnings, strain_tensor, strain_tensor_eigenvalues, volumetric_strain, volume_dilation, isochoric_strain, isochoric_strain_eigenvalues, eigenstrain_norm, eigenstrain_rms, von_mises_strain, disorientation_angle, disorientation_axis, rotation)[source]#

Supplementary information about a mapping of a structure onto a reference structure, as returned by map_structure_to_reference.

The quantities are reached as attributes, which is the intended form. The class is at the same time a read-only mapping, so result['drmax'], result.get('drmax'), result.keys(), result.items(), result.values(), 'drmax' in result and len(result) all behave as they do for a dictionary. That is there so that code written against the dictionary that this class replaces continues to work.

Nothing can be assigned to, and neither can the contents of what is handed out: the arrays are not writeable and the sequences are tuples. Instances are therefore not hashable, as a mapping is not.

drmax#: Largest distance in Ångstrom between a position of the input structure and the site it was assigned to.

dravg#: Average of that distance over the atoms. Both are taken over the atoms that are present rather than over the sites, so vacancies do not dilute them.

transformation_matrix#: Integer matrix relating the cell of the reference structure to that of the returned supercell.

translation#: Rigid offset in Ångstrom that was removed from the input structure, which vanishes if none was removed.

ambiguous_sites#: Sites whose assignment is ambiguous, either because the atom was assigned past a closer site or because it is almost equally far from two sites.

warnings#: Tags of the warnings that were triggered, reported whether or not they were printed.

strain_tensor#: Biot strain tensor of the input structure relative to the returned supercell.

strain_tensor_eigenvalues#: Its three eigenvalues in ascending order.

volumetric_strain#: Relative change of the volume.

volume_dilation#: Ratio of the volumes, i.e. the volumetric strain plus one.

isochoric_strain#: Biot strain tensor of the volume-preserving part of the deformation, which describes the change of shape alone.

isochoric_strain_eigenvalues#: Its three eigenvalues in ascending order.

eigenstrain_norm#: Norm of those eigenvalues, a single measure of the change of shape.

eigenstrain_rms#: Root mean square of those eigenvalues.

von_mises_strain#: Von Mises strain of the isochoric strain, the conventional measure of the change of shape. For small strains it is the eigenstrain norm divided by the square root of three halves.

disorientation_angle#: Angle of the rotation in degrees, reduced over the symmetry of the reference structure. A supercell related to the returned one by a symmetry operation describes the same crystal, so this is the smallest angle among the descriptions that the structure cannot tell apart, and it is the rotation that is meaningful on its own.

disorientation_axis#: Unit vector in Cartesian coordinates that the disorientation is taken about. It belongs to the same description as disorientation_angle, so the two together give that rotation, and it is not the axis of rotation. It vanishes if the angle does, since a vanishing rotation has no axis, and it is not determined by the crystal if the misorientation lies on a symmetry element, in which case one of the equivalent axes is reported.

rotation#: Rotation relating the returned supercell to the cell of the input structure. It describes the cells and is defined only up to a symmetry operation of the reference structure, so its angle can be much larger than disorientation_angle; a misorientation of twenty degrees can appear here as one hundred and seventy two.

Structure enumeration#

icet.tools.enumerate_structures(structure, sizes, chemical_symbols, concentration_restrictions=None, niggli_reduce=None, symprec=1e-05, position_tolerance=None)[source]#

Yields a sequence of enumerated structures. The function generates all inequivalent structures that are permissible given a certain lattice. Using the chemical_symbols and concentration_restrictions keyword arguments it is possible to specify which chemical_symbols are to be included on which site and in which concentration range.

The function is sensitive to the boundary conditions of the input structure. An enumeration of, for example, a surface can thus be performed by setting structure.pbc = [True, True, False].

The algorithm implemented here was developed by Gus L. W. Hart and Rodney W. Forcade in Phys. Rev. B 77, 224115 (2008) [HarFor08] and Phys. Rev. B 80, 014120 (2009) [HarFor09].

Parameters:

structure (Atoms) – Primitive structure from which derivative superstructures should be generated.
sizes (list[int] | range) – Number of sites (included in enumeration).
chemical_symbols (list) – Chemical species with which to decorate the structure, e.g., ['Au', 'Ag']; see below for more examples.
concentration_restrictions (dict | None) – Allowed concentration range for one or more element in chemical_symbols, e.g., {'Au': (0, 0.2)} will only enumerate structures in which the Au content is between 0 and 20 %. Here, concentration is always defined as the number of atoms of the specified kind divided by the number of all atoms.
niggli_reduce (bool | None) – If True perform a Niggli reduction with spglib for each structure. The default is True if structure is periodic in all directions, False otherwise.
symprec (float) – Tolerance imposed when analyzing the symmetry using spglib.
position_tolerance (float | None) – Tolerance applied when comparing positions in Cartesian coordinates; by default this value is set equal to symprec.

Return type:

Atoms

Examples

The following code snippet illustrates how to enumerate structures with up to 6 atoms in the unit cell for a binary alloy without any constraints:

>>> from ase.build import bulk
>>> prim = bulk('Ag')
>>> for structure in enumerate_structures(structure=prim,
...                                       sizes=range(1, 5),
...                                       chemical_symbols=['Ag', 'Au']):
...     pass # Do something with the structure

To limit the concentration range to 10 to 40% Au the code should be modified as follows:

>>> conc_restr = {'Au': (0.1, 0.4)}
>>> for structure in enumerate_structures(
...         structure=prim, sizes=range(1, 5),
...         chemical_symbols=['Ag', 'Au'],
...         concentration_restrictions=conc_restr):
...     pass # Do something with the structure

Often one would like to consider mixing on only one sublattice. This can be achieved as illustrated for a Ga(1-x)Al(x)As alloy as follows:

>>> prim = bulk('GaAs', crystalstructure='zincblende', a=5.65)
>>> for structure in enumerate_structures(
...         structure=prim, sizes=range(1, 9),
...         chemical_symbols=[['Ga', 'Al'], ['As']]):
...     pass # Do something with the structure

icet.tools.enumerate_supercells(structure, sizes, niggli_reduce=None, symprec=1e-05, position_tolerance=None)[source]#

Yields a sequence of enumerated supercells. The function generates all inequivalent supercells that are permissible given a certain lattice. Any supercell can be reduced to one of the supercells generated.

The function is sensitive to the boundary conditions of the input structure. An enumeration of, for example, a surface can thus be performed by setting structure.pbc = [True, True, False].

The algorithm is based on Gus L. W. Hart and Rodney W. Forcade in Phys. Rev. B 77, 224115 (2008) [HarFor08] and Phys. Rev. B 80, 014120 (2009) [HarFor09].

Parameters:

structure (Atoms) – Primitive structure from which supercells should be generated.
sizes (list[int] | range) – Number of sites (included in enumeration).
niggli_reduce (bool | None) – If True perform a Niggli reduction with spglib for each supercell. The default is True if structure is periodic in all directions, False otherwise.
symprec (float) – Tolerance imposed when analyzing the symmetry using spglib.
position_tolerance (float | None) – Tolerance applied when comparing positions in Cartesian coordinates. By default this value is set equal to symprec.

Return type:

Atoms

Examples

The following code snippet illustrates how to enumerate supercells with up to 6 atoms in the unit cell:

>>> from ase.build import bulk
>>> prim = bulk('Ag')
>>> for supercell in enumerate_supercells(structure=prim, sizes=range(1, 7)):
...     pass # Do something with the supercell

Generation of training structures#

icet.tools.training_set_generation.structure_selection_annealing(cluster_space, monte_carlo_structures, n_structures_to_add, n_steps, base_structures=None, cooling_start=5, cooling_stop=0.001, cooling_function='exponential', initial_indices=None)[source]#

Given a cluster space, a base pool of structures, and a new pool of structures, this function uses a Monte Carlo inspired annealing method to find a good structure pool for training.

Return type:

tuple[list[int], list[float]]

Returns:

A tuple comprising the indices of the optimal structures in the monte_carlo_structures pool and a list of accepted metric values.

Parameters:

cluster_space (ClusterSpace) – A cluster space defining the lattice to be occupied.
monte_carlo_structures (list[Atoms]) – A list of candidate training structures.
n_structures_to_add (int) – How many of the structures in the monte_carlo_structures pool that should be kept for training.
n_steps (int) – Number of steps in the annealing algorithm.
base_structures (list[Atoms] | None) – A list of structures that is already in your training pool; can be None if you do not have any structures yet.
cooling_start (float) – Initial value of the cooling_function.
cooling_stop (float) – Last value of the cooling_function.
cooling_function (str | Callable) – Artificial number that rescales the difference between the metric value between two iterations. Available options are 'linear' and 'exponential'.
initial_indices (list[int] | None) – Picks out the starting structure from the monte_carlo_structures pool. Can be used if you want to continue from an old run for example.

Example

The following snippet demonstrates the use of this function for generating an optimized structure pool. Here, we first set up a pool of candidate structures by randomly occupying a FCC supercell with Au and Pd:

>>> from ase.build import bulk
>>> from icet.tools.structure_generation import occupy_structure_randomly

>>> prim = bulk('Au', a=4.0)
>>> cs = ClusterSpace(prim, [6.0], [['Au', 'Pd']])
>>> structure_pool = []
>>> for _ in range(500):
>>>     # Create random supercell.
>>>     supercell = np.random.randint(1, 4, size=3)
>>>     structure = prim.repeat(supercell)
>>>
>>>     # Randomize concentrations in the supercell
>>>     n_atoms = len(structure)
>>>     n_Au = np.random.randint(0, n_atoms)
>>>     n_Pd = n_atoms - n_Au
>>>     concentration = {'Au': n_Au / n_atoms, 'Pd': n_Pd / n_atoms}
>>>
>>>     # Occupy the structure randomly and store it.
>>>     occupy_structure_randomly(structure, cs, concentration)
>>>     structure_pool.append(structure)
>>> start_inds = [f for f in range(10)]

Now we can use the structure_selection_annealing() function to find an optimized structure pool:

>>> inds, cond = structure_selection_annealing(cs,
>>>                                            structure_pool,
>>>                                            n_structures_to_add=10,
>>>                                            n_steps=100)
>>> training_structures = [structure_pool[ind] for ind in inds]
>>> print(training_structures)

Generation of special structures#

icet.tools.structure_generation.generate_sqs(cluster_space, max_size, target_concentrations, include_smaller_cells=True, pbc=None, T_start=5.0, T_stop=0.001, n_steps=None, optimality_weight=1.0, random_seed=None, tol=1e-05)[source]#

Given a cluster_space, generate a special quasirandom structure (SQS), i.e., a structure that for a given supercell size provides the best possible approximation to a random alloy [ZunWeiFer90].

In the present case, this means that the generated structure will have a cluster vector that as closely as possible matches the cluster vector of an infintely large randomly occupied supercell. Internally the function uses a simulated annealing algorithm and the difference between two cluster vectors is calculated with the measure suggested by A. van de Walle et al. in Calphad 42, 13-18 (2013) [WalTiwJon13] (for more information, see mchammer.calculators.TargetVectorCalculator).

Parameters:

cluster_space (ClusterSpace) – Cluster space defining the lattice to be occupied.
max_size (int) – Maximum supercell size.
target_concentrations (dict) – Concentration of each species in the target structure, per sublattice (for example {'Au': 0.5, 'Pd': 0.5} for a single sublattice Au-Pd structure, or {'A': {'Au': 0.5, 'Pd': 0.5}, 'B': {'H': 0.25, 'X': 0.75}} for a system with two sublattices. The symbols defining sublattices (‘A’, ‘B’ etc) can be found by printing the cluster_space.
include_smaller_cells (bool) – If True, search among all supercell sizes including max_size, else search only among those exactly matching max_size
pbc (tuple[bool, bool, bool] | tuple[int, int, int] | None) – Periodic boundary conditions for each direction, e.g., (True, True, False). The axes are defined by the cell of cluster_space.primitive_structure. Default is periodic boundary in all directions.
T_start (float) – Artificial temperature at which the simulated annealing starts.
T_stop (float) – Artifical temperature at which the simulated annealing stops.
n_steps (int | None) – Total number of Monte Carlo steps in the simulation.
optimality_weight (float) – Controls weighting \(L\) of perfect correlations, see mchammer.calculators.TargetVectorCalculator.
random_seed (int | None) – Seed for the random number generator used in the Monte Carlo simulation.
tol (float) – Numerical tolerance.

Return type:

Atoms

icet.tools.structure_generation.generate_sqs_by_enumeration(cluster_space, max_size, target_concentrations, include_smaller_cells=True, pbc=None, optimality_weight=1.0, tol=1e-05)[source]#

Given a cluster_space, generate a special quasirandom structure (SQS), i.e., a structure that for a given supercell size provides the best possible approximation to a random alloy [ZunWeiFer90].

In the present case, this means that the generated structure will have a cluster vector that as closely as possible matches the cluster vector of an infintely large randomly occupied supercell. Internally the function uses a simulated annealing algorithm and the difference between two cluster vectors is calculated with the measure suggested by A. van de Walle et al. in Calphad 42, 13-18 (2013) [WalTiwJon13] (for more information, see mchammer.calculators.TargetVectorCalculator).

This functions generates SQS cells by exhaustive enumeration, which means that the generated SQS cell is guaranteed to be optimal with regard to the specified measure and cell size.

Parameters:

cluster_space (ClusterSpace) – Cluster space defining the lattice to be occupied.
max_size (int) – Maximum supercell size.
target_concentrations (dict) – Concentration of each species in the target structure, per sublattice (for example {'Au': 0.5, 'Pd': 0.5} for a single sublattice Au-Pd structure, or {'A': {'Au': 0.5, 'Pd': 0.5}, 'B': {'H': 0.25, 'X': 0.75}} for a system with two sublattices. The symbols defining sublattices (‘A’, ‘B’ etc) can be found by printing the cluster_space.
include_smaller_cells (bool) – if True search among all supercell sizes including max_size, else search only among those exactly matching max_size.
pbc (tuple[bool, bool, bool] | tuple[int, int, int] | None) – Periodic boundary conditions for each direction, e.g., (True, True, False). The axes are defined by the cell of cluster_space.primitive_structure. Default is periodic boundary in all directions.
optimality_weight (float) – Controls weighting \(L\) of perfect correlations, see mchammer.calculators.TargetVectorCalculator.
tol (float) – Numerical tolerance.

Return type:

Atoms

icet.tools.structure_generation.generate_sqs_from_supercells(cluster_space, supercells, target_concentrations, T_start=5.0, T_stop=0.001, n_steps=None, optimality_weight=1.0, random_seed=None, random_start=True, tol=1e-05)[source]#

Given a cluster_space` and one or more supercells, generate a special quasirandom structure (SQS), i.e., a structure that for the provided supercells size provides the best possible approximation to a random alloy [ZunWeiFer90].

In the present case, this means that the generated structure will have a cluster vector that as closely as possible matches the cluster vector of an infintely large randomly occupied supercell. Internally the function uses a simulated annealing algorithm and the difference between two cluster vectors is calculated with the measure suggested by A. van de Walle et al. in Calphad 42, 13-18 (2013) [WalTiwJon13] (for more information, see mchammer.calculators.TargetVectorCalculator).

Parameters:

cluster_space (ClusterSpace) – Cluster space defining the lattice to be occupied.
supercells (list[Atoms]) – List of one or more supercells among which an optimal structure will be searched for.
target_concentrations (dict) – Concentration of each species in the target structure, per sublattice (for example {'Au': 0.5, 'Pd': 0.5} for a single sublattice Au-Pd structure, or {'A': {'Au': 0.5, 'Pd': 0.5}, 'B': {'H': 0.25, 'X': 0.75}} for a system with two sublattices. The symbols defining sublattices (‘A’, ‘B’ etc) can be found by printing the cluster_space.
T_start (float) – Artificial temperature at which the simulated annealing starts.
T_stop (float) – Artifical temperature at which the simulated annealing stops.
n_steps (int | None) – Total number of Monte Carlo steps in the simulation.
optimality_weight (float) – Controls weighting \(L\) of perfect correlations, see mchammer.calculators.TargetVectorCalculator.
random_seed (int | None) – Seed for the random number generator used in the Monte Carlo simulation and used for initializing the occupation of the supercells if random_start is True.
random_start (bool) – Randomly occupy starting structure, can be disabled if the user prefers to pass an initial structure.
tol (float) – Numerical tolerance.

Return type:

Atoms

icet.tools.structure_generation.generate_target_structure(cluster_space, max_size, target_concentrations, target_cluster_vector, include_smaller_cells=True, pbc=None, T_start=5.0, T_stop=0.001, n_steps=None, optimality_weight=1.0, random_seed=None, tol=1e-05)[source]#

Given a cluster_space and a target_cluster_vector, generate a structure that as closely as possible matches that cluster vector. The search is performed among all inequivalent supercells shapes up to a certain size.

Internally the function uses a simulated annealing algorithm and the difference between two cluster vectors is calculated with the measure suggested by A. van de Walle et al. in Calphad 42, 13-18 (2013) [WalTiwJon13] (for more information, see mchammer.calculators.TargetVectorCalculator).

Parameters:

cluster_space (ClusterSpace) – Cluster space defining the lattice to be occupied.
max_size (int) – Maximum supercell size.
target_concentrations (dict) – Concentration of each species in the target structure, per sublattice (for example {'Au': 0.5, 'Pd': 0.5} for a single sublattice Au-Pd structure, or {'A': {'Au': 0.5, 'Pd': 0.5}, 'B': {'H': 0.25, 'X': 0.75}} for a system with two sublattices. The symbols defining sublattices (‘A’, ‘B’ etc) can be found by printing the cluster_space.
target_cluster_vector (list[float]) – Cluster vector that the generated structure should match as closely as possible.
include_smaller_cells (bool) – If True, search among all supercell sizes including max_size, else search only among those exactly matching max_size
pbc (tuple[bool, bool, bool] | tuple[int, int, int] | None) – Periodic boundary conditions for each direction, e.g., (True, True, False). The axes are defined by the cell of cluster_space.primitive_structure`. Default is periodic boundary in all directions.
T_start (float) – Artificial temperature at which the simulated annealing starts.
T_stop (float) – Artifical temperature at which the simulated annealing stops.
n_steps (int | None) – Total number of Monte Carlo steps in the simulation.
optimality_weight (float) – Controls weighting \(L\) of perfect correlations, see mchammer.calculators.TargetVectorCalculator.
random_seed (int | None) – Seed for the random number generator used in the Monte Carlo simulation.
tol (float) – Numerical tolerance.

Return type:

Atoms

icet.tools.structure_generation.generate_target_structure_from_supercells(cluster_space, supercells, target_concentrations, target_cluster_vector, T_start=5.0, T_stop=0.001, n_steps=None, optimality_weight=1.0, random_seed=None, random_start=True, tol=1e-05)[source]#

Given a cluster_space and a target_cluster_vector and one or more supercells, generate a structure that as closely as possible matches that cluster vector.

Internally the function uses a simulated annealing algorithm and the difference between two cluster vectors is calculated with the measure suggested by A. van de Walle et al. in Calphad 42, 13-18 (2013) [WalTiwJon13] (for more information, see mchammer.calculators.TargetVectorCalculator).

Parameters:

cluster_space (ClusterSpace) – A cluster space defining the lattice to be occupied.
supercells (list[Atoms]) – List of one or more supercells among which an optimal structure will be searched for.
target_concentrations (dict) – Concentration of each species in the target structure, per sublattice (for example {'Au': 0.5, 'Pd': 0.5} for a single sublattice Au-Pd structure, or {'A': {'Au': 0.5, 'Pd': 0.5}, 'B': {'H': 0.25, 'X': 0.75}} for a system with two sublattices. The symbols defining sublattices (‘A’, ‘B’ etc) can be found by printing the cluster_space.
target_cluster_vector (list[float]) – Cluster vector that the generated structure should match as closely as possible.
T_start (float) – Artificial temperature at which the simulated annealing starts.
T_stop (float) – Artifical temperature at which the simulated annealing stops.
n_steps (int | None) – Total number of Monte Carlo steps in the simulation.
optimality_weight (float) – Controls weighting \(L\) of perfect correlations, see mchammer.calculators.TargetVectorCalculator.
random_seed (int | None) – Seed for the random number generator used in the Monte Carlo simulation and used for initializing the occupation of the supercells if random_start is True.
random_start (bool) – Randomly occupy starting structure, can be disabled if the user prefers to pass an initial structure.
tol (float) – Numerical tolerance.

Return type:

Atoms

icet.tools.structure_generation.occupy_structure_randomly(structure, cluster_space, target_concentrations, random_seed=None)[source]#

Occupy a structure with quasirandom order but fulfilling target_concentrations.

Parameters:

structure (Atoms) – ASE Atoms object that will be occupied randomly.
cluster_space (ClusterSpace) – Cluster space (needed as it carries information about sublattices).
target_concentrations (dict) – Concentration of each species in the target structure, per sublattice (for example {'Au': 0.5, 'Pd': 0.5} for a single sublattice Au-Pd structure, or {'A': {'Au': 0.5, 'Pd': 0.5}, 'B': {'H': 0.25, 'X': 0.75}} for a system with two sublattices. The symbols defining sublattices (‘A’, ‘B’ etc) can be found by printing the cluster_space.
random_seed (int | None) – Seed for the random number generator.

Return type:

None

Ground state finder#

class icet.tools.ground_state_finder.GroundStateFinder(cluster_expansion, structure)[source]#

This class provides functionality for determining the ground states using a binary cluster expansion. This is efficiently achieved through the use of mixed integer programming (MIP) as developed by Larsen et al. in Phys. Rev. Lett. 120, 256101 (2018).

This class relies on the HiGHS package.

Note

The current implementation only works for binary systems.

Parameters:

cluster_expansion (ClusterExpansion) – Cluster expansion for which to find ground states.
structure (Atoms) – Atomic configuration.

Example

The following snippet illustrates how to determine the ground state for a Au-Ag alloy. Here, the parameters of the cluster expansion are set to emulate a simple Ising model in order to obtain an example that can be run without modification. In practice, one should of course use a proper cluster expansion:

>>> from ase.build import bulk
>>> from icet import ClusterExpansion, ClusterSpace

>>> # prepare cluster expansion
>>> # the setup emulates a second nearest-neighbor (NN) Ising model
>>> # (zerolet and singlet parameters are zero; only first and second
>>> # neighbor pairs are included)
>>> prim = bulk('Au')
>>> chemical_symbols = ['Ag', 'Au']
>>> cs = ClusterSpace(prim, cutoffs=[4.3], chemical_symbols=chemical_symbols)
>>> ce = ClusterExpansion(cs, [0, 0, 0.1, -0.02])

>>> # prepare initial configuration
>>> structure = prim.repeat(3)

>>> # set up the ground state finder and calculate the ground state energy
>>> gsf = GroundStateFinder(ce, structure)
>>> ground_state = gsf.get_ground_state({'Ag': 5})
>>> print('Ground state energy:', ce.predict(ground_state))

property constraints: dict[str, highs_cons]#: Dictionary with highs_cons objects and the names as keys.

get_ground_state(species_count=None, max_seconds=inf, threads=0)[source]#

Finds the ground state for a given structure and species count. If species_count is not provided when initializing the instance of this class the first species in the list of chemical symbols for the active sublattice will be used.

Parameters:

species_count (dict[str, int] | None) – Dictionary with count for one of the species on each active sublattice. If no count is provided for a sublattice, the concentration is allowed to vary.
max_seconds (float) – Maximum runtime in seconds.
threads (int) – Number of threads to be used when solving the problem, given that a positive integer has been provided. If set to \(0\) the solver default configuration is used while \(-1\) corresponds to all available processing cores.

Return type:

Atoms

property model: Highs#: HiGHS model.

property optimization_status: HighsModelStatus#: Optimization status.

Convex hull construction#

class icet.tools.ConvexHull(concentrations, energies)[source]#

This class provides functionality for extracting the convex hull of the (free) energy of mixing. It is based on the convex hull calculator in SciPy.

Only the lower part of the hull is retained, i.e., the part that minimizes rather than maximizes the energy. It is obtained directly from the facets of the hull, namely those whose outward normal points downward along the energy axis.

Parameters:

concentrations (list[float] | list[list[float]]) – Concentrations for each structure listed as [[c1, c2], [c1, c2], ...]; for binaries, in which case there is only one independent concentration, the format [c1, c2, c3, ...] works as well.
energies (list[float]) – Energy (or energy of mixing) for each structure.

concentrations#

Concentrations of the structures on the convex hull.

Type:: np.ndarray

energies#

Energies of the structures on the convex hull.

Type:: np.ndarray

dimensions#

Number of independent concentrations needed to specify a point in concentration space (1 for binaries, 2 for ternaries and so on).

Type:: int

structures#

Indices of structures that constitute the convex hull (indices are defined by the order of their concentrations and energies are fed when initializing the ConvexHull object).

Type:: list[int]

Examples

A ConvexHull object is easily initialized by providing lists of concentrations and energies:

>>> data = {'concentration': [0,    0.2,  0.2,  0.3,  0.4,  0.5,  0.8,  1.0],
...         'mixing_energy': [0.1, -0.2, -0.1, -0.2,  0.2, -0.4, -0.2, -0.1]}
>>> hull = ConvexHull(data['concentration'], data['mixing_energy'])

Now one can for example access the points along the convex hull directly:

>>> for c, e in zip(hull.concentrations, hull.energies):
...     print(c, e)
0.0 0.1
0.2 -0.2
0.5 -0.4
1.0 -0.1

or plot the convex hull along with the original data using e.g., matplotlib:

>>> import matplotlib.pyplot as plt
>>> plt.scatter(data['concentration'], data['mixing_energy'], color='darkred')
>>> plt.plot(hull.concentrations, hull.energies)
>>> plt.show(block=False)

It is also possible to extract structures at or close to the convex hull:

>>> low_energy_structures = hull.extract_low_energy_structures(
...     data['concentration'], data['mixing_energy'],
...     energy_tolerance=0.005)

A complete example can be found in the basic tutorial.

property concentration_labels: list[tuple[str, str]] | None#

The species that each independent concentration refers to, as (sublattice symbol, chemical symbol) pairs in the order of the columns of concentrations, of the target concentrations of the various methods, and of the columns returned by get_chemical_potentials().

This is None unless the object was created with from_sublattice_concentrations().

extract_low_energy_structures(concentrations, energies, energy_tolerance)[source]#

Returns the indices of energies that lie within a certain tolerance of the convex hull.

Parameters:

concentrations (list[float] | list[list[float]]) –
Concentrations of candidate structures.

If there is one independent concentration, a list of floats is sufficient. Otherwise, the concentrations must be provided as a list of lists, such as [[0.1, 0.2], [0.3, 0.1], ...].
energies (list[float]) – Energies of candidate structures.
energy_tolerance (float) – Include structures with an energy that is at most this far from the convex hull.

Return type:

list[int]

Returns:

The indices of the structures that lie within the tolerance of the
convex hull.

classmethod from_sublattice_concentrations(concentrations, energies, site_fractions=None, cluster_space=None)[source]#

Constructs a convex hull from concentrations that are resolved by sublattice.

Parameters:

concentrations (list[dict]) – One entry per structure, in the form used throughout icet, namely {'A': {'Ag': 0.3, 'Pd': 0.7}, 'B': {'H': 0.2, 'X': 0.8}}, where the keys of the outer dictionary are the symbols of the sublattices of a ClusterSpace and the concentrations sum to one on each sublattice. If there is only one sublattice the outer dictionary can be omitted, as in {'Ag': 0.3, 'Pd': 0.7}. Such concentrations are conveniently obtained with get_sublattice_concentrations().
energies (list[float]) – Energy (or energy of mixing) for each structure.
site_fractions (dict[str, float] | None) – Fraction of the sites of the lattice that each sublattice occupies, as in {'A': 0.5, 'B': 0.5}, conveniently obtained with get_sublattice_site_fractions(). This is needed by get_chemical_potentials() and by get_species_chemical_potentials(), which cannot recover it from the concentrations. The geometry of the hull, including get_facet_gradients() and get_decomposition(), does not use it.
cluster_space (ClusterSpace | None) – Cluster space that the structures belong to. When it is given, the concentrations are checked against it, so that a sublattice which allows more than one species and is left out is reported here instead of through a missing chemical potential during a simulation.

Return type:

ConvexHull

Returns:

A convex hull over the independent concentrations that the supplied
sublattices define.

Notes

A species whose concentration is the same in every structure carries no information and is left out, as is a sublattice all of whose species are. Of the species that remain on a sublattice, the alphabetically last is dropped, since the others determine it. The concentrations that result are ordered by sublattice symbol and then by chemical symbol, and concentration_labels records that order.

A ternary cluster space fed with data along one binary edge therefore gives concentrations for two species, not three, and a hull built while one sublattice is held fixed gives none for that sublattice.

All structures must share the same underlying lattice, so that the number of sites of each sublattice stands in the same ratio in all of them. This holds for structures that belong to one cluster space, and it is what makes the concentration of a mixture of phases the average of their concentrations on every sublattice at once.

Every sublattice that allows more than one species has to appear in the concentrations, at a fixed composition if it does not vary. A sublattice that is left out is neither described by the hull nor visible to it, and a simulation that samples it looks up a chemical potential that was never formed. get_sublattice_concentrations() covers every sublattice on its own, and passing cluster_space here checks it for concentrations that were assembled by hand.

Examples

>>> from icet.tools import ConvexHull
>>> concentrations = [{'A': {'Ag': 1.0, 'Pd': 0.0}, 'B': {'H': 0.0, 'X': 1.0}},
...                   {'A': {'Ag': 0.0, 'Pd': 1.0}, 'B': {'H': 0.0, 'X': 1.0}},
...                   {'A': {'Ag': 0.5, 'Pd': 0.5}, 'B': {'H': 0.0, 'X': 1.0}},
...                   {'A': {'Ag': 0.5, 'Pd': 0.5}, 'B': {'H': 1.0, 'X': 0.0}}]
>>> hull = ConvexHull.from_sublattice_concentrations(
...     concentrations, [0.0, 0.0, -0.1, 0.3])
>>> hull.concentration_labels
[('A', 'Ag'), ('B', 'H')]

get_chemical_potentials()[source]#

Returns the chemical potentials at which the structures on each face of the lower convex hull coexist, as an array of shape (n_faces, dimensions).

Entry \(i\) of a row is the difference between the chemical potential of the species that the \(i\)-th concentration counts and that of the species of the same sublattice whose concentration is fixed by the others. The rows are ordered as the faces returned by get_facets().

These differences are resolved by sublattice. A simulation in the semi-grand canonical ensemble takes one chemical potential per species instead, which get_species_chemical_potentials() assembles from them. The two carry the same information only when no species occurs on more than one sublattice.

A concentration counts the species of one sublattice relative to the sites of that sublattice, while the energy of a cluster expansion is given per site of the whole lattice. Changing the occupation of a single site therefore changes the concentration by \(1/N_s\) rather than by \(1/N\), and the chemical potential is the gradient of get_facet_gradients() divided by the fraction of the sites that the sublattice occupies.

For a hull that was built from an array of concentrations that fraction is taken to be one, i.e., the concentrations are assumed to refer to all sites. For a hull built with from_sublattice_concentrations() the fractions have to be supplied there, since they cannot be recovered from the concentrations.

Return type:

ndarray

Returns:

The chemical potential differences of each face, of shape
(n_faces, dimensions).

Examples

>>> from icet.tools import ConvexHull
>>> hull = ConvexHull([0.0, 0.5, 1.0], [0.0, -1.0, 0.0])
>>> hull.get_chemical_potentials()
array([[-2.],
       [ 2.]])

get_decomposition(target_concentrations)[source]#

Returns the decomposition of each target concentration into the structures on the convex hull.

At every concentration the hull is spanned by one of its facets. The structures that define this facet are the phases that a sample of the given overall concentration separates into, and the barycentric coordinates of the concentration within the facet are the fractions in which they occur.

Parameters:

target_concentrations (list[float] | list[list[float]]) –

Concentrations at target points.

If there is one independent concentration, a list of floats is sufficient. Otherwise, the concentrations ought to be provided as a list of lists, such as [[0.1, 0.2], [0.3, 0.1], ...].

Return type:

list[dict]

Returns:

One dictionary per target concentration with the keys structures,
holding the indices of the structures the concentration decomposes
into, fractions, holding the fraction of each of them, and
energy, holding the resulting energy. For a concentration outside
the allowed range, and for one that is not a number, the lists are
empty and the energy is NaN.
The fractions reproduce the target concentration. Right at the rim of
the sampled range a target can be accepted as lying inside while
falling outside every simplex by a small amount, in which case the
fractions describe the closest concentration that is inside.
A structure can carry a fraction of zero, which happens for a target on
the boundary between two simplices of one face. The structures reported
are therefore those of a simplex that contains the target, and which
simplex that is depends on the order of the input.

Examples

>>> from icet.tools import ConvexHull
>>> hull = ConvexHull([0.0, 0.5, 1.0], [0.0, -1.0, 0.0])
>>> decomposition = hull.get_decomposition([0.25])
>>> decomposition[0]['structures']
[0, 1]
>>> [round(f, 6) for f in decomposition[0]['fractions']]
[0.5, 0.5]
>>> round(decomposition[0]['energy'], 6)
-0.5

get_energy_above_convex_hull(concentrations, energies)[source]#

Returns the energy of each structure relative to the convex hull, commonly referred to as the energy above the hull. The value is zero for structures on the hull and positive otherwise. If a concentration lies outside the allowed range, NaN is returned.

Parameters:

concentrations (list[float] | list[list[float]]) –
Concentrations of the structures.

If there is one independent concentration, a list of floats is sufficient. Otherwise, the concentrations must be provided as a list of lists, such as [[0.1, 0.2], [0.3, 0.1], ...].
energies (list[float]) – Energies of the structures.

Return type:

ndarray

Returns:

The energy of each structure relative to the convex hull, and NaN
where the concentration lies outside the sampled range.

Examples

>>> from icet.tools import ConvexHull
>>> hull = ConvexHull([0.0, 0.5, 1.0], [0.0, -1.0, 0.0])
>>> hull.get_energy_above_convex_hull([0.25, 0.5], [0.0, -1.0])
array([0.5, 0. ])

get_energy_at_convex_hull(target_concentrations)[source]#

Returns the energy of the convex hull at specified concentrations. If any concentration is outside the allowed range, NaN is returned.

Parameters:

target_concentrations (list[float] | list[list[float]]) –

Concentrations at target points.

If there is one independent concentration, a list of floats is sufficient. Otherwise, the concentrations ought to be provided as a list of lists, such as [[0.1, 0.2], [0.3, 0.1], ...].

Return type:

ndarray

Returns:

The energy of the convex hull at each target concentration, and NaN
where the concentration lies outside the sampled range.

Examples

>>> from icet.tools import ConvexHull
>>> hull = ConvexHull([0.0, 0.5, 1.0], [0.0, -1.0, 0.0])
>>> hull.get_energy_at_convex_hull([0.25, 0.75])
array([-0.5, -0.5])

get_facet_gradients()[source]#

Returns the gradient of the energy along each face of the lower convex hull as an array of shape (n_faces, dimensions).

Entry \(i\) of a row is \(\partial e/\partial c_i\), the derivative of the energy with respect to the \(i\)-th independent concentration along that face. The rows are ordered as the faces returned by get_facets().

This is a property of the geometry of the hull alone. See get_chemical_potentials() for the related quantity that a semi-grand canonical simulation takes, which differs by the fraction of the sites that the concentration refers to.

Return type:

ndarray

Returns:

The gradient of the energy along each face, of shape
(n_faces, dimensions).

Examples

>>> from icet.tools import ConvexHull
>>> hull = ConvexHull([0.0, 0.5, 1.0], [0.0, -1.0, 0.0])
>>> hull.get_facet_gradients()
array([[-2.],
       [ 2.]])

get_facets()[source]#

Returns the faces of the lower convex hull, each as the list of indices of the structures that span it.

In a binary these are the tie lines, in a ternary the three-phase triangles, and so on. Together they enumerate the coexistence regions of the system.

A face can have more vertices than dimensions + 1, in which case the structures on it are degenerate in energy over the whole face. get_decomposition() still reports dimensions + 1 of them, since a given concentration separates into that many phases.

Return type:

list[list[int]]

Returns:

The indices of the structures that span each face of the lower convex
hull, one list per face.

Examples

>>> from icet.tools import ConvexHull
>>> hull = ConvexHull([0.0, 0.5, 1.0], [0.0, -1.0, 0.0])
>>> hull.get_facets()
[[0, 1], [1, 2]]

get_species_chemical_potentials(strict=True)[source]#

Returns the chemical potential of each species at which the structures on each face of the lower convex hull coexist, as one dictionary per face in the order of get_facets().

This is the chemical_potentials argument of SemiGrandCanonicalEnsemble. Setting it to the values of a face makes the structures on that face degenerate, so that the simulation sits at the corresponding phase boundary. Only differences between species that share a sublattice matter, since those are the exchanges a simulation makes, and the species that comes first alphabetically is given a value of zero. A difference between species of different sublattices is whatever the solution happened to give and carries no meaning.

A chemical potential belongs to a species and not to a sublattice, while get_chemical_potentials() gives one difference per concentration and therefore per sublattice. When a species occurs on more than one sublattice those differences can constrain the same pair of species twice, and there is then no set of chemical potentials that satisfies all of them at once. Such a face is not a state that a semi-grand canonical simulation can reach, whatever chemical potentials are chosen, and this method raises rather than return values that do not describe it.

The result covers the species whose concentration varies across the structures the hull was built from, and only those. A species held at one composition throughout, and every species of a sublattice held that way, carry no information for the hull to report and are absent from the dictionary.

A simulation therefore has to leave them where the hull assumed them. SemiGrandCanonicalEnsemble holds a whole sublattice fixed through its sublattice_probabilities argument, and HybridEnsemble holds individual species fixed through the allowed_symbols of a semi-grand canonical step. A simulation that samples a species the dictionary does not name changes a composition the hull took as given, and no longer sits at the phase boundary it reports.

Parameters:

strict (bool) – Raise for a face that no set of chemical potentials describes. When it is False such a face yields None in its place and the remaining faces are returned as usual.

Return type:

list[dict[str, float] | None]

Returns:

The chemical potential of each species at which the structures on each
face coexist, as one dictionary per face in the order of
get_facets(), with None for a face that no set of chemical
potentials describes when strict is False.

Raises:

ValueError – If the hull was not built from concentrations that are resolved by sublattice, if the site fractions were not supplied, or if a face cannot be described by one chemical potential per species while strict is True.

Examples

>>> from icet.tools import ConvexHull
>>> concentrations = [{'A': {'Ag': 1.0, 'Pd': 0.0}},
...                   {'A': {'Ag': 0.0, 'Pd': 1.0}},
...                   {'A': {'Ag': 0.5, 'Pd': 0.5}}]
>>> hull = ConvexHull.from_sublattice_concentrations(
...     concentrations, [0.0, 0.0, -1.0], site_fractions={'A': 1.0})
>>> [{k: round(v, 6) for k, v in mu.items()}
...  for mu in hull.get_species_chemical_potentials()]
[{'Ag': 0.0, 'Pd': -2.0}, {'Ag': 0.0, 'Pd': 2.0}]

is_on_convex_hull(concentrations, energies, energy_tolerance=1e-06)[source]#

Returns a boolean array that is True for every structure that lies on the convex hull, which is convenient for filtering a table of results.

A structure whose energy lies below the hull is not on it, and the result is False. That situation means that the hull was constructed without that structure and no longer describes the system. Use extract_low_energy_structures() to collect the structures that are at or below the hull, which is what a search for new ground states needs.

Parameters:

concentrations (list[float] | list[list[float]]) –
Concentrations of the structures.

If there is one independent concentration, a list of floats is sufficient. Otherwise, the concentrations must be provided as a list of lists, such as [[0.1, 0.2], [0.3, 0.1], ...].
energies (list[float]) – Energies of the structures.
energy_tolerance (float) – Consider a structure to be on the hull if its energy is at most this far above it.

Return type:

ndarray

Returns:

A boolean for each structure that is True where it lies on the
convex hull.

Examples

>>> from icet.tools import ConvexHull
>>> hull = ConvexHull([0.0, 0.5, 1.0], [0.0, -1.0, 0.0])
>>> hull.is_on_convex_hull([0.25, 0.5], [0.0, -1.0])
array([False,  True])

icet.tools.get_sublattice_concentrations(structures, cluster_space)[source]#

Returns the concentrations of each structure resolved by sublattice, in the form expected by ConvexHull.from_sublattice_concentrations.

Parameters:

structures (list[Atoms]) – Atomic configurations, each of which must be a supercell of the primitive structure of the cluster space.
cluster_space (ClusterSpace) – Cluster space that defines the sublattices and the species allowed on them.

Return type:

list[dict]

Returns:

One dictionary per structure of the form ``{‘A’ ({‘Ag’: 0.3, ‘Pd’: 0.7},)
’B’ ({‘H’: 0.2, ‘X’: 0.8}}``. Every species that the cluster space allows)
on a sublattice appears, including those that the structure does not
contain.

Examples

>>> from ase.build import bulk
>>> from icet import ClusterSpace
>>> from icet.tools import get_sublattice_concentrations
>>> primitive_structure = bulk('Ag', 'fcc', a=4.0)
>>> cluster_space = ClusterSpace(primitive_structure, [4.0], ['Ag', 'Pd'])
>>> structure = primitive_structure.repeat(2)
>>> structure.symbols[:4] = 'Pd'
>>> get_sublattice_concentrations([structure], cluster_space)
[{'A': {'Ag': 0.5, 'Pd': 0.5}}]

icet.tools.get_sublattice_site_fractions(structure, cluster_space)[source]#

Returns the fraction of the sites of the lattice that each sublattice occupies, in the form expected by ConvexHull.from_sublattice_concentrations.

A concentration counts the species of one sublattice relative to the sites of that sublattice, while the energy of a cluster expansion is given per site of the whole lattice. These fractions relate the two and are what turns the gradient of the energy along a face of the hull into a chemical potential.

Parameters:

structure (Atoms) – Atomic configuration, which must be a supercell of the primitive structure of the cluster space. The fractions are a property of the lattice rather than of the occupation, so any such supercell will do.
cluster_space (ClusterSpace) – Cluster space that defines the sublattices.

Return type:

dict[str, float]

Returns:

The fraction of the sites that each sublattice occupies, keyed by the
symbol of the sublattice.

Examples

>>> from ase.build import bulk
>>> from icet import ClusterSpace
>>> from icet.tools import get_sublattice_site_fractions
>>> primitive_structure = bulk('NaCl', 'rocksalt', a=5.6)
>>> cluster_space = ClusterSpace(primitive_structure, [6.0],
...                              [['Na', 'K'], ['Cl', 'X']])
>>> get_sublattice_site_fractions(primitive_structure, cluster_space)
{'A': 0.5, 'B': 0.5}

Fitting with constraints#

class icet.tools.constraints.Constraints(n_params)[source]#

Class for handling linear constraints with right-hand-side equal to zero.

Parameters:: n_params (int) – Number of parameters in model.

Example

The following example demonstrates fitting of a cluster expansion under the constraint that parameter 2 and parameter 4 should be equal:

>>> import numpy as np
>>> from icet.tools import Constraints
>>> from trainstation import Optimizer

>>> # Set up random sensing matrix and target "energies"
>>> n_params = 10
>>> n_energies = 20
>>> A = np.random.random((n_energies, n_params))
>>> y = np.random.random(n_energies)

>>> # Define constraints
>>> c = Constraints(n_params=n_params)
>>> M = np.zeros((1, n_params))
>>> M[0, 2] = 1
>>> M[0, 4] = -1
>>> c.add_constraint(M)

>>> # Do the actual fit and finally extract parameters
>>> A_constrained = c.transform(A)
>>> opt = Optimizer((A_constrained, y), fit_method='ridge')
>>> opt.train()
>>> parameters = c.inverse_transform(opt.parameters)

add_constraint(M)[source]#

Add a constraint matrix and resolve for the constraint space.

Parameters:: M (ndarray) – Constraint matrix with each constraint as a row. Can (but need not be) cluster vectors.
Return type:: None

inverse_transform(A)[source]#

Inverse transform array from constrained parameter space to unconstrained space.

Parameters:: A (ndarray) – Array to be inverse transformed.
Return type:: ndarray

transform(A)[source]#

Transform array to constrained parameter space.

Parameters:: A (ndarray) – Array to be transformed.
Return type:: ndarray

icet.tools.constraints.get_mixing_energy_constraints(cluster_space)[source]#

A cluster expansion of the mixing energy should ideally predict zero energy for concentrations 0 and 1. This function constructs a Constraints object that enforces that condition during training.

Parameters:: cluster_space – Cluster space corresponding to cluster expansion for which constraints should be imposed.
Return type:: Constraints

Example

This example demonstrates how to constrain the mixing energy to zero at the pure phases in a toy example with random cluster vectors and random target energies:

>>> import numpy as np
>>> from ase.build import bulk
>>> from icet import ClusterSpace
>>> from icet.tools import get_mixing_energy_constraints
>>> from trainstation import Optimizer

>>> # Set up cluster space along with random sensing matrix
>>> # and target "energies"
>>> prim = bulk('Au')
>>> cs = ClusterSpace(prim, cutoffs=[6.0, 5.0],
...                   chemical_symbols=['Au', 'Ag'])
>>> n_params = len(cs)
>>> n_energies = 20
>>> A = np.random.random((n_energies, n_params))
>>> y = np.random.random(n_energies)

>>> # Define constraints
>>> c = get_mixing_energy_constraints(cs)

>>> # Do the actual fit and finally extract parameters
>>> A_constrained = c.transform(A)
>>> opt = Optimizer((A_constrained, y), fit_method='ridge')
>>> opt.train()
>>> parameters = c.inverse_transform(opt.parameters)

Warning

Constraining the energy of one structure is always done at the expense of the fit quality of the others. Always expect that your cross-validation scores will increase somewhat when using this function.

Constituent strain#

class icet.tools.ConstituentStrain(supercell, primitive_structure, chemical_symbols, concentration_symbol, strain_energy_function, k_to_parameter_function=None, damping=1.0, tol=1e-06)[source]#

Class for handling constituent strain in cluster expansions (see Laks et al., Phys. Rev. B 46, 12587 (1992) [LakFerFro92]). This makes it possible to use cluster expansions to describe systems with strain due to, for example, coherent phase separation. For an extensive example on how to use this module, please see this example.

Parameters:

supercell (Atoms) – Defines supercell that will be used when calculating constituent strain.
primitive_structure (Atoms) – Primitive structure the supercell is based on.
chemical_symbols (list[str]) – List with chemical symbols involved, such as ['Ag', 'Cu'].
concentration_symbol (str) – Chemical symbol used to define concentration, such as 'Ag'.
strain_energy_function (Callable[[float, list[float]], float]) – A function that takes two arguments, a list of parameters and concentration (e.g., [0.5, 0.5, 0.5] and 0.3), and returns the corresponding strain energy. The parameters are in turn determined by k_to_parameter_function (see below). If k_to_parameter_function is None, the parameters list will be the k-point. For more information, see this example.
k_to_parameter_function (Callable[[list[float]], list[float]] | None) – A function that takes a k-point as a list of three floats and returns a parameter vector that will be fed into strain_energy_function (see above). If None, the k-point itself will be the parameter vector to strain_energy_function. The purpose of this function is to be able to precompute any factor in the strain energy that depends on the k-point but not the concentration. For more information, see this example.
damping (float) – Damping factor \(\eta\) used to suppress impact of large-magnitude k-points by multiplying strain with \(\exp(-(\eta \mathbf{k})^2)\) (unit Ångstrom).
tol (float) – Numerical tolerance when comparing k-points (units of inverse Ångstrom).

accept_change()[source]#

Update structure factor for each kpoint to the value in structure_factor_after. This makes it possible to efficiently calculate changes in constituent strain with the get_constituent_strain_change() function; this function should be called if the last occupations used to call get_constituent_strain_change() should be the starting point for the next call of get_constituent_strain_change(). This is taken care of automatically by the Monte Carlo simulations in mchammer.

Return type:: None

get_concentration(occupations)[source]#

Calculate current concentration.

Parameters:: occupations (ndarray) – Current occupations.
Return type:: float

get_constituent_strain(occupations, update_structure_factors=True)[source]#

Calculate total constituent strain.

Parameters:

occupations (list[int]) – Current occupations.
update_structure_factors (bool) – If True the structure factor stored for each k-point is replaced by the one for occupations, which is what get_constituent_strain_change() reads as the value before a change. Pass False to evaluate a configuration without moving that state, which a caller scoring a hypothetical configuration, or observing one alongside a simulation, has to do.

Return type:

float

get_constituent_strain_change(occupations, atom_index)[source]#

Calculate change in constituent strain upon change of the occupation of one site.

Warning

This function is dependent on the internal state of the ConstituentStrain object and should typically only be used internally by mchammer. Specifically, the structure factor is saved internally to speed up computation. The first time this function is called, occupations must be the same array as was used to initialize the ConstituentStrain object, or the same as was last used when get_constituent_strain() was called. After the present function has been called, the same occupations vector need to be used the next time as well, unless accept_change() has been called, in which case occupations should incorporate the changes implied by the previous call to the function.

Parameters:

occupations (ndarray) – Occupations before change.
atom_index (int) – Index of site the occupation of which is to be changed.

Return type:

float

Other structure tools#

icet.tools.get_primitive_structure(structure, no_idealize=True, to_primitive=True, symprec=1e-05)[source]#

Returns the primitive structure using spglib.

Parameters:

structure (Atoms) – Input atomic structure.
no_idealize (bool) – If True lengths and angles are not idealized.
to_primitive (bool) – If True convert to primitive structure.
symprec (float) – Tolerance imposed when analyzing the symmetry using spglib.

Return type:

Atoms

icet.tools.get_wyckoff_sites(structure, map_occupations=None, symprec=1e-05, include_representative_atom_index=False)[source]#

Returns the Wyckoff symbols of the input structure. The Wyckoff sites are of general interest for symmetry analysis but can be especially useful when setting up, e.g., a SiteOccupancyObserver. The Wyckoff labels can be conveniently attached as an array to the structure object as demonstrated in the examples section below.

By default the occupation of the sites is part of the symmetry analysis. If a chemically disordered structure is provided this will usually reduce the symmetry substantially. If one is interested in the symmetry of the underlying structure one can control how occupations are handled. To this end, one can provide the map_occupations keyword argument. The latter must be a list, each entry of which is a list of species that should be treated as indistinguishable. As a shortcut, if all species should be treated as indistinguishable one can provide an empty list. Examples that illustrate the usage of the keyword are given below.

Parameters:

structure (Atoms) – Input structure. Note that the occupation of the sites is included in the symmetry analysis.
map_occupations (list[list[str]] | None) – Each sublist in this list specifies a group of chemical species that shall be treated as indistinguishable for the purpose of the symmetry analysis.
symprec (float) – Tolerance imposed when analyzing the symmetry using spglib.
include_representative_atom_index (bool) – If True the index of the first atom in the structure that is representative of the Wyckoff site is included in the symbol. This is in particular useful in cases when there are multiple Wyckoff sites sites with the same Wyckoff letter.

Return type:

list[str]

Examples

Wyckoff sites of a hexagonal-close packed structure:

>>> from ase.build import bulk
>>> structure = bulk('Ti')
>>> wyckoff_sites = get_wyckoff_sites(structure)
>>> print(wyckoff_sites)
['2d', '2d']

The Wyckoff labels can also be attached as an array to the structure, in which case the information is also included when storing the Atoms object:

>>> from ase.io import write
>>> structure.new_array('wyckoff_sites', wyckoff_sites, str)
>>> write('structure.xyz', structure)  # xdoctest: +SKIP

The function can also be applied to supercells:

>>> structure = bulk('GaAs', crystalstructure='zincblende', a=3.0).repeat(2)
>>> wyckoff_sites = get_wyckoff_sites(structure)
>>> print(wyckoff_sites)
['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c',
 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c']

Now assume that one is given a supercell of a (Ga,Al)As alloy. Applying the function directly yields much lower symmetry since the symmetry of the original structure is broken:

>>> structure.set_chemical_symbols(
...        ['Ga', 'As', 'Al', 'As', 'Ga', 'As', 'Al', 'As',
...         'Ga', 'As', 'Ga', 'As', 'Al', 'As', 'Ga', 'As'])
>>> print(get_wyckoff_sites(structure))
['8g', '8i', '4e', '8i', '8g', '8i', '2c', '8i',
 '2d', '8i', '8g', '8i', '4e', '8i', '8g', '8i']

Since Ga and Al occupy the same sublattice, they should, however, be treated as indistinguishable for the purpose of the symmetry analysis, which can be achieved via the map_occupations keyword:

>>> print(get_wyckoff_sites(structure,
...       map_occupations=[['Ga', 'Al'], ['As']]))
['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c',
 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c']

If occupations are to ignored entirely, one can simply provide an empty list. In the present case, this turns the zincblende lattice into a diamond lattice, on which case there is only one Wyckoff site:

>>> print(get_wyckoff_sites(structure, map_occupations=[]))
['8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a',
 '8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a']