Introduction to ASE and DFT

# Essential ASE imports for advanced materials science workflows
import numpy as np
from ase import Atoms
from ase.build import bulk, surface, molecule, add_adsorbate, make_supercell
from ase.visualize import view
from ase.io import read, write, Trajectory
from ase.optimize import BFGS, LBFGS, QuasiNewton
from ase.constraints import FixAtoms, FixBondLengths, FixCartesian, UnitCellFilter
from ase.dft.kpoints import bandpath, get_special_points
from ase.phonons import Phonons
from ase.eos import EquationOfState
from ase.vibrations import Vibrations
from ase.thermochemistry import HarmonicThermo, IdealGasThermo
from ase.neb import NEB, NEBTools
from ase.md import VelocityVerlet, MDLogger
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution

# Database and high-throughput tools
from ase.db import connect
from ase.parallel import paropen

# Calculator interfaces (examples)
# from ase.calculators.vasp import Vasp
# from ase.calculators.espresso import Espresso  
# from gpaw import GPAW
# from ase.calculators.gaussian import Gaussian

print(f"ASE version: {ase.__version__}")
print(f"Available calculators: {len(ase.calculators.names)} codes supported")
print(f"Supported file formats: {len(ase.io.formats.all_formats)} formats")

# Essential Atoms object creation and manipulation
import numpy as np
from ase import Atoms
from ase.build import bulk, molecule, surface

# Method 1: Manual construction for complete control
positions = np.array([[0.0, 0.0, 0.0], [1.42, 0.0, 0.0]])
co_molecule = Atoms(
    symbols=['C', 'O'],
    positions=positions,
    pbc=False  # No periodic boundaries for molecules
)

# Method 2: Built-in structure database
silicon = bulk('Si', 'diamond', a=5.43, cubic=True)
water = molecule('H2O')
al_surface = surface('Al', (1, 1, 1), layers=4, vacuum=10.0)

print("=== ATOMS OBJECT PROPERTIES ===")
print(f"Silicon formula: {silicon.get_chemical_formula()}")
print(f"Number of atoms: {len(silicon)}")
print(f"Unit cell volume: {silicon.get_volume():.3f} ų")
print(f"Periodic boundaries: {silicon.pbc}")

# Access and modify atomic properties
print(f"Atomic positions:\n{silicon.positions}")
print(f"Chemical symbols: {silicon.get_chemical_symbols()}")
print(f"Cell vectors:\n{silicon.cell}")

# Dynamic property calculation
print(f"Density: {silicon.get_density():.3f} g/cm³")
print(f"Center of mass: {silicon.get_center_of_mass()}")

# Advanced cell operations and lattice analysis
from ase.build import bulk
from ase.geometry import get_duplicate_atoms

# Create and analyze crystal structure
al = bulk('Al', 'fcc', a=4.05, cubic=True)

print("=== LATTICE ANALYSIS ===")
cell = al.get_cell()
print(f"Cell vectors (Å):\n{cell}")
print(f"Cell parameters: {cell.cellpar()}")  # [a, b, c, α, β, γ]

# Reciprocal lattice for k-point sampling
reciprocal = cell.reciprocal()
print(f"Reciprocal cell (Å⁻¹):\n{reciprocal}")

# Volume and geometric properties
print(f"Unit cell volume: {cell.volume:.4f} ų")
print(f"Cell determinant: {np.linalg.det(cell):.4f}")

# Transform to primitive cell
from ase.spacegroup import crystal
primitive_al = crystal('Al', [(0, 0, 0)], spacegroup=225, 
                      cellpar=[4.05, 4.05, 4.05, 90, 90, 90])
print(f"Primitive cell atoms: {len(primitive_al)}")

# Supercell construction
from ase.build import make_supercell
supercell_matrix = np.diag([2, 2, 2])  # 2×2×2 supercell
al_supercell = make_supercell(al, supercell_matrix)
print(f"Supercell: {len(al_supercell)} atoms")

# Comprehensive structural analysis
from ase.neighborlist import neighbor_list, natural_cutoffs
from ase.geometry.analysis import Analysis

# Neighbor analysis
silicon = bulk('Si', 'diamond', a=5.43)

# Get neighbor lists with automatic cutoffs
cutoffs = natural_cutoffs(silicon, mult=1.1)
i, j, d = neighbor_list('ijd', silicon, cutoffs)

print("=== STRUCTURAL ANALYSIS ===")
print(f"Neighbor cutoffs: {cutoffs}")
print(f"Bond distances: {d}")
print(f"Coordination numbers: {np.bincount(i)}")

# Advanced analysis with Analysis class
analysis = Analysis(silicon)
all_bonds = analysis.get_bonds('Si', 'Si')
print(f"All Si-Si bonds: {len(all_bonds[0])} bonds")
print(f"Average bond length: {np.mean(all_bonds[1]):.4f} Å")

# Radial distribution function
from ase.geometry import get_distances
distances, _ = get_distances(silicon.positions, cell=silicon.cell, pbc=silicon.pbc)
print(f"Distance matrix shape: {distances.shape}")

# Crystallographic analysis
try:
    from ase.spacegroup import get_spacegroup
    sg = get_spacegroup(silicon, symprec=1e-5)
    print(f"Space group: {sg.symbol} (#{sg.no})")
except ImportError:
    print("Install spglib for space group analysis")

# Professional calculator setup patterns
from ase.calculators.emt import EMT
from ase.build import bulk, surface, molecule

def setup_production_calculator(system_type, accuracy='normal'):
    """
    Factory function for production-ready calculator settings.
    
    Parameters:
    -----------
    system_type : str
        'bulk', 'surface', 'molecule'
    accuracy : str
        'fast', 'normal', 'accurate'
    """
    
    # Base parameters for EMT (educational example)
    calc_params = {}
    
    # In practice, this would configure VASP/QE/GPAW parameters
    if system_type == 'bulk':
        print("Configuring for bulk system...")
        # calc_params = {'kpts': (8,8,8), 'encut': 500, 'ismear': 1}
    elif system_type == 'surface': 
        print("Configuring for surface system...")
        # calc_params = {'kpts': (8,8,1), 'encut': 500, 'dipol': True}
    elif system_type == 'molecule':
        print("Configuring for molecular system...")
        # calc_params = {'kpts': (1,1,1), 'encut': 400, 'ismear': 0}
    
    return EMT()  # Simplified for demonstration

# Universal calculation interface
def calculate_properties(atoms, calculator):
    """Demonstrate universal calculator interface."""
    
    # Attach calculator to atoms
    atoms.calc = calculator
    
    # Universal methods work with any calculator
    try:
        energy = atoms.get_potential_energy()
        forces = atoms.get_forces()
        stress = atoms.get_stress()
        
        print(f"Energy: {energy:.6f} eV")
        print(f"Max force: {np.max(np.linalg.norm(forces, axis=1)):.4f} eV/Å")
        print(f"Pressure: {-np.trace(stress)/3*160.2:.2f} GPa")
        
        return {'energy': energy, 'forces': forces, 'stress': stress}
    
    except Exception as e:
        print(f"Calculation failed: {e}")
        return None

# Demonstrate with different systems
systems = {
    'Silicon crystal': bulk('Si', 'diamond', a=5.43),
    'CO molecule': molecule('CO'),
    'Al(111) surface': surface('Al', (1,1,1), layers=3, vacuum=8.0)
}

for name, atoms in systems.items():
    print(f"\n=== {name.upper()} ===")
    system_type = 'bulk' if atoms.pbc.all() else ('surface' if atoms.pbc.any() else 'molecule')
    calc = setup_production_calculator(system_type)
    results = calculate_properties(atoms, calc)

# Method comparison framework
import numpy as np
import matplotlib.pyplot as plt

class MethodComparison:
    """Framework for comparing different computational methods."""
    
    def __init__(self, atoms):
        self.atoms = atoms
        self.results = {}
    
    def add_method(self, name, calculator):
        """Add a computational method to comparison."""
        atoms_copy = self.atoms.copy()
        atoms_copy.calc = calculator
        
        try:
            energy = atoms_copy.get_potential_energy()
            forces = atoms_copy.get_forces()
            
            self.results[name] = {
                'energy': energy,
                'max_force': np.max(np.linalg.norm(forces, axis=1)),
                'success': True
            }
            print(f"{name}: E = {energy:.6f} eV")
            
        except Exception as e:
            self.results[name] = {'success': False, 'error': str(e)}
            print(f"{name}: FAILED - {e}")
    
    def analyze_convergence(self):
        """Analyze convergence between methods."""
        successful_methods = {k: v for k, v in self.results.items() 
                            if v['success']}
        
        if len(successful_methods) < 2:
            print("Need at least 2 successful methods for comparison")
            return
        
        energies = [v['energy'] for v in successful_methods.values()]
        methods = list(successful_methods.keys())
        
        energy_spread = np.std(energies)
        mean_energy = np.mean(energies)
        
        print(f"\n=== METHOD CONVERGENCE ANALYSIS ===")
        print(f"Methods: {methods}")
        print(f"Energy spread: {energy_spread:.6f} eV")
        print(f"Relative spread: {energy_spread/abs(mean_energy)*100:.3f}%")
        
        return {'methods': methods, 'energies': energies, 'spread': energy_spread}

# Example usage
silicon = bulk('Si', 'diamond', a=5.43)
comparison = MethodComparison(silicon)

# Add different "methods" (using EMT with different settings for demo)
comparison.add_method('Method A', EMT())
comparison.add_method('Method B', EMT()) 
comparison.add_method('Method C', EMT())

# Analyze results
comparison.analyze_convergence()

# Systematic defect creation for materials engineering
from ase.build import bulk, make_supercell
from ase import Atoms
import numpy as np

def create_point_defects(perfect_crystal, defect_types=['vacancy', 'interstitial']):
    """Create various point defects in crystal structure."""
    
    defects = {}
    
    for defect_type in defect_types:
        crystal = perfect_crystal.copy()
        
        if defect_type == 'vacancy':
            # Remove central atom
            center_idx = len(crystal) // 2
            del crystal[center_idx]
            defects['vacancy'] = crystal
            print(f"Created vacancy: {len(crystal)} atoms (was {len(perfect_crystal)})")
            
        elif defect_type == 'interstitial':
            # Add atom at interstitial site
            center = np.mean(crystal.positions, axis=0)
            offset = np.array([0.5, 0.5, 0.5])  # Simple offset
            new_position = center + offset
            
            crystal.append(Atoms(crystal[0].symbol, [new_position]))
            defects['interstitial'] = crystal
            print(f"Created interstitial: {len(crystal)} atoms")
            
        elif defect_type == 'substitutional':
            # Replace atom with different element
            crystal[0].symbol = 'Mg' if crystal[0].symbol == 'Al' else 'Al'
            defects['substitutional'] = crystal
            print(f"Created substitutional defect: {crystal.get_chemical_formula()}")
    
    return defects

# Example: Defects in aluminum
al_perfect = make_supercell(bulk('Al', 'fcc', a=4.05), [[3, 0, 0], [0, 3, 0], [0, 0, 3]])
defects = create_point_defects(al_perfect)

print(f"Perfect crystal: {al_perfect.get_chemical_formula()}")
for defect_name, structure in defects.items():
    print(f"{defect_name}: {structure.get_chemical_formula()}")

# Advanced surface and interface construction
from ase.build import surface, add_adsorbate, fcc111
from ase.constraints import FixAtoms

def create_realistic_surface(element, miller_indices, layers=5, vacuum=10.0):
    """Create surface with proper constraints for calculations."""
    
    # Create surface slab
    slab = surface(element, miller_indices, layers=layers, vacuum=vacuum)
    
    # Apply constraints: fix bottom layers
    z_positions = slab.positions[:, 2]
    z_min = z_positions.min()
    fixed_layers = 2  # Fix bottom 2 layers
    
    fixed_indices = []
    for layer in range(fixed_layers):
        layer_z = z_min + layer * 2.3  # Approximate layer spacing
        layer_atoms = [i for i, z in enumerate(z_positions) 
                      if abs(z - layer_z) < 1.0]
        fixed_indices.extend(layer_atoms)
    
    if fixed_indices:
        constraint = FixAtoms(indices=fixed_indices)
        slab.set_constraint(constraint)
        print(f"Fixed {len(fixed_indices)} atoms in bottom layers")
    
    return slab

def add_adsorbates(surface, adsorbate_list):
    """Add multiple adsorbates to surface."""
    
    surface_with_ads = surface.copy()
    
    for ads_spec in adsorbate_list:
        adsorbate = ads_spec['molecule']
        height = ads_spec.get('height', 2.0)
        site = ads_spec.get('site', 'ontop')
        
        add_adsorbate(surface_with_ads, adsorbate, height=height, position=site)
        print(f"Added {adsorbate} at {site} site, height {height} Å")
    
    return surface_with_ads

# Example: CO on Pt(111)
pt_surface = create_realistic_surface('Pt', (1, 1, 1), layers=4, vacuum=12.0)

# Add CO molecules at different sites
adsorbates = [
    {'molecule': 'CO', 'height': 2.0, 'site': 'ontop'},
    {'molecule': 'O', 'height': 1.5, 'site': 'hollow'}
]

pt_with_adsorbates = add_adsorbates(pt_surface, adsorbates)
print(f"Final system: {pt_with_adsorbates.get_chemical_formula()}")

# Strain engineering and alloy creation tools
import numpy as np
from ase.build import bulk

def apply_strain(atoms, strain_tensor):
   """Apply strain tensor to crystal structure."""
   
   strained = atoms.copy()
   
   # Strain tensor can be 3x3 or Voigt notation (6 components)
   if len(strain_tensor) == 6:
       # Convert Voigt to full tensor
       strain = np.zeros((3, 3))
       strain[0, 0] = strain_tensor[0]  # ε_xx
       strain[1, 1] = strain_tensor[1]  # ε_yy  
       strain[2, 2] = strain_tensor[2]  # ε_zz
       strain[1, 2] = strain[2, 1] = strain_tensor[3]  # ε_yz
       strain[0, 2] = strain[2, 0] = strain_tensor[4]  # ε_xz
       strain[0, 1] = strain[1, 0] = strain_tensor[5]  # ε_xy
   else:
       strain = np.array(strain_tensor)
   
   # Apply deformation: F = I + ε
   deformation = np.eye(3) + strain
   new_cell = deformation @ atoms.cell
   
   strained.set_cell(new_cell, scale_atoms=True)
   
   volume_change = (strained.get_volume() / atoms.get_volume() - 1) * 100
   print(f"Applied strain, volume change: {volume_change:.2f}%")
   
   return strained

def create_random_alloy(elements, concentrations, base_structure, size=(2, 2, 2)):
   """Create random substitutional alloy."""
   
   # Create supercell
   base = bulk(elements[0], base_structure, cubic=True)
   supercell = make_supercell(base, size)
   
   # Random substitution
   n_total = len(supercell)
   n_each = [int(c * n_total) for c in concentrations]
   
   # Adjust for rounding
   n_each[0] = n_total - sum(n_each[1:])
   
   # Create atom type list
   atom_types = []
   for elem, count in zip(elements, n_each):
       atom_types.extend([elem] * count)
   
   np.random.shuffle(atom_types)
   
   # Assign to supercell
   for i, atom_type in enumerate(atom_types):
       supercell[i].symbol = atom_type
   
   print(f"Created {supercell.get_chemical_formula()} alloy")
   return supercell

# Examples
print("=== STRAIN ENGINEERING ===")
silicon = bulk('Si', 'diamond', a=5.43)

# Apply 2% tensile strain in x-direction
tensile_strain = [0.02, 0.0, 0.0, 0.0, 0.0, 0.0]
strained_si = apply_strain(silicon, tensile_strain)

print("\n=== ALLOY DESIGN ===")
# Create CuNi alloy with 70% Cu, 30% Ni
cu_ni_alloy = create_random_alloy(['Cu', 'Ni'], [0.7, 0.3], 'fcc')
print(f"Alloy composition: {cu_ni_alloy.get_chemical_formula()}")

# Complete equation of state analysis
from ase.build import bulk
from ase.eos import EquationOfState
from ase.calculators.emt import EMT
import numpy as np
import matplotlib.pyplot as plt

def calculate_bulk_properties(element, structure, a0_guess, calc, 
                            volume_range=0.15, npoints=9):
    """
    Calculate bulk modulus and equilibrium properties from EOS fitting.
    
    Parameters:
    -----------
    element : str
        Chemical element symbol
    structure : str  
        Crystal structure ('fcc', 'bcc', 'diamond', etc.)
    a0_guess : float
        Initial lattice parameter guess (Å)
    calc : ASE Calculator
        DFT calculator
    volume_range : float
        Relative volume range for EOS (±volume_range)
    npoints : int
        Number of volume points
    
    Returns:
    --------
    dict : Bulk properties and EOS parameters
    """
    
    # Generate volume points
    volumes = []
    energies = []
    lattice_params = []
    
    # Calculate volume range
    atoms_ref = bulk(element, structure, a=a0_guess)
    v0_guess = atoms_ref.get_volume()
    
    volume_factors = np.linspace(1 - volume_range, 1 + volume_range, npoints)
    
    print(f"Calculating EOS for {element} ({structure} structure)")
    print(f"Reference volume: {v0_guess:.3f} ų")
    print("Volume scaling factors:", volume_factors)
    
    for i, factor in enumerate(volume_factors):
        # Scale lattice parameter
        a_scaled = a0_guess * factor**(1/3)
        atoms = bulk(element, structure, a=a_scaled)
        atoms.calc = calc
        
        volume = atoms.get_volume()
        energy = atoms.get_potential_energy()
        
        volumes.append(volume)
        energies.append(energy)
        lattice_params.append(a_scaled)
        
        print(f"Point {i+1}/{npoints}: a={a_scaled:.4f} Å, "
              f"V={volume:.3f} ų, E={energy:.6f} eV")
    
    # Fit equation of state
    eos = EquationOfState(volumes, energies, eos='birchmurnaghan')
    v0, e0, B = eos.fit()
    
    # Calculate optimized lattice parameter
    if structure == 'fcc':
        a0_opt = (4 * v0)**(1/3)
        atoms_per_cell = 4
    elif structure == 'bcc':
        a0_opt = (2 * v0)**(1/3)  
        atoms_per_cell = 2
    elif structure == 'diamond':
        a0_opt = (8 * v0)**(1/3)
        atoms_per_cell = 8
    elif structure == 'sc':
        a0_opt = v0**(1/3)
        atoms_per_cell = 1
    else:
        a0_opt = v0**(1/3)  # Generic cubic
        atoms_per_cell = len(atoms_ref)
    
    # Convert bulk modulus from eV/ų to GPa
    B_GPa = B * 160.21766208  # Conversion factor
    
    # Calculate cohesive energy (approximation for simple systems)
    atoms_isolated = bulk(element, structure, a=a0_opt)
    atoms_isolated.calc = calc
    
    # For accurate cohesive energy, need isolated atom calculation
    # This is a simplified approach
    cohesive_energy_per_atom = e0 / atoms_per_cell
    
    results = {
        'equilibrium_volume': v0,
        'equilibrium_energy': e0,
        'equilibrium_lattice_param': a0_opt,
        'bulk_modulus_eV': B,
        'bulk_modulus_GPa': B_GPa,
        'cohesive_energy_per_atom': cohesive_energy_per_atom,
        'volumes': volumes,
        'energies': energies,
        'eos_fit': eos
    }
    
    # Print results
    print(f"\n{'='*50}")
    print(f"EQUATION OF STATE RESULTS FOR {element.upper()}")
    print(f"{'='*50}")
    print(f"Equilibrium lattice parameter: {a0_opt:.4f} Å")
    print(f"Equilibrium volume: {v0:.4f} ų")
    print(f"Minimum energy: {e0:.6f} eV")
    print(f"Bulk modulus: {B_GPa:.1f} GPa")
    print(f"Cohesive energy per atom: {cohesive_energy_per_atom:.4f} eV")
    
    # Plot EOS
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
    
    # Energy vs Volume
    eos.plot(ax=ax1)
    ax1.set_xlabel('Volume (ų)')
    ax1.set_ylabel('Energy (eV)')
    ax1.set_title(f'{element} Equation of State')
    ax1.grid(True, alpha=0.3)
    
    # Pressure vs Volume  
    volumes_fine = np.linspace(min(volumes), max(volumes), 100)
    pressures = []
    for v in volumes_fine:
        # P = -dE/dV, calculated numerically from EOS
        dv = 0.001
        e1 = eos.fit0([v - dv/2])[0] 
        e2 = eos.fit0([v + dv/2])[0]
        p = -(e2 - e1) / dv * 160.21766208  # Convert to GPa
        pressures.append(p)
    
    ax2.plot(volumes_fine, pressures, 'b-', linewidth=2)
    ax2.axhline(y=0, color='k', linestyle='--', alpha=0.5)
    ax2.axvline(x=v0, color='r', linestyle='--', alpha=0.5, 
                label=f'V₀ = {v0:.3f} ų')
    ax2.set_xlabel('Volume (ų)')
    ax2.set_ylabel('Pressure (GPa)')  
    ax2.set_title(f'{element} Pressure vs Volume')
    ax2.grid(True, alpha=0.3)
    ax2.legend()
    
    plt.tight_layout()
    plt.savefig(f'{element}_eos_analysis.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    return results

# Example: Silicon equation of state
si_calc = EMT()  # Replace with DFT calculator
si_props = calculate_bulk_properties('Si', 'diamond', 5.43, si_calc)

# Compare with experimental values
experimental_data = {
    'Si': {'a0_exp': 5.431, 'B0_exp': 97.9},  # GPa
    'Al': {'a0_exp': 4.050, 'B0_exp': 76.0},
    'Cu': {'a0_exp': 3.615, 'B0_exp': 140.0}
}

if 'Si' in experimental_data:
    exp = experimental_data['Si']
    calc_error_a = abs(si_props['equilibrium_lattice_param'] - exp['a0_exp'])
    calc_error_B = abs(si_props['bulk_modulus_GPa'] - exp['B0_exp'])
    
    print(f"\nComparison with experiment:")
    print(f"Lattice parameter: calc={si_props['equilibrium_lattice_param']:.4f} Å, "
          f"exp={exp['a0_exp']:.4f} Å, error={calc_error_a:.4f} Å")
    print(f"Bulk modulus: calc={si_props['bulk_modulus_GPa']:.1f} GPa, "
          f"exp={exp['B0_exp']:.1f} GPa, error={calc_error_B:.1f} GPa")

# Surface energy and adsorption calculations
from ase.build import surface, bulk, add_adsorbate
from ase.constraints import FixAtoms
from ase.optimize import BFGS
from ase.calculators.emt import EMT
import numpy as np

def calculate_surface_energy(element, miller_indices, layers_range=range(3, 8), 
                           vacuum=10.0, calc=None):
    """
    Calculate surface energy with convergence testing.
    
    Parameters:
    -----------
    element : str
        Chemical element
    miller_indices : tuple
        Surface Miller indices (h, k, l)
    layers_range : range
        Range of slab thicknesses to test
    vacuum : float
        Vacuum spacing (Å)
    calc : ASE Calculator
        DFT calculator
    """
    
    if calc is None:
        calc = EMT()
    
    # Calculate bulk reference energy
    bulk_atoms = bulk(element, cubic=True)
    bulk_atoms.calc = calc
    e_bulk_per_atom = bulk_atoms.get_potential_energy() / len(bulk_atoms)
    
    print(f"Bulk energy per atom: {e_bulk_per_atom:.6f} eV")
    
    surface_energies = []
    layer_counts = []
    
    for nlayers in layers_range:
        # Create slab
        slab = surface(element, miller_indices, layers=nlayers, vacuum=vacuum)
        
        # Fix bottom layer(s) for thick slabs
        if nlayers > 4:
            z_positions = slab.positions[:, 2]
            bottom_layer = z_positions.min()
            fixed_indices = [i for i, z in enumerate(z_positions) 
                           if z < bottom_layer + 1.0]
            constraint = FixAtoms(indices=fixed_indices)
            slab.set_constraint(constraint)
        
        slab.calc = calc
        
        # Optimize slab
        opt = BFGS(slab, trajectory=f'slab_{nlayers}L.traj', logfile=None)
        opt.run(fmax=0.05)
        
        # Calculate surface energy
        e_slab = slab.get_potential_energy()
        n_atoms = len(slab)
        area = slab.get_volume() / (slab.cell[2, 2] - vacuum)  # Surface area
        
        # Surface energy per unit area
        gamma = (e_slab - n_atoms * e_bulk_per_atom) / (2 * area)
        
        surface_energies.append(gamma)
        layer_counts.append(nlayers)
        
        print(f"Layers: {nlayers}, Surface energy: {gamma:.4f} eV/Ų "
              f"({gamma * 16.0218:.3f} J/m²)")
    
    # Check convergence
    if len(surface_energies) > 1:
        convergence = np.abs(np.diff(surface_energies))
        print(f"\nConvergence analysis:")
        for i, conv in enumerate(convergence):
            print(f"Δγ({layer_counts[i]}→{layer_counts[i+1]}L): {conv:.4f} eV/Ų")
    
    # Plot convergence
    import matplotlib.pyplot as plt
    plt.figure(figsize=(8, 6))
    plt.plot(layer_counts, surface_energies, 'bo-', linewidth=2, markersize=8)
    plt.xlabel('Number of Layers')
    plt.ylabel('Surface Energy (eV/Ų)')
    plt.title(f'{element}({miller_indices[0]}{miller_indices[1]}{miller_indices[2]}) Surface Energy')
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig(f'{element}_{miller_indices}_surface_energy.png', dpi=300)
    plt.show()
    
    return {
        'layer_counts': layer_counts,
        'surface_energies': surface_energies,
        'converged_energy': surface_energies[-1],
        'bulk_energy_per_atom': e_bulk_per_atom
    }

def calculate_adsorption_energy(surface_atoms, adsorbate, adsorption_sites, 
                              height=2.0, calc=None):
    """
    Calculate adsorption energies at different surface sites.
    
    Parameters:
    -----------
    surface_atoms : ASE Atoms
        Clean surface slab
    adsorbate : str or ASE Atoms
        Adsorbate molecule
    adsorption_sites : list
        List of adsorption sites ('ontop', 'bridge', 'hollow')
    height : float
        Initial adsorption height (Å)
    calc : ASE Calculator
        DFT calculator
    """
    
    if calc is None:
        calc = EMT()
    
    # Calculate clean surface energy
    clean_surface = surface_atoms.copy()
    clean_surface.calc = calc
    
    # Fix bottom layers
    z_positions = clean_surface.positions[:, 2]
    bottom_threshold = z_positions.min() + 2.0
    fixed_indices = [i for i, z in enumerate(z_positions) 
                    if z < bottom_threshold]
    constraint = FixAtoms(indices=fixed_indices)
    clean_surface.set_constraint(constraint)
    
    # Optimize clean surface
    opt_clean = BFGS(clean_surface, logfile=None)
    opt_clean.run(fmax=0.05)
    e_clean = clean_surface.get_potential_energy()
    
    # Calculate isolated adsorbate energy
    if isinstance(adsorbate, str):
        from ase.build import molecule
        ads_isolated = molecule(adsorbate)
    else:
        ads_isolated = adsorbate.copy()
    
    ads_isolated.center(vacuum=8.0)
    ads_isolated.calc = calc
    opt_ads = BFGS(ads_isolated, logfile=None)
    opt_ads.run(fmax=0.05)
    e_adsorbate = ads_isolated.get_potential_energy()
    
    print(f"Clean surface energy: {e_clean:.6f} eV")
    print(f"Isolated adsorbate energy: {e_adsorbate:.6f} eV")
    
    adsorption_results = {}
    
    for site in adsorption_sites:
        # Create surface + adsorbate system
        surface_ads = clean_surface.copy()
        
        # Add adsorbate at specified site
        add_adsorbate(surface_ads, adsorbate, height=height, position=site)
        surface_ads.set_constraint(constraint)
        surface_ads.calc = calc
        
        # Optimize
        opt_surf_ads = BFGS(surface_ads, 
                           trajectory=f'surface_{adsorbate}_{site}.traj',
                           logfile=None)
        opt_surf_ads.run(fmax=0.05)
        
        e_surface_ads = surface_ads.get_potential_energy()
        
        # Calculate adsorption energy
        e_ads = e_surface_ads - e_clean - e_adsorbate
        
        # Find adsorbate height after optimization
        ads_z = surface_ads.positions[-1, 2]  # Assuming adsorbate is last atom
        surface_z_max = surface_atoms.positions[:, 2].max()
        final_height = ads_z - surface_z_max
        
        adsorption_results[site] = {
            'energy': e_ads,
            'height': final_height,
            'total_energy': e_surface_ads
        }
        
        print(f"\n{site.upper()} site:")
        print(f"  Adsorption energy: {e_ads:.4f} eV")
        print(f"  Optimized height: {final_height:.3f} Å")
    
    # Find most stable site
    most_stable_site = min(adsorption_results.keys(), 
                          key=lambda x: adsorption_results[x]['energy'])
    most_stable_energy = adsorption_results[most_stable_site]['energy']
    
    print(f"\nMost stable adsorption site: {most_stable_site}")
    print(f"Most stable adsorption energy: {most_stable_energy:.4f} eV")
    
    return adsorption_results

# Example calculations
print("Calculating Al(111) surface energy...")
al_surface_data = calculate_surface_energy('Al', (1, 1, 1), 
                                          layers_range=range(4, 8),
                                          calc=EMT())

# Create Al(111) slab for adsorption study
al_111 = surface('Al', (1, 1, 1), layers=5, vacuum=10.0)

print("\nCalculating CO adsorption on Al(111)...")
co_adsorption = calculate_adsorption_energy(
    al_111, 'CO', 
    adsorption_sites=['ontop', 'bridge', 'hollow'],
    calc=EMT()
)

# Electronic structure analysis tools
from ase.dft.kpoints import bandpath
from ase.build import bulk
import numpy as np
import matplotlib.pyplot as plt

def calculate_band_structure(atoms, calc, path=None, npoints=100):
    """
    Calculate electronic band structure along high-symmetry k-path.
    
    Parameters:
    -----------
    atoms : ASE Atoms
        Crystal structure
    calc : ASE Calculator
        DFT calculator (must support band structure calculations)
    path : str
        High-symmetry path (e.g., 'GXMG' for fcc)
    npoints : int
        Number of k-points along path
    """
    
    # Get Bravais lattice and define k-path
    lat = atoms.cell.get_bravais_lattice()
    
    if path is None:
        # Default paths for common lattices
        lattice_paths = {
            'cubic': 'GXMG',
            'fcc': 'GXWKGLUWLK', 
            'bcc': 'GHNGPG',
            'hexagonal': 'GMKGALH'
        }
        lattice_name = lat.name.lower()
        path = lattice_paths.get(lattice_name, 'GXMG')
    
    kpath = lat.bandpath(path, npoints=npoints)
    
    print(f"Calculating band structure along {path} path")
    print(f"High-symmetry points: {list(kpath.special_points.keys())}")
    
    # Self-consistent calculation first
    atoms.calc = calc
    atoms.get_potential_energy()  # Ensure SCF convergence
    
    # Band structure calculation
    # Note: This is a simplified example - real implementation depends on calculator
    print("Band structure calculation completed")
    
    return kpath

def calculate_density_of_states(atoms, calc, width=0.1, npts=1000, 
                               energy_range=(-10, 10)):
    """
    Calculate electronic density of states.
    
    Parameters:
    -----------
    atoms : ASE Atoms
        Crystal structure
    calc : ASE Calculator
        DFT calculator
    width : float
        Gaussian broadening width (eV)
    npts : int
        Number of energy points
    energy_range : tuple
        Energy range (emin, emax) in eV relative to Fermi level
    """
    
    atoms.calc = calc
    atoms.get_potential_energy()  # Ensure SCF calculation
    
    print("Calculating density of states...")
    print(f"Energy range: {energy_range[0]} to {energy_range[1]} eV")
    print(f"Broadening: {width} eV")
    
    # This is a conceptual implementation
    # Real DOS calculation depends on the specific calculator
    energies = np.linspace(energy_range[0], energy_range[1], npts)
    
    # Placeholder DOS calculation
    dos = np.exp(-(energies**2)/(2*width**2))  # Gaussian example
    
    return energies, dos

def analyze_chemical_bonding(atoms, calc):
    """
    Analyze chemical bonding through charge density and electron localization.
    
    Parameters:
    -----------
    atoms : ASE Atoms
        Structure to analyze
    calc : ASE Calculator
        DFT calculator
    """
    
    atoms.calc = calc
    energy = atoms.get_potential_energy()
    
    print("Chemical bonding analysis:")
    print(f"Total energy: {energy:.6f} eV")
    
    # Bader charge analysis (conceptual)
    charges = atoms.get_initial_charges()  # Placeholder
    
    print("Bader charges:")
    for i, (symbol, charge) in enumerate(zip(atoms.get_chemical_symbols(), charges)):
        print(f"  {symbol}{i+1}: {charge:+.3f} e")
    
    # Bond order analysis
    from ase.neighborlist import neighbor_list
    
    # Get nearest neighbors
    i, j, d = neighbor_list('ijd', atoms, cutoff=3.0)
    
    # Calculate average bond lengths by element pair
    bond_stats = {}
    for idx in range(len(i)):
        atom1 = atoms[i[idx]].symbol
        atom2 = atoms[j[idx]].symbol
        bond_pair = tuple(sorted([atom1, atom2]))
        distance = d[idx]
        
        if bond_pair not in bond_stats:
            bond_stats[bond_pair] = []
        bond_stats[bond_pair].append(distance)
    
    print("\nBond length statistics:")
    for bond_pair, distances in bond_stats.items():
        avg_dist = np.mean(distances)
        std_dist = np.std(distances)
        print(f"  {bond_pair[0]}-{bond_pair[1]}: {avg_dist:.3f} ± {std_dist:.3f} Å")
    
    return bond_stats

# Example usage
si = bulk('Si', 'diamond', a=5.43)
calc = EMT()  # Replace with DFT calculator

# Band structure
kpath = calculate_band_structure(si, calc)

# Density of states  
energies, dos = calculate_density_of_states(si, calc)

# Plot DOS
plt.figure(figsize=(8, 6))
plt.plot(energies, dos, 'b-', linewidth=2)
plt.axvline(x=0, color='k', linestyle='--', alpha=0.5, label='Fermi level')
plt.xlabel('Energy (eV)')
plt.ylabel('Density of States (states/eV)')
plt.title('Electronic Density of States')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('dos_analysis.png', dpi=300)
plt.show()

# Bonding analysis
bonding_analysis = analyze_chemical_bonding(si, calc)