# Copyright (c) 2022-2026, Yongchao Wu in Aalto University
# This file is from the mdapy project, released under the BSD 3-Clause License.
"""
Embedded Atom Method (EAM) Potential Calculator
================================================
This module implements the Embedded Atom Method (EAM) for static simulations.
The EAM potential is widely used for metallic systems and consists of two main components:
1. **Embedding energy**: :math:`F_i(\\rho_i)` - energy to embed atom i into the electron density :math:`\\rho_i`
2. **Pair interaction**: :math:`\\phi_{ij}(r_{ij})` - pair potential between atoms i and j
The total energy is given by:
.. math::
E_{tot} = \\sum_i F_i(\\rho_i) + \\frac{1}{2} \\sum_{i,j \\neq i} \\phi_{ij}(r_{ij})
where the electron density at atom i is:
.. math::
\\rho_i = \\sum_{j \\neq i} \\rho_j(r_{ij})
The forces on atom i are computed as:
.. math::
\\mathbf{F}_i = -\\sum_{j \\neq i} \\left[ \\left(\\frac{dF_i}{d\\rho_i} + \\frac{dF_j}{d\\rho_j}\\right) \\frac{d\\rho_j}{dr_{ij}} + \\frac{d\\phi_{ij}}{dr_{ij}} \\right] \\frac{\\mathbf{r}_{ij}}{r_{ij}}
The stress tensor is computed from the virial theorem:
.. math::
\\sigma_{\\alpha\\beta} = -\\frac{1}{V} \\sum_i \\sum_{j \\neq i} r_{ij,\\alpha} F_{ij,\\beta}
References
----------
Daw, M. S., & Baskes, M. I. (1984). Embedded-atom method: Derivation and
application to impurities, surfaces, and other defects in metals.
Physical Review B, 29(12), 6443.
Foiles, S. M., Baskes, M. I., & Daw, M. S. (1986). Embedded-atom-method
functions for the fcc metals Cu, Ag, Au, Ni, Pd, Pt, and their alloys.
Physical Review B, 33(12), 7983.
"""
from __future__ import annotations
import numpy as np
from numpy.typing import NDArray
from mdapy.neighbor import Neighbor
import polars as pl
from mdapy.box import Box
from mdapy.calculator import CalculatorMP
from mdapy import _eam
from mdapy.parallel import get_num_threads
import mdapy.tool_function as tool
from mdapy.data import atomic_masses, atomic_numbers
from typing import TYPE_CHECKING, Optional, Tuple, List, Dict
import numpy.typing as npt
import datetime
if TYPE_CHECKING:
from matplotlib.figure import Figure
from matplotlib.axes import Axes
[docs]
class EAM(CalculatorMP):
"""
Embedded Atom Method (EAM) potential calculator.
This class implements the EAM potential for metallic systems, supporting both
single-element and multi-element alloy systems. It reads EAM potential files
in the LAMMPS alloy format and provides methods to calculate energies, forces,
stresses, and virials.
Parameters
----------
filename : str
Path to the EAM potential file in LAMMPS alloy format.
mass_list : List[float], optional
Atomic masses in the same order as ``elements_list``. When ``None``
(default) masses are looked up from the built-in periodic table via
the element symbols read from the potential file.
Notes
-----
After parsing, the instance exposes the potential as these attributes:
* ``filename`` — path to the EAM file.
* ``header`` — the first 3 lines of the file.
* ``Nelements``, ``elements_list`` — element count and symbols.
* ``nrho`` / ``drho`` — size and spacing of the electron-density grid.
* ``nr`` / ``dr`` — size and spacing of the distance grid.
* ``rc`` — cutoff radius, in Å.
* ``F_rho`` — embedding function F(ρ), shape ``(Nelements, nrho)``.
* ``rho_r`` — electron density ρ(r), shape ``(Nelements, nr)``.
* ``phi_r`` — pair potential φ(r), shape ``(Nelements, Nelements, nr)``.
* ``r``, ``rho`` — grid points, shape ``(nr,)`` and ``(nrho,)``.
* ``results`` — dict storing computed energies, forces, stresses, virials.
"""
def __init__(self, filename: str, mass_list: Optional[List[float]] = None) -> None:
self.filename: str = filename
self.header: List[str] = []
self.Nelements: int = 0
self.elements_list: List[str] = []
self.nrho: int = 0
self.drho: float = 0.0
self.nr: int = 0
self.dr: float = 0.0
self.rc: float = 0.0
self.F_rho: NDArray[np.float64] = np.array([])
self.rho_r: NDArray[np.float64] = np.array([])
self.phi_r: NDArray[np.float64] = np.array([])
self._rphi_r: NDArray[np.float64] = np.array([])
self.r: NDArray[np.float64] = np.array([])
self.rho: NDArray[np.float64] = np.array([])
self.results: Dict[str, NDArray] = {}
self._read_eam_alloy()
if mass_list is None:
mass_list = []
for i in range(self.Nelements):
species = self.elements_list[i]
if species in atomic_numbers.keys():
aindex = atomic_numbers[species]
amass = atomic_masses[aindex]
else:
amass = 0.0
mass_list.append(amass)
else:
assert len(mass_list) == len(self.elements_list)
self.mass_list = mass_list
self._eam = _eam.EAM(
self.rc, self.dr, self.drho, self.F_rho, self.rho_r, self._rphi_r
)
def _read_eam_alloy(self) -> None:
"""
Read EAM potential file in LAMMPS alloy format.
Each data section (F(rho), rho(r), r*phi(r)) is parsed starting on its
own line. Any trailing tokens on the line where a section ends are
discarded — matching LAMMPS' ``pair_eam_alloy`` reader and GPUMD's
``eam_alloy`` reader. This avoids the off-by-one shift that would occur
if values from the padding past the end of one section were fed into
the next section.
"""
with open(self.filename) as f:
lines = f.readlines()
self.header = lines[:3]
line4 = lines[3].split()
self.Nelements = int(line4[0])
self.elements_list = line4[1 : 1 + self.Nelements]
# nrho, drho, nr, dr, rc
line5 = lines[4].split()
self.nrho = int(line5[0])
self.drho = float(line5[1])
self.nr = int(line5[2])
self.dr = float(line5[3])
self.rc = float(line5[4])
# Initialize storage arrays
self.F_rho = np.zeros((self.Nelements, self.nrho))
self.rho_r = np.zeros((self.Nelements, self.nr))
self.phi_r = np.zeros((self.Nelements, self.Nelements, self.nr))
self.r = np.arange(0, self.nr) * self.dr
self.rho = np.arange(0, self.nrho) * self.drho
line_idx = [5]
def read_section(n: int) -> NDArray[np.float64]:
"""Read ``n`` floats starting at ``lines[line_idx[0]]``; discard
any leftover tokens on the last touched line."""
out = np.empty(n, dtype=np.float64)
count = 0
while count < n and line_idx[0] < len(lines):
toks = lines[line_idx[0]].split("#")[0].split()
for t in toks:
if count >= n:
break
try:
out[count] = float(t)
count += 1
except ValueError:
break
line_idx[0] += 1
if count < n:
raise ValueError(
f"EAM file '{self.filename}' ended while expecting {n} "
f"values (got {count})."
)
return out
for element in range(self.Nelements):
# skip the per-element info line (atomic number, mass, ...)
line_idx[0] += 1
self.F_rho[element] = read_section(self.nrho)
self.rho_r[element] = read_section(self.nr)
# r*phi(r) stored in the lower triangle (j <= i); mirror to j > i.
self._rphi_r = np.zeros_like(self.phi_r)
for element_i in range(self.Nelements):
for element_j in range(element_i + 1):
self._rphi_r[element_i, element_j] = read_section(self.nr)
if element_i != element_j:
self._rphi_r[element_j, element_i] = self._rphi_r[
element_i, element_j
]
# Convert r*phi to phi (avoid division by zero at r=0)
self.phi_r[:, :, 1:] = self._rphi_r[:, :, 1:] / self.r[1:]
self.phi_r[:, :, 0] = self.phi_r[:, :, 1]
[docs]
def write_eam_alloy(self, output_name: Optional[str] = None):
"""Write to an eam.alloy file.
Args:
output_name (Optional[str], optional): output filename, it will automatically generated if not given.
"""
if output_name is None:
output_name = ""
for i in self.elements_list:
output_name += i
output_name += ".new.eam.alloy"
with open(output_name, "w") as op:
header = self.header.copy()
h0 = f" eam/alloy {self.Nelements}"
for i in self.elements_list:
h0 += f" {i}"
h0 += "\n"
header[0] = h0
header[1] = " Generated by mdapy!\n"
header[2] = "\n"
op.write("".join(header))
op.write(f" {self.Nelements} ")
for i in self.elements_list:
op.write(f"{i} ")
op.write("\n")
op.write(f" {self.nrho} {self.drho} {self.nr} {self.dr} {self.rc}\n")
num = 1
colnum = 5
for i in range(self.Nelements):
op.write(f" 0 {self.mass_list[i]} 0 dummy\n ")
num = 1
for j in range(self.nrho):
op.write(f"{self.F_rho[i, j]:.16E} ")
if num >= colnum:
op.write("\n ")
num = 0
num += 1
if num > 1:
op.write("\n ")
num = 1
for j in range(self.nr):
op.write(f"{self.rho_r[i, j]:.16E} ")
if num >= colnum:
op.write("\n ")
num = 0
num += 1
if num > 1:
op.write("\n ")
num = 1
for i1 in range(self.Nelements):
for i2 in range(i1 + 1):
for j in range(self.nr):
op.write(f"{self._rphi_r[i1, i2, j]:.16E} ")
if num >= colnum:
op.write("\n ")
num = 0
num += 1
if num > 1:
op.write("\n")
[docs]
def get_energies(self, data: pl.DataFrame, box: Box) -> NDArray[np.float64]:
"""
Calculate per-atom potential energies.
Parameters
----------
data : pl.DataFrame
Atomic configuration containing columns: 'x', 'y', 'z', 'element'.
box : Box
Simulation box object containing boundary conditions and dimensions.
Returns
-------
NDArray[np.float64]
Per-atom potential energies in eV. Shape: (N,) where N is the number of atoms.
"""
if "energies" not in self.results.keys():
self._calculate(data, box)
return self.results["energies"]
[docs]
def get_energy(self, data: pl.DataFrame, box: Box) -> float:
"""
Calculate total system potential energy.
Parameters
----------
data : pl.DataFrame
Atomic configuration containing columns: 'x', 'y', 'z', 'element'.
box : Box
Simulation box object containing boundary conditions and dimensions.
Returns
-------
float
Total potential energy in eV.
"""
return self.get_energies(data, box).sum()
[docs]
def get_forces(self, data: pl.DataFrame, box: Box) -> NDArray[np.float64]:
"""
Calculate atomic forces.
Parameters
----------
data : pl.DataFrame
Atomic configuration containing columns: 'x', 'y', 'z', 'element'.
box : Box
Simulation box object containing boundary conditions and dimensions.
Returns
-------
NDArray[np.float64]
Atomic forces in eV/Å. Shape: (N, 3) where N is the number of atoms.
Columns represent [Fx, Fy, Fz].
"""
if "forces" not in self.results.keys():
self._calculate(data, box)
return self.results["forces"]
[docs]
def get_stress(self, data: pl.DataFrame, box: Box) -> NDArray[np.float64]:
"""
Calculate stress tensor in Voigt notation.
Parameters
----------
data : pl.DataFrame
Atomic configuration containing columns: 'x', 'y', 'z', 'element'.
box : Box
Simulation box object containing boundary conditions and dimensions.
Returns
-------
NDArray[np.float64]
Stress tensor in Voigt notation [σ_xx, σ_yy, σ_zz, σ_yz, σ_xz, σ_xy].
Units: eV/ų or GPa depending on unit conversion.
Shape: (6,).
"""
if "stress" not in self.results.keys():
self._calculate(data, box)
return self.results["stress"]
[docs]
def get_virials(self, data: pl.DataFrame, box: Box) -> NDArray[np.float64]:
"""
Calculate per-atom virial tensors.
Parameters
----------
data : pl.DataFrame
Atomic configuration containing columns: 'x', 'y', 'z', 'element'.
box : Box
Simulation box object containing boundary conditions and dimensions.
Returns
-------
NDArray[np.float64]
Per-atom virial tensors. Shape: (N, 9) where N is the number of atoms.
Components: [W_xx, W_xy, W_xz, W_yx, W_yy, W_yz, W_zx, W_zy, W_zz].
Units: eV.
"""
if "virials" not in self.results.keys():
self._calculate(data, box)
return self.results["virials"]
def _calculate(self, data: pl.DataFrame, box: Box) -> None:
"""
Main calculation method for energies, forces, stress, and virials.
This method performs the complete EAM calculation:
1. Validates input data
2. Constructs neighbor lists
3. Computes electron densities at each atom
4. Evaluates embedding energies
5. Calculates pair interactions
6. Computes forces and virials
7. Derives stress tensor
Parameters
----------
data : pl.DataFrame
Atomic configuration with required columns: 'x', 'y', 'z', 'element'.
box : Box
Simulation box with boundary conditions.
Raises
------
AssertionError
If required columns are missing from data.
AssertionError
If element types in data are not supported by the potential.
Notes
-----
Results are stored in `self.results` dictionary with keys:
- 'energies': Per-atom potential energies
- 'forces': Atomic forces
- 'stress': Global stress tensor (Voigt notation)
- 'virials': Per-atom virial tensors
"""
for i in ["x", "y", "z", "element"]:
assert i in data.columns, f"Required column '{i}' not found in data."
for i in data["element"].unique():
assert (
i in self.elements_list
), f"{i} not in current EAM potential supported elements: {self.elements_list}."
old_N = data.shape[0]
ele2type = {j: i for i, j in enumerate(self.elements_list)}
data = data.with_columns(
pl.col("element")
.replace_strict(ele2type, return_dtype=pl.Int32)
.rechunk()
.alias("type")
)
repeat = box.check_small_box(self.rc)
if sum(repeat) != 3:
data, box = tool.replicate(data, box, *repeat)
type_list = data["type"].to_numpy(allow_copy=False)
neigh = Neighbor(self.rc, box, data, max_neigh=200)
neigh.compute()
N = data.shape[0] # Number of atoms
# Allocate arrays for outputs
potential = np.zeros(N, float) # Per-atom energies
force = np.zeros((N, 3), float) # Per-atom forces [fx, fy, fz]
virial = np.zeros((N, 9), float) # Per-atom virials (9 components)
self._eam.calculate(
data["x"].to_numpy(allow_copy=False),
data["y"].to_numpy(allow_copy=False),
data["z"].to_numpy(allow_copy=False),
type_list,
box.box,
box.origin,
box.boundary,
neigh.verlet_list,
neigh.distance_list,
neigh.neighbor_number,
force,
virial,
potential,
get_num_threads(),
)
# Store results
self.results["energies"] = potential[:old_N]
self.results["forces"] = force[:old_N]
self.results["virials"] = virial[:old_N]
# Calculate stress tensor from virials
v = virial.sum(axis=0) # Sum virials over all atoms
# Reshape to 3×3 matrix: v_xx, v_xy, v_xz, v_yx, v_yy, v_yz, v_zx, v_zy, v_zz
v = v.reshape(3, 3)
# Stress = -(virial + virial^T) / (2 * volume)
stress = (-0.5 * (v + v.T) / box.volume).ravel()
# Convert to Voigt notation: [σ_xx, σ_yy, σ_zz, σ_yz, σ_xz, σ_xy]
stress = stress[[0, 4, 8, 5, 2, 1]]
self.results["stress"] = stress
[docs]
def plot(
self,
fig: Optional[Figure] = None,
ax: Optional[List[Axes]] = None,
) -> Tuple[Figure, List[Axes]]:
"""
Plot EAM potential functions: embedding function, electron density, and pair potential.
Creates three subplots showing:
1. Embedding function F(ρ) vs electron density ρ
2. Electron density function ρ(r) vs distance r
3. Pair potential φ(r) vs distance r
Parameters
----------
fig : Figure, optional
Matplotlib figure object. If None, creates a new figure.
ax : List[Axes], optional
List of three matplotlib axes objects. If None, creates new axes.
Returns
-------
fig : Figure
Matplotlib figure object containing the plots.
ax : List[Axes]
List of three matplotlib axes objects [ax_F, ax_rho, ax_phi].
"""
from mdapy.plotset import set_figure
if fig is None and ax is None:
fig, ax = set_figure(ncol=3, figsize=(18, 6))
else:
assert len(ax) == 3, "ax must be a list of 3 axes objects"
# Plot data for each element
for i, element in enumerate(self.elements_list):
ax[0].plot(self.rho, self.F_rho[i], label=element)
ax[1].plot(self.r, self.rho_r[i], label=element)
ax[2].plot(self.r, self.phi_r[i, i], label=f"{element}-{element}")
# Set labels
ax[0].set_xlabel("Electron Density ρ")
ax[0].set_ylabel("Embedding Function F(ρ) (eV)")
ax[1].set_xlabel("Distance r (Å)")
ax[1].set_ylabel("Electron Density ρ(r)")
ax[2].set_xlabel("Distance r (Å)")
ax[2].set_ylabel("Pair Potential φ(r) (eV)")
ax[0].set_xlim(0, self.rho[-1])
ax[1].set_xlim(0, self.rc)
ax[2].set_xlim(0, self.rc)
for i in range(3):
ax[i].legend()
return fig, ax
[docs]
class EAMAverage(EAM):
"""
Generate average EAM (A-atom) potential for alloy investigation.
This class creates an average EAM potential, which is useful in alloy investigation.
The A-atom potential has a similar formula to the original EAM potential:
.. math:: E_{i}=\\sum_{i} F^{A}\\left(\\bar{\\rho}_{i}\\right)+\\frac{1}{2} \\sum_{i, j \\neq i} \\phi_{ij}^{A A},
.. math:: F^{A}\\left(\\bar{\\rho}_{i}\\right)=\\sum_{\\alpha} c_{\\alpha} F^{\\alpha}\\left(\\bar{\\rho}_{i}\\right),
.. math:: \\phi_{ij}^{A A}=\\sum_{\\alpha, \\beta } c_{\\alpha} c_{\\beta} \\phi_{ij}^{\\alpha \\beta},
.. math:: \\quad \\bar{\\rho}_{i}=\\sum_{j \\neq i} \\sum_{\\alpha} c_{\\alpha} \\rho_{ij}^{\\alpha},
where :math:`A` denotes an average-atom.
.. note:: If you use this module in publication, you should also cite the original paper.
`Average-atom interatomic potential for random alloys <https://doi.org/10.1103/PhysRevB.93.104201>`_
Parameters
----------
filename : str
Filename of eam.alloy file.
concentration : List[float]
Atomic ratio list, such as [0.5, 0.5]. The summation should equal 1.
output_name : str, optional
Filename of generated average EAM potential. If None, auto-generates name.
Examples
--------
>>> potential = EAMAverage("CuNiAl.eam.alloy", [0.25, 0.25, 0.5])
"""
def __init__(
self,
filename: str,
concentration: List[float],
output_name: Optional[str] = None,
) -> None:
super().__init__(filename)
self.concentration = concentration
assert (
len(self.concentration) == self.Nelements
), f"Number of concentration list should be equal to {self.Nelements}."
assert np.isclose(
np.sum(concentration), 1.0
), "Concentration summation should be equal to 1."
# Perform averaging
self._average()
# Determine output filename
if output_name is None:
self.output_name = ""
for i in self.elements_list[:-1]: # Exclude 'A' element
self.output_name += i
self.output_name += ".average.eam.alloy"
else:
self.output_name = output_name
# Write the averaged potential
self.write_eam_alloy(self.output_name)
def _average(self) -> None:
"""
Average the EAM potential arrays according to concentration.
This method:
1. Adds an 'A' element to represent the average atom
2. Averages embedding and electron density functions
3. Averages pair potentials with proper weighting
"""
# Calculate average mass for 'A' element
avg_mass = np.average(self.mass_list, weights=self.concentration)
# Update element count and list
self.Nelements += 1
self.elements_list.append("A")
self.mass_list.append(avg_mass)
# Average embedded_data and elec_density_data
new_F_rho = np.r_[
self.F_rho,
np.zeros((1, self.F_rho.shape[1]), dtype=self.F_rho.dtype),
]
new_rho_r = np.r_[
self.rho_r,
np.zeros((1, self.rho_r.shape[1]), dtype=self.rho_r.dtype),
]
new_F_rho[-1, :] = np.average(self.F_rho, axis=0, weights=self.concentration)
new_rho_r[-1, :] = np.average(self.rho_r, axis=0, weights=self.concentration)
# Average rphi_data
new_rphi_r = np.concatenate(
(
self._rphi_r,
np.zeros(
(self._rphi_r.shape[0], 1, self._rphi_r.shape[2]),
dtype=self._rphi_r.dtype,
),
),
axis=1,
)
new_rphi_r = np.concatenate(
(
new_rphi_r,
np.zeros(
(1, new_rphi_r.shape[1], new_rphi_r.shape[2]),
dtype=new_rphi_r.dtype,
),
),
axis=0,
)
# Average for A-X interactions
new_rphi_r[-1, :-1, :] = np.average(
self._rphi_r, axis=0, weights=self.concentration
)
new_rphi_r[:-1, -1, :] = new_rphi_r[-1, :-1, :]
# Average for A-A interaction
column = new_rphi_r[:-1, -1, :]
new_rphi_r[-1, -1, :] = np.average(column, axis=0, weights=self.concentration)
# Convert rphi to phi
new_phi_r = np.zeros_like(new_rphi_r)
new_phi_r[:, :, 1:] = new_rphi_r[:, :, 1:] / self.r[1:]
new_phi_r[:, :, 0] = new_phi_r[:, :, 1]
# Update attributes
self.F_rho = new_F_rho
self.rho_r = new_rho_r
self._rphi_r = new_rphi_r
self.phi_r = new_phi_r
[docs]
class EAMGenerator:
"""
Generator for EAM (Embedded Atom Method) potential files.
This class creates EAM potential files in the eam.alloy format for single
or multi-element systems. The implementation is based on the Zhou-Johnson-Wadley
parameterization.
Parameters
----------
elements_list : List[str]
List of element symbols (e.g., ['Cu', 'Ni', 'Al']).
output_filename : str, optional
Name of the output file. If None, generates name from elements.
nr : int, optional
Number of points in distance grid (default: 2000).
nrho : int, optional
Number of points in electron density grid (default: 2000).
rst : float, optional
Minimum distance for pair potential computation (default: 0.5 Å).
Raises
------
ValueError
If any element is not in the supported elements list.
References
----------
Zhou, X. W., Johnson, R. A., & Wadley, H. N. G. (2004).
Misfit-energy-increasing dislocations in vapor-deposited CoFe/NiFe multilayers.
Physical Review B, 69(14), 144113.
https://doi.org/10.1103/PhysRevB.69.144113
Examples
--------
>>> gen = EAMGenerator(["Cu", "Ni"])
>>> # This automatically generates a 'CuNi.eam.alloy' file
"""
# Class constants
SUPPORTED_ELEMENTS: Tuple[str, ...] = (
"Cu",
"Ag",
"Au",
"Ni",
"Pd",
"Pt",
"Al",
"Pb",
"Fe",
"Mo",
"Ta",
"W",
"Mg",
"Co",
"Ti",
"Zr",
)
DEFAULT_NR: int = 2000
DEFAULT_NRHO: int = 2000
DEFAULT_RST: float = 0.5
# Raw parameter data (Zhou et al. parameters)
# 27 parameters per element in this order:
# re, fe, rhoe, rhos, alpha, beta, A, B, kappa, lambda,
# Fm0, Fm1, Fm2, Fm3, Fm4, Fi0, Fi1, Fi2, Fi3, Fe,
# atomic_number, atomic_mass, eta, beta1, lambda1, rhol, rhoh
_RAW_PARAMETERS: Tuple[str, ...] = (
"Cu",
"2.556162",
"1.554485",
"21.175871",
"21.175395",
"8.127620",
"4.334731",
"0.396620",
"0.548085",
"0.308782",
"0.756515",
"-2.170269",
"-0.263788",
"1.088878",
"-0.817603",
"-2.19",
"0.00",
"0.561830",
"-2.100595",
"0.310490",
"-2.186568",
"29",
"63.546",
"-2.100595",
"4.334731",
"0.756515",
"0.85",
"1.15",
"Ag",
"2.891814",
"1.106232",
"14.604100",
"14.604144",
"9.132010",
"4.870405",
"0.277758",
"0.419611",
"0.339710",
"0.750758",
"-1.729364",
"-0.255882",
"0.912050",
"-0.561432",
"-1.75",
"0.00",
"0.744561",
"-1.150650",
"0.783924",
"-1.748423",
"47",
"107.8682",
"-1.150650",
"4.870405",
"0.750758",
"0.85",
"1.15",
"Au",
"2.885034",
"1.529021",
"19.991632",
"19.991509",
"9.516052",
"5.075228",
"0.229762",
"0.356666",
"0.356570",
"0.748798",
"-2.937772",
"-0.500288",
"1.601954",
"-0.835530",
"-2.98",
"0.00",
"1.706587",
"-1.134778",
"1.021095",
"-2.978815",
"79",
"196.96654",
"-1.134778",
"5.075228",
"0.748798",
"0.85",
"1.15",
"Ni",
"2.488746",
"2.007018",
"27.562015",
"27.562031",
"8.383453",
"4.471175",
"0.429046",
"0.633531",
"0.443599",
"0.820658",
"-2.693513",
"-0.076445",
"0.241442",
"-2.375626",
"-2.70",
"0.00",
"0.265390",
"-0.152856",
"0.445470",
"-2.7",
"28",
"58.6934",
"-0.152856",
"4.471175",
"0.820658",
"0.85",
"1.15",
"Pd",
"2.750897",
"1.595417",
"21.335246",
"21.940073",
"8.697397",
"4.638612",
"0.406763",
"0.598880",
"0.397263",
"0.754799",
"-2.321006",
"-0.473983",
"1.615343",
"-0.231681",
"-2.36",
"0.00",
"1.481742",
"-1.675615",
"1.130000",
"-2.352753",
"46",
"106.42",
"-1.675615",
"4.638612",
"0.754799",
"0.85",
"1.15",
"Pt",
"2.771916",
"2.336509",
"33.367564",
"35.205357",
"7.105782",
"3.789750",
"0.556398",
"0.696037",
"0.385255",
"0.770510",
"-1.455568",
"-2.149952",
"0.528491",
"1.222875",
"-4.17",
"0.00",
"3.010561",
"-2.420128",
"1.450000",
"-4.145597",
"78",
"195.08",
"-2.420128",
"3.789750",
"0.770510",
"0.25",
"1.15",
"Al",
"2.863924",
"1.403115",
"20.418205",
"23.195740",
"6.613165",
"3.527021",
"0.314873",
"0.365551",
"0.379846",
"0.759692",
"-2.807602",
"-0.301435",
"1.258562",
"-1.247604",
"-2.83",
"0.00",
"0.622245",
"-2.488244",
"0.785902",
"-2.824528",
"13",
"26.981539",
"-2.488244",
"3.527021",
"0.759692",
"0.85",
"1.15",
"Pb",
"3.499723",
"0.647872",
"8.450154",
"8.450063",
"9.121799",
"5.212457",
"0.161219",
"0.236884",
"0.250805",
"0.764955",
"-1.422370",
"-0.210107",
"0.682886",
"-0.529378",
"-1.44",
"0.00",
"0.702726",
"-0.538766",
"0.935380",
"-1.439436",
"82",
"207.2",
"-0.538766",
"5.212457",
"0.764955",
"0.85",
"1.15",
"Fe",
"2.481987",
"1.885957",
"20.041463",
"20.041463",
"9.818270",
"5.236411",
"0.392811",
"0.646243",
"0.170306",
"0.340613",
"-2.534992",
"-0.059605",
"0.193065",
"-2.282322",
"-2.54",
"0.00",
"0.200269",
"-0.148770",
"0.391750",
"-2.539945",
"26",
"55.847",
"-0.148770",
"5.236411",
"0.340613",
"0.85",
"1.15",
"Mo",
"2.728100",
"2.723710",
"29.354065",
"29.354065",
"8.393531",
"4.476550",
"0.708787",
"1.120373",
"0.137640",
"0.275280",
"-3.692913",
"-0.178812",
"0.380450",
"-3.133650",
"-3.71",
"0.00",
"0.875874",
"0.776222",
"0.790879",
"-3.712093",
"42",
"95.94",
"0.776222",
"4.476550",
"0.275280",
"0.85",
"1.15",
"Ta",
"2.860082",
"3.086341",
"33.787168",
"33.787168",
"8.489528",
"4.527748",
"0.611679",
"1.032101",
"0.176977",
"0.353954",
"-5.103845",
"-0.405524",
"1.112997",
"-3.585325",
"-5.14",
"0.00",
"1.640098",
"0.221375",
"0.848843",
"-5.141526",
"73",
"180.9479",
"0.221375",
"4.527748",
"0.353954",
"0.85",
"1.15",
"W",
"2.740840",
"3.487340",
"37.234847",
"37.234847",
"8.900114",
"4.746728",
"0.882435",
"1.394592",
"0.139209",
"0.278417",
"-4.946281",
"-0.148818",
"0.365057",
"-4.432406",
"-4.96",
"0.00",
"0.661935",
"0.348147",
"0.582714",
"-4.961306",
"74",
"183.84",
"0.348147",
"4.746728",
"0.278417",
"0.85",
"1.15",
"Mg",
"3.196291",
"0.544323",
"7.132600",
"7.132600",
"10.228708",
"5.455311",
"0.137518",
"0.225930",
"0.5",
"1.0",
"-0.896473",
"-0.044291",
"0.162232",
"-0.689950",
"-0.90",
"0.00",
"0.122838",
"-0.226010",
"0.431425",
"-0.899702",
"12",
"24.305",
"-0.226010",
"5.455311",
"1.0",
"0.85",
"1.15",
"Co",
"2.505979",
"1.975299",
"27.206789",
"27.206789",
"8.679625",
"4.629134",
"0.421378",
"0.640107",
"0.5",
"1.0",
"-2.541799",
"-0.219415",
"0.733381",
"-1.589003",
"-2.56",
"0.00",
"0.705845",
"-0.687140",
"0.694608",
"-2.559307",
"27",
"58.9332",
"-0.687140",
"4.629134",
"1.0",
"0.85",
"1.15",
"Ti",
"2.933872",
"1.863200",
"25.565138",
"25.565138",
"8.775431",
"4.680230",
"0.373601",
"0.570968",
"0.5",
"1.0",
"-3.203773",
"-0.198262",
"0.683779",
"-2.321732",
"-3.22",
"0.00",
"0.608587",
"-0.750710",
"0.558572",
"-3.219176",
"22",
"47.88",
"-0.750710",
"4.680230",
"1.0",
"0.85",
"1.15",
"Zr",
"3.199978",
"2.230909",
"30.879991",
"30.879991",
"8.559190",
"4.564902",
"0.424667",
"0.640054",
"0.5",
"1.0",
"-4.485793",
"-0.293129",
"0.990148",
"-3.202516",
"-4.51",
"0.00",
"0.928602",
"-0.981870",
"0.597133",
"-4.509025",
"40",
"91.224",
"-0.981870",
"4.564902",
"1.0",
"0.85",
"1.15",
)
def __init__(
self,
elements_list: List[str],
output_filename: Optional[str] = None,
nr: int = DEFAULT_NR,
nrho: int = DEFAULT_NRHO,
rst: float = DEFAULT_RST,
) -> None:
# Validate input elements
for element in elements_list:
if element not in self.SUPPORTED_ELEMENTS:
raise ValueError(
f"Element '{element}' is not supported. "
f"Supported elements: {', '.join(self.SUPPORTED_ELEMENTS)}"
)
self.elements_list = elements_list
self.n_elements = len(elements_list)
self.nr = nr
self.nrho = nrho
self.rst = rst
# Output filename
if output_filename is None:
self.output_filename = self._generate_filename()
else:
self.output_filename = output_filename
# Initialize parameter storage
self._get_eam_parameters()
# Write to file
self._write_eam_file()
def _diedai(self):
for i1 in range(self.ntypes):
for i2 in range(i1 + 1):
if i1 == i2:
for i in range(self.nr):
r = i * self.dr
if r < self.rst:
r = self.rst
fvalue = self._prof(i1, r)
if self.fmax < fvalue:
self.fmax = fvalue
self.rho[i, i1] = fvalue
psi = self._pair(i1, i2, r)
self.rphi[i, i1, i2] = r * psi
else:
for i in range(self.nr):
r = i * self.dr
if r < self.rst:
r = self.rst
psi = self._pair(i1, i2, r)
self.rphi[i, i1, i2] = r * psi
self.rphi[i, i2, i1] = self.rphi[i, i1, i2]
rhom = self.fmax
if rhom < 2.0 * self.rhoemax:
rhom = 2.0 * self.rhoemax
if rhom < 100.0:
rhom = 100.0
self.drho = rhom / (self.nrho - 1.0)
for it in range(self.ntypes):
for i in range(self.nrho):
rhoF = i * self.drho
self.embedding_function[i, it] = self._embed(it, rhoF)
def _prof(self, it, r):
f = self.fe[it] * np.exp(-self.beta1[it] * (r / self.re[it] - 1.0))
f = f / (1.0 + (r / self.re[it] - self.ramda1[it]) ** 20)
return f
def _pair(self, it1, it2, r):
if it1 == it2:
psi1 = self.A[it1] * np.exp(-self.alpha[it1] * (r / self.re[it1] - 1.0))
psi1 = psi1 / (1.0 + (r / self.re[it1] - self.cai[it1]) ** 20)
psi2 = self.B[it1] * np.exp(-self.beta[it1] * (r / self.re[it1] - 1.0))
psi2 = psi2 / (1.0 + (r / self.re[it1] - self.ramda[it1]) ** 20)
psi = psi1 - psi2
else:
psiab, fab = [], []
for it in [it1, it2]:
psi1 = self.A[it] * np.exp(-self.alpha[it] * (r / self.re[it] - 1.0))
psi1 = psi1 / (1.0 + (r / self.re[it] - self.cai[it]) ** 20)
psi2 = self.B[it] * np.exp(-self.beta[it] * (r / self.re[it] - 1.0))
psi2 = psi2 / (1.0 + (r / self.re[it] - self.ramda[it]) ** 20)
psiab.append(psi1 - psi2)
fab.append(self._prof(it, r))
psi = 0.5 * (fab[1] / fab[0] * psiab[0] + fab[0] / fab[1] * psiab[1])
return psi
def _embed(self, it, rho):
if rho < self.rhoe[it]:
Fm33 = self.Fm3[it]
else:
Fm33 = self.Fm4[it]
if rho < self.rhoin[it]:
emb = (
self.Fi0[it]
+ self.Fi1[it] * (rho / self.rhoin[it] - 1.0)
+ self.Fi2[it] * (rho / self.rhoin[it] - 1.0) ** 2
+ self.Fi3[it] * (rho / self.rhoin[it] - 1.0) ** 3
)
elif rho < self.rhoout[it]:
emb = (
self.Fm0[it]
+ self.Fm1[it] * (rho / self.rhoe[it] - 1.0)
+ self.Fm2[it] * (rho / self.rhoe[it] - 1.0) ** 2
+ Fm33 * (rho / self.rhoe[it] - 1.0) ** 3
)
else:
emb = (
self.Fn[it]
* (1.0 - self.fnn[it] * np.log(rho / self.rhos[it]))
* (rho / self.rhos[it]) ** self.fnn[it]
)
return emb
def _get_eam_parameters(self):
name, data = [], []
a = 0
for _ in range(16):
name.append(self._RAW_PARAMETERS[a].strip())
data.append(self._RAW_PARAMETERS[a + 1 : a + 28])
a += 28
self.data = pl.from_numpy(np.array(data, dtype=np.float64).T, schema=name)
(
self.re,
self.fe,
self.rhoe,
self.rhos,
self.alpha,
self.beta,
self.A,
self.B,
self.cai,
self.ramda,
self.Fi0,
self.Fi1,
self.Fi2,
self.Fi3,
self.Fm0,
self.Fm1,
self.Fm2,
self.Fm3,
self.fnn,
self.Fn,
self.atomic_number,
self.atomic_mass,
self.Fm4,
self.beta1,
self.ramda1,
self.rhol,
self.rhoh,
) = self.data.select(self.elements_list).to_numpy()
self.atomic_number = np.array(self.atomic_number, np.int32)
self.lattice_constant = np.sqrt(2.0) * self.re
self.rhoin = self.rhol * self.rhoe
self.rhoout = self.rhoh * self.rhoe
self.alatmax = self.lattice_constant.max()
self.rhoemax = self.rhoe.max()
self.ntypes = len(self.elements_list)
self.rc = np.sqrt(10.0) / 2.0 * self.alatmax
self.dr = self.rc / (self.nr - 1.0)
self.fmax = -1.0
self.rho = np.zeros((self.nrho, len(self.elements_list)))
self.rphi = np.zeros(
(self.nr, len(self.elements_list), len(self.elements_list))
)
self.embedding_function = np.zeros((self.nr, len(self.elements_list)))
self._diedai()
def _generate_filename(self) -> str:
"""Generate output filename from element list."""
element_string = "".join(self.elements_list)
return f"{element_string}.eam.alloy"
def _write_eam_file(self) -> None:
"""Write the EAM potential to file in eam.alloy format."""
crystal_structure = "fcc"
with open(self.output_filename, "w") as file:
# Header section
self._write_header(file)
# Grid parameters
file.write(
f" {self.nrho} {self.drho:.16E} "
f"{self.nr} {self.dr:.16E} {self.rc:.16E}\n"
)
# Element-specific data
for i in range(self.n_elements):
self._write_element_data(file, i, crystal_structure)
# Pair potentials
self._write_pair_potentials(file)
def _write_header(self, file) -> None:
"""Write file header with metadata."""
# First line with element count and symbols
file.write(f" eam/alloy {self.n_elements}")
for element in self.elements_list:
file.write(f" {element}")
file.write("\n")
# Citation and generation info
file.write(
f" Python version generated by EAMGenerator! Based on the Fortran version by Zhou et al. Generated: {datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S')}\n"
)
file.write(
" CITATION: X. W. Zhou, R. A. Johnson, H. N. G. Wadley, "
"Phys. Rev. B, 69, 144113 (2004)\n"
)
# Element list (repeated)
file.write(f" {self.n_elements} ")
for element in self.elements_list:
file.write(f"{element} ")
file.write("\n")
def _write_element_data(
self, file, element_idx: int, crystal_structure: str
) -> None:
"""Write embedding function and electron density for one element."""
file.write(
f" {self.atomic_number[element_idx]} "
f"{self.atomic_mass[element_idx]:.10f} "
f"{self.lattice_constant[element_idx]:.6f} "
f"{crystal_structure}\n"
)
self._write_array_data(file, self.embedding_function[:, element_idx])
self._write_array_data(file, self.rho[:, element_idx])
def _write_pair_potentials(self, file) -> None:
"""Write pair potential data for all element pairs."""
for i in range(self.n_elements):
for j in range(i + 1):
self._write_array_data(file, self.rphi[:, i, j])
def _write_array_data(
self, file, data: npt.NDArray[np.float64], values_per_line: int = 5
) -> None:
"""Write array data to file with specified formatting."""
for idx, value in enumerate(data):
if idx % values_per_line == 0:
if idx > 0:
file.write("\n")
file.write(" ")
file.write(f"{value:.16E} ")
# Ensure we end on a new line
file.write("\n")
if __name__ == "__main__":
# import os
# EAMGenerator(
# ["Co", "Ni", "Fe", "Al", "Cu"], "tests/input_files/CoNiFeAlCu.eam.alloy"
# )
# EAMAverage("CuNiAl.eam.alloy", [0.25, 0.25, 0.5])
# eam = EAM("CuNiAl.average.eam.alloy")
eam = EAM(r"tests/input_files/NiCoCr.lammps.eam")
eam.plot()
import matplotlib.pyplot as plt
plt.show()
eam2 = EAMAverage(r"tests/input_files/NiCoCr.lammps.eam", [0.2, 0.2, 0.6])
eam2.plot()
plt.show()
eam1 = EAM(r"tests/input_files/FeNiCrCoTi-heamix.setfl")
eam1.plot()
plt.show()
# os.system("rm CuNiAl*.eam.alloy")