Merge pull request #5 from openkim-hackathons/ilia-updates-May2026

Ilia updates may2026

Author: ilia-nikiforov-umn

.gitattributes

@@ -1,3 +1,4 @@
 output export-ignore
 local-props export-ignore
 run.py export-ignore
+job.txt export-ignore

.gitignore

@@ -1,4 +1,7 @@
-output/*
+output*
+log.*
+*.log
+*.out
 
 # Byte-compiled / optimized / DLL files
 __pycache__/

README.rst

@@ -8,7 +8,7 @@ This test driver repeats the unit cell to build a supercell and then runs a mole
 NPT ensemble using Lammps.
 
 This test driver uses kim-convergence to detect an equilibrated molecular-dynamics simulation. It checks
-convergence of the volume, temperature, enthalpy and cell shape parameters every 10000 timesteps.
+convergence of the volume, temperature, enthalpy and cell shape parameters every 10000 timesteps (default).
 
 During the equilibrated part of the simulation, the test driver averages the cell parameters and atomic positions to
 obtain the equilibrium crystal structure. This includes an average over time, and an average over the replicated unit

output/.gitkeep

@@ -1,9 +1,19 @@
-import argparse
-import subprocess
-from kim_tools.test_driver.core import query_crystal_structures
 from test_driver.test_driver import TestDriver
+from ase.build import bulk
 
+atoms = bulk("Al")
+td=TestDriver("EAM_Dynamo_ErcolessiAdams_1994_Al__MO_123629422045_006")
+print(td(
+        atoms,
+        temperature_K=300,
+        target_size=500,
+        lammps_command= "lmp",
+))
 
+
+
+
+"""
 if __name__ == '__main__':
     # Argument parsing
     parser = argparse.ArgumentParser(description='Pass arguments to the KIM test driver')
@@ -35,3 +45,4 @@
         except Exception as e:
             print(f"Got exception {repr(e)}")
     test_driver.write_property_instances_to_file()
+"""

run.py

@@ -19,8 +19,15 @@ if temperature_K != "":
 else:
     print("No temperature given")
     temperature_K = 0
+pressure_eV_angstrom3 = input("pressure (eV/angstrom^3)?\n")
+if pressure_eV_angstrom3 != "":
+    pressure_eV_angstrom3 = float(pressure_eV_angstrom3)
+    print(pressure_eV_angstrom3)
+else:
+    print("No pressure given")
+    pressure_eV_angstrom3 = None
 cell_cauchy_stress_eV_angstrom3 = input(
-    "Cauchy stress (literal list of floats, Voigt order xx,yy,zz,yz,xz,xy, eV/A^3)?\n"
+    "Cauchy stress (literal list of floats, Voigt order xx,yy,zz,yz,xz,xy, eV/angstrom^3)?\n"
 ).split()
 if cell_cauchy_stress_eV_angstrom3 != []:
     cell_cauchy_stress_eV_angstrom3 = [
@@ -29,7 +36,7 @@ if cell_cauchy_stress_eV_angstrom3 != []:
     print(cell_cauchy_stress_eV_angstrom3)
 else:
     print("No stress given")
-    cell_cauchy_stress_eV_angstrom3 = [0, 0, 0, 0, 0, 0]
+    cell_cauchy_stress_eV_angstrom3 = None
 runtime_args_text = input("Runtime arguments (literal dictonary)?\n")
 if runtime_args_text != "":
     runtime_args = literal_eval(runtime_args_text)
@@ -42,24 +49,40 @@ query_result = literal_eval(
 )
 print(json.dumps(query_result, indent=2))
 test = TestDriver(model_name)
+successful_run = False
 for structure in query_result:
-    species_match = (
-        structure["stoichiometric-species"]["source-value"]
-        == query_result[0]["stoichiometric-species"]["source-value"]
-    )
-    prototype_match = (
-        structure["prototype-label"]["source-value"]
-        == query_result[0]["prototype-label"]["source-value"]
-    )
-    assert species_match and prototype_match, (
-        "The input structures must have identical "
-        "prototype-label and stoichiometric-species"
-    )
-    test(
-        structure,
-        temperature_K=temperature_K,
-        cell_cauchy_stress_eV_angstrom3=cell_cauchy_stress_eV_angstrom3,
-        **runtime_args
-    )
+    try:
+        species_match = (
+            structure["stoichiometric-species"]["source-value"]
+            == query_result[0]["stoichiometric-species"]["source-value"]
+        )
+        prototype_match = (
+            structure["prototype-label"]["source-value"]
+            == query_result[0]["prototype-label"]["source-value"]
+        )
+        assert species_match and prototype_match, (
+            "The input structures must have identical "
+            "prototype-label and stoichiometric-species"
+        )
+        test(
+            structure,
+            temperature_K=temperature_K,
+            pressure_eV_angstrom3=pressure_eV_angstrom3,
+            cell_cauchy_stress_eV_angstrom3=cell_cauchy_stress_eV_angstrom3,
+            **runtime_args
+        )
+        successful_run = True
+    except Exception as e:
+        print(
+            "The following exception prevented one of the parameter sets from "
+            f"returning a result: \n{repr(e)}"
+        )
+        exc = e
+if not successful_run:
+    raise RuntimeError(
+        "No parameter sets successfully ran. Last exception is attached, "
+        "see output file for others, if any."
+    ) from exc
+test.deduplicate_property_instances(allow_rotation=False)
 test.write_property_instances_to_file()
 copy("kim-tools.log", "output/kim-tools.log")

runner

@@ -9,9 +9,10 @@
 
 
 def run_lammps(modelname: str, temperature_K: float, pressure_bar: float, timestep_ps: float,
-               number_sampling_timesteps: int, species: List[str],
+               thermo_sampling_period: int, species: List[str],
                msd_threshold_angstrom_squared_per_sampling_timesteps: float, number_msd_timesteps: int,
-               lammps_command: str, random_seed: int) -> Tuple[str, str, str, str]:
+               rlc_run_length: int, rlc_n_every: int, output_dir: str, equilibration_plots: bool, lammps_command: str,
+               random_seed: int) -> Tuple[str, str, str, str, str]:
     """
     Run LAMMPS NPT simulation with the given parameters.
 
@@ -33,9 +34,9 @@ def run_lammps(modelname: str, temperature_K: float, pressure_bar: float, timest
     :param timestep_ps:
         Timestep in picoseconds.
     :type timestep_ps: float
-    :param number_sampling_timesteps:
+    :param thermo_sampling_period:
         Number of timesteps for sampling thermodynamic quantities.
-    :type number_sampling_timesteps: int
+    :type thermo_sampling_period: int
     :param species:
         List of chemical species in the system.
     :type species: List[str]
@@ -47,6 +48,19 @@ def run_lammps(modelname: str, temperature_K: float, pressure_bar: float, timest
         Before the mean-squared displacement is monitored, the system will be equilibrated for the same number of
         timesteps.
     :type number_msd_timesteps: int
+    :param rlc_run_length:
+        Number of timesteps after which kim-convergence will check for convergence.
+        This is also the timestep interval for generated trajectories.
+    :type rlc_run_length: int
+    :param rlc_n_every:
+        Number of timesteps between storage of values for the run-length control in kim-convergence.
+    :type rlc_n_every: int
+    :param output_dir:
+        Directory to store the output files.
+    :type output_dir: str
+    :param equilibration_plots:
+        Whether to plot the equilibration plots.
+    :type equilibration_plots: bool
     :param lammps_command:
         Command to run LAMMPS (e.g., "mpirun -np 4 lmp_mpi" or "lmp").
     :type lammps_command: str
@@ -55,15 +69,19 @@ def run_lammps(modelname: str, temperature_K: float, pressure_bar: float, timest
     :type random_seed: int
 
     :return:
-        A tuple containing paths to the LAMMPS log file, restart file, full average position file, and full average cell
-        file.
-    :rtype: Tuple[str, str, str, str]
+        A tuple containing paths to the LAMMPS log file, restart file, full average position file, full average cell
+        file, and melted crystal dump file (only exists if crystal melted).
+    :rtype: Tuple[str, str, str, str, str]
     """
-    pdamp = timestep_ps * 100.0
-    tdamp = timestep_ps * 1000.0
-
-    log_filename = "output/lammps.log"
-    restart_filename = "output/final_configuration.restart"
+    pdamp = timestep_ps * 1000.0
+    tdamp = timestep_ps * 100.0
+
+    # Lammps will be run directly in output_dir so all paths are with respect to that directory.
+    log_filename = "lammps.log"
+    restart_filename = "final_configuration.restart"
+    melted_crystal_filename = "melted_crystal.dump"
+    average_position_filename = "average_position.dump"
+    average_cell_filename = "average_cell.dump"
     variables = {
         "modelname": modelname,
         "temperature": temperature_K,
@@ -72,43 +90,48 @@ def run_lammps(modelname: str, temperature_K: float, pressure_bar: float, timest
         "pressure": pressure_bar,
         "pressure_damping": pdamp,
         "timestep": timestep_ps,
-        "number_sampling_timesteps": number_sampling_timesteps,
+        "thermo_sampling_period": thermo_sampling_period,
         "species": " ".join(species),
-        "average_position_filename": "output/average_position.dump.*",
-        "average_cell_filename": "output/average_cell.dump",
+        "zero_temperature_crystal_filename": "zero_temperature_crystal.lmp",
+        "average_position_filename": f"{average_position_filename}.*",
+        "average_cell_filename": average_cell_filename,
         "write_restart_filename": restart_filename,
-        "trajectory_filename": "output/trajectory.lammpstrj",
-        "msd_trajectory_filename": "output/msd_trajectory.lammpstrj",
+        "trajectory_filename": "trajectory.lammpstrj",
+        "msd_trajectory_filename": "msd_trajectory.lammpstrj",
         "msd_threshold": msd_threshold_angstrom_squared_per_sampling_timesteps,
         "msd_timesteps": number_msd_timesteps,
-	    "melted_crystal_output": "output/melted_crystal.dump"
+        "rlc_run_length": rlc_run_length,
+        "rlc_n_every": rlc_n_every,
+        "melted_crystal_output": melted_crystal_filename
     }
 
     command = (
             f"{lammps_command} "
             + " ".join(f"-var {key} '{item}'" for key, item in variables.items())
             + f" -log {log_filename}"
-            + " -in npt.lammps")
+            + f" -in npt.lammps")
 
-    subprocess.run(command, check=True, shell=True)
+    subprocess.run(command, check=True, shell=True, cwd=output_dir)
 
-    plot_property_from_lammps_log(log_filename, ("v_vol_metal", "v_temp_metal", "v_enthalpy_metal"))
+    if equilibration_plots:
+        plot_property_from_lammps_log(f"{output_dir}/{log_filename}",
+                                      ("v_vol_metal", "v_temp_metal", "v_enthalpy_metal"))
 
-    # 10000 offset from MSD detection during which kim_convergence was not used.
-    equilibration_time = extract_equilibration_step_from_logfile(log_filename) + 10000
-    # Round to next multiple of 10000.
-    equilibration_time = int(ceil(equilibration_time / 10000.0)) * 10000
+    equilibration_time = extract_equilibration_step_from_logfile(f"{output_dir}/{log_filename}")
+    # Round to next multiple of rlc_run_length.
+    equilibration_time = int(ceil(equilibration_time / float(rlc_run_length))) * rlc_run_length
 
-    full_average_position_file = "output/average_position.dump.full"
-    compute_average_positions_from_lammps_dump("output",
-                                               "average_position.dump",
+    full_average_position_file = f"{output_dir}/{average_position_filename}.full"
+    compute_average_positions_from_lammps_dump(output_dir,
+                                               average_position_filename,
                                                full_average_position_file, equilibration_time)
 
-    full_average_cell_file = "output/average_cell.dump.full"
-    compute_average_cell_from_lammps_dump("output/average_cell.dump",
+    full_average_cell_file = f"{output_dir}/{average_cell_filename}.full"
+    compute_average_cell_from_lammps_dump(f"{output_dir}/{average_cell_filename}",
                                           full_average_cell_file, equilibration_time)
 
-    return log_filename, restart_filename, full_average_position_file, full_average_cell_file
+    return (f"{output_dir}/{log_filename}", f"{output_dir}/{restart_filename}", full_average_position_file,
+            full_average_cell_file, f"{output_dir}/{melted_crystal_filename}")
 
 
 def plot_property_from_lammps_log(in_file_path: str, property_names: Iterable[str]) -> None:

test_driver/helper_functions.py

@@ -8,7 +8,7 @@ boundary p p p
 neigh_modify delay 0 every 1 check yes one 4000
 
 # Read initial zero-temperature crystal.
-read_data output/zero_temperature_crystal.lmp
+read_data ${zero_temperature_crystal_filename}
 
 # Change to triclinic box.
 change_box all triclinic
@@ -66,23 +66,23 @@ run 0
 velocity all scale $(v_temperature*v__u_temperature)
 
 # Write trajectory.
-dump msd_traj all atom 10000 ${msd_trajectory_filename}
+dump msd_traj all atom ${rlc_run_length} ${msd_trajectory_filename}
 
 # Compute mean squared displacement 
 compute msd all msd com yes
 
 # Perform a short run to equilibrate MSD
 run ${msd_timesteps}
-        
+
 # Calculate slope of MSD to detect diffusion
-fix msd_vector all vector ${number_sampling_timesteps} c_msd[4]
+fix msd_vector all vector ${thermo_sampling_period} c_msd[4]
 variable msd_slope equal slope(f_msd_vector)
 
 # Thermodynamic output
 thermo_style custom step v_pe_metal v_ke_metal v_temp_metal v_press_metal v_vol_metal &
                          v_xlo_metal v_xhi_metal v_ylo_metal v_yhi_metal v_zlo_metal v_zhi_metal &
                          v_xy_metal v_xz_metal v_yz_metal v_enthalpy_metal c_msd[4] v_msd_slope
-thermo ${number_sampling_timesteps}
+thermo ${thermo_sampling_period}
 
 # Perform a short run and decide whether or not to quit
 run ${msd_timesteps}
@@ -96,13 +96,13 @@ undump msd_traj
 reset_timestep 0
 
 # Write trajectory.
-dump traj all atom 10000 ${trajectory_filename}
+dump traj all atom ${rlc_run_length} ${trajectory_filename}
 
 # Set up logging of thermodynamic information.
 thermo_style custom step v_pe_metal v_ke_metal v_temp_metal v_press_metal v_vol_metal &
                          v_xlo_metal v_xhi_metal v_ylo_metal v_yhi_metal v_zlo_metal v_zhi_metal &
                          v_xy_metal v_xz_metal v_yz_metal v_enthalpy_metal
-thermo ${number_sampling_timesteps}
+thermo ${thermo_sampling_period}
 
 # Store unwrapped coordinates for every atom. 
 # Note that Lammps cannot store scaled unwrapped coordinates (it can only dump them).
@@ -125,34 +125,33 @@ variable ysu atom "(c_up[2]-ylo)/(yhi-ylo) - (yz*(c_up[3]-zlo))/((yhi-ylo)*(zhi-
 variable zsu atom "(c_up[3]-zlo)/(zhi-zlo)"
 
 # Since Lammps can dump scaled unwrapped coordinates, one can test above variables by including the following two lines:
-# dump test all custom 10000 output/test.dump.* id xsu ysu zsu v_xsu v_ysu v_zsu
-# dump_modify test delay 10000
+# dump test all custom ${rlc_run_length} output/test.dump.* id xsu ysu zsu v_xsu v_ysu v_zsu
+# dump_modify test delay ${rlc_run_length}
 
 # Average the scaled unwrapped positions.
-# Note that kim-convergence interrupts simulations after 10000 timesteps.
-# We write out the average scaled unwrapped positions after every 10000 timesteps.
-fix avePos all ave/atom 1 10000 10000 v_xsu v_ysu v_zsu
-dump avePosDump all custom 10000 ${average_position_filename} id f_avePos[1] f_avePos[2] f_avePos[3]
+# Note that kim-convergence interrupts simulations after ${rlc_run_length} timesteps.
+# We write out the average scaled unwrapped positions after every ${rlc_run_length} timesteps.
+fix avePos all ave/atom 1 ${rlc_run_length} ${rlc_run_length} v_xsu v_ysu v_zsu
+dump avePosDump all custom ${rlc_run_length} ${average_position_filename} id f_avePos[1] f_avePos[2] f_avePos[3]
 
 # Prevent dump at timestep 0.
-dump_modify avePosDump delay 10000
+dump_modify avePosDump delay ${rlc_run_length}
 
 # Average the cell information.
 variable lx_metal equal "v_xhi_metal - v_xlo_metal"
 variable ly_metal equal "v_yhi_metal - v_ylo_metal"
 variable lz_metal equal "v_zhi_metal - v_zlo_metal"
-fix aveCell all ave/time 1 10000 10000 v_lx_metal v_ly_metal v_lz_metal v_xy_metal v_xz_metal v_yz_metal file ${average_cell_filename}
+fix aveCell all ave/time 1 ${rlc_run_length} ${rlc_run_length} v_lx_metal v_ly_metal v_lz_metal v_xy_metal v_xz_metal v_yz_metal file ${average_cell_filename}
 
 # Set up convergence check with kim-convergence.
-# python run_length_control input 8 SELF 1 variable vol_metal variable temp_metal variable enthalpy_metal format pissssss file run_length_control.py
-python run_length_control input 20 SELF 1 variable vol_metal variable temp_metal variable enthalpy_metal variable lx_metal variable ly_metal variable lz_metal variable xy_metal variable xz_metal variable yz_metal format pissssssssssssssssss file run_length_control.py
+python run_length_control input 20 SELF ${rlc_n_every} variable vol_metal variable temp_metal variable enthalpy_metal variable lx_metal variable ly_metal variable lz_metal variable xy_metal variable xz_metal variable yz_metal format pissssssssssssssssss file run_length_control.py
 
 # Run until converged.
 python run_length_control invoke
 
-# Run 10000 steps more in the converged regime. 
+# Run ${rlc_run_length} steps more in the converged regime.
 # This makes sure that we have average fractional coordinates and cell information that were obtained entirely in equilibrium.
-run 10000
+run ${rlc_run_length}
 
 # Reset.
 unfix NPT