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,
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:
Deduce them analytically from the stiffness tensors of the two phases, or
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 [RahLofErh22], 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,
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)
from trainstation import 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.
Presently it can therefore not be used together with the canonical ensemble.
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)
from trainstation import 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)