Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
45 changes: 26 additions & 19 deletions .github/workflows/python-package-conda.yml
Original file line number Diff line number Diff line change
Expand Up @@ -5,30 +5,37 @@ 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

- name: Set up Miniconda
uses: conda-incubator/setup-miniconda@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
auto-update-conda: true
environment-file: environment.yml
activate-environment: base
auto-activate-base: true

- 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
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 install pytest
pytest -q -m "not slow"
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 iro2_test_1.py
ls -R .

- name: Upload simulation outputs
uses: actions/upload-artifact@v4
with:
name: irO2-simulation-outputs
path: outputs/
4 changes: 2 additions & 2 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -17,8 +17,8 @@ build/
.ipynb_checkpoints/

# Common simulation outputs (adjust as needed)
*.traj
#*.traj
*.xyz
*.cif
*.log
#*.log
*.out
206 changes: 103 additions & 103 deletions iro2_test_1.py
Original file line number Diff line number Diff line change
@@ -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()
Loading
Loading