From 0bb362445064aaaf03e63ef8104aca2b76f31339 Mon Sep 17 00:00:00 2001 From: Prabhath C Date: Mon, 16 Mar 2026 00:11:56 +0100 Subject: [PATCH 1/3] Change potential to dataframe, add few type hints --- src/atomistics/calculators/lammps/melting.py | 49 ++++++++++++-------- 1 file changed, 29 insertions(+), 20 deletions(-) diff --git a/src/atomistics/calculators/lammps/melting.py b/src/atomistics/calculators/lammps/melting.py index 460411ec..2cd36b6a 100644 --- a/src/atomistics/calculators/lammps/melting.py +++ b/src/atomistics/calculators/lammps/melting.py @@ -10,6 +10,7 @@ calc_molecular_dynamics_npt_with_lammpslib, ) import numpy as np +import pandas as pd from structuretoolkit.analyse import ( get_adaptive_cna_descriptors, get_diamond_structure_descriptors, @@ -39,7 +40,6 @@ def _check_diamond(structure: Atoms) -> bool: > dia_dict["IdentifyDiamond.counts.OTHER"] ) - def _analyse_structure( structure: Atoms, mode: str = "total", diamond: bool = False ) -> dict: @@ -63,7 +63,6 @@ def _analyse_structure( structure=structure, mode=mode, ovito_compatibility=True ) - def _analyse_minimized_structure(structure: Atoms) -> tuple: """ @@ -93,18 +92,25 @@ def _analyse_minimized_structure(structure: Atoms) -> tuple: final_structure_dict, ) - -def _next_calc(structure, potential, temperature, seed, run_time_steps=10000): +def _next_calc( + structure: Atoms, + potential_dataframe: pd.DataFrame, + temperature: float, + seed: int, + run_time_steps: int = 10000 +) -> Atoms: """ Calculate NPT ensemble at a given temperature using the job defined in the project parameters: - job_type: Type of Simulation code to be used - project: Project object used to create the job - - potential: Interatomic Potential + - potential_dataframe: Interatomic Potential dataframe - queue (optional): HPC Job queue to be used Args: structure (pyiron_atomistics.structure.atoms.Atoms): Atomistic Structure object to be set to the job as input sturcture + potential_dataframe (pd.DataFrame): Interatomic Potential dataframe temperature (float): Temperature of the Molecular dynamics calculation + seed (int): Random seed for the simulation run_time_steps (int): Number of Molecular dynamics steps Returns: @@ -112,7 +118,7 @@ def _next_calc(structure, potential, temperature, seed, run_time_steps=10000): """ output_md_dict = calc_molecular_dynamics_npt_with_lammpslib( structure=structure, - potential_dataframe=potential, + potential_dataframe=potential_dataframe, Tstart=temperature, Tstop=temperature, Tdamp=0.1, @@ -131,15 +137,15 @@ def _next_calc(structure, potential, temperature, seed, run_time_steps=10000): structure_md = structure.copy() structure_md.set_positions(output_md_dict["positions"][-1]) structure_md.set_cell(output_md_dict["cell"][-1]) + return structure_md - def _next_step_funct( number_of_atoms, key_max, structure_left, structure_right, - potential, + potential_dataframe: pd.DataFrame, temperature_left, temperature_right, distribution_initial_half, @@ -160,19 +166,22 @@ def _next_step_funct( distribution_initial_half: structure_after_minimization: run_time_steps: + seed: Returns: """ + diamond_flag = crystalstructure.lower() == "diamond" + structure_left_dict = _analyse_structure( structure=structure_left, mode="total", - diamond=crystalstructure.lower() == "diamond", + diamond=diamond_flag, ) structure_right_dict = _analyse_structure( structure=structure_right, mode="total", - diamond=crystalstructure.lower() == "diamond", + diamond=diamond_flag, ) temperature_diff = temperature_right - temperature_left if ( @@ -185,7 +194,7 @@ def _next_step_funct( structure_right = _next_calc( structure=structure_after_minimization, temperature=temperature_right, - potential=potential, + potential_dataframe=potential_dataframe, seed=seed, run_time_steps=run_time_steps, ) @@ -199,7 +208,7 @@ def _next_step_funct( structure_left = _next_calc( structure=structure_after_minimization, temperature=temperature_left, - potential=potential, + potential_dataframe=potential_dataframe, seed=seed, run_time_steps=run_time_steps, ) @@ -214,7 +223,7 @@ def _next_step_funct( structure_left = _next_calc( structure=structure_after_minimization, temperature=temperature_left, - potential=potential, + potential_dataframe=potential_dataframe, seed=seed, run_time_steps=run_time_steps, ) @@ -222,10 +231,9 @@ def _next_step_funct( raise ValueError("We should never reach this point!") return structure_left, structure_right, temperature_left, temperature_right - -def estimate_melting_temperature( +def estimate_melting_temperature_using_bisection_CNA( + potential_dataframe: pd.DataFrame, element, - potential, strain_run_time_steps=1000, temperature_left=0, temperature_right=1000, @@ -246,7 +254,7 @@ def estimate_melting_temperature( structure_opt = optimize_positions_and_volume_with_lammpslib( structure=basis, - potential_dataframe=potential, + potential_dataframe=potential_dataframe, min_style="cg", etol=0.0, ftol=0.0001, @@ -269,7 +277,7 @@ def estimate_melting_temperature( structure=structure_after_minimization, temperature=temperature_right, seed=seed, - potential=potential, + potential_dataframe=potential_dataframe, run_time_steps=strain_run_time_steps, ) temperature_step = temperature_right - temperature_left @@ -285,7 +293,7 @@ def estimate_melting_temperature( key_max=key_max, structure_left=structure_left, structure_right=structure_right, - potential=potential, + potential_dataframe=potential_dataframe, temperature_left=temperature_left, temperature_right=temperature_right, distribution_initial_half=distribution_initial_half, @@ -295,4 +303,5 @@ def estimate_melting_temperature( crystalstructure=crystalstructure, ) temperature_step = temperature_right - temperature_left - return int(round(temperature_left)) + + return int(round(temperature_left)) \ No newline at end of file From 563f3ed39eafbfa6b401824a9cd49825ee294d15 Mon Sep 17 00:00:00 2001 From: Prabhath C Date: Mon, 16 Mar 2026 00:22:20 +0100 Subject: [PATCH 2/3] Take Atoms as input for bisection, add helper function for basis size --- src/atomistics/calculators/lammps/melting.py | 64 ++++++++++++++------ 1 file changed, 45 insertions(+), 19 deletions(-) diff --git a/src/atomistics/calculators/lammps/melting.py b/src/atomistics/calculators/lammps/melting.py index 2cd36b6a..fa54dfba 100644 --- a/src/atomistics/calculators/lammps/melting.py +++ b/src/atomistics/calculators/lammps/melting.py @@ -1,5 +1,6 @@ import operator import random +import re from ase import Atoms from ase.build import bulk @@ -63,7 +64,10 @@ def _analyse_structure( structure=structure, mode=mode, ovito_compatibility=True ) -def _analyse_minimized_structure(structure: Atoms) -> tuple: +def _analyse_minimized_structure( + structure: Atoms, + diamond_flag: bool +) -> tuple: """ Args: @@ -75,7 +79,6 @@ def _analyse_minimized_structure(structure: Atoms) -> tuple: and the final structure dictionary. """ - diamond_flag = _check_diamond(structure=structure) final_structure_dict = _analyse_structure( structure=structure, mode="total", diamond=diamond_flag ) @@ -92,6 +95,30 @@ def _analyse_minimized_structure(structure: Atoms) -> tuple: final_structure_dict, ) +def _get_repeated_structure( + structure: Atoms, target_number_of_atoms: int +) -> Atoms: + """ + Get a repeated structure that is as close as possible to the target number of atoms. + + Args: + structure (Atoms): The input structure to be repeated. + target_number_of_atoms (int): The target number of atoms for the simulation cell. + + Returns: + Atoms: The repeated structure. + """ + r_est = (target_number_of_atoms / len(structure)) ** (1 / 3) + candidates = np.array( + [max(1, int(np.floor(r_est))), int(np.round(r_est)), int(np.ceil(r_est))] + ) + basis_lst = [structure.repeat([i, i, i]) for i in candidates] + basis = basis_lst[ + np.argmin([np.abs(len(b) - target_number_of_atoms) for b in basis_lst]) + ] + + return basis + def _next_calc( structure: Atoms, potential_dataframe: pd.DataFrame, @@ -151,7 +178,7 @@ def _next_step_funct( distribution_initial_half, structure_after_minimization, run_time_steps, - crystalstructure, + diamond_flag, seed, ): """ @@ -171,8 +198,6 @@ def _next_step_funct( Returns: """ - diamond_flag = crystalstructure.lower() == "diamond" - structure_left_dict = _analyse_structure( structure=structure_left, mode="total", @@ -232,28 +257,26 @@ def _next_step_funct( return structure_left, structure_right, temperature_left, temperature_right def estimate_melting_temperature_using_bisection_CNA( + structure: Atoms, potential_dataframe: pd.DataFrame, - element, + target_number_of_atoms: int, strain_run_time_steps=1000, temperature_left=0, temperature_right=1000, number_of_atoms=8000, seed=None, ): + if seed is None: seed = random.randint(0, 99999) - crystalstructure = reference_states[atomic_numbers[element]]["symmetry"] - if crystalstructure == "hcp": - basis = bulk(name=element, orthorhombic=True) - else: - basis = bulk(name=element, cubic=True) - basis_lst = [basis.repeat([i, i, i]) for i in range(5, 30)] - basis = basis_lst[ - np.argmin([np.abs(len(b) - number_of_atoms / 2) for b in basis_lst]) - ] + + diamond_flag = _check_diamond(structure=structure) + repeated_structure = _get_repeated_structure( + structure=structure, target_number_of_atoms=target_number_of_atoms + ) - structure_opt = optimize_positions_and_volume_with_lammpslib( - structure=basis, + position_and_volume_optimized_structure = optimize_positions_and_volume_with_lammpslib( + structure=repeated_structure, potential_dataframe=potential_dataframe, min_style="cg", etol=0.0, @@ -270,7 +293,10 @@ def estimate_melting_temperature_using_bisection_CNA( number_of_atoms, distribution_initial_half, _, - ) = _analyse_minimized_structure(structure=structure_opt) + ) = _analyse_minimized_structure( + structure=position_and_volume_optimized_structure, + diamond_flag=diamond_flag + ) structure_left = structure_after_minimization structure_right = _next_calc( @@ -300,7 +326,7 @@ def estimate_melting_temperature_using_bisection_CNA( structure_after_minimization=structure_after_minimization, run_time_steps=strain_run_time_steps, seed=seed, - crystalstructure=crystalstructure, + diamond_flag=diamond_flag, ) temperature_step = temperature_right - temperature_left From 74899d3416439ce5ca05f6c285348ae68441771c Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sun, 15 Mar 2026 23:23:56 +0000 Subject: [PATCH 3/3] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/atomistics/calculators/lammps/melting.py | 58 ++++++++++---------- 1 file changed, 30 insertions(+), 28 deletions(-) diff --git a/src/atomistics/calculators/lammps/melting.py b/src/atomistics/calculators/lammps/melting.py index fa54dfba..990c45b2 100644 --- a/src/atomistics/calculators/lammps/melting.py +++ b/src/atomistics/calculators/lammps/melting.py @@ -41,6 +41,7 @@ def _check_diamond(structure: Atoms) -> bool: > dia_dict["IdentifyDiamond.counts.OTHER"] ) + def _analyse_structure( structure: Atoms, mode: str = "total", diamond: bool = False ) -> dict: @@ -64,10 +65,8 @@ def _analyse_structure( structure=structure, mode=mode, ovito_compatibility=True ) -def _analyse_minimized_structure( - structure: Atoms, - diamond_flag: bool -) -> tuple: + +def _analyse_minimized_structure(structure: Atoms, diamond_flag: bool) -> tuple: """ Args: @@ -95,9 +94,8 @@ def _analyse_minimized_structure( final_structure_dict, ) -def _get_repeated_structure( - structure: Atoms, target_number_of_atoms: int -) -> Atoms: + +def _get_repeated_structure(structure: Atoms, target_number_of_atoms: int) -> Atoms: """ Get a repeated structure that is as close as possible to the target number of atoms. @@ -119,12 +117,13 @@ def _get_repeated_structure( return basis + def _next_calc( - structure: Atoms, - potential_dataframe: pd.DataFrame, - temperature: float, - seed: int, - run_time_steps: int = 10000 + structure: Atoms, + potential_dataframe: pd.DataFrame, + temperature: float, + seed: int, + run_time_steps: int = 10000, ) -> Atoms: """ Calculate NPT ensemble at a given temperature using the job defined in the project parameters: @@ -164,9 +163,10 @@ def _next_calc( structure_md = structure.copy() structure_md.set_positions(output_md_dict["positions"][-1]) structure_md.set_cell(output_md_dict["cell"][-1]) - + return structure_md + def _next_step_funct( number_of_atoms, key_max, @@ -256,6 +256,7 @@ def _next_step_funct( raise ValueError("We should never reach this point!") return structure_left, structure_right, temperature_left, temperature_right + def estimate_melting_temperature_using_bisection_CNA( structure: Atoms, potential_dataframe: pd.DataFrame, @@ -266,25 +267,27 @@ def estimate_melting_temperature_using_bisection_CNA( number_of_atoms=8000, seed=None, ): - + if seed is None: seed = random.randint(0, 99999) - + diamond_flag = _check_diamond(structure=structure) repeated_structure = _get_repeated_structure( structure=structure, target_number_of_atoms=target_number_of_atoms ) - position_and_volume_optimized_structure = optimize_positions_and_volume_with_lammpslib( - structure=repeated_structure, - potential_dataframe=potential_dataframe, - min_style="cg", - etol=0.0, - ftol=0.0001, - maxiter=100000, - maxeval=10000000, - thermo=10, - lmp=None, + position_and_volume_optimized_structure = ( + optimize_positions_and_volume_with_lammpslib( + structure=repeated_structure, + potential_dataframe=potential_dataframe, + min_style="cg", + etol=0.0, + ftol=0.0001, + maxiter=100000, + maxeval=10000000, + thermo=10, + lmp=None, + ) ) ( @@ -294,8 +297,7 @@ def estimate_melting_temperature_using_bisection_CNA( distribution_initial_half, _, ) = _analyse_minimized_structure( - structure=position_and_volume_optimized_structure, - diamond_flag=diamond_flag + structure=position_and_volume_optimized_structure, diamond_flag=diamond_flag ) structure_left = structure_after_minimization @@ -330,4 +332,4 @@ def estimate_melting_temperature_using_bisection_CNA( ) temperature_step = temperature_right - temperature_left - return int(round(temperature_left)) \ No newline at end of file + return int(round(temperature_left))