Constituent strain calculations

This example demonstrates how strain can be incorporated in cluster expansions in icet using the constituent strain formalism originally introduced by Laks et al. [LakFerFro92].

Warning

Since this module has only been tested for FCC crystals, great caution should be exercised when using it, in particular for non-cubic systems.

Warning

Due to the requirement to calculate \(\Delta E_\text{CS}(\hat{\mathbf{k}}, c)\) separately (see below), this module is only recommended to advanced users. We can only provide very limited support.

In a regular cluster expansion, strain is typically not captured properly. Consider, for example, a phase separated system with phases 1 and 2. If the interface is coherent (i.e., the atomic planes line up without defects), the two phases will have to adopt the same in-plane lattice parameter. This leads to strain, unless the two phases happen to have the same equilibrium lattice parameter. A cluster expansion, however, normally has fairly short cutoffs, and will not be able to describe this situation, since the strain extends into the phases far from the interface, well beyond the cutoffs of the cluster expansion. Furthermore, the strain energy is dependent not only on the composition of the phases, but also the orientation of the interface, since the stiffness tensors are not isotropic.

To overcome this limitation of cluster expansions, the constituent strain formalism adds a term to the cluster expansion that is a sum in reciprocal space,

\[\text{strain energy} = \sum_\mathbf{k} \Delta E_\text{CS} (\hat{\mathbf{k}}, c) F(\mathbf{k}, \sigma)\]

Here, the sum runs over all \(\mathbf{k}\) points in the first Brillouin zone. The term \(\Delta E_\text{CS}(\hat{\mathbf{k}}, c)\) expresses the dependency of strain energy on crystallographic direction and overall concentration \(c\), while \(F(\mathbf{k}, \sigma)\) quantifies the variation of concentration with periodicity \(\mathbf{k}\) for an atomic structure with occupation vector \(\sigma\).

icet provides functionality for efficiently evaluating the above sum when fitting cluster expansions and running Monte Carlo simulations. It does, however, only provide a framework for evaluating \(\Delta E_\text{CS}(\hat{\mathbf{k}}, c)\), not explicit functionality for its definition for a specific material. For more information, including an example, see below.

Here, we will demonstrate the usage of the constituent strain module using the Ag-Cu alloy as an example. Ag-Cu is immiscible and has a large size mismatch. This means that it will phase separate, and if it does so coherently there will be significant strain energy constributions. In practice, it is likely that dislocations will form that make the interface bewteen Ag and Cu incoherent. If that happens, the strain energy will be relaxed, but for the purpose of demonstration we will ignore this aspect here.

Define \(\Delta E_\text{CS} (\hat{\mathbf{k}}, c)\)

The first step in the constituent strain calculation is to define \(\Delta E_\text{CS}(\hat{\mathbf{k}}, c)\). This module does not provide explicit functionality to calculate this term, but an example for how it can be done is provided below.

Ideally, we need to be able to evaluate this term for any \(\mathbf{k}\) point (or rather for any direction \(\hat{\mathbf{k}}\) in reciprocal space) and for any overall concentration \(c\). When we later on create instances of the ConstituentStrain class, we will thus supply it with functions (and not floats or arrays) that evaluate \(\Delta E_\text{CS} (\hat{\mathbf{k}}, c)\) for given \(\hat{\mathbf{k}}\) and \(c\). There are at least two approaches for constructing such functions:

  1. Deduce them analytically from the stiffness tensors of the two phases, or

  2. Do a long series of DFT calculations with strain in different orientations and interpolate the results for remaining \(\hat{\mathbf{k}}\) and \(c\).

The first of these approaches is described in [LakFerFro92]. The second was outlined in [OzoWolZun98], and is demonstrated in the cubic case in this Jupyter notebook. The notebook demonstrates a procedure for obtaining the necessary energies in the form of expansion coefficients for Redlich-Kister polynomials in six crystallographic directions. To handle arbitrary directions, these expansion coefficients are interpolated using the functions in this module. Functions from this module (which is specific for the cubic case and the choice of crystallographic directions in the notebook) will be supplied as input arguments to ConstituentStrain. For further details, please see [REFERENCE TO OUR PAPER WHEN PUBLISHED], in particular its Supporting Information.

Fit a cluster expansion with strain

Fitting of a cluster expansion with strain is not very different from fitting a regular cluster expansion. The only thing we need to keep in mind is that the total energy of a configuration \(\sigma\) consists of two terms,

\[E(\sigma) = E_\text{CE}(\sigma) + \Delta E_\text{CS} (\hat{\mathbf{k}}, c).\]

Here, the first term is described by a regular cluster expansion and the second term is added afterwards. When fitting the cluster expansion to energies calculated with DFT (or any other method), we therefore first need to subtract the second term from this energy, and fit only the remainder (\(E_\text{CE}(\sigma)\)).

The following example illustrates the process. For the sake of simplicity, we will use the EMT calculator in ASE instead of DFT to calculate energies. We begin by some imports (including our custom functions from the previous section), we define a functon that will allow us to calculate energies, and we calculate the energies of the pure phases (Ag and Cu) to be able to define mixing energy properly.

from custom_functions import (custom_k_to_parameter_function,
                              custom_strain_energy_function)
from icet import (ClusterSpace,
                  StructureContainer,
                  ClusterExpansion,
                  Optimizer)
from icet.tools import (ConstituentStrain,
                        enumerate_structures,
                        get_mixing_energy_constraints)
from ase.build import bulk
from ase.calculators.emt import EMT
from ase.optimize import QuasiNewton
from ase.constraints import UnitCellFilter
import numpy as np


def get_relaxed_energy(structure):
    calc = EMT()
    structure = structure.copy()  # Since we will relax it
    structure.set_calculator(calc)
    ucf = UnitCellFilter(structure)
    qn = QuasiNewton(ucf)
    qn.run(fmax=0.05)
    return structure.get_potential_energy()


# Calculate energy of elements so we can define mixing energy properly
a = 4.0
prim = bulk('Ag', a=a)
elemental_energies = {}
for element in ['Ag', 'Cu']:
    structure = prim.copy()
    structure[0].symbol = element
    elemental_energies[element] = get_relaxed_energy(structure)

We then define a cluster space just as we normally would. We use short cutoffs here for the sake of demonstration.

cutoffs = a * np.array([1.2, 0.8])
cs = ClusterSpace(prim, cutoffs=cutoffs, chemical_symbols=[['Ag', 'Cu']])

We then start to fill a structure container with data using enumeration. Here, we need to calculate strain by initiating an instance of ConstituentStrain for each structure. To this end, we use the functions defined in custom_functions.

sc = StructureContainer(cluster_space=cs)
for structure in enumerate_structures(prim, range(6), ['Ag', 'Cu']):
    conc = structure.get_chemical_symbols().count('Cu') / len(structure)

    energy = get_relaxed_energy(structure)
    mix_energy = energy / len(structure) \
        - conc * elemental_energies['Cu'] - (1 - conc) * elemental_energies['Ag']

    strain = ConstituentStrain(supercell=structure,
                               primitive_structure=prim,
                               chemical_symbols=['Ag', 'Cu'],
                               concentration_symbol='Cu',
                               strain_energy_function=custom_strain_energy_function,
                               k_to_parameter_function=custom_k_to_parameter_function,
                               damping=10.0)
    strain_energy = strain.get_constituent_strain(structure.get_atomic_numbers())

    sc.add_structure(structure, properties={'concentration': conc,
                                            'mix_energy': mix_energy,
                                            'mix_energy_wo_strain': mix_energy - strain_energy})

Finally, we fit the cluster expansion. Here, we use get_mixing_energy_constraints to reproduce the pure phases exactly.

A, y = sc.get_fit_data(key='mix_energy_wo_strain')
constr = get_mixing_energy_constraints(cs)
A_constrained = constr.transform(A)

opt = Optimizer((A_constrained, y))
opt.train()
parameters = constr.inverse_transform(opt.parameters)
ce = ClusterExpansion(cluster_space=cs, parameters=parameters, metadata=opt.summary)
ce.write('mixing_energy_wo_strain.ce')

Run Monte Carlo with strain cluster expansions

The mchammer module includes custom calculators and observers for running Monte Carlo simulations with cluster expansions with strain. If these are used, the procedure is not very different from Monte Carlo simulations with a regular cluster expansion. Instead of creating a ClusterExpansionCalculator, we create a ConstituentStrainCalculator with a ConstituentStrain instance as input. The latter is constructed in the same way as demonstrated above. If we want to be able to extract the strain energy separately, we also need to attach an observer to the ensemble object.

Note

The ConstituentStrainCalculator is currently only compatible with ensembles that flip one site at a time. It can therefore not be used together with the canonical ensemble currently.

from custom_functions import (custom_k_to_parameter_function,
                              custom_strain_energy_function)
from icet import ClusterExpansion
from icet.tools import ConstituentStrain
from mchammer.ensembles import VCSGCEnsemble
from mchammer.calculators import ConstituentStrainCalculator
from ase.build import make_supercell
import numpy as np
import os

# Read cluster expansion and define supercell
ce = ClusterExpansion.read('mixing_energy_wo_strain.ce')
prim = ce.get_cluster_space_copy().primitive_structure
P = np.array([[-1, 1, 1],
              [1, -1, 1],
              [1, 1, -1]])
structure = make_supercell(prim, P)
structure = structure.repeat((20, 3, 3))

output_directory = 'monte_carlo_data'
try:
    os.mkdir(output_directory)
except FileExistsError:
    pass

T = 300
for phi in np.arange(0.1, -2.1, -0.1):
    # The constituent strain object can be used standalone
    # but we will use it together with a calucalator
    strain = ConstituentStrain(supercell=structure,
                               primitive_structure=prim,
                               chemical_symbols=['Ag', 'Cu'],
                               concentration_symbol='Cu',
                               strain_energy_function=custom_strain_energy_function,
                               k_to_parameter_function=custom_k_to_parameter_function,
                               damping=10.0)

    # Here we define a calculator to be used with mchammer
    calc = ConstituentStrainCalculator(constituent_strain=strain,
                                       cluster_expansion=ce)

    # Now we can run Monte Carlo in regular fashion
    ens = VCSGCEnsemble(calculator=calc,
                        structure=structure,
                        temperature=T,
                        dc_filename=f'{output_directory}/vcsgc-T{T}-phi{phi:+.3f}.dc',
                        phis={'Cu': phi},
                        kappa=200)
    ens.run(number_of_trial_steps=len(structure) * 10)
    structure = ens.structure  # Use last structure as starting point at next phi
    print(phi, structure)

Note

Due to the non-locality of the structure factor, Monte Carlo simulations with constituent strain are typically significantly more time-consuming than if strain is not included.

Source code

The complete source code is available in examples/advanced_topics/constituent_strain/1_fit_cluster_expansion.py

from custom_functions import (custom_k_to_parameter_function,
                              custom_strain_energy_function)
from icet import (ClusterSpace,
                  StructureContainer,
                  ClusterExpansion,
                  Optimizer)
from icet.tools import (ConstituentStrain,
                        enumerate_structures,
                        get_mixing_energy_constraints)
from ase.build import bulk
from ase.calculators.emt import EMT
from ase.optimize import QuasiNewton
from ase.constraints import UnitCellFilter
import numpy as np


def get_relaxed_energy(structure):
    calc = EMT()
    structure = structure.copy()  # Since we will relax it
    structure.set_calculator(calc)
    ucf = UnitCellFilter(structure)
    qn = QuasiNewton(ucf)
    qn.run(fmax=0.05)
    return structure.get_potential_energy()


# Calculate energy of elements so we can define mixing energy properly
a = 4.0
prim = bulk('Ag', a=a)
elemental_energies = {}
for element in ['Ag', 'Cu']:
    structure = prim.copy()
    structure[0].symbol = element
    elemental_energies[element] = get_relaxed_energy(structure)

# Define a cluster space
cutoffs = a * np.array([1.2, 0.8])
cs = ClusterSpace(prim, cutoffs=cutoffs, chemical_symbols=[['Ag', 'Cu']])

# Fill a structure container with data
sc = StructureContainer(cluster_space=cs)
for structure in enumerate_structures(prim, range(6), ['Ag', 'Cu']):
    conc = structure.get_chemical_symbols().count('Cu') / len(structure)

    energy = get_relaxed_energy(structure)
    mix_energy = energy / len(structure) \
        - conc * elemental_energies['Cu'] - (1 - conc) * elemental_energies['Ag']

    strain = ConstituentStrain(supercell=structure,
                               primitive_structure=prim,
                               chemical_symbols=['Ag', 'Cu'],
                               concentration_symbol='Cu',
                               strain_energy_function=custom_strain_energy_function,
                               k_to_parameter_function=custom_k_to_parameter_function,
                               damping=10.0)
    strain_energy = strain.get_constituent_strain(structure.get_atomic_numbers())

    sc.add_structure(structure, properties={'concentration': conc,
                                            'mix_energy': mix_energy,
                                            'mix_energy_wo_strain': mix_energy - strain_energy})

# Constrain sensing matrix so as to reproduce c=0 and c=1 exactly,
# then fit a cluster expansion
A, y = sc.get_fit_data(key='mix_energy_wo_strain')
constr = get_mixing_energy_constraints(cs)
A_constrained = constr.transform(A)

opt = Optimizer((A_constrained, y))
opt.train()
parameters = constr.inverse_transform(opt.parameters)
ce = ClusterExpansion(cluster_space=cs, parameters=parameters, metadata=opt.summary)
ce.write('mixing_energy_wo_strain.ce')

The complete source code is available in examples/advanced_topics/constituent_strain/2_run_monte_carlo.py

from custom_functions import (custom_k_to_parameter_function,
                              custom_strain_energy_function)
from icet import ClusterExpansion
from icet.tools import ConstituentStrain
from mchammer.ensembles import VCSGCEnsemble
from mchammer.calculators import ConstituentStrainCalculator
from ase.build import make_supercell
import numpy as np
import os

# Read cluster expansion and define supercell
ce = ClusterExpansion.read('mixing_energy_wo_strain.ce')
prim = ce.get_cluster_space_copy().primitive_structure
P = np.array([[-1, 1, 1],
              [1, -1, 1],
              [1, 1, -1]])
structure = make_supercell(prim, P)
structure = structure.repeat((20, 3, 3))

output_directory = 'monte_carlo_data'
try:
    os.mkdir(output_directory)
except FileExistsError:
    pass

T = 300
for phi in np.arange(0.1, -2.1, -0.1):
    # The constituent strain object can be used standalone
    # but we will use it together with a calucalator
    strain = ConstituentStrain(supercell=structure,
                               primitive_structure=prim,
                               chemical_symbols=['Ag', 'Cu'],
                               concentration_symbol='Cu',
                               strain_energy_function=custom_strain_energy_function,
                               k_to_parameter_function=custom_k_to_parameter_function,
                               damping=10.0)

    # Here we define a calculator to be used with mchammer
    calc = ConstituentStrainCalculator(constituent_strain=strain,
                                       cluster_expansion=ce)

    # Now we can run Monte Carlo in regular fashion
    ens = VCSGCEnsemble(calculator=calc,
                        structure=structure,
                        temperature=T,
                        dc_filename=f'{output_directory}/vcsgc-T{T}-phi{phi:+.3f}.dc',
                        phis={'Cu': phi},
                        kappa=200)
    ens.run(number_of_trial_steps=len(structure) * 10)
    structure = ens.structure  # Use last structure as starting point at next phi
    print(phi, structure)

The source code for interpolation of \(\Delta E_\text{CS} (\hat{\mathbf{k}}, c)\) in the cubic case is available in examples/advanced_topics/constituent_strain/custom_functions.py

from icet.tools import redlich_kister
import numpy as np
import itertools


def _fit_cs_rk_parameters(cs_rk_parameters):
    """
    Make a polynomial fit of Redlich-Kister parameters on
    the stereographic projection of the 1BZ of FCC.
    """
    A = []
    y = []
    for orientation, parameters in cs_rk_parameters.items():
        projection = _get_projection(orientation)
        A.append(_get_xy_vector(*projection))
        y.append(parameters)

    rk_fit, _, _, _ = np.linalg.lstsq(A, y, rcond=None)
    return rk_fit


def _get_projection(orientation):
    """
    Stereographic projection of the orientation of an interface.
    Assumes cubic symmetry.
    """
    orientation = np.array(sorted(np.abs(orientation)))
    orientation = orientation / np.linalg.norm(orientation)
    return orientation[:2]


def _get_xy_vector(x, y, deg=2):
    """
    2D polynomial of degree `deg`
    """
    vector = []
    for i, j in itertools.product(range(deg + 1), repeat=2):
        if i + j > deg:
            continue
        vector.append(x**i * y**j)
    return vector


def _k_to_parameter_function(k, cs_fitted_rk_parameters):
    """
    Calculate strain energy at a specific concentration and k point
    using constitutent_strain_functions as fitted with
    _fit_cs_parameters.
    """
    projection = _get_projection(k)
    vector = _get_xy_vector(*projection)
    parameters = np.dot(vector, cs_fitted_rk_parameters)
    return parameters


def custom_k_to_parameter_function(k):
    """
    Create function that precomputes Redlich-Kister parameters
    for a specific k point.
    """
    # Read Redlich-Kister parameters for constituent strain
    cs_data = np.loadtxt('constituent-strain-RK-parameters.data')
    cs_rk_parameters = {}
    for row in cs_data:
        cs_rk_parameters[tuple(int(i) for i in row[:3])] = row[3:]
    cs_fitted_rk_parameters = _fit_cs_rk_parameters(cs_rk_parameters)

    # Define function
    f = _k_to_parameter_function(k, cs_fitted_rk_parameters)
    return f


def custom_strain_energy_function(parameters, c):
    """
    Create function that takes Redlich-Kister parameters and
    concentration and returns corresponding strain energy.
    """
    return redlich_kister(c, *parameters)