import numpy as np
from ase import Atoms
from icet.input_output.logging_tools import logger
from .structure_mapping_support.alignment import (DEFAULT_TOL_CELL, DEFAULT_TOL_POSITIONS,
get_reference_supercell)
from .structure_mapping_support.candidate_search import get_minimum_site_distance
from .structure_mapping_support.deformation import CellDeformation
from .structure_mapping_support.diagnostics import (DEFAULT_TRIGGER_LEVELS,
DEFAULT_TRIGGER_LEVELS_NO_CELL_RELAXATION,
collect_warnings,
get_displacement_metrics,
get_trigger_levels)
from .structure_mapping_support.mapping_result import StructureMapping
from .structure_mapping_support.orientation import (calculate_disorientation,
get_point_group_operations)
from .structure_mapping_support.site_matching import match_positions
from .structure_mapping_support.strain import (calculate_strain_eigenvalues,
calculate_strain_tensor,
calculate_volumetric_strain,
get_periodic_frame,
transfer_positions)
__all__ = ['DEFAULT_TRIGGER_LEVELS', 'DEFAULT_TRIGGER_LEVELS_NO_CELL_RELAXATION',
'StructureMapping', 'map_structure_to_reference']
[docs]
def map_structure_to_reference(structure: Atoms,
reference: Atoms,
*,
inert_species: list[str] | None = None,
tol_positions: float = DEFAULT_TOL_POSITIONS,
tol_cell: float = DEFAULT_TOL_CELL,
symprec: float = 1e-5,
find_translation: bool = True,
trigger_levels: dict[str, float] | None = None,
suppress_warnings: bool = False,
assume_no_cell_relaxation: bool = False) \
-> tuple[Atoms, StructureMapping]:
"""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 :class:`StructureMapping
<icet.tools.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 instead of 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 :attr:`structure` and :attr:`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 :attr:`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 :program:`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 instead of over the sites.
This is also the reason for :attr:`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
:class:`ValueError`.
There is no mechanism for introducing a site that the reference structure
does not define.
Parameters
----------
structure
Input structure, typically a relaxed structure.
reference
Reference structure, which can but need not be the primitive structure.
inert_species
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
Tolerance factor applied when scanning for overlapping positions in
Ångstrom (forwarded to :func:`ase.build.make_supercell`).
tol_cell
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
Tolerance imposed when analyzing the symmetry of the reference
structure using :program:`spglib`.
The symmetry is used to discard candidate transformation matrices that
describe the same mapping.
By default 1e-5.
find_translation
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 improves the agreement with
the sites, so that a structure that is already aligned is left alone.
The offset that was removed is reported as ``translation``.
By default the offset is removed.
trigger_levels
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 :class:`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 :attr:`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
If ``True`` no warnings are printed.
They are still reported in the supplementary information, so that they
can be inspected programmatically.
By default the warnings are printed.
assume_no_cell_relaxation
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.
By default the volume and the cell metric are rescaled.
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 :attr:`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 :attr:`trigger_levels` contains a key that is not a known trigger level.
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
"""
# Obtain supercell of reference structure that is compatible
# with relaxed structure
reference_supercell, P, reduced_translation = get_reference_supercell(
structure, reference, inert_species=inert_species,
tol_positions=tol_positions, tol_cell=tol_cell, symprec=symprec,
assume_no_cell_relaxation=assume_no_cell_relaxation)
# The strain and its decomposition describe the periodic subspace, the cell vectors
# along the other directions being replaced by an orthonormal basis of its complement.
# Reporting them for the cells as they stand would let the amount and the direction of
# vacuum enter, and a slab whose cell is skewed would appear strained by that alone.
reference_frame = get_periodic_frame(reference_supercell.cell, reference.pbc)
structure_frame = get_periodic_frame(structure.cell, reference.pbc)
# Calculate strain tensor
strain_tensor = calculate_strain_tensor(reference_frame, structure_frame)
# Symmetric matrix has real eigenvalues
eigenvalues = calculate_strain_eigenvalues(strain_tensor)
volumetric_strain = calculate_volumetric_strain(strain_tensor)
logger.debug('Eigenvalues of the strain tensor: {}, volumetric strain: {:.6f}'.format(
np.round(eigenvalues, 6), volumetric_strain))
# Separate the deformation into its change of volume, change of shape and rotation.
# The rotation is reduced over the symmetry of the reference structure, since a
# supercell related to the chosen one by a symmetry operation describes the same
# crystal and would otherwise give a different angle.
deformation = CellDeformation(reference_frame, structure_frame)
# The rotation comes from the periodic cell vectors rather than from the frames above.
# A frame completes the periodic vectors with a basis of what is left, and that basis is
# only determined where one direction is left over, so for a structure periodic along a
# single direction the frames of two cells can differ by a rotation that neither structure
# has undergone, or agree although one has been turned.
operations = get_point_group_operations(reference, symprec)
disorientation_angle, disorientation_axis = calculate_disorientation(
deformation.rotation, operations)
logger.debug('Volume dilation {:.6f}, isochoric eigenstrains {}, disorientation {:.4f} '
'degrees about {} (from {:.4f} degrees before reducing over {} '
'operations)'.format(
deformation.volume_dilation,
np.round(deformation.isochoric_strain_eigenvalues, 6),
disorientation_angle, np.round(disorientation_axis, 4),
deformation.rotation_angle, len(operations)))
# Express the input structure in the cell of the supercell it is mapped onto. The
# positions follow the periodic cell vectors and keep their Cartesian component along the
# remaining directions, since the extent of the cell along those is arbitrary and scaling
# by it would move the atoms.
structure_scaled = structure.copy()
structure_scaled.positions = transfer_positions(
structure.positions, structure_frame, reference_frame)
structure_scaled.set_cell(reference_supercell.cell, scale_atoms=False)
# Remove a rigid offset of the input structure relative to the reference
# structure. The offset is only non-zero if it improves the agreement with the
# reference sites, so that a structure that is already aligned is left alone
# rather than being displaced by a spurious offset.
translation = np.zeros(3)
if find_translation and np.abs(reduced_translation).max() > 1e-12:
# the offset comes back in the coordinates it was sought in, which are those of the
# frame of the reference structure rather than of its cell
translation = np.dot(reduced_translation,
get_periodic_frame(reference.cell, reference.pbc))
logger.debug('Removing a rigid offset of {} Angstrom, which is {} in reduced '
'coordinates of the reference cell'.format(np.round(translation, 5),
np.round(reduced_translation, 5)))
structure_scaled.positions -= translation
structure_scaled.wrap()
# Match positions
mapped_structure = match_positions(structure_scaled, reference_supercell)
drmax, dravg = get_displacement_metrics(mapped_structure)
# Collect the warnings unconditionally, so that they are reported in the
# supplementary information even when they are not printed
levels = get_trigger_levels(assume_no_cell_relaxation, trigger_levels)
warnings, ambiguous_sites = collect_warnings(mapped_structure, volumetric_strain,
eigenvalues, drmax, dravg, levels,
get_minimum_site_distance(reference))
n_vacancies = mapped_structure.get_chemical_symbols().count('X')
logger.debug('Assigned {} atoms to {} sites, leaving {} vacant; maximum relaxation '
'distance {:.5f} Angstrom, average {:.5f} Angstrom, {} ambiguous '
'site(s)'.format(len(structure), len(mapped_structure), n_vacancies,
drmax, dravg, len(ambiguous_sites)))
warning_tags = []
for tag, message in warnings:
if tag not in warning_tags:
warning_tags.append(tag)
if not suppress_warnings:
logger.warning(message)
# Populate dictionary with supplementary information
result = StructureMapping(
drmax=drmax,
dravg=dravg,
transformation_matrix=P,
translation=translation,
ambiguous_sites=ambiguous_sites,
warnings=warning_tags,
strain_tensor=strain_tensor,
strain_tensor_eigenvalues=eigenvalues,
volumetric_strain=volumetric_strain,
volume_dilation=deformation.volume_dilation,
isochoric_strain=deformation.isochoric_strain,
isochoric_strain_eigenvalues=deformation.isochoric_strain_eigenvalues,
eigenstrain_norm=deformation.eigenstrain_norm,
eigenstrain_rms=deformation.eigenstrain_rms,
von_mises_strain=deformation.von_mises_strain,
disorientation_angle=disorientation_angle,
disorientation_axis=disorientation_axis,
rotation=deformation.rotation)
return mapped_structure, result
