#!/usr/bin/env python3 """ Phase diagram generator using Materials Project data. This script generates phase diagrams for chemical systems using data from the Materials Project database via pymatgen's MPRester. Usage: python phase_diagram_generator.py chemical_system [options] Examples: python phase_diagram_generator.py Li-Fe-O python phase_diagram_generator.py Li-Fe-O --output li_fe_o_pd.png python phase_diagram_generator.py Fe-O --show python phase_diagram_generator.py Li-Fe-O --analyze "LiFeO2" """ import argparse import os import sys from pathlib import Path try: from pymatgen.core import Composition from pymatgen.analysis.phase_diagram import PhaseDiagram, PDPlotter except ImportError: print("Error: pymatgen is not installed. Install with: pip install pymatgen") sys.exit(1) try: from mp_api.client import MPRester except ImportError: print("Error: mp-api is not installed. Install with: pip install mp-api") sys.exit(1) def get_api_key() -> str: """Get Materials Project API key from environment.""" api_key = os.environ.get("MP_API_KEY") if not api_key: print("Error: MP_API_KEY environment variable not set.") print("Get your API key from https://next-gen.materialsproject.org/") print("Then set it with: export MP_API_KEY='your_key_here'") sys.exit(1) return api_key def generate_phase_diagram(chemsys: str, args): """ Generate and analyze phase diagram for a chemical system. Args: chemsys: Chemical system (e.g., "Li-Fe-O") args: Command line arguments """ api_key = get_api_key() print(f"\n{'='*60}") print(f"PHASE DIAGRAM: {chemsys}") print(f"{'='*60}\n") # Get entries from Materials Project print("Fetching data from Materials Project...") with MPRester(api_key) as mpr: entries = mpr.get_entries_in_chemsys(chemsys) print(f"✓ Retrieved {len(entries)} entries") if len(entries) == 0: print(f"Error: No entries found for chemical system {chemsys}") sys.exit(1) # Build phase diagram print("Building phase diagram...") pd = PhaseDiagram(entries) # Get stable entries stable_entries = pd.stable_entries print(f"✓ Phase diagram constructed with {len(stable_entries)} stable phases") # Print stable phases print("\n--- STABLE PHASES ---") for entry in stable_entries: formula = entry.composition.reduced_formula energy = entry.energy_per_atom print(f" {formula:<20} E = {energy:.4f} eV/atom") # Analyze specific composition if requested if args.analyze: print(f"\n--- STABILITY ANALYSIS: {args.analyze} ---") try: comp = Composition(args.analyze) # Find closest entry closest_entry = None min_distance = float('inf') for entry in entries: if entry.composition.reduced_formula == comp.reduced_formula: closest_entry = entry break if closest_entry: # Calculate energy above hull e_above_hull = pd.get_e_above_hull(closest_entry) print(f"Energy above hull: {e_above_hull:.4f} eV/atom") if e_above_hull < 0.001: print(f"Status: STABLE (on convex hull)") elif e_above_hull < 0.05: print(f"Status: METASTABLE (nearly stable)") else: print(f"Status: UNSTABLE") # Get decomposition decomp = pd.get_decomposition(comp) print(f"\nDecomposes to:") for entry, fraction in decomp.items(): formula = entry.composition.reduced_formula print(f" {fraction:.3f} × {formula}") # Get reaction energy rxn_energy = pd.get_equilibrium_reaction_energy(closest_entry) print(f"\nDecomposition energy: {rxn_energy:.4f} eV/atom") else: print(f"No entry found for composition {args.analyze}") print("Checking stability of hypothetical composition...") # Analyze hypothetical composition decomp = pd.get_decomposition(comp) print(f"\nWould decompose to:") for entry, fraction in decomp.items(): formula = entry.composition.reduced_formula print(f" {fraction:.3f} × {formula}") except Exception as e: print(f"Error analyzing composition: {e}") # Get chemical potentials if args.chemical_potentials: print("\n--- CHEMICAL POTENTIALS ---") print("(at stability regions)") try: chempots = pd.get_all_chempots() for element, potentials in chempots.items(): print(f"\n{element}:") for potential in potentials[:5]: # Show first 5 print(f" {potential:.4f} eV") except Exception as e: print(f"Could not calculate chemical potentials: {e}") # Plot phase diagram print("\n--- GENERATING PLOT ---") plotter = PDPlotter(pd, show_unstable=args.show_unstable) if args.output: output_path = Path(args.output) plotter.write_image(str(output_path), image_format=output_path.suffix[1:]) print(f"✓ Phase diagram saved to {output_path}") if args.show: print("Opening interactive plot...") plotter.show() print(f"\n{'='*60}\n") def main(): parser = argparse.ArgumentParser( description="Generate phase diagrams using Materials Project data", formatter_class=argparse.RawDescriptionHelpFormatter, epilog=""" Requirements: - Materials Project API key (set MP_API_KEY environment variable) - mp-api package: pip install mp-api Examples: %(prog)s Li-Fe-O %(prog)s Li-Fe-O --output li_fe_o_phase_diagram.png %(prog)s Fe-O --show --analyze "Fe2O3" %(prog)s Li-Fe-O --analyze "LiFeO2" --show-unstable """ ) parser.add_argument( "chemsys", help="Chemical system (e.g., Li-Fe-O, Fe-O)" ) parser.add_argument( "--output", "-o", help="Output file for phase diagram plot (PNG, PDF, SVG)" ) parser.add_argument( "--show", "-s", action="store_true", help="Show interactive plot" ) parser.add_argument( "--analyze", "-a", help="Analyze stability of specific composition (e.g., LiFeO2)" ) parser.add_argument( "--show-unstable", action="store_true", help="Include unstable phases in plot" ) parser.add_argument( "--chemical-potentials", action="store_true", help="Calculate chemical potentials" ) args = parser.parse_args() # Validate chemical system format elements = args.chemsys.split("-") if len(elements) < 2: print("Error: Chemical system must contain at least 2 elements") print("Example: Li-Fe-O") sys.exit(1) # Generate phase diagram generate_phase_diagram(args.chemsys, args) if __name__ == "__main__": main()