From 50028fa80b2bd932930414b76b8ebc6c0d0d2e25 Mon Sep 17 00:00:00 2001 From: abhijeetgangan Date: Fri, 28 Feb 2025 22:25:51 -0800 Subject: [PATCH] Update docs --- a2c_ase/runner.py | 121 +++++++++------- a2c_ase/utils.py | 358 +++++++++++++++++++++------------------------- example/Si64.py | 43 +++--- 3 files changed, 257 insertions(+), 265 deletions(-) diff --git a/a2c_ase/runner.py b/a2c_ase/runner.py index 0df86c9..aef72e6 100644 --- a/a2c_ase/runner.py +++ b/a2c_ase/runner.py @@ -29,31 +29,48 @@ def melt_quench_md( log_interval: int = 100, ): """ - Run a melt-quench molecular dynamics simulation to generate amorphous structures. - - This function performs a three-stage MD simulation commonly used to create amorphous materials: - 1. High-temperature equilibration (melting) - 2. Controlled cooling (quenching) - 3. Low-temperature equilibration (structural relaxation) - - Args: - atoms: ASE Atoms object containing initial atomic structure - calculator: ASE calculator to compute forces and energies - equi_steps: Number of steps for high-temperature equilibration - cool_steps: Number of steps for cooling phase - final_steps: Number of steps for final low-temperature equilibration - T_high: High temperature for melting phase (Kelvin) - T_low: Low temperature for final phase (Kelvin) - time_step: MD time step in femtoseconds - friction: Langevin friction parameter in atomic units - trajectory_file: Path to save MD trajectory (None = no saving) - seed: Random seed for reproducibility - verbose: If True, print progress information - log_interval: Number of steps between logging progress - - Returns: - atoms: ASE Atoms object with final amorphous structure configuration - log_data: Dictionary containing temperature and energy trajectories + Run melt-quench molecular dynamics to generate amorphous structures. + + Performs a three-stage MD simulation: + 1. High-temperature equilibration to melt the structure + 2. Controlled cooling to quench the liquid + 3. Low-temperature equilibration to relax the structure + + Parameters + ---------- + atoms : Atoms + ASE Atoms object with initial atomic structure + calculator : Calculator + ASE calculator for computing forces and energies + equi_steps : int, default=2500 + Number of steps for high-temperature equilibration + cool_steps : int, default=2500 + Number of steps for cooling phase + final_steps : int, default=2500 + Number of steps for final low-temperature equilibration + T_high : float, default=2000.0 + High temperature for melting phase (K) + T_low : float, default=300.0 + Low temperature for final phase (K) + time_step : float, default=2.0 + MD time step (fs) + friction : float, default=0.01 + Langevin friction parameter (atomic units) + trajectory_file : str, optional + Path to save MD trajectory + seed : int, default=42 + Random seed for reproducibility + verbose : bool, default=True + Whether to print progress information + log_interval : int, default=100 + Number of steps between progress logs + + Returns + ------- + atoms : Atoms + ASE Atoms with final amorphous structure + log_data : dict + Dictionary with temperature and energy trajectories """ # Set the random seed for reproducibility np.random.seed(seed) @@ -147,29 +164,37 @@ def relax_unit_cell( trajectory_file: str | None = None, verbose: bool = True, ) -> tuple[Atoms, dict]: - """Relax both atomic positions and cell parameters using the FIRE algorithm. - - This function performs geometry optimization of both atomic positions and unit cell - parameters simultaneously using ASE's implementation of the Fast Inertial Relaxation - Engine (FIRE) algorithm. The cell relaxation is handled by ASE's ExpCellFilter, which - properly handles the unit cell degrees of freedom. - - Args: - atoms: ASE Atoms object containing atomic positions, species and lattice information - calculator: ASE calculator object that can compute energies, forces and stresses - max_iter: Maximum number of FIRE iterations. Defaults to 200. - fmax: Maximum force criterion for convergence in eV/Å. Defaults to 0.01. - trajectory_file: Optional path to save optimization trajectory. Defaults to None. - verbose: If True, prints optimization progress. Defaults to True. - - Returns: - tuple containing: - - Atoms: Relaxed ASE Atoms object with optimized positions and cell - - dict: Dictionary with energy, forces, stress, and volume trajectories - - Note: - The calculator must be able to compute stresses for cell optimization to work. - The Atoms object must have periodic boundary conditions set appropriately. + """Relax atomic positions and cell parameters using FIRE optimization. + + Performs simultaneous optimization of atomic positions and unit cell parameters using + the Fast Inertial Relaxation Engine (FIRE) algorithm. The cell optimization is handled + through ASE's FrechetCellFilter. + + Parameters + ---------- + atoms : Atoms + ASE Atoms object containing atomic positions, species and cell information + calculator : Calculator + ASE calculator for computing energies, forces and stresses + max_iter : int, optional + Maximum FIRE iterations, by default 200 + fmax : float, optional + Force convergence criterion in eV/Å, by default 0.01 + trajectory_file : str, optional + Path to save optimization trajectory, by default None + verbose : bool, optional + Print optimization progress, by default True + + Returns + ------- + tuple[Atoms, dict] + - Relaxed Atoms object with optimized positions and cell + - Dictionary containing energy, forces, stress and volume trajectories + + Notes + ----- + The calculator must support stress tensor calculations. + Periodic boundary conditions must be enabled on the Atoms object. """ # Ensure atoms has a calculator atoms.calc = calculator diff --git a/a2c_ase/utils.py b/a2c_ase/utils.py index 1912fe6..840b2d4 100644 --- a/a2c_ase/utils.py +++ b/a2c_ase/utils.py @@ -14,44 +14,29 @@ def get_diameter(composition: Composition) -> float: - """Determine the characteristic atomic diameter for a chemical composition. + """Calculate characteristic atomic diameter for a chemical composition. - Estimates the minimum distance between atoms based on their elemental properties. - For compositions with multiple elements, calculates the smallest possible interatomic - distance by examining all possible element pairs and summing their ionic radii. - For single-element compositions, applies element-specific logic: using metallic - radius for metals and atomic/ionic radius for non-metals. + Args: + composition (Composition): Pymatgen Composition object specifying the chemical formula. - This diameter value serves as a practical approximation for the minimum separation - between atoms, which can be valuable for setting up simulation parameters or - initializing atomic structures. + Returns: + float: Estimated minimum atomic diameter in Angstroms. - Parameters: - composition (Composition): Pymatgen Composition object specifying the chemical - formula of interest. May contain one or more elements. + For multi-element compositions: + - Uses ionic radii to find minimum possible interatomic distance + - Examines all element pairs and sums their ionic radii + - Returns smallest pair distance as diameter - Returns: - float: Estimated minimum atomic diameter in Angstroms. For multi-component - systems, represents the smallest expected distance between any atom pair. - For pure elements, equals twice the appropriate atomic radius. + For single-element compositions: + - Uses metallic radius for metals + - Uses atomic radius for non-metals (falls back to ionic radius if needed) + - Returns twice the appropriate radius as diameter Examples: - >>> from pymatgen.core.composition import Composition - >>> # For a compound: Fe2O3 >>> comp = Composition("Fe2O3") - >>> min_distance = get_diameter(comp) # Minimum separation in iron oxide - >>> print(f"{min_distance:.2f}") - - >>> # For a pure element: Cu + >>> get_diameter(comp) # Minimum separation in iron oxide >>> comp = Composition("Cu") - >>> min_distance = get_diameter(comp) # Twice the metallic radius - >>> print(f"{min_distance:.2f}") - - Technical details: - - Multi-element compositions use ionic radii for consistent handling - - Pure metals use metallic radius to better reflect metallic bonding - - Pure non-metals use atomic radius when available, ionic radius as fallback - - All radius data comes from pymatgen's element database + >>> get_diameter(comp) # Twice Cu metallic radius """ elements = composition.elements @@ -82,26 +67,22 @@ def min_distance( structure: Structure, distance_tolerance: float = 0.0001, ) -> float: - """Calculate the minimum distance between any pair of atoms in a periodic structure. - - This function computes all pairwise distances between atoms in a periodic system and - returns the smallest non-zero distance. Self-interactions (an atom's distance to - itself) are excluded using the distance_tolerance parameter. - - Args: - structure: A pymatgen Structure object containing atomic positions, species, - and lattice information. - distance_tolerance: Minimum distance threshold used to exclude self-interactions. - Distances smaller than this value are masked out. Defaults to 0.0001 Å. - - Returns: - float: The minimum distance between any two different atoms in the - structure, considering periodic boundary conditions. - - Note: - The function uses periodic boundary conditions by default. This means atoms near - the cell boundaries are properly connected to atoms on the opposite side of the - cell. + """Calculate minimum interatomic distance in a periodic structure. + + Computes the smallest non-zero distance between any pair of atoms, accounting for + periodic boundary conditions. Self-interactions are excluded via a distance tolerance. + + Parameters + ---------- + structure : Structure + Pymatgen Structure object containing atomic positions and cell information + distance_tolerance : float, default=0.0001 + Distances below this value (in Å) are considered self-interactions and ignored + + Returns + ------- + float + Minimum distance between any two different atoms in Å """ # Get the distance matrix from pymatgen distance_matrix = structure.distance_matrix @@ -130,35 +111,39 @@ def random_packed_structure( trajectory_file: Optional[str] = None, verbose: bool = True, ) -> tuple[Atoms, Optional[list[dict]]]: - """Generates a random packed atomic structure and minimizes atomic overlaps. - - This function creates a random atomic structure within a given cell and optionally - minimizes atomic overlaps using a soft-sphere potential and FIRE optimization. - - Args: - composition: A pymatgen Composition object specifying the atomic composition - (e.g. Fe80B20). The numbers indicate actual atom counts. - cell: A 3x3 numpy array defining the triclinic simulation box in Angstroms. - Does not need to be cubic. - seed: Random seed for reproducible structure generation. If None, uses - random initialization. - diameter: The minimum allowed interatomic distance. Atoms closer than this - distance are considered overlapping. Used for soft-sphere potential. - auto_diameter: If True, automatically calculates appropriate diameter based - on atomic/ionic radii from pymatgen. - max_iter: Maximum number of FIRE optimization steps to minimize overlaps. - Stops early if minimum distance criterion is met. - fmax: Maximum force criterion for convergence in eV/Å. - distance_tolerance: Threshold below which atoms are considered at the same - position when computing minimum distances. - trajectory_file: Path to save the optimization trajectory. If None, no trajectory - is saved. Supports .traj format. - verbose: If True, prints progress information during optimization. - - Returns: - tuple: A tuple containing: - - Atoms: The ASE Atoms object with optimized positions - - list or None: The optimization log if verbose is True, otherwise None + """Generate a random packed atomic structure with minimal atomic overlaps. + + Parameters + ---------- + composition : Composition + Pymatgen Composition object specifying atomic composition (e.g. Fe80B20). + Numbers indicate actual atom counts. + cell : np.ndarray + 3x3 array defining triclinic simulation box in Angstroms. + seed : int, optional + Random seed for reproducible generation, by default 42. + If None, uses random initialization. + diameter : float, optional + Minimum allowed interatomic distance for overlap detection. + Used for soft-sphere potential, by default None. + auto_diameter : bool, optional + If True, automatically calculate diameter from atomic radii, by default False. + max_iter : int, optional + Maximum FIRE optimization steps to minimize overlaps, by default 100. + fmax : float, optional + Maximum force criterion for convergence in eV/Å, by default 0.01. + distance_tolerance : float, optional + Distance threshold for considering atoms at same position, by default 0.0001. + trajectory_file : str, optional + Path to save optimization trajectory (.traj format), by default None. + verbose : bool, optional + Print progress information during optimization, by default True. + + Returns + ------- + tuple[Atoms, Optional[list[dict]]] + - ASE Atoms object with optimized positions + - Optimization log if verbose=True, otherwise None """ # Extract number of atoms for each element from composition element_counts = [int(i) for i in composition.as_dict().values()] @@ -310,12 +295,6 @@ def valid_subcell( Returns: bool: True if the structure passes all validation checks, False otherwise. - - Notes: - - The function uses periodic boundary conditions when checking atomic distances - - Formation energies are typically negative but not extremely negative - - The optimization should reduce energy unless already at a minimum - - Atomic fusion (distances < ~1.5 Å) indicates an unphysical structure """ # Check if formation energy is unphysically negative if final_energy < fe_lower_limit: @@ -356,29 +335,24 @@ def subcells_to_structures( ) -> list[tuple[np.ndarray, np.ndarray, list[str]]]: """Convert subcell candidates to structure tuples. - Transforms subcell regions defined by atom IDs and boundary coordinates into - normalized crystallographic structures. Each subcell is rescaled to have fractional - coordinates in the [0,1] range with an appropriately scaled cell. - - Args: - candidates: List of (ids, lower_bound, upper_bound) tuples representing subcell - candidates. Each tuple contains NumPy arrays of atom indices and the lower/upper - boundaries of the subcell in fractional coordinates. - fractional_positions: NumPy array of shape [n_atoms, 3] containing fractional - coordinates of all atoms in the parent structure. - cell: 3x3 NumPy array representing the unit cell of the parent structure. - species: List of atomic species symbols for all atoms in the parent structure. - - Returns: - List of (fractional_positions, cell, species) tuples for each subcell, where: - - fractional_positions: NumPy array of normalized [0,1] fractional coordinates - - cell: 3x3 NumPy array representing the scaled subcell - - species: List of chemical symbols for the atoms in the subcell - - Notes: - - The function normalizes coordinates to start from origin [0,0,0] - - Cell vectors are scaled according to the size of the subcell region - - No filtering of valid subcells is performed in this function + Parameters + ---------- + candidates : list[tuple[np.ndarray, np.ndarray, np.ndarray]] + List of (atom_indices, lower_bound, upper_bound) tuples defining subcell regions + fractional_positions : np.ndarray + Fractional coordinates of all atoms in parent structure, shape [n_atoms, 3] + cell : np.ndarray + 3x3 array representing parent unit cell + species : list[str] + Chemical symbols for all atoms in parent structure + + Returns + ------- + list[tuple[np.ndarray, np.ndarray, list[str]]] + List of (frac_coords, cell, species) tuples for each subcell: + - frac_coords: Normalized [0,1] fractional coordinates + - cell: Scaled 3x3 subcell + - species: Chemical symbols for subcell atoms """ list_subcells = [] for ids, l, h in candidates: # noqa: E741 @@ -405,38 +379,30 @@ def subcells_to_structures( def get_target_temperature( step: int, equi_steps: int, cool_steps: int, T_high: float, T_low: float ) -> float: - """Calculate target temperature for each step in a melt-quench-equilibrate simulation. - - Implements a three-stage temperature profile commonly used in materials simulations: - 1. Equilibration: Maintains a constant high temperature (T_high) for 'equi_steps' - 2. Cooling: Linearly decreases temperature from T_high to T_low over 'cool_steps' - 3. Final equilibration: Maintains the low temperature (T_low) for all remaining steps - - This temperature profile simulates the physical process of melting a material at high - temperature, then cooling it to create amorphous or crystalline structures. + """Calculate target temperature for a melt-quench-equilibrate simulation step. Args: - step: Current simulation step (zero-indexed). - equi_steps: Number of initial equilibration steps at high temperature. - cool_steps: Number of steps over which the temperature is linearly decreased. - T_high: Initial high temperature in Kelvin for equilibration phase. - T_low: Final low temperature in Kelvin for the quenched state. + step (int): Current simulation step number (0-indexed) + equi_steps (int): Number of steps for initial high-temperature equilibration + cool_steps (int): Number of steps for linear cooling + T_high (float): Initial high temperature in Kelvin + T_low (float): Final low temperature in Kelvin Returns: - float: Target temperature in Kelvin for the current simulation step. + float: Target temperature in Kelvin for the current step - Examples: - >>> # During high-temperature equilibration - >>> temp = get_target_temperature(10, 100, 200, 2000.0, 300.0) - >>> print(f"{temp:.1f} K") # Returns 2000.0 K + The temperature profile has three phases: + 1. Hold at T_high for equi_steps + 2. Cool linearly from T_high to T_low over cool_steps + 3. Hold at T_low for remaining steps - >>> # During cooling phase (halfway through cooling) - >>> temp = get_target_temperature(200, 100, 200, 2000.0, 300.0) - >>> print(f"{temp:.1f} K") # Returns ~1150.0 K - - >>> # After cooling complete - >>> temp = get_target_temperature(350, 100, 200, 2000.0, 300.0) - >>> print(f"{temp:.1f} K") # Returns 300.0 K + Examples: + >>> get_target_temperature(10, 100, 200, 2000.0, 300.0) # During equilibration + 2000.0 + >>> get_target_temperature(200, 100, 200, 2000.0, 300.0) # During cooling + 1150.0 + >>> get_target_temperature(350, 100, 200, 2000.0, 300.0) # After cooling + 300.0 """ # Initial high-temperature equilibration phase if step < equi_steps: @@ -464,41 +430,30 @@ def get_subcells_to_crystallize( ) -> list[tuple[np.ndarray, np.ndarray, np.ndarray]]: """Extract subcell structures from a larger structure for crystallization. - This function systematically divides a large structure (typically amorphous) into - smaller subcells that can be independently relaxed to find stable crystal structures. - The approach creates a grid in fractional coordinates and identifies groups of atoms - within each grid cell that meet size and composition criteria. + Divides a structure into overlapping subcells that can be relaxed to find stable crystal + structures. Uses a grid in fractional coordinates to identify atom groups meeting size and + composition criteria. Args: - fractional_positions: NumPy array of shape [n_atoms, 3] containing fractional - coordinates of all atoms in the parent structure. - species: List of chemical element symbols corresponding to each atom. - d_frac: Grid spacing in fractional coordinates used to define subcell boundaries. - Smaller values create more overlap between subcells. Defaults to 0.05. - n_min: Minimum number of atoms required in a subcell. Defaults to 1. - n_max: Maximum number of atoms allowed in a subcell (controls computational - cost during relaxation). Defaults to 48. - restrict_to_compositions: Optional list of chemical formulas (e.g. ["AB", "AB2"]) - to filter subcells by composition. Only subcells with matching reduced - formulas will be returned. Defaults to None (no filtering). - max_coeff: Maximum stoichiometric coefficient to consider when automatically - generating composition restrictions. For example, max_coeff=2 would include - formulas like AB2 but not AB3. Defaults to None. - elements: List of elements to consider when generating stoichiometries. - Required if max_coeff is provided. Defaults to None. + fractional_positions (np.ndarray): Fractional coordinates of atoms, shape [n_atoms, 3]. + species (list[str]): Chemical element symbols for each atom. + d_frac (float, optional): Grid spacing in fractional coordinates. Smaller values create + more overlap. Defaults to 0.05. + n_min (int, optional): Minimum atoms per subcell. Defaults to 1. + n_max (int, optional): Maximum atoms per subcell. Defaults to 48. + restrict_to_compositions (Sequence[str], optional): Chemical formulas to filter subcells by + (e.g. ["AB", "AB2"]). Only matching compositions are returned. Defaults to None. + max_coeff (int, optional): Maximum stoichiometric coefficient for auto-generating + restrictions. E.g. max_coeff=2 allows AB2 but not AB3. Defaults to None. + elements (Sequence[str], optional): Elements for generating stoichiometries. Required if + max_coeff provided. Defaults to None. Returns: - List of tuples, where each tuple contains: - - indices: NumPy array of atom indices included in the subcell - - lower: NumPy array of lower bounds for subcell in fractional coords [3] - - upper: NumPy array of upper bounds for subcell in fractional coords [3] - - Notes: - - Overlapping subcells provide redundancy in structure discovery - - Using composition restrictions can significantly reduce the number of - candidates and focus on chemically relevant structures - - For multi-component systems, setting appropriate n_min and n_max values - helps ensure reasonable stoichiometries + list[tuple[np.ndarray, np.ndarray, np.ndarray]]: List of (indices, lower_bounds, + upper_bounds) where: + - indices: Atom indices in subcell + - lower_bounds: Lower bounds in fractional coords [3] + - upper_bounds: Upper bounds in fractional coords [3] """ # Convert species list to numpy array for easier composition handling species_array = np.array(species) @@ -562,15 +517,21 @@ def default_subcell_filter( cubic_only: bool = True, allowed_atom_counts: Optional[list[int]] = None, ) -> bool: - """Default filter function for subcells. - - Args: - subcell: A tuple containing (indices, lower_bound, upper_bound) - cubic_only: If True, only accept subcells where all dimensions are equal - allowed_atom_counts: If provided, only accept subcells with these atom counts - - Returns: - bool: True if subcell passes all filters, False otherwise + """Filter subcells based on shape and size criteria. + + Parameters + ---------- + subcell : tuple[np.ndarray, np.ndarray, np.ndarray] + Tuple containing (atom_indices, lower_bound, upper_bound) arrays that define the subcell + cubic_only : bool, default=True + If True, only accept subcells with equal dimensions in x, y, and z + allowed_atom_counts : list[int], optional + List of allowed atom counts. If provided, only accept subcells with specified atom counts + + Returns + ------- + bool + True if subcell meets all criteria, False otherwise """ indices, lower_bound, upper_bound = subcell cell_dimensions = upper_bound - lower_bound @@ -596,28 +557,33 @@ def extract_crystallizable_subcells( allowed_atom_counts: Optional[list[int]] = None, filter_function: Optional[Callable] = None, ) -> list[Atoms]: - """ - Extract and filter subcells from an amorphous structure for crystallization analysis. - - This function divides an amorphous structure into overlapping subcells, then filters them - based on shape criteria (e.g., cubic cells) and atom count to identify potential - crystalline structural motifs. - - Args: - atoms: ASE Atoms object representing the amorphous structure - d_frac: Grid spacing in fractional coordinates for subcell division. Defaults to 0.2. - n_min: Minimum number of atoms required in a subcell. Defaults to 2. - n_max: Maximum number of atoms allowed in a subcell. Defaults to 8. - cubic_only: If True, only keep subcells where all dimensions are equal. Defaults to True. - allowed_atom_counts: If provided, only keep subcells with these atom counts. - Defaults to None (no filtering by atom count). - filter_function: Optional callable that takes a subcell tuple and returns a boolean. - If provided, this function is used to filter subcells instead of the default filter. - The function should accept a tuple (indices, lower_bound, upper_bound) and return - True if the subcell should be kept, False otherwise. - - Returns: - list[Atoms]: List of ASE Atoms objects representing the filtered subcells + """Extract and filter subcells from an amorphous structure for crystallization analysis. + + Divides an amorphous structure into overlapping subcells and filters them based on shape + criteria and atom count to identify potential crystalline structural motifs. + + Parameters + ---------- + atoms : Atoms + ASE Atoms object representing the amorphous structure + d_frac : float, optional + Grid spacing in fractional coordinates for subcell division, by default 0.2 + n_min : int, optional + Minimum number of atoms required in a subcell, by default 2 + n_max : int, optional + Maximum number of atoms allowed in a subcell, by default 8 + cubic_only : bool, optional + If True, only keep subcells where all dimensions are equal, by default True + allowed_atom_counts : list[int], optional + If provided, only keep subcells with these atom counts, by default None + filter_function : callable, optional + Custom filter function that takes a subcell tuple (indices, lower_bound, upper_bound) + and returns a boolean. If provided, used instead of default filter, by default None + + Returns + ------- + list[Atoms] + List of ASE Atoms objects representing the filtered subcells """ # Get fractional positions from the Atoms object cell = atoms.get_cell() diff --git a/example/Si64.py b/example/Si64.py index 89038ca..f6e729c 100644 --- a/example/Si64.py +++ b/example/Si64.py @@ -17,32 +17,32 @@ from ase.build import bulk from mace.calculators.foundations_models import mace_mp -# Define system +# System configuration comp = Composition("Si64") cell = np.array([[11.1, 0.0, 0.0], [0.0, 11.1, 0.0], [0.0, 0.0, 11.1]]) -# Relaxation parameters +# Optimization parameters global_seed = 42 -fmax = 0.01 +fmax = 0.01 # Force convergence criterion in eV/Å max_iter = 100 -# MD parameters +# Molecular dynamics parameters md_log_interval = 50 -md_equi_steps = 2500 -md_cool_steps = 2500 -md_final_steps = 2500 -md_T_high = 2000.0 -md_T_low = 300.0 -md_time_step = 2.0 -md_friction = 0.01 - -# Define model +md_equi_steps = 2500 # High temperature equilibration steps +md_cool_steps = 2500 # Cooling steps +md_final_steps = 2500 # Low temperature equilibration steps +md_T_high = 2000.0 # Initial melting temperature (K) +md_T_low = 300.0 # Final temperature (K) +md_time_step = 2.0 # fs +md_friction = 0.01 # Langevin friction + +# Initialize MACE calculator mace_checkpoint_url = ( "https://github.com/ACEsuit/mace-mp/releases/download/mace_omat_0/mace-omat-0-medium.model" ) calculator = mace_mp(model=mace_checkpoint_url, device="cuda", dtype="float32", enable_cueq=False) -# Generate random packed structure +# Generate initial random structure packed_atoms, log_data = random_packed_structure( composition=comp, cell=cell, @@ -73,13 +73,12 @@ ) print(f"Final amorphous structure is ready {amorphous_atoms}") -# After running melt-quench MD to create amorphous structure -# Extract crystallizable subcells +# Extract potential crystalline regions crystallizable_cells = extract_crystallizable_subcells( atoms=amorphous_atoms, - d_frac=0.2, - n_min=2, - n_max=8, + d_frac=0.2, # Size of subcell as fraction of original cell + n_min=2, # Min atoms per subcell + n_max=8, # Max atoms per subcell cubic_only=True, allowed_atom_counts=[2, 4, 8], ) @@ -101,7 +100,7 @@ # Store the relaxed structure and its properties relaxed_atoms_list.append((relaxed_atoms, logger, final_energy, final_pressure)) -# Find the structure with the lowest energy per atom +# Find lowest energy structure lowest_e_candidate = min(relaxed_atoms_list, key=lambda x: x[-2] / len(x[0])) lowest_e_atoms, lowest_e_logger, lowest_e_energy, lowest_e_pressure = lowest_e_candidate @@ -121,6 +120,7 @@ ) print(f"Final pressure: {lowest_e_pressure:.6f} eV/ų") +# Count frequency of space groups across all candidates spg_counter = defaultdict(lambda: 0) for s in relaxed_atoms_list: @@ -138,6 +138,7 @@ print("All space groups encountered:", dict(spg_counter)) +# Compare to reference diamond structure si_diamond = bulk("Si", "diamond", a=5.43) pymatgen_ref_struct = Structure( lattice=si_diamond.get_cell(), @@ -147,6 +148,6 @@ ) print("Prediction matches diamond-cubic Si?", StructureMatcher().fit(pymatgen_struct, pymatgen_ref_struct)) -# Save the best structure +# Save the lowest energy structure lowest_e_atoms.write("final_crystal_structure.cif") print("Lowest energy structure saved to final_crystal_structure.cif")