From cebcb9bbc00ff743ec0eea86f24e2c92632938cf Mon Sep 17 00:00:00 2001 From: sapphire10-k Date: Tue, 23 Dec 2025 10:02:23 +0000 Subject: [PATCH 01/15] Update .gitignore to comment out .log file Need to check output in log files --- .gitignore | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index c291fe1..f819342 100644 --- a/.gitignore +++ b/.gitignore @@ -20,5 +20,5 @@ build/ *.traj *.xyz *.cif -*.log +#*.log *.out From 8aab7b5bee3971c541188ce9b5da289722babc35 Mon Sep 17 00:00:00 2001 From: sapphire10-k Date: Tue, 23 Dec 2025 10:03:13 +0000 Subject: [PATCH 02/15] Update .gitignore to comment out .traj Removing .traj from the list to check outputs of script --- .gitignore | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index f819342..65ebff9 100644 --- a/.gitignore +++ b/.gitignore @@ -17,7 +17,7 @@ build/ .ipynb_checkpoints/ # Common simulation outputs (adjust as needed) -*.traj +#*.traj *.xyz *.cif #*.log From c6bf2eb703891f0be63e06a160681f9392ff07fc Mon Sep 17 00:00:00 2001 From: sapphire10-k Date: Tue, 23 Dec 2025 10:15:59 +0000 Subject: [PATCH 03/15] Create iro2_slab_setup.py Split original script into 2 parts. This one sets up the slab, indexes top layer O atoms and adds hydrogen atom --- scripts/iro2_slab_setup.py | 133 +++++++++++++++++++++++++++++++++++++ 1 file changed, 133 insertions(+) create mode 100644 scripts/iro2_slab_setup.py diff --git a/scripts/iro2_slab_setup.py b/scripts/iro2_slab_setup.py new file mode 100644 index 0000000..e0ff263 --- /dev/null +++ b/scripts/iro2_slab_setup.py @@ -0,0 +1,133 @@ +from pathlib import Path +from ase.io import read, write +from ase import Atom +from ase.data import covalent_radii +from ase.constraints import FixAtoms +from ase.neighborlist import neighbor_list +import numpy as np +import json + +def setup_structure( + input_file: str = "inputs/slab_clean_2x2.in", + o_index: int = 20, + oh_distance: float = 1.0, + z_freeze: float = 20.0, + neighbor_cutoff: float = 1.5, +): + # Create output directory + Path("structures").mkdir(exist_ok=True) + + # Adjust H covalent radius + covalent_radii[1] = 0.6 + + print(f"Reading {input_file}...") + slab = read(slab_clean_2x2.in, format="espresso-in") + + # Disable PBC for xTB + slab.set_pbc((False, False, False)) + + print(f"Initial structure: {len(slab)} atoms") + print(f"Cell: {slab.cell.cellpar()}") + + # Find surface oxygen atoms + positions = slab.get_positions() + symbols = slab.get_chemical_symbols() + z_top = positions[:, 2].max() + + surface_O = [ + i for i, s in enumerate(symbols) + if s == "O" and (z_top - positions[i, 2]) < 1.5 + ] + + print(f"\nSurface O atoms (within 1.5 Å of top): {surface_O}") + + if o_index not in surface_O: + print(f" Warning: O[{o_index}] is not in surface list!") + print(f" O[{o_index}] is at z={positions[o_index, 2]:.2f} Å") + print(f" Top of slab: z={z_top:.2f} Å") + + # Add hydrogen atom + O_pos = slab.positions[o_index] + H_pos = O_pos + np.array([0.0, 0.0, oh_distance]) + slab.append(Atom("H", position=H_pos)) + h_index = len(slab) - 1 + + print(f"\nAdded H atom:") + print(f" Index: {h_index}") + print(f" Position: {H_pos}") + print(f" Total atoms: {len(slab)}") + + # Freeze bottom layers + z = slab.positions[:, 2] + freeze_mask = z < z_freeze + slab.set_constraint(FixAtoms(mask=freeze_mask)) + n_frozen = int(freeze_mask.sum()) + n_mobile = len(slab) - n_frozen + + print(f"\nFreezing atoms below z={z_freeze:.1f} Å:") + print(f" Frozen: {n_frozen} atoms") + print(f" Mobile: {n_mobile} atoms") + + # Check neighbors around H + i, j, d = neighbor_list("ijd", slab, cutoff=[neighbor_cutoff] * len(slab)) + neighbors = [(int(jj), float(dd)) for ii, jj, dd in zip(i, j, d) if ii == h_index] + + print(f"\nNeighbors within {neighbor_cutoff} Å of H:") + for jj, dd in neighbors: + sym = slab[jj].symbol + pos = slab[jj].position + print(f" {sym}[{jj}]: {dd:.3f} Å at ({pos[0]:.2f}, {pos[1]:.2f}, {pos[2]:.2f})") + + # Save prepared structure + output_file = f"structures/slab_H_o{o_index}_ready.xyz" + write(output_file, slab) + + # Save metadata + metadata = { + "input_file": input_file, + "o_index": o_index, + "h_index": h_index, + "oh_distance": oh_distance, + "z_freeze": z_freeze, + "n_atoms": len(slab), + "n_frozen": n_frozen, + "n_mobile": n_mobile, + "neighbors": neighbors, + "h_position": H_pos.tolist(), + "o_position": O_pos.tolist(), + } + + meta_file = f"structures/metadata_o{o_index}.json" + with open(meta_file, 'w') as f: + json.dump(metadata, f, indent=2) + + print(f"\n{'='*60}") + print(f" Structure prepared successfully!") + print(f"{'='*60}") + print(f"Output files:") + print(f" Structure: {output_file}") + print(f" Metadata: {meta_file}") + print(f"\nNext step:") + print(f" python run_optimization.py {output_file}") + print(f"{'='*60}") + + return slab, metadata + + +if __name__ == "__main__": + import argparse + + parser = argparse.ArgumentParser(description="Prepare H* adsorption structure") + parser.add_argument("--input", default="inputs/slab_clean_2x2.in", help="Input structure file") + parser.add_argument("--o-index", type=int, default=20, help="Oxygen atom index for H placement") + parser.add_argument("--oh-dist", type=float, default=1.0, help="O-H distance (Å)") + parser.add_argument("--z-freeze", type=float, default=20.0, help="Freeze atoms below this z (Å)") + + args = parser.parse_args() + + slab, meta = setup_structure( + input_file=args.input, + o_index=args.o_index, + oh_distance=args.oh_dist, + z_freeze=args.z_freeze, + ) From 15da3b4ecbe9c54881bd5807feb9296b31670b04 Mon Sep 17 00:00:00 2001 From: sapphire10-k Date: Tue, 23 Dec 2025 10:38:07 +0000 Subject: [PATCH 04/15] Update iro2_slab_setup.py Fixed filename error --- scripts/iro2_slab_setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/iro2_slab_setup.py b/scripts/iro2_slab_setup.py index e0ff263..4045d32 100644 --- a/scripts/iro2_slab_setup.py +++ b/scripts/iro2_slab_setup.py @@ -21,7 +21,7 @@ def setup_structure( covalent_radii[1] = 0.6 print(f"Reading {input_file}...") - slab = read(slab_clean_2x2.in, format="espresso-in") + slab = read(slab_clean_2x2, format="espresso-in") # Disable PBC for xTB slab.set_pbc((False, False, False)) From 7906c1db3ddb62601090a91bbdb0836f69d9122a Mon Sep 17 00:00:00 2001 From: sapphire10-k Date: Tue, 23 Dec 2025 10:41:15 +0000 Subject: [PATCH 05/15] Create optimization_iro2_test.py Coped script from main branch --- scripts/optimization_iro2_test.py | 189 ++++++++++++++++++++++++++++++ 1 file changed, 189 insertions(+) create mode 100644 scripts/optimization_iro2_test.py diff --git a/scripts/optimization_iro2_test.py b/scripts/optimization_iro2_test.py new file mode 100644 index 0000000..f223fe1 --- /dev/null +++ b/scripts/optimization_iro2_test.py @@ -0,0 +1,189 @@ +from pathlib import Path +from ase.io import read, write +from xtb.ase.calculator import XTB +from ase.optimize import BFGS +import numpy as np +import json +import time + +def run_optimization( + structure_file: str slab_clean_2x2.in, + fmax: float = 0.05, + max_steps: int = 200, + method: str = "GFN2-xTB", + output_dir: str = "results", +): + # Create output directory + Path(output_dir).mkdir(exist_ok=True) + + # Start timing + start_time = time.time() + + print(f"{'='*70}") + print(f"XTB GEOMETRY OPTIMIZATION") + print(f"{'='*70}") + print(f"Structure: {structure_file}") + print(f"Method: {method}") + print(f"Convergence: fmax < {fmax} eV/Å") + print(f"Max steps: {max_steps}") + print(f"{'='*70}\n") + + # Load structure + print("Loading structure...") + slab = read(structure_file) + + # Check constraints + if slab.constraints: + n_constrained = sum(c.get_indices().size for c in slab.constraints) + print(f" Constraints loaded: {n_constrained} frozen atoms") + else: + print("WARNING: No constraints found!") + + print(f"Total atoms: {len(slab)}") + print(f"Mobile atoms: {len(slab) - n_constrained}") + + # Setup xTB calculator + print(f"\nSetting up {method} calculator...") + slab.calc = XTB(method=method) + + # Initial energy and forces + print("\nCalculating initial state...") + e_initial = slab.get_potential_energy() + forces = slab.get_forces() + fmax_initial = np.sqrt((forces**2).sum(axis=1)).max() + + print(f" Initial energy: {e_initial:.6f} eV") + print(f" Initial fmax: {fmax_initial:.6f} eV/Å") + + if fmax_initial < fmax: + print(f"\n Already converged! No optimization needed.") + final_file = f"{output_dir}/{Path(structure_file).stem}_final.xyz" + write(final_file, slab) + return slab, { + "converged": True, + "initial": True, + "e_initial": e_initial, + "e_final": e_initial, + "fmax_initial": fmax_initial, + "fmax_final": fmax_initial, + "n_steps": 0, + "time_seconds": time.time() - start_time, + } + + # Setup output files + base_name = Path(structure_file).stem + traj_file = f"{output_dir}/{base_name}_opt.traj" + log_file = f"{output_dir}/{base_name}_opt.log" + + print(f"\n{'='*70}") + print(f"Starting optimization...") + print(f"{'='*70}") + print(f"Trajectory: {traj_file}") + print(f"Log file: {log_file}") + print(f"{'='*70}\n") + + # Run optimization + opt = BFGS(slab, trajectory=traj_file, logfile=log_file) + + try: + opt.run(fmax=fmax, steps=max_steps) + converged = True + except Exception as e: + print(f"\n Optimization stopped: {e}") + converged = False + + # Get final results + e_final = slab.get_potential_energy() + forces_final = slab.get_forces() + fmax_final = np.sqrt((forces_final**2).sum(axis=1)).max() + n_steps = opt.get_number_of_steps() + elapsed_time = time.time() - start_time + + # Print results + print(f"\n{'='*70}") + print(f"OPTIMIZATION COMPLETE") + print(f"{'='*70}") + print(f"Steps taken: {n_steps}") + print(f"Time elapsed: {elapsed_time:.1f} seconds ({elapsed_time/60:.1f} min)") + print(f"Initial energy: {e_initial:.6f} eV") + print(f"Final energy: {e_final:.6f} eV") + print(f"Energy change: {e_final - e_initial:.6f} eV") + print(f"Initial fmax: {fmax_initial:.6f} eV/Å") + print(f"Final fmax: {fmax_final:.6f} eV/Å") + print(f"Converged: {'YES' if fmax_final < fmax else 'NO'}") + print(f"{'='*70}") + + # Save final structure + final_file = f"{output_dir}/{base_name}_final.xyz" + write(final_file, slab) + print(f"\nFinal structure saved to: {final_file}") + + # Save results summary + results = { + "structure_file": structure_file, + "method": method, + "fmax_target": fmax, + "max_steps": max_steps, + "n_steps": n_steps, + "converged": fmax_final < fmax, + "e_initial": float(e_initial), + "e_final": float(e_final), + "delta_e": float(e_final - e_initial), + "fmax_initial": float(fmax_initial), + "fmax_final": float(fmax_final), + "time_seconds": elapsed_time, + "time_minutes": elapsed_time / 60, + } + + results_file = f"{output_dir}/{base_name}_results.json" + with open(results_file, 'w') as f: + json.dump(results, f, indent=2) + + print(f"Results summary saved to: {results_file}") + print(f"{'='*70}\n") + + return slab, results + + +#Quick test with limited steps + +def quick_test(structure_file: str, steps: int = 5): + + print(f"\n{'='*70}") + print(f"QUICK TEST MODE - {steps} STEPS ONLY") + print(f"{'='*70}\n") + + return run_optimization( + structure_file, + fmax=0.05, + max_steps=steps, + output_dir="test_results" + ) + + +if __name__ == "__main__": + import argparse + + parser = argparse.ArgumentParser(description="Run xTB optimization") + parser.add_argument("structure", help="Prepared structure file (.xyz)") + parser.add_argument("--fmax", type=float, default=0.05, help="Force convergence (eV/Å)") + parser.add_argument("--steps", type=int, default=200, help="Maximum steps") + parser.add_argument("--method", default="GFN2-xTB", choices=["GFN1-xTB", "GFN2-xTB"], + help="xTB method (GFN1 is faster)") + parser.add_argument("--test", action="store_true", help="Quick test mode (5 steps)") + parser.add_argument("--output", default="results", help="Output directory") + + args = parser.parse_args() + + if args.test: + slab, results = quick_test(args.structure, steps=5) + else: + slab, results = run_optimization( + args.structure, + fmax=args.fmax, + max_steps=args.steps, + method=args.method, + output_dir=args.output, + ) + + print("\n Done!") From 5421e46775b85b57ab2916d6fe6c69cf24797fb1 Mon Sep 17 00:00:00 2001 From: sapphire10-k Date: Tue, 23 Dec 2025 10:42:51 +0000 Subject: [PATCH 06/15] Update iro2_slab_setup.py Corrected filename repetition --- scripts/iro2_slab_setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/iro2_slab_setup.py b/scripts/iro2_slab_setup.py index 4045d32..e0fa275 100644 --- a/scripts/iro2_slab_setup.py +++ b/scripts/iro2_slab_setup.py @@ -21,7 +21,7 @@ def setup_structure( covalent_radii[1] = 0.6 print(f"Reading {input_file}...") - slab = read(slab_clean_2x2, format="espresso-in") + slab = read(input_file, format="espresso-in") # Disable PBC for xTB slab.set_pbc((False, False, False)) From 7a089ed2bcfce772e117209791169d06a0fda85b Mon Sep 17 00:00:00 2001 From: sapphire10-k Date: Tue, 23 Dec 2025 10:47:03 +0000 Subject: [PATCH 07/15] Update optimization_iro2_test.py Testing structure_file declaration issue --- scripts/optimization_iro2_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/optimization_iro2_test.py b/scripts/optimization_iro2_test.py index f223fe1..9aa5957 100644 --- a/scripts/optimization_iro2_test.py +++ b/scripts/optimization_iro2_test.py @@ -7,7 +7,7 @@ import time def run_optimization( - structure_file: str slab_clean_2x2.in, + structure_file: str, fmax: float = 0.05, max_steps: int = 200, method: str = "GFN2-xTB", From c5315fdc887f7733f816babd0cbaf8c9145395a9 Mon Sep 17 00:00:00 2001 From: sapphire10-k Date: Tue, 23 Dec 2025 10:53:11 +0000 Subject: [PATCH 08/15] Create optimization_iro2_3.py Created a script for faster optimization using xTB --- scripts/optimization_iro2_3.py | 106 +++++++++++++++++++++++++++++++++ 1 file changed, 106 insertions(+) create mode 100644 scripts/optimization_iro2_3.py diff --git a/scripts/optimization_iro2_3.py b/scripts/optimization_iro2_3.py new file mode 100644 index 0000000..20fa13e --- /dev/null +++ b/scripts/optimization_iro2_3.py @@ -0,0 +1,106 @@ +import json +from pathlib import Path +from ase.io import read +import numpy as np + +def analyze_result(result_dir: str = "results"): + #Analyze optimization results and print summary + result_path = Path(result_dir) + + # Find all result JSON files + result_files = list(result_path.glob("*_results.json")) + + if not result_files: + print(f"No results found in {result_dir}/") + return + + print(f"\n{'='*70}") + print(f"OPTIMIZATION RESULTS SUMMARY") + print(f"{'='*70}\n") + + for result_file in sorted(result_files): + with open(result_file) as f: + data = json.load(f) + + print(f"File: {result_file.name}") + print(f" Structure: {Path(data['structure_file']).name}") + print(f" Method: {data['method']}") + print(f" Steps: {data['n_steps']}/{data['max_steps']}") + print(f" Time: {data['time_minutes']:.1f} min") + print(f" Energy: {data['e_initial']:.4f} → {data['e_final']:.4f} eV") + print(f" ΔE: {data['delta_e']:.4f} eV") + print(f" fmax: {data['fmax_initial']:.4f} → {data['fmax_final']:.4f} eV/Å") + print(f" Converged: {'YES' if data['converged'] else 'NO'}") + print() + + +def compare_structures(initial_xyz: str, final_xyz: str): + """ + Compare initial and final structures + """ + initial = read(initial_xyz) + final = read(final_xyz) + + # Find H atom (last atom) + h_idx = len(initial) - 1 + + h_init = initial[h_idx].position + h_final = final[h_idx].position + h_displacement = np.linalg.norm(h_final - h_init) + + print(f"\n{'='*70}") + print(f"STRUCTURE COMPARISON") + print(f"{'='*70}") + print(f"Initial: {initial_xyz}") + print(f"Final: {final_xyz}") + print(f"\nH atom displacement:") + print(f" Initial position: ({h_init[0]:.3f}, {h_init[1]:.3f}, {h_init[2]:.3f})") + print(f" Final position: ({h_final[0]:.3f}, {h_final[1]:.3f}, {h_final[2]:.3f})") + print(f" Total movement: {h_displacement:.3f} Å") + print(f"{'='*70}\n") + + # Calculate RMSD for all atoms + rmsd = np.sqrt(((initial.positions - final.positions)**2).sum(axis=1).mean()) + print(f"Overall RMSD: {rmsd:.4f} Å") + + +def visualize_trajectory(traj_file: str): + """ + Quick stats from trajectory + """ + from ase.io import read as read_traj + + traj = read_traj(traj_file, index=":") + + energies = [atoms.get_potential_energy() for atoms in traj] + + print(f"\n{'='*70}") + print(f"TRAJECTORY ANALYSIS") + print(f"{'='*70}") + print(f"File: {traj_file}") + print(f"Steps: {len(traj)}") + print(f"Energy range: {min(energies):.4f} to {max(energies):.4f} eV") + print(f"Total ΔE: {energies[-1] - energies[0]:.4f} eV") + print(f"{'='*70}\n") + + +if __name__ == "__main__": + import argparse + + parser = argparse.ArgumentParser(description="Analyze optimization results") + parser.add_argument("--summary", action="store_true", help="Show all results summary") + parser.add_argument("--compare", nargs=2, metavar=("INITIAL", "FINAL"), + help="Compare two structures") + parser.add_argument("--traj", help="Analyze trajectory file") + parser.add_argument("--dir", default="results", help="Results directory") + + args = parser.parse_args() + + if args.summary or (not args.compare and not args.traj): + analyze_result(args.dir) + + if args.compare: + compare_structures(args.compare[0], args.compare[1]) + + if args.traj: + visualize_trajectory(args.traj) From 2a638c8a067db3fa549b3aaa21563b2a67281723 Mon Sep 17 00:00:00 2001 From: sapphire10-k Date: Tue, 23 Dec 2025 10:54:23 +0000 Subject: [PATCH 09/15] Update iro2_test_1.py Commenting out all lines to test split scripts --- iro2_test_1.py | 206 ++++++++++++++++++++++++------------------------- 1 file changed, 103 insertions(+), 103 deletions(-) diff --git a/iro2_test_1.py b/iro2_test_1.py index 7fa1b1b..d2f73e7 100644 --- a/iro2_test_1.py +++ b/iro2_test_1.py @@ -1,119 +1,119 @@ -# pip install ase -# conda config --add channels conda-forge -# conda install xtb-python -# conda search xtb-python --channel conda-forge -# conda install qcelemental ase -from pathlib import Path -from ase.io import read -import numpy as np -from ase.constraints import FixAtoms -from ase import Atom -from ase.visualize import view -from ase.io import write -from ase.visualize import view -from ase.data import covalent_radii -from ase.constraints import FixAtoms -from ase.neighborlist import neighbor_list -from xtb.ase.calculator import XTB -from ase.optimize import BFGS - -Path("outputs").mkdir(exist_ok=True) - - -def run_relaxation( - qe_input: str = "slab_clean_2x2.in", - traj: str = "relax_H.traj", - log: str = "relax_H.log", - fmax: float = 0.05, - cut: float = 1.3, -): - - covalent_radii[1] = 0.6 - - slab = read("slab_clean_2x2.in", format="espresso-in") - # xTB does not support PBC multipoles in this mode → treat slab as cluster - slab.set_pbc((False, False, False)) - print(slab) - print("Atoms:", len(slab)) - print("Cell (Å):\n", slab.cell) - print("PBC:", slab.pbc) - - # Get atomic symbols and positions - symbols = slab.get_chemical_symbols() - positions = slab.get_positions() - - # Find top of slab - z_positions = positions[:, 2] - z_top = z_positions.max() - - print(f"Top of slab z = {z_top:.3f} Å") - - # Find surface O atoms within 1.5 Å of top - surface_O_indices = [ - i for i, s in enumerate(symbols) - if s == "O" and (z_top - positions[i, 2]) < 1.5 - ] - - print("Candidate surface O atoms:") - for i in surface_O_indices: - x, y, z = positions[i] - print(f"Index {i}: position = ({x:.3f}, {y:.3f}, {z:.3f})") - - #Adding H* on top of Oxygen 20 - - o_idx = 20 - O_pos = slab.positions[o_idx] - - # Typical O–H bond length - OH_distance = 1.0 # Å - - # Place H above O along +z - H_pos = O_pos + np.array([0.0, 0.0, OH_distance]) - - # Add H atom - slab.append(Atom("H", position=H_pos)) - - #for atom in slab: - #if atom.symbol == "H": - #atom.tag = 1 +# # pip install ase +# # conda config --add channels conda-forge +# # conda install xtb-python +# # conda search xtb-python --channel conda-forge +# # conda install qcelemental ase +# #from pathlib import Path +# from ase.io import read +# import numpy as np +# from ase.constraints import FixAtoms +# from ase import Atom +# from ase.visualize import view +# from ase.io import write +# from ase.visualize import view +# from ase.data import covalent_radii +# from ase.constraints import FixAtoms +# from ase.neighborlist import neighbor_list +# from xtb.ase.calculator import XTB +# from ase.optimize import BFGS + +# Path("outputs").mkdir(exist_ok=True) + + +# def run_relaxation( +# qe_input: str = "slab_clean_2x2.in", +# traj: str = "relax_H.traj", +# log: str = "relax_H.log", +# fmax: float = 0.05, +# cut: float = 1.3, +# ): + +# covalent_radii[1] = 0.6 + +# slab = read("slab_clean_2x2.in", format="espresso-in") +# # xTB does not support PBC multipoles in this mode → treat slab as cluster +# slab.set_pbc((False, False, False)) +# print(slab) +# print("Atoms:", len(slab)) +# print("Cell (Å):\n", slab.cell) +# print("PBC:", slab.pbc) + +# # Get atomic symbols and positions +# symbols = slab.get_chemical_symbols() +# positions = slab.get_positions() + +# # Find top of slab +# z_positions = positions[:, 2] +# z_top = z_positions.max() + +# print(f"Top of slab z = {z_top:.3f} Å") + +# # Find surface O atoms within 1.5 Å of top +# surface_O_indices = [ +# i for i, s in enumerate(symbols) +# if s == "O" and (z_top - positions[i, 2]) < 1.5 +# ] + +# print("Candidate surface O atoms:") +# for i in surface_O_indices: +# x, y, z = positions[i] +# print(f"Index {i}: position = ({x:.3f}, {y:.3f}, {z:.3f})") + +# #Adding H* on top of Oxygen 20 + +# o_idx = 20 +# O_pos = slab.positions[o_idx] + +# # Typical O–H bond length +# OH_distance = 1.0 # Å + +# # Place H above O along +z +# H_pos = O_pos + np.array([0.0, 0.0, OH_distance]) + +# # Add H atom +# slab.append(Atom("H", position=H_pos)) + +# #for atom in slab: +# #if atom.symbol == "H": +# #atom.tag = 1 - #set_highlight_atoms([i for i,a in enumerate(slab) if a.symbol=="H"]) +# #set_highlight_atoms([i for i,a in enumerate(slab) if a.symbol=="H"]) - print("Total atoms:", len(slab)) - print("Last atom symbol:", slab[-1].symbol) - print("Last atom position:", slab[-1].position) +# print("Total atoms:", len(slab)) +# print("Last atom symbol:", slab[-1].symbol) +# print("Last atom position:", slab[-1].position) - print(f"Added H at position {H_pos}") - print("Total atoms now:", len(slab)) +# print(f"Added H at position {H_pos}") +# print("Total atoms now:", len(slab)) - write("slab_with_H.xyz", slab) +# write("slab_with_H.xyz", slab) - view(slab) +# view(slab) - z = slab.positions[:, 2] - z_freeze = 20.0 # Å (tune if needed) +# z = slab.positions[:, 2] +# z_freeze = 20.0 # Å (tune if needed) - freeze_mask = z < z_freeze - slab.set_constraint(FixAtoms(mask=freeze_mask)) +# freeze_mask = z < z_freeze +# slab.set_constraint(FixAtoms(mask=freeze_mask)) - print("Frozen atoms:", int(freeze_mask.sum()), "/", len(slab)) +# print("Frozen atoms:", int(freeze_mask.sum()), "/", len(slab)) - iH = len(slab) - 1 # last atom is H - cut = 1.3 # Å check radius around H +# iH = len(slab) - 1 # last atom is H +# cut = 1.3 # Å check radius around H - i, j, d = neighbor_list("ijd", slab, cutoff=[cut]*len(slab)) - close = [(jj, dd) for ii, jj, dd in zip(i, j, d) if ii == iH] +# i, j, d = neighbor_list("ijd", slab, cutoff=[cut]*len(slab)) +# close = [(jj, dd) for ii, jj, dd in zip(i, j, d) if ii == iH] - print("Neighbors within 1.3 Å of H:", close) +# print("Neighbors within 1.3 Å of H:", close) - # Optimize structure with xTB - slab.calc = XTB(method="GFN2-xTB") - opt = BFGS(slab, trajectory=traj, logfile=log) - opt.run(fmax=fmax) +# # Optimize structure with xTB +# slab.calc = XTB(method="GFN2-xTB") +# opt = BFGS(slab, trajectory=traj, logfile=log) +# opt.run(fmax=fmax) - return slab, close +# return slab, close -if __name__ == "__main__": - run_relaxation() +# if __name__ == "__main__": +# run_relaxation() From d339562cd3eec452952ebad9cf89bbb8e5afd60c Mon Sep 17 00:00:00 2001 From: sapphire10-k Date: Tue, 23 Dec 2025 11:05:43 +0000 Subject: [PATCH 10/15] Update python-package-conda.yml Added miniconda install to check outputs --- .github/workflows/python-package-conda.yml | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/.github/workflows/python-package-conda.yml b/.github/workflows/python-package-conda.yml index 63d6113..af4fd3f 100644 --- a/.github/workflows/python-package-conda.yml +++ b/.github/workflows/python-package-conda.yml @@ -32,3 +32,11 @@ jobs: run: | conda install pytest pytest -q -m "not slow" + - name: Set up Miniconda + uses: conda-incubator/setup-miniconda@v3 + with: + auto-update-conda: true + python-version: "3.10" + environment-file: environment.yml + activate-environment: base + auto-activate-base: true From a74d00bf149b79abb80ceb7f8a3aec13332d9245 Mon Sep 17 00:00:00 2001 From: sapphire10-k Date: Tue, 23 Dec 2025 11:11:14 +0000 Subject: [PATCH 11/15] Update python-package-conda.yml Added steps for uploading simulation outputs as artifacts --- .github/workflows/python-package-conda.yml | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/.github/workflows/python-package-conda.yml b/.github/workflows/python-package-conda.yml index af4fd3f..043de35 100644 --- a/.github/workflows/python-package-conda.yml +++ b/.github/workflows/python-package-conda.yml @@ -33,10 +33,20 @@ jobs: conda install pytest pytest -q -m "not slow" - name: Set up Miniconda - uses: conda-incubator/setup-miniconda@v3 + - uses: conda-incubator/setup-miniconda@v3 with: auto-update-conda: true python-version: "3.10" environment-file: environment.yml activate-environment: base auto-activate-base: true + - name: Run simulation + run: | + mkdir -p outputs + conda run -n base python path/to/your_driver_script.py + ls -R outputs + - name: Upload simulation outputs + - uses: actions/upload-artifact@v4 + with: + name: irO2-simulation-outputs + path: outputs/ From 6abe26bf5de0caddc7c1d56653eb7fbbf9e5587b Mon Sep 17 00:00:00 2001 From: sapphire10-k Date: Tue, 23 Dec 2025 11:13:07 +0000 Subject: [PATCH 12/15] Update python-package-conda.yml Fixed indentation errors --- .github/workflows/python-package-conda.yml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/workflows/python-package-conda.yml b/.github/workflows/python-package-conda.yml index 043de35..f8d1775 100644 --- a/.github/workflows/python-package-conda.yml +++ b/.github/workflows/python-package-conda.yml @@ -33,7 +33,7 @@ jobs: conda install pytest pytest -q -m "not slow" - name: Set up Miniconda - - uses: conda-incubator/setup-miniconda@v3 + uses: conda-incubator/setup-miniconda@v3 with: auto-update-conda: true python-version: "3.10" @@ -46,7 +46,7 @@ jobs: conda run -n base python path/to/your_driver_script.py ls -R outputs - name: Upload simulation outputs - - uses: actions/upload-artifact@v4 - with: - name: irO2-simulation-outputs - path: outputs/ + uses: actions/upload-artifact@v4 + with: + name: irO2-simulation-outputs + path: outputs/ From b0d595865a14a89d8301ef4ee44d974008e03705 Mon Sep 17 00:00:00 2001 From: sapphire10-k Date: Tue, 23 Dec 2025 14:23:49 +0000 Subject: [PATCH 13/15] Updated python-package-conda.yml Removed extraneous additions of packages --- .github/workflows/python-package-conda.yml | 42 +++++++++------------- 1 file changed, 16 insertions(+), 26 deletions(-) diff --git a/.github/workflows/python-package-conda.yml b/.github/workflows/python-package-conda.yml index f8d1775..9d197dd 100644 --- a/.github/workflows/python-package-conda.yml +++ b/.github/workflows/python-package-conda.yml @@ -5,33 +5,10 @@ on: [push] jobs: build-linux: runs-on: ubuntu-latest - strategy: - max-parallel: 5 steps: - uses: actions/checkout@v4 - - name: Set up Python 3.10 - uses: actions/setup-python@v3 - with: - python-version: '3.10' - - name: Add conda to system path - run: | - # $CONDA is an environment variable pointing to the root of the miniconda directory - echo $CONDA/bin >> $GITHUB_PATH - - name: Install dependencies - run: | - conda env update --file environment.yml --name base - - name: Lint with flake8 - run: | - conda install flake8 - # stop the build if there are Python syntax errors or undefined names - flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics - # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide - flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics - - name: Test with pytest - run: | - conda install pytest - pytest -q -m "not slow" + - name: Set up Miniconda uses: conda-incubator/setup-miniconda@v3 with: @@ -40,11 +17,24 @@ jobs: environment-file: environment.yml activate-environment: base auto-activate-base: true + + - name: Lint with flake8 + run: | + conda run -n base python -m pip install flake8 + conda run -n base flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics + conda run -n base flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics + + - name: Test with pytest + run: | + conda run -n base python -m pip install pytest + conda run -n base pytest -q -m "not slow" + - name: Run simulation run: | mkdir -p outputs - conda run -n base python path/to/your_driver_script.py - ls -R outputs + conda run -n base python iro2_test_1.py + ls -R . + - name: Upload simulation outputs uses: actions/upload-artifact@v4 with: From a7ee30df8f9fb85a86d7cb4c408545d01b52087b Mon Sep 17 00:00:00 2001 From: sapphire10-k Date: Tue, 23 Dec 2025 14:29:37 +0000 Subject: [PATCH 14/15] Update python-package-conda.yml Removed conflicting python version tag from miniconda --- .github/workflows/python-package-conda.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/python-package-conda.yml b/.github/workflows/python-package-conda.yml index 9d197dd..071b3a5 100644 --- a/.github/workflows/python-package-conda.yml +++ b/.github/workflows/python-package-conda.yml @@ -13,7 +13,6 @@ jobs: uses: conda-incubator/setup-miniconda@v3 with: auto-update-conda: true - python-version: "3.10" environment-file: environment.yml activate-environment: base auto-activate-base: true From 92ad2eb86f323c9424280476cec2808d60df050e Mon Sep 17 00:00:00 2001 From: sapphire10-k Date: Tue, 23 Dec 2025 14:55:50 +0000 Subject: [PATCH 15/15] Updated iro2_slab_setup.py Revised code due to no artifact generation --- scripts/iro2_slab_setup.py | 210 +++++++++++++++++-------------------- 1 file changed, 99 insertions(+), 111 deletions(-) diff --git a/scripts/iro2_slab_setup.py b/scripts/iro2_slab_setup.py index e0fa275..b5c204e 100644 --- a/scripts/iro2_slab_setup.py +++ b/scripts/iro2_slab_setup.py @@ -1,11 +1,16 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import json from pathlib import Path -from ase.io import read, write + +import numpy as np from ase import Atom -from ase.data import covalent_radii from ase.constraints import FixAtoms +from ase.io import read, write from ase.neighborlist import neighbor_list -import numpy as np -import json + def setup_structure( input_file: str = "inputs/slab_clean_2x2.in", @@ -13,121 +18,104 @@ def setup_structure( oh_distance: float = 1.0, z_freeze: float = 20.0, neighbor_cutoff: float = 1.5, + outputs_dir: str = "outputs", ): - # Create output directory - Path("structures").mkdir(exist_ok=True) - - # Adjust H covalent radius - covalent_radii[1] = 0.6 - - print(f"Reading {input_file}...") - slab = read(input_file, format="espresso-in") - - # Disable PBC for xTB - slab.set_pbc((False, False, False)) - - print(f"Initial structure: {len(slab)} atoms") - print(f"Cell: {slab.cell.cellpar()}") - - # Find surface oxygen atoms - positions = slab.get_positions() - symbols = slab.get_chemical_symbols() - z_top = positions[:, 2].max() - - surface_O = [ - i for i, s in enumerate(symbols) - if s == "O" and (z_top - positions[i, 2]) < 1.5 - ] - - print(f"\nSurface O atoms (within 1.5 Å of top): {surface_O}") - - if o_index not in surface_O: - print(f" Warning: O[{o_index}] is not in surface list!") - print(f" O[{o_index}] is at z={positions[o_index, 2]:.2f} Å") - print(f" Top of slab: z={z_top:.2f} Å") - - # Add hydrogen atom - O_pos = slab.positions[o_index] - H_pos = O_pos + np.array([0.0, 0.0, oh_distance]) - slab.append(Atom("H", position=H_pos)) + """ + Read slab, place H above a chosen O, freeze atoms below z_freeze, + and write outputs for the next step. + """ + input_path = Path(input_file) + if not input_path.exists(): + raise FileNotFoundError(f"Input file not found: {input_path}") + + slab = read(str(input_path)) + + n0 = len(slab) + if not (0 <= o_index < n0): + raise IndexError(f"o_index={o_index} out of range for {n0} atoms (0..{n0-1})") + + # Slab z-extent (for diagnostics) + zmax = slab.positions[:, 2].max() + + # Position of selected O + o_pos = slab.positions[o_index].copy() + + # Place H above O along +z + h_pos = o_pos + np.array([0.0, 0.0, float(oh_distance)]) + + # Sanity check: ensure H is above the slab top + if h_pos[2] <= zmax: + print( + f"[warn] H z={h_pos[2]:.3f} Å is not above slab top zmax={zmax:.3f} Å. " + "Raising H to be above slab." + ) + h_pos[2] = zmax + float(oh_distance) + + slab.append(Atom("H", position=h_pos)) + h_index = len(slab) - 1 + + slab.append(Atom("H", position=h_pos)) h_index = len(slab) - 1 - - print(f"\nAdded H atom:") - print(f" Index: {h_index}") - print(f" Position: {H_pos}") - print(f" Total atoms: {len(slab)}") - - # Freeze bottom layers - z = slab.positions[:, 2] - freeze_mask = z < z_freeze + + # Freeze atoms below z_freeze + freeze_mask = slab.positions[:, 2] < float(z_freeze) slab.set_constraint(FixAtoms(mask=freeze_mask)) - n_frozen = int(freeze_mask.sum()) - n_mobile = len(slab) - n_frozen - - print(f"\nFreezing atoms below z={z_freeze:.1f} Å:") - print(f" Frozen: {n_frozen} atoms") - print(f" Mobile: {n_mobile} atoms") - - # Check neighbors around H - i, j, d = neighbor_list("ijd", slab, cutoff=[neighbor_cutoff] * len(slab)) - neighbors = [(int(jj), float(dd)) for ii, jj, dd in zip(i, j, d) if ii == h_index] - - print(f"\nNeighbors within {neighbor_cutoff} Å of H:") - for jj, dd in neighbors: - sym = slab[jj].symbol - pos = slab[jj].position - print(f" {sym}[{jj}]: {dd:.3f} Å at ({pos[0]:.2f}, {pos[1]:.2f}, {pos[2]:.2f})") - - # Save prepared structure - output_file = f"structures/slab_H_o{o_index}_ready.xyz" - write(output_file, slab) - - # Save metadata - metadata = { - "input_file": input_file, - "o_index": o_index, - "h_index": h_index, - "oh_distance": oh_distance, - "z_freeze": z_freeze, - "n_atoms": len(slab), - "n_frozen": n_frozen, - "n_mobile": n_mobile, - "neighbors": neighbors, - "h_position": H_pos.tolist(), - "o_position": O_pos.tolist(), + + # Neighbors of H (for sanity check) + i, j, d = neighbor_list("ijd", slab, cutoff=float(neighbor_cutoff)) + neigh = [(int(jj), float(dd)) for ii, jj, dd in zip(i, j, d) if int(ii) == h_index] + neigh_sorted = sorted(neigh, key=lambda x: x[1]) + + outdir = Path(outputs_dir) + (outdir / "results").mkdir(parents=True, exist_ok=True) + + # Write structures + write(str(outdir / "slab_with_H.traj"), slab) + write(str(outdir / "slab_with_H.xyz"), slab) + + meta = { + "input_file": str(input_path), + "num_atoms_initial": int(n0), + "num_atoms_total": int(len(slab)), + "o_index": int(o_index), + "h_index": int(h_index), + "oh_distance_A": float(oh_distance), + "z_freeze_A": float(z_freeze), + "neighbor_cutoff_A": float(neighbor_cutoff), + "o_position_A": [float(x) for x in o_pos], + "h_position_A": [float(x) for x in h_pos], + "num_frozen": int(np.count_nonzero(freeze_mask)), + "neighbors_of_H": [{"index": idx, "distance_A": dist} for idx, dist in neigh_sorted], } - - meta_file = f"structures/metadata_o{o_index}.json" - with open(meta_file, 'w') as f: - json.dump(metadata, f, indent=2) - - print(f"\n{'='*60}") - print(f" Structure prepared successfully!") - print(f"{'='*60}") - print(f"Output files:") - print(f" Structure: {output_file}") - print(f" Metadata: {meta_file}") - print(f"\nNext step:") - print(f" python run_optimization.py {output_file}") - print(f"{'='*60}") - - return slab, metadata + (outdir / "results" / "slab_setup_metadata.json").write_text(json.dumps(meta, indent=2)) + print(f"[slab_setup] Loaded {input_path} (atoms={n0})") + print(f"[slab_setup] Placed H at index {h_index} above O index {o_index} by {oh_distance} Å") + print(f"[slab_setup] Frozen atoms: {meta['num_frozen']} (z < {z_freeze} Å)") + print(f"[slab_setup] H neighbors within {neighbor_cutoff} Å: {neigh_sorted}") -if __name__ == "__main__": - import argparse - - parser = argparse.ArgumentParser(description="Prepare H* adsorption structure") - parser.add_argument("--input", default="inputs/slab_clean_2x2.in", help="Input structure file") - parser.add_argument("--o-index", type=int, default=20, help="Oxygen atom index for H placement") - parser.add_argument("--oh-dist", type=float, default=1.0, help="O-H distance (Å)") - parser.add_argument("--z-freeze", type=float, default=20.0, help="Freeze atoms below this z (Å)") - - args = parser.parse_args() - - slab, meta = setup_structure( + return slab, meta + + +def main(): + p = argparse.ArgumentParser(description="IrO2 slab setup: place H above O, freeze bottom layers, write outputs/") + p.add_argument("--input", default="inputs/slab_clean_2x2.in") + p.add_argument("--o-index", type=int, default=20, help="0-based index of O atom to adsorb H onto") + p.add_argument("--oh-dist", type=float, default=1.0) + p.add_argument("--z-freeze", type=float, default=20.0) + p.add_argument("--neighbor-cutoff", type=float, default=1.5) + p.add_argument("--outputs", default="outputs") + args = p.parse_args() + + setup_structure( input_file=args.input, o_index=args.o_index, oh_distance=args.oh_dist, z_freeze=args.z_freeze, + neighbor_cutoff=args.neighbor_cutoff, + outputs_dir=args.outputs, ) + + +if __name__ == "__main__": + main()