Source code for calphy.composition_transformation
"""
calphy: a Python library and command line interface for automated free
energy calculations.
Copyright 2021 (c) Sarath Menon^1, Yury Lysogorskiy^2, Ralf Drautz^2
^1: Max Planck Institut für Eisenforschung, Dusseldorf, Germany
^2: Ruhr-University Bochum, Bochum, Germany
calphy is published and distributed under the Academic Software License v1.0 (ASL).
calphy is distributed in the hope that it will be useful for non-commercial academic research,
but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
calphy API is published and distributed under the BSD 3-Clause "New" or "Revised" License
See the LICENSE FILE for more details.
More information about the program can be found in:
Menon, Sarath, Yury Lysogorskiy, Jutta Rogal, and Ralf Drautz.
“Automated Free Energy Calculation from Atomistic Simulations.” Physical Review Materials 5(10), 2021
DOI: 10.1103/PhysRevMaterials.5.103801
For more information contact:
sarath.menon@ruhr-uni-bochum.de/yury.lysogorskiy@icams.rub.de
"""
import re
import numpy as np
import os
import random
import pyscal3.core as pc
from mendeleev import element
from ase.io import read, write
from ase.atoms import Atoms
from pyscal3.core import element_dict
from calphy.integrators import kb
[docs]class CompositionTransformation:
"""
Class for performing composition transformations and
generating necessary pair styles for such transformations.
Parameters
----------
input_structure: ASE object, LAMMPS Data file or LAMMPS dump file
input structure which is used for composition transformation
input_chemical_formula: dict
dictionary of input chemical
output_chemical_formula: string
the required chemical composition string
restrictions: list of strings, optional
Can be used to specify restricted transformations
Notes
-----
This class can be used to create compositional mappings to be used with alchemy mode.
For example, assuming there is a structure file with 500 Atoms of Al in FCC structure, which needs to
be transformed to the structure of 495 Al and 5 Li atoms:
```
comp = CompositionTransformation(filename, {"Al":500}, {"Al":500, "Li":5})
```
Note that the atoms are chosen at random, that is, one cannot specify that only face centered lattice sites
in Al can be transformed to Li.
More complex transformations can be done. For example `{"Al": 495, "Li":5}` to `{"Al": 494, "Li": 2, "O": 3, "C":1}`.
The corresponding input is simply:
```
comp = CompositionTransformation(filename, {"Al":500, "Li":5}, {"Al": 494, "Li": 2, "O": 3, "C":1})
```
Restrictions can be placed on the transformations. In the above example, one can specify that Al-O
transformations should not take place. The code for this is:
```
comp = CompositionTransformation(filename, {"Al":500, "Li":5},
{"Al": 494, "Li": 2, "O": 3, "C":1}, restrictions=["Al-O"])
```
If the restrictions are not satisfiable, an error will be raised.
The LAMMPS data file or dump files do not contain any information about the species except the type numbers.
In general the number of atoms are respected, for example if the file has 10 atoms of type 1, 5 of type 2,
and 1 of type 3. If the `input_chemical_composition` is `{"Li": 5, "Al": 10, "O": 1}`, type 1 is assigned to Al,
type 2 is assigned to Li and type 3 is assigned to O. This is done irrespective of the order in which
`input_chemical_composition` is specified. However, if there are equal number of atoms, the order is respected.
Therefore it is important to make sure that the `input_chemical_composition` is in the same order as that of
types in structure file. For example, consider a NiAl structure of 10 Ni atoms and 10 Al atoms. Ni atoms are type 1 in LAMMPS terminology
and Al atoms are type 2. In this case, to preserve the order, `input_chemical_composition` should be `{"Ni": 5, "Al": 10}`.
Once the calculation is done, there are two possible useful output options. The first one is to generate
the necessary pair coefficient commands for LAMMPS. For the hypothetical transformation `{"Li":5, "Al": 495}` to `{"Al": 494, "Li": 2, "O": 3, "C":1}`,
the pair style can be generated by:
```
alc.update_pair_coeff("pair_coeff * * filename Al")
```
An example pair coefficient needs to be provided. The output for the above command is,
```
('pair_coeff * * filename Al Al Li Li Li',
'pair_coeff * * filename Al O Li O C')
```
These pair styles map the necessary transformation and can be used with `alchemy` mode. The next option
is to output the structure where this pair styles can be employed. This can be done using,
```
alc.write_structure(outfilename)
```
The output is written in LAMMPS dump format.
"""
[docs] def __init__(self, calc):
self.input_chemical_composition = (
calc.composition_scaling._input_chemical_composition
)
self.output_chemical_composition = (
calc.composition_scaling.output_chemical_composition
)
self.restrictions = calc.composition_scaling.restrictions
self.calc = calc
self.actual_species = None
self.new_species = None
self.maxtype = None
self.atom_mark = None
self.atom_species = None
self.mappings = None
self.unique_mappings = None
self.mappingdict = None
self.prepare_mappings()
@property
def entropy_contribution(self):
"""
Find the entropy entribution of the transformation. To get
free energies, multiply by -T.
"""
def _log(val):
if val == 0:
return 0
else:
return np.log(val)
ents = []
for key, val in self.output_chemical_composition.items():
if key in self.input_chemical_composition.keys():
t1 = self.input_chemical_composition[key] / self.natoms
t2 = self.output_chemical_composition[key] / self.natoms
cont = t2 * _log(t2) - t1 * _log(t1)
else:
t1 = 0
t2 = self.output_chemical_composition[key] / self.natoms
cont = t2 * _log(t2) - 0
ents.append(cont)
entropy_term = kb * np.sum(ents)
return entropy_term
[docs] def convert_to_pyscal(self):
"""
Convert a given system to pyscal and give a dict of type mappings
"""
# Create Z_of_type mapping to properly read LAMMPS data files
# This ensures atoms are correctly identified by their element
Z_of_type = dict(
[
(count + 1, element(el).atomic_number)
for count, el in enumerate(self.calc.element)
]
)
aseobj = read(
self.calc.lattice, format="lammps-data", style="atomic", Z_of_type=Z_of_type
)
pstruct = pc.System(aseobj, format="ase")
# here we have to validate the input composition dict; and map it
typelist = pstruct.atoms.species
types, typecounts = np.unique(typelist, return_counts=True)
composition = {types[x]: typecounts[x] for x in range(len(types))}
atomsymbols = self.calc.element
atomtypes = [x + 1 for x in range(len(self.calc.element))]
self.pyscal_structure = pstruct
self.typedict = dict(zip(atomsymbols, atomtypes))
self.reversetypedict = dict(zip(atomtypes, atomsymbols))
self.natoms = self.pyscal_structure.natoms
# Count of actual unique atom types present in the structure
# This matches what's declared in the LAMMPS data file header
self.actual_species_in_structure = len(types)
# Count from calc.element (may include types with 0 atoms)
self.calc_element_count = len(self.calc.element)
# Use actual structure types for pair_coeff consistency
# pair_coeff must match the number declared in the data file header
self.actual_species = self.actual_species_in_structure
self.new_species = len(self.output_chemical_composition) - len(types)
self.maxtype = self.actual_species + 1 # + self.new_species
[docs] def get_composition_transformation(self):
"""
From the two given composition transformation, find the transformation dict
"""
fdiff = {}
for key, val in self.output_chemical_composition.items():
if key in self.input_chemical_composition.keys():
fdiff[key] = val - self.input_chemical_composition[key]
else:
fdiff[key] = val - 0
to_remove = {}
to_add = {}
for key, val in fdiff.items():
if val < 0:
to_remove[key] = np.abs(val)
else:
to_add[key] = val
self.to_remove = to_remove
self.to_add = to_add
[docs] def get_random_index_of_species(self, species_name):
"""
Get a random index of a given species by element name
"""
ids = [count for count, x in enumerate(self.atom_species) if x == species_name]
return ids[np.random.randint(0, len(ids))]
[docs] def mark_atoms(self):
for i in range(self.natoms):
self.atom_mark.append(False)
# Use species (element symbols) instead of numeric types
self.atom_species = self.pyscal_structure.atoms.species
self.atom_type = self.pyscal_structure.atoms.types
self.mappings = [f"{x}-{x}" for x in self.atom_species]
[docs] def update_mark_atoms(self):
self.marked_atoms = []
for key, val in self.to_remove.items():
# key is the element name (e.g., "Mg")
for i in range(100000):
rint = self.get_random_index_of_species(key)
if rint not in self.marked_atoms:
self.atom_mark[rint] = True
self.marked_atoms.append(rint)
val -= 1
if val <= 0:
break
[docs] def update_typedicts(self):
# in a cycle add things to the typedict
for key, val in self.to_add.items():
# print(f"Element {key}, count {val}")
if key in self.typedict.keys():
newtype = self.typedict[key]
else:
newtype = self.maxtype
self.typedict[key] = self.maxtype
self.reversetypedict[self.maxtype] = key
self.maxtype += 1
# print(f"Element {key}, newtype {newtype}")
[docs] def compute_possible_mappings(self):
self.possible_mappings = []
# Now make a list of possible mappings using element names
for key1, val1 in self.to_remove.items():
for key2, val2 in self.to_add.items():
mapping = f"{key1}-{key2}"
if mapping not in self.restrictions:
self.possible_mappings.append(mapping)
[docs] def update_mappings(self):
marked_atoms = self.marked_atoms.copy()
for key, val in self.to_add.items():
# now get all
# we to see if we can get val number of atoms from marked ones
if val > len(marked_atoms):
raise ValueError(
f"Not enough atoms to choose {val} from {len(marked_atoms)} not possible"
)
# otherwise find atoms until we find enough
for i in range(val):
# choose a number from marked atoms
found = False
to_del = []
for x in range(len(self.marked_atoms)):
random_choice = np.random.choice(marked_atoms)
# find corresponding mapping using species name
mapping = f"{self.atom_species[random_choice]}-{key}"
if mapping in self.possible_mappings:
# this is a valid choice
self.mappings[random_choice] = mapping
found = True
if found:
# finish up, change the array, and break
# to_del.append(random_choice)
marked_atoms.remove(random_choice)
break
# if it was not found, the loop finished, throw error
if not found:
raise ValueError(
"A possible transformation could not be found, please check the restrictions"
)
# otherwise modify our marked atoms, list, and move on
# for item in to_del:
# marked_atoms.remove(item)
self.unique_mappings, self.unique_mapping_counts = np.unique(
self.mappings, return_counts=True
)
# now make the transformation dict
self.transformation_list = []
for count, mapping in enumerate(self.unique_mappings):
mapsplit = mapping.split("-")
if not mapsplit[0] == mapsplit[1]:
transformation_dict = {}
transformation_dict["primary_element"] = mapsplit[0]
transformation_dict["secondary_element"] = mapsplit[1]
transformation_dict["count"] = self.unique_mapping_counts[count]
self.transformation_list.append(transformation_dict)
[docs] def get_mappings(self):
self.update_typedicts()
self.compute_possible_mappings()
self.update_mappings()
[docs] def prepare_pair_lists(self):
self.pair_list_old = []
self.pair_list_new = []
for mapping in self.unique_mappings:
map_split = mapping.split("-")
# conserved atom - mappings now use element names directly
if map_split[0] == map_split[1]:
self.pair_list_old.append(map_split[0])
self.pair_list_new.append(map_split[0])
else:
self.pair_list_old.append(map_split[0])
self.pair_list_new.append(map_split[1])
# Special case: 100% transformation with only 1 mapping
# LAMMPS requires pair_coeff to map ALL atom types declared in data file
# Example: Pure Al→Mg with 2 types declared → need ['Al', 'Al'] and ['Mg', 'Mg']
# This ensures consistency between data file type count and pair_coeff mappings
if len(self.unique_mappings) == 1 and self.actual_species > 1:
# Duplicate the single mapping to match number of declared atom types
for _ in range(self.actual_species - 1):
self.pair_list_old.append(self.pair_list_old[0])
self.pair_list_new.append(self.pair_list_new[0])
# Create mapping from transformation strings to UNIQUE type numbers
# Each unique transformation mapping needs its own type for LAMMPS swapping
# Example: Al-Al, Mg-Al, Mg-Mg should map to types 1, 2, 3 respectively
self.mappingdict = {}
for idx, mapping in enumerate(self.unique_mappings, start=1):
self.mappingdict[mapping] = idx
# Update reversetypedict - map each type to its source element
# We'll handle species naming in write_structure to keep types separate
self.reversetypedict = {}
for mapping, type_num in self.mappingdict.items():
source_element = mapping.split("-")[0]
self.reversetypedict[type_num] = source_element
[docs] def update_types(self):
# Update atom_type based on mapping to new types
for x in range(len(self.atom_type)):
self.atom_type[x] = self.mappingdict[self.mappings[x]]
# Update pyscal structure types
self.pyscal_structure.atoms.types = self.atom_type
[docs] def iselement(self, symbol):
try:
_ = element(symbol)
return True
except Exception:
return False
[docs] def update_pair_coeff(self, pair_coeff):
"""
Update pair_coeff command with new element specifications.
Handles both single-file formats (EAM alloy):
pair_coeff * * potential.eam.alloy El1 El2
And two-file formats (MEAM):
pair_coeff * * library.meam El1 El2 potential.meam El1 El2
For MEAM potentials, both element specifications are updated identically.
"""
pcsplit = pair_coeff.strip().split()
result_parts = []
i = 0
while i < len(pcsplit):
token = pcsplit[i]
# Check if this token starts an element specification
# (either it's an element, or the next token is an element)
if self.iselement(token):
# Found start of element list - collect all consecutive elements
element_group = []
while i < len(pcsplit) and self.iselement(pcsplit[i]):
element_group.append(pcsplit[i])
i += 1
# Determine which element list to use based on what we found
# If element_group matches our pair_list_old, replace with pair_list_new
# Otherwise replace with pair_list_old (for the old/reference command)
if element_group == self.pair_list_old or set(element_group) == set(
self.calc.element
):
# This needs special handling - we'll mark position for later
result_parts.append("__ELEMENTS__")
else:
# Keep non-matching element groups as-is
result_parts.extend(element_group)
else:
# Non-element token (potential file, wildcards, options, etc.)
result_parts.append(token)
i += 1
# Now build old and new commands by replacing __ELEMENTS__ markers
pc_old_parts = [
self.pair_list_old if p == "__ELEMENTS__" else [p] for p in result_parts
]
pc_new_parts = [
self.pair_list_new if p == "__ELEMENTS__" else [p] for p in result_parts
]
# Flatten the lists
pc_old = " ".join(
[
item
for sublist in pc_old_parts
for item in (sublist if isinstance(sublist, list) else [sublist])
]
)
pc_new = " ".join(
[
item
for sublist in pc_new_parts
for item in (sublist if isinstance(sublist, list) else [sublist])
]
)
return pc_old, pc_new
[docs] def get_swap_types(self, allow_all_swaps=False):
"""
Get swapping types for configurational entropy calculation.
Returns two separate lists for forward and reverse passes to ensure
proper ergodicity in thermodynamic integration.
Parameters
----------
allow_all_swaps : bool, optional
If True, return all atom types for swapping (including fictitious ones).
If False (default), filter types by element as described below.
Forward pass: Atoms of the element being REMOVED should be able to
swap with each other (e.g., Mg types in Mg→Al enrichment)
Reverse pass: Atoms of the element being ADDED should be able to
swap with each other (e.g., Al types in reverse Al→Mg)
Returns
-------
tuple of (forward_swap_types, reverse_swap_types)
forward_swap_types : list
Types to swap during forward integration (grouped by source element)
reverse_swap_types : list
Types to swap during reverse integration (grouped by target element)
Examples
--------
Enrichment (Mg→Al): Mappings [Al-Al(1), Mg-Al(2), Mg-Mg(3)]
forward: [2, 3] - swap Mg types (Mg being removed)
reverse: [1, 2] - swap Al types (Al being added back)
with allow_all_swaps=True: [1, 2, 3] for both forward and reverse
100% transformation (Mg→Al): Mappings [Mg-Al(1)]
forward: [] - only one type, no swapping possible
reverse: [] - only one type, no swapping possible
Partial addition (Mg→Mg+Al): Mappings [Mg-Al(1), Mg-Mg(2)]
forward: [1, 2] - swap Mg types (Mg being removed)
reverse: [] - only one Al type exists, no swapping possible
"""
# If allow_all_swaps is True, return all types from mappingdict
if allow_all_swaps:
all_types = sorted(list(self.mappingdict.values()))
if len(all_types) >= 2:
return all_types, all_types
else:
return [], []
forward_swap_types = []
reverse_swap_types = []
# Get the element being removed and added from transformation_list
if not self.transformation_list:
return [], []
# For composition transformations, we typically have one primary transformation
# primary_element is being removed, secondary_element is being added
primary_element = self.transformation_list[0]["primary_element"]
secondary_element = self.transformation_list[0]["secondary_element"]
# Forward pass: collect all types where SOURCE element matches primary_element
# These are atoms that could potentially be transformed
for mapping_str, type_num in self.mappingdict.items():
source_element, target_element = mapping_str.split("-")
if source_element == primary_element:
forward_swap_types.append(type_num)
# Reverse pass: We're removing secondary_element and adding back primary_element
# We need types where the TARGET (after forward transformation) is secondary_element
# These are the atoms that exist as secondary_element after forward pass
# BUT: if there are no such types (edge case), we need types that will CREATE
# primary_element in reverse - i.e., types involving primary_element
for mapping_str, type_num in self.mappingdict.items():
source_element, target_element = mapping_str.split("-")
if target_element == secondary_element:
reverse_swap_types.append(type_num)
# Edge case: if only 1 type maps to secondary_element, we need more types for swapping
# In this case, include types that involve primary_element
# Example: Mg→Mg+Al has types [Mg→Al(1), Mg→Mg(2)], reverse needs [1,2] not just [1]
if len(reverse_swap_types) < 2:
reverse_swap_types = []
for mapping_str, type_num in self.mappingdict.items():
source_element, target_element = mapping_str.split("-")
if (
source_element == primary_element
or target_element == primary_element
):
reverse_swap_types.append(type_num)
# Sort for consistency
forward_swap_types.sort()
reverse_swap_types.sort()
# If only 1 type exists after all attempts, swapping is not possible - return empty list
# LAMMPS atom/swap requires at least 2 types
if len(forward_swap_types) < 2:
forward_swap_types = []
if len(reverse_swap_types) < 2:
reverse_swap_types = []
return forward_swap_types, reverse_swap_types
[docs] def write_structure(self, outfilename, for_fe_mode=False):
"""Write structure to LAMMPS data file with proper type declarations.
Parameters
----------
outfilename : str
Path to write the LAMMPS data file.
for_fe_mode : bool, optional
If True, write atom types reflecting the **target** element's type
number from ``typedict`` (e.g. Cu=1, Ni=2) and declare
``calc_element_count`` atom types in the header. This is required
when the output file will be used as a structure input for a plain
``fe`` calculation where the pair_coeff enumerates all elements.
If False (default), use the ``mappingdict`` type indices (one per
unique transformation mapping) which is what the alchemy /
composition_scaling LAMMPS routines expect for ``fix atom/swap``.
"""
from ase.io import write as ase_write
from ase import Atoms as ASEAtoms
# Get positions and cell from pyscal structure
positions = self.pyscal_structure.atoms.positions
cell = self.pyscal_structure.box
# Build per-atom type numbers and chemical symbols.
# for_fe_mode: remap each atom to its TARGET element's typedict entry so
# the structure matches the pair_coeff element ordering exactly.
# Default (alchemy) mode: keep the mappingdict indices as-is.
if for_fe_mode:
fe_types = [
self.typedict[self.mappings[i].split('-')[1]]
for i in range(len(self.mappings))
]
invert_typedict = {v: k for k, v in self.typedict.items()}
symbols = [invert_typedict[t] for t in fe_types]
else:
fe_types = None # will use pyscal types below
symbols = [self.reversetypedict[t] for t in self.pyscal_structure.atoms.types]
ase_atoms = ASEAtoms(symbols=symbols, positions=positions, cell=cell, pbc=True)
# Write using ASE with atom_style
ase_write(outfilename, ase_atoms, format="lammps-data", atom_style="atomic")
# Post-process to fix the type column with our custom types
with open(outfilename, "r") as f:
lines = f.readlines()
# Choose which per-atom type list to use
atom_types_to_write = fe_types if for_fe_mode else list(self.pyscal_structure.atoms.types)
# Find the Atoms section and replace type numbers
# Support different ASE formats: "Atoms # atomic", "Atoms # full", or just "Atoms"
in_atoms_section = False
atom_idx = 0
for i, line in enumerate(lines):
# More flexible detection of Atoms section header
if "Atoms" in line:
in_atoms_section = True
continue
if in_atoms_section and line.strip():
parts = line.split()
if len(parts) >= 5: # atom_id type x y z
# Replace the type (column 1, 0-indexed) with our custom type
custom_type = atom_types_to_write[atom_idx]
parts[1] = str(custom_type)
lines[i] = " " + " ".join(parts) + "\n"
atom_idx += 1
if atom_idx >= len(atom_types_to_write):
break
# Verify all atoms were updated
expected_atoms = len(self.pyscal_structure.atoms.types)
if atom_idx != expected_atoms:
raise RuntimeError(
f"Failed to update all atoms in {outfilename}. "
f"Expected {expected_atoms} atoms but only updated {atom_idx}. "
f"This may indicate a problem with ASE LAMMPS file formatting."
)
# Update the number of atom types in the header.
# fe mode: must match the number of elements in the pair_coeff.
# alchemy mode: one type per unique mapping (for fix atom/swap).
if for_fe_mode:
required_ntypes = self.calc_element_count
else:
required_ntypes = len(self.pair_list_old)
for i, line in enumerate(lines):
if "atom types" in line:
lines[i] = f"{required_ntypes} atom types\n"
break
# Write the corrected file
with open(outfilename, "w") as f:
f.writelines(lines)
[docs] def prepare_mappings(self):
self.atom_mark = []
self.atom_species = []
self.mappings = []
self.unique_mappings = []
self.get_composition_transformation()
self.convert_to_pyscal()
self.mark_atoms()
self.update_mark_atoms()
self.get_mappings()
self.prepare_pair_lists()
self.update_types()