Hongyu-yu/materials-properties

Materials Properties Calculation Skill

You are an expert assistant for calculating materials properties from first-principles using the Atomic Simulation Environment (ASE) and specialized packages. Help users perform structure relaxations, compute ground state properties, and calculate advanced materials properties using scientifically rigorous methods with proper citations.

Overview

This skill covers comprehensive materials property calculations including:

Core Capabilities:

• Structure relaxation (geometry optimization, stress relaxation)
• Ground state properties (lattice constants, space groups, crystal structure)
• Calculator setup (EMT, GPAW, VASP, Quantum ESPRESSO)

Advanced Properties:

• Surface energies
• Adsorption energies
• Reaction barriers (NEB)
• Vibrational analysis
• Phonon calculations (phonopy)
• Elastic constants (elastic)
• Cluster expansions (icet)
• CALPHAD integration (pycalphad)
• Defect formation energies
• Interface/grain boundary energies
• Magnetic properties
• Thermal expansion
• Electronic structure

Installation

# Core packages
pip install ase spglib matplotlib

# Specialized packages
pip install phonopy elastic icet pycalphad

# Optional DFT calculators
pip install gpaw  # Real DFT (requires compilation)
# VASP requires separate license and installation
# Quantum ESPRESSO via ase.calculators.espresso

Core Workflow: Structure Relaxation

1. Basic Geometry Optimization

from ase import Atoms
from ase.optimize import BFGS
from ase.calculators.emt import EMT

# Create or load structure
atoms = Atoms('Cu', positions=[[0, 0, 0]], cell=[2.5, 2.5, 2.5], pbc=True)

# Set calculator
atoms.calc = EMT()

# Optimize geometry
opt = BFGS(atoms, trajectory='opt.traj')
opt.run(fmax=0.01)  # Force convergence criterion

# Get optimized energy
E_opt = atoms.get_potential_energy()
print(f"Optimized energy: {E_opt:.3f} eV")

Key Parameters:

• fmax: Maximum force (eV/Å) - typical values: 0.01-0.05
• steps: Maximum optimization steps
• Optimizers: BFGS, LBFGS, FIRE, GPMin

2. Unit Cell Optimization (Stress Relaxation)

from ase.optimize import BFGS
from ase.constraints import ExpCellFilter

# Relax both positions and cell
ecf = ExpCellFilter(atoms)
opt = BFGS(ecf, trajectory='cell_opt.traj')
opt.run(fmax=0.01)

# Get optimized lattice
a = atoms.cell.cellpar()[0]
print(f"Optimized lattice constant: {a:.3f} Å")

Applications:

• Lattice constant determination
• Equation of state calculations
• Pressure-dependent structures

3. Ground State Properties

import spglib

# Get space group
cell = (atoms.cell, atoms.get_scaled_positions(), atoms.get_atomic_numbers())
spacegroup = spglib.get_spacegroup(cell, symprec=1e-5)
print(f"Space group: {spacegroup}")

# Get lattice parameters
a, b, c, alpha, beta, gamma = atoms.cell.cellpar()
print(f"Lattice: a={a:.3f}, b={b:.3f}, c={c:.3f} Å")
print(f"Angles: α={alpha:.1f}, β={beta:.1f}, γ={gamma:.1f}°")

# Volume
V = atoms.get_volume()
print(f"Volume: {V:.3f} ų")

# Crystal system
dataset = spglib.get_symmetry_dataset(cell)
print(f"Crystal system: {dataset['international']}")

Calculator Setup

EMT Calculator (Fast, for Testing)

from ase.calculators.emt import EMT

atoms.calc = EMT()

Advantages:

• Very fast
• No installation issues
• Good for workflow development

Limitations:

• Only works for a few elements (Cu, Ag, Au, Ni, Pd, Pt, Al, Pb, Fe)
• Approximate potential

GPAW Calculator (Real DFT)

from gpaw import GPAW, PW

atoms.calc = GPAW(mode=PW(500),  # Plane-wave cutoff (eV)
                 xc='PBE',       # Exchange-correlation functional
                 kpts=(8, 8, 8), # k-point sampling
                 txt='gpaw.txt') # Output file

Parameters:

• mode: PW(cutoff) for plane waves, grid-based for real-space
• xc: PBE, LDA, RPBE, BEEF-vdW, etc.
• kpts: k-point mesh or specific k-points
• convergence: Energy convergence criterion

VASP Calculator

from ase.calculators.vasp import Vasp

atoms.calc = Vasp(xc='PBE',
                 encut=500,     # Cutoff energy (eV)
                 kpts=(8,8,8),
                 ismear=1,      # Smearing method
                 sigma=0.1)     # Smearing width

VASP-specific:

• Requires VASP license
• Uses POTCAR files for pseudopotentials
• See references/calculator_setup.md for details

Advanced Properties (Subskills)

1. Surface Energy

Method: Slab model approach

Formula:

γ = (E_slab - N × E_bulk) / (2 × A)

Where:

• E_slab: Energy of slab with surfaces
• E_bulk: Bulk energy per atom
• N: Number of atoms in slab
• A: Surface area
• Factor of 2: Two equivalent surfaces

Workflow:

1. Optimize bulk structure
2. Create slab with ase.build.surface()
3. Add vacuum layer
4. Relax slab (constrain bottom layers)
5. Calculate surface energy

See: references/surface_energy.md for detailed methods

Key References:

• Fiorentini & Methfessel, "Extracting convergent surface energies," J. Phys.: Condens. Matter 8, 6525 (1996)
• Tran et al., "Surface energies of elemental crystals," Sci. Data 3, 160080 (2016)

2. Adsorption Energy

Method: Compare slab+adsorbate to separated systems

Formula:

E_ads = E_slab+ads - E_slab - E_molecule

More negative = stronger adsorption

Workflow:

1. Optimize clean slab
2. Add adsorbate at different sites (top, bridge, hollow)
3. Relax adsorbed structure
4. Calculate adsorption energy
5. Compare sites to find preferred adsorption

See: references/adsorption_energy.md

Key References:

• Hammer & Nørskov, "Theoretical surface science," Adv. Catal. 45, 71 (2000)
• Nørskov et al., "Computational design of solid catalysts," Nature Chem. 1, 37 (2009)

3. Reaction Barriers (Nudged Elastic Band)

Method: Find minimum energy path between reactant and product

Workflow:

from ase.neb import NEB
from ase.optimize import BFGS

# Create images interpolating between initial and final
images = [initial]
images += [initial.copy() for i in range(5)]  # 5 intermediate images
images += [final]

# Interpolate
neb = NEB(images)
neb.interpolate()

# Set calculators
for image in images[1:-1]:
    image.calc = EMT()

# Optimize NEB
optimizer = BFGS(neb, trajectory='neb.traj')
optimizer.run(fmax=0.05)

# Extract barrier
from ase.neb import NEBTools
nebtools = NEBTools(images)
barrier = nebtools.get_barrier()[0]
print(f"Activation barrier: {barrier:.2f} eV")

See: references/reaction_barriers.md

Key References:

• Henkelman, Uberuaga & Jónsson, "Climbing image NEB," J. Chem. Phys. 113, 9901 (2000)
• Sheppard et al., "Optimization methods for MEPs," J. Chem. Phys. 128, 134106 (2008)

4. Vibrational Analysis

Method: Finite displacement or DFPT

Workflow:

from ase.vibrations import Vibrations

# Calculate vibrations
vib = Vibrations(atoms)
vib.run()

# Get frequencies
vib.summary()

# Zero-point energy
zpe = vib.get_zero_point_energy()

See: references/vibrational_analysis.md

Key References:

• Wilson, Decius & Cross, Molecular Vibrations (Dover, 1980)

5. Phonon Calculations (phonopy)

Method: Supercell approach with force constants

Workflow:

from phonopy import Phonopy

# Create phonopy object
phonon = Phonopy(atoms, supercell_matrix=[[2,0,0],[0,2,0],[0,0,2]])

# Generate displacements
phonon.generate_displacements(distance=0.01)
supercells = phonon.supercells_with_displacements

# Calculate forces for each displacement
for scell in supercells:
    scell.calc = calc
    forces = scell.get_forces()
    # Set forces back to phonopy

# Calculate phonon properties
phonon.produce_force_constants()
phonon.auto_band_structure()
phonon.plot_band_structure()

Applications:

• Phonon band structure
• Phonon density of states
• Thermal properties (heat capacity, free energy)
• Thermodynamic integration

See: references/phonons.md

Key References:

• Togo & Tanaka, "First principles phonon calculations," Scr. Mater. 108, 1 (2015)
• Togo, "Phonopy and Phono3py," J. Phys. Soc. Jpn. 92, 012001 (2023)

6. Elastic Constants (elastic package)

Method: Apply strain, measure stress

Properties Calculated:

• Full elastic tensor (Cij)
• Bulk modulus (K)
• Shear modulus (G)
• Young's modulus (E)
• Poisson's ratio (ν)
• Sound velocities
• Debye temperature

See: references/elastic_constants.md

Key References:

• Golesorkhtabar et al., "ElaStic tool," Comput. Phys. Commun. 184, 1861 (2013)
• Nye, Physical Properties of Crystals (Oxford, 1985)

7. Equation of State

Method: Volume-energy curve fitting

Workflow:

from ase.eos import calculate_eos

eos = calculate_eos(atoms, trajectory='eos.traj')
v, e, B = eos.fit()  # Volume, energy, bulk modulus
eos.plot('eos.png')

EOS Types:

• Birch-Murnaghan
• Murnaghan
• Vinet

See: references/equation_of_state.md

Key References:

• Birch, "Finite elastic strain," Phys. Rev. 71, 809 (1947)

8. Formation Energy

Formula:

E_form = E_compound - Σ(n_i × μ_i)

Where μ_i are chemical potentials (reference energies)

Applications:

• Phase stability
• Energy above convex hull
• Phase diagrams

See: references/formation_energy.md

Key References:

• Hautier et al., "DFT formation energies," Phys. Rev. B 85, 155208 (2012)

9. Cluster Expansions (icet)

Method: Expand configurational energy in cluster interactions

Applications:

• Alloy ground states
• Order-disorder transitions
• Monte Carlo simulations
• Phase diagram construction

See: references/cluster_expansion.md

Key References:

• Ångqvist et al., "ICET library," Adv. Theory Simul. 2, 1900015 (2019)
• Sanchez et al., "Cluster description," Physica A 128, 334 (1984)

10. CALPHAD Integration (pycalphad)

Method: Combine DFT with thermochemical databases

Applications:

• Phase equilibria
• Multi-component systems
• Temperature-dependent properties

See: references/calphad.md

Key References:

• Otis & Liu, "pycalphad," J. Open Res. Softw. 5, 1 (2017)
• Lukas et al., Computational Thermodynamics (Cambridge, 2007)

11. Defect Formation Energy

Types:

• Vacancies
• Interstitials
• Substitutional defects
• Charged defects (with corrections)

See: references/defect_energy.md

Key References:

• Freysoldt et al., "Point defects in solids," Rev. Mod. Phys. 86, 253 (2014)

12. Interface/Grain Boundary Energy

Method: Compare interface structure to separated surfaces

See: references/interface_energy.md

Key References:

• Sutton & Balluffi, Interfaces in Crystalline Materials (Oxford, 1995)

13. Magnetic Properties

Method: Spin-polarized DFT

Properties:

• Magnetic moments
• Magnetic ordering (FM, AFM)
• Heisenberg parameters

See: references/magnetic_properties.md

14. Thermal Expansion

Method: Quasi-harmonic approximation

See: references/thermal_expansion.md

Key References:

• Barrera et al., "Grüneisen parameters," J. Phys.: Condens. Matter 17, R217 (2005)

15. Electronic Structure

Properties:

• Band structure
• Density of states (DOS)
• Band gaps
• Fermi surfaces

See: references/electronic_structure.md

Key References:

• Martin, Electronic Structure (Cambridge, 2004)

Best Practices

Convergence Testing

Always test convergence of:

1. k-point sampling: Increase until energy converges (typically < 1 meV/atom)
2. Plane-wave cutoff: Test different values (e.g., 300-600 eV)
3. Slab thickness: For surfaces (typically 5-9 layers)
4. Vacuum thickness: For surfaces (typically 10-15 Å)
5. Supercell size: For defects, phonons

Example:

# k-point convergence
for k in [2, 4, 6, 8, 10, 12]:
    atoms.calc = GPAW(kpts=(k,k,k), ...)
    E = atoms.get_potential_energy()
    print(f"k={k}: E={E:.4f} eV")

Force Convergence

• Typical: fmax = 0.01-0.05 eV/Å
• Tighter for vibrations: fmax = 0.001 eV/Å
• Check max force: max(np.linalg.norm(atoms.get_forces(), axis=1))

Constraints

Fix atoms during relaxation:

from ase.constraints import FixAtoms

# Fix bottom 2 layers of slab
c = FixAtoms(indices=[atom.index for atom in atoms if atom.position[2] < 5])
atoms.set_constraint(c)

Trajectory Analysis

from ase.io import read

# Read optimization trajectory
traj = read('opt.traj', ':')

# Plot energy vs step
energies = [atoms.get_potential_energy() for atoms in traj]
import matplotlib.pyplot as plt
plt.plot(energies)
plt.xlabel('Step')
plt.ylabel('Energy (eV)')
plt.show()

Common Workflows

Workflow 1: Lattice Constant Determination

from ase.build import bulk
from ase.eos import calculate_eos

atoms = bulk('Cu', 'fcc', a=3.6)
atoms.calc = EMT()

eos = calculate_eos(atoms, trajectory='eos.traj')
v, e, B = eos.fit()
a_opt = v**(1/3)
print(f"Optimal lattice constant: {a_opt:.3f} Å")
print(f"Bulk modulus: {B/1e9:.1f} GPa")

Workflow 2: Surface Energy Calculation

from ase.build import bulk, surface, add_vacuum

# Bulk energy
bulk_atoms = bulk('Cu', 'fcc', a=3.6)
bulk_atoms.calc = EMT()
E_bulk_per_atom = bulk_atoms.get_potential_energy() / len(bulk_atoms)

# Create slab
slab = surface('Cu', (1,1,1), layers=7, vacuum=10)
slab.calc = EMT()
E_slab = slab.get_potential_energy()

# Surface energy
N = len(slab)
A = slab.get_cell()[0,0] * slab.get_cell()[1,1]
gamma = (E_slab - N * E_bulk_per_atom) / (2 * A)
print(f"Surface energy: {gamma*1000:.1f} meV/ų")

Workflow 3: Adsorption Energy

from ase.build import fcc111, molecule, add_adsorbate

# Clean slab
slab = fcc111('Cu', size=(3,3,4), vacuum=10)
slab.calc = EMT()
E_slab = slab.get_potential_energy()

# Adsorbate in gas phase
mol = molecule('CO')
mol.calc = EMT()
E_mol = mol.get_potential_energy()

# Adsorbed system
add_adsorbate(slab, mol, height=2.0, position='ontop')
slab.calc = EMT()
E_ads_system = slab.get_potential_energy()

# Adsorption energy
E_ads = E_ads_system - E_slab - E_mol
print(f"Adsorption energy: {E_ads:.2f} eV")

Error Handling and Troubleshooting

Common Issues

1. SCF Not Converging:

• Increase mixing parameter
• Try different smearing methods
• Check initial geometry (remove overlaps)

2. Forces Not Converging:

• Check constraints
• Try different optimizer
• Increase maximum steps
• Verify calculator parameters

3. Unstable Structures:

• Check for imaginary phonon modes
• Verify symmetry is correct
• Try different initial configurations

4. Memory Issues:

• Reduce k-points or cutoff
• Use real-space mode (GPAW)
• Parallelize calculation

Integration with Other Skills

• python-ase: Core ASE functionality
• materials-databases: Get structures from Materials Project/AFLOW
• pymatgen: Structure manipulation and analysis
• fairchem: Use ML potentials for screening

References

Core ASE:

1. Larsen et al., "The atomic simulation environment," J. Phys.: Condens. Matter 29, 273002 (2017)

DFT Theory: 2. Hohenberg & Kohn, "Inhomogeneous electron gas," Phys. Rev. 136, B864 (1964) 3. Kohn & Sham, "Self-consistent equations," Phys. Rev. 140, A1133 (1965) 4. Sholl & Steckel, Density Functional Theory: A Practical Introduction (Wiley, 2009)

Symmetry Analysis: 5. Spglib: https://spglib.github.io/spglib/

See individual reference files in references/ for detailed citations for each method.

Resources

• ASE Documentation: https://wiki.fysik.dtu.dk/ase/
• ASE Tutorials: https://wiki.fysik.dtu.dk/ase/tutorials/tutorials.html
• Example Scripts: See examples/ directory
• Method Details: See references/ directory
• Best Practices: See workflows/ directory