#!/usr/bin/env python3

import numpy as np
from ase import build
from ase.visualize import view
from ase.io import write

# 1) Build a graphene slab and add vacuum along the z-axis
#    - formula='C2' with lattice constant a=2.46 Å
#    - 5×5×1 supercell
#    - center() adds ~10 Å of vacuum in z; the CIF remains 3D with vacuum
slab = build.graphene(formula='C2', a=2.46, size=(5, 5, 1))
slab.center(vacuum=10, axis=2)

# 2) Option A: substitute by the atom index seen in your viewer (e.g., #24)
slab[24].symbol = 'N'           # alternatively: slab[24].number = 7

# --- If you prefer selecting by a point shown in the viewer, use Option B below
#     and comment out the line above. Adjust the coordinates to match your viewer.
# target = np.array([2.460, 4.261, 20.000])  # example from your screenshot (Å)
# i = np.linalg.norm(slab.get_positions() - target, axis=1).argmin()
# slab[i].symbol = 'N'

# 3) Visualize (optional)
# view(slab)

# 4) Write the CIF file (ASE infers the format from the extension)
write('ndg.cif', slab)
print('File written: ndg.cif')
