Mapping structures#

A cluster vector calculation requires all atoms to reside on a fixed lattice. Properties of interest, on the other hand, are typically calculated for a structure in which cell metric and atoms have been allowed to relax. Unless the ideal structures have been saved prior to relaxation, one is therefore faced with the task of mapping back the relaxed structure onto the ideal one. In some cases, in particular involving vacancies, relaxation can also lead to atoms moving between sites, in which case remapping is mandatory.

This is the purpose of the function map_structure_to_reference(). The function is also useful to analyze whether the relaxation has gone too far for the cluster expansion to be viable, i.e., whether the ideal structure from which the relaxation started is not a valid representation of the structure for which the property has been obtained.

Import modules#

The map_structure_to_reference() function needs to be imported together with some additional functionality from ASE.

from icet.tools import map_structure_to_reference
from ase.build import bulk
from ase import Atom

Prepare dummy structures#

First, for the sake of demonstration, a reference structure defining the ideal lattice is created, and a supercell thereof is scaled and rattled to simulate relaxation in an energy minimixation.

reference = bulk('Au', a=4.00)

supercell = reference.repeat(3)
supercell.rattle(0.1, seed=42)
supercell.set_cell(1.05 * supercell.cell, scale_atoms=True)

# Switch some atoms to Pd
for i in [0, 1, 5, 8, 10]:
    supercell[i].symbol = 'Pd'

Map relaxed structure onto ideal structure#

The structure can now be mapped onto a structure in which all atoms reside on ideal lattice sites. The function returns that structure together with a StructureMapping object carrying the supplementary information. Among other things it reports the maximum and the average distance between the relaxed positions and the sites they were assigned to (drmax and dravg), the strain and its separation into a change of volume and a change of shape, the rotation, the transformation matrix, and the tags of any warnings that were triggered. Printing it, or displaying it in a notebook, gives a summary of the mapping.

The rotation is reported as disorientation_angle and disorientation_axis, which belong together and describe one rotation. Rotating a structure several times in succession does not give several rotations to report, since any number of rotations compose into a single rotation about a single axis, and that is what the mapping recovers. Both are reduced over the symmetry of the reference structure, so the angle is the smallest one among the descriptions that the structure cannot tell apart. For a cubic crystal it therefore never exceeds about 62.8 degrees.

ideal_structure, info = map_structure_to_reference(supercell, reference)
print('Maximum displacement: {:.3f} Angstrom'.format(info.drmax))
print('Average displacement: {:.3f} Angstrom'.format(info.dravg))

Structures with vacancies#

For cluster expansions with vacancies, one typically wants to map the relaxed structure onto an ideal lattice that explicitly contains vacant sites. In that case, if the volume of the cell has changed during relaxation, it can be tricky to determine the size of the ideal supercell. To help the function with this task, an additional keyword, inert_species, can be specified, which is a list of species that reside on sublattices without vacancies.

In the example below, a Au-Pd-H-vacancy system is created. The system of choice consists of two sublattices, one occupied by Au and Pd and another occupied by H and vacancies. Since Au and Pd belong to a sublattice in which we do not allow vacancies, we may set inert_species = ['Au', 'Pd'].

reference = bulk('Au', a=4.00)
reference.append(Atom(('H'), (2.0, 2.0, 2.0)))

supercell = reference.repeat(3)
supercell.rattle(0.1, seed=42)
supercell.set_cell(1.05 * supercell.cell, scale_atoms=True)

# Switch some Au to Pd and delete some H (to create vacancies)
for i in [0, 4, 6, 2, 7, 3, 17]:
    if supercell[i].symbol == 'Au':
        supercell[i].symbol = 'Pd'
    elif supercell[i].symbol == 'H':
        del supercell[i]

ideal_structure, info = map_structure_to_reference(supercell, reference,
                                                   inert_species=['Au', 'Pd'])
print('Maximum displacement: {:.3f} Angstrom'.format(info.drmax))
print('Average displacement: {:.3f} Angstrom'.format(info.dravg))

The mapped structure will contain atoms of type X, which represent vacancies. Which sites end up vacant is decided by the mapping itself, since leaving a site vacant carries no cost in the underlying Hungarian algorithm. The vacancies are therefore placed wherever they leave the remaining atoms closest to their sites.

The reference structure in this example illustrates at the same time how interstitials are treated. The interstitial sites are simply made part of the reference structure, and the sites that are not occupied in a given relaxed structure come out as vacancies. The converse does not work: a relaxed structure cannot contain atoms for which the reference structure provides no site, and attempting it raises a ValueError. If a relaxed structure contains interstitial atoms, the interstitial sites therefore have to be added to the reference structure.

If there is no sublattice without vacancies, one typically has to set the keyword argument assume_no_cell_relaxation to True. The volume of the relaxed structure is then not rescaled before the size of the ideal supercell is determined, which is what the number of inert sites would otherwise be needed for. Note that the ideal supercell always has the cell metric of the reference structure; the keyword argument affects only how its size is determined, and it implies that the relaxed structure is expected to be obtainable from the reference structure by an integer transformation.

Inspecting a mapping#

The mapping reports what it did at the DEBUG level, which is not printed by default. This includes how the volume was rescaled, whether the transformation matrix was obtained by rounding or by a search, the matrix itself, the eigenvalues of the strain tensor, any rigid offset that was removed, and how many atoms were assigned to how many sites. It is therefore the first thing to look at if a mapping fails or does not give the expected result. The output is switched on as follows:

from icet.input_output.logging_tools import set_log_config
set_log_config(level='DEBUG')

Source code#

The complete source code is available in examples/advanced_topics/map_structure_to_reference.py

"""
This example demonstrates how to map a structure in which the cell has
been scaled and/or the atoms displaced onto an ideal (primitive) structure
"""

# Import modules
from icet.tools import map_structure_to_reference
from ase.build import bulk
from ase import Atom
# End import

# Begin by creating a reference structure, in this case fcc Au.
# Then create a supercell structure, scale the cell and displace the atoms to
# simulate a relaxed structure.
reference = bulk('Au', a=4.00)

supercell = reference.repeat(3)
supercell.rattle(0.1, seed=42)
supercell.set_cell(1.05 * supercell.cell, scale_atoms=True)

# Switch some atoms to Pd
for i in [0, 1, 5, 8, 10]:
    supercell[i].symbol = 'Pd'

# Map the "relaxed" structure onto an ideal supercell
ideal_structure, info = map_structure_to_reference(supercell, reference)
print('Maximum displacement: {:.3f} Angstrom'.format(info.drmax))
print('Average displacement: {:.3f} Angstrom'.format(info.dravg))

# Map a structure that contains vacancies, in this case Pd-Au-H-Vac, in which
# Pd and Au share one sublattice and Pd and H another.
reference = bulk('Au', a=4.00)
reference.append(Atom(('H'), (2.0, 2.0, 2.0)))

supercell = reference.repeat(3)
supercell.rattle(0.1, seed=42)
supercell.set_cell(1.05 * supercell.cell, scale_atoms=True)

# Switch some Au to Pd and delete some H (to create vacancies)
for i in [0, 4, 6, 2, 7, 3, 17]:
    if supercell[i].symbol == 'Au':
        supercell[i].symbol = 'Pd'
    elif supercell[i].symbol == 'H':
        del supercell[i]

ideal_structure, info = map_structure_to_reference(supercell, reference,
                                                   inert_species=['Au', 'Pd'])
print('Maximum displacement: {:.3f} Angstrom'.format(info.drmax))
print('Average displacement: {:.3f} Angstrom'.format(info.dravg))