This tutorial targets beginners to intermediate users of SIESTA. This is the second part and it covers nudged elestic band (NEB) calculations. The first part covers basic runs, k-points convergence, bulk and slab geometry optimization, surface reconstructions, density of states, and workfucntion. The tutorial utilizes the atomistic simulation enviroment (ASE) python wrapper for SIESTA. All cells are supposed to run separately, but in sequence, there are no global variables. The tutorial requires SIESTA 5.2.4 and ASE (3.29 not strict). A few python packages should be installed like matplotlib, pymatgen, etc. You will need a proper instalation and linking between ASE and SIESTA, described in the ASE manuals. You will need to download and adjust the path to the pseudopotential files, the tutorial used pseudodojo SIESTA PBE pseudopotentials. Basic understanding of python is required. Once you go through it, you will enjoy coding computational materials science!
About the author: Dr. Aleksandar Staykov, associate professor at Kyushu Universtity, Japan. E-mail: alex@i2cner.kyushu-u.ac.jp
for more information on theory, read my book which covers DFT fundamentals, wavefunctions, methods and applicatoins:
Ab Initio Simulations in Materials Science: Hands-On Introduction to Electronic Structure Modeling with VASP
https://www.amazon.com/Initio-Simulations-Materials-Science-Hands/dp/B0FGSNL5QR
The Nudged Elastic Band (NEB) method is a computational technique used to determine the minimum‑energy path (MEP) between two known states of a system—typically an initial configuration and a final configuration connected by some physical process such as diffusion, adsorption, desorption, or a chemical reaction. Instead of guessing the transition pathway, NEB constructs a series of intermediate “images” of the system and links them together like an elastic band. During optimization, forces perpendicular to the path push the images toward the true MEP, while spring forces along the path maintain a smooth connection between them. The result is a physically meaningful reaction pathway and an accurate estimate of the activation energy barrier, which is essential for understanding kinetics, transition mechanisms, and rate‑limiting steps in materials and surface processes.
An optimized initial geometry is essential because all subsequent calculations—forces, stresses, electronic structure, and transition pathways—depend sensitively on how close the system is to its true equilibrium configuration. If the starting structure contains unrealistic bond lengths, strained angles, or residual forces, the SCF cycle may become unstable, the optimizer may take inefficient or incorrect steps, and the system may relax toward an artificial local minimum rather than the physically meaningful one. Poor initial geometries also distort quantities such as surface dipoles, work functions, reaction barriers, and NEB pathways, since these properties are defined relative to the equilibrium atomic arrangement. By ensuring that the initial and final geometries are fully relaxed with respect to the chosen method, basis set, and pseudopotentials, we eliminate spurious stresses and guarantee that all subsequent simulations begin from a consistent, physically accurate reference state.
We will use Ni atom migrating over graphene surface. First, we will converge k-points for graphene and optimize the starting and end points with adsorbed Ni. Then, we will perform the NEB calculations
from ase.build import graphene
from ase.visualize import view
from ase.calculators.siesta import Siesta
import os
from ase.units import Ry
import numpy as np
import matplotlib.pyplot as plt
os.mkdir("gr_siesta_kpts")
os.chdir("gr_siesta_kpts")
# Build graphene with total 8 Å vacuum along z
atoms = graphene(a=2.46, vacuum=8.0)
# Optional: center the sheet exactly in the cell along z
atoms.center(axis=2)
kk_array = []
energy_array = []
for kk in range(1, 15, 1): # loop over k-points from 1x1x1 to 8x8x8
atoms.calc = Siesta(label='c2',
xc='PBE',
pseudo_path='/Users/staykov/pseudodojo',
pseudo_qualifier='gga',
symlink_pseudos=True,
mesh_cutoff=200 * Ry,
energy_shift=0.01 * Ry,
basis_set='DZP',
spin='non-polarized',
kpts=(kk, kk, 1),
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 100},
)
energy = atoms.get_total_energy()
print("k-points:", kk)
kk_array.append(kk)
print("Energy:", energy)
energy_array.append(energy)
# Plot k-point energy convergence
plt.figure(figsize=(10, 6)) # Set the figure size (optional)
plt.plot(kk_array, energy_array, marker='o', linestyle='-', color='b') # Plot the data
plt.title('Energy vs k-points') # Add a title to the plot
plt.xlabel('k-points') # Label for the x-axis
plt.ylabel('Energy') # Label for the y-axis
plt.grid(True) # Add a grid (optional)
plt.show() # Display the plot
os.chdir("../")
view(atoms, viewer='x3d')
Gamma-point calculation with interaction between periodic images Some features might not work optimally: e.g. DM initialization from atomic data Job completed
k-points: 1 Energy: -320.285938
Job completed
k-points: 2 Energy: -325.10686
Job completed
k-points: 3 Energy: -326.622313
Job completed
k-points: 4 Energy: -327.237976
Job completed
k-points: 5 Energy: -327.182288
Job completed
k-points: 6 Energy: -327.148775
Job completed
k-points: 7 Energy: -327.219756
Job completed
k-points: 8 Energy: -327.20646
Job completed
k-points: 9 Energy: -327.189635
Job completed
k-points: 10 Energy: -327.211477
Job completed
k-points: 11 Energy: -327.207561
Job completed
k-points: 12 Energy: -327.199597
Job completed
k-points: 13 Energy: -327.20903 k-points: 14 Energy: -327.207521
Job completed
Initial geometry of Ni on graphene
from ase.build import graphene
from ase.visualize import view
from ase.calculators.siesta import Siesta
import os
from ase.optimize.bfgs import BFGS
from ase.units import Ry
from ase import Atom
from ase.io import write
# Build graphene with total 8 Å vacuum along z
atoms = graphene(a=2.46, vacuum=8.0)
# Optional: center the sheet exactly in the cell along z
atoms.center(axis=2)
atoms = atoms.repeat((3,3,1))
atoms.append(Atom('Ni', (2.5, 1.5, 10.0)))
os.mkdir("gr1_siesta_opt")
os.chdir("gr1_siesta_opt")
atoms.calc = Siesta(label='Graphene_Ni',
xc='PBE',
pseudo_path='/Users/staykov/pseudodojo',
pseudo_qualifier='gga',
symlink_pseudos=True,
mesh_cutoff=230 * Ry,
energy_shift=0.01 * Ry,
basis_set='DZP',
spin='non-polarized',
kpts=(3, 3, 1),
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 200},
)
energy = atoms.get_total_energy()
opt = BFGS(atoms, trajectory='atoms.traj')
opt.run(fmax=0.03)
os.chdir("../")
write('gr1_siesta.traj', atoms)
view(atoms, viewer='x3d')
Job completed
Step Time Energy fmax BFGS: 0 02:03:42 -7503.991448 3.959180
Job completed
BFGS: 1 02:03:54 -7504.260850 4.002311
Job completed
BFGS: 2 02:04:11 -7505.074641 2.860394
Job completed
BFGS: 3 02:04:25 -7505.376135 0.514493
Job completed
BFGS: 4 02:04:35 -7505.394248 0.453758
Job completed
BFGS: 5 02:04:47 -7505.405263 0.429110
Job completed
BFGS: 6 02:05:01 -7505.424424 0.264263
Job completed
BFGS: 7 02:05:14 -7505.431695 0.309149
Job completed
BFGS: 8 02:05:28 -7505.449402 0.234121
Job completed
BFGS: 9 02:05:40 -7505.456623 0.148337
Job completed
BFGS: 10 02:05:51 -7505.460228 0.153918
Job completed
BFGS: 11 02:06:01 -7505.462354 0.154895
Job completed
BFGS: 12 02:06:13 -7505.464785 0.121334
Job completed
BFGS: 13 02:06:25 -7505.466130 0.094136
Job completed
BFGS: 14 02:06:36 -7505.466934 0.073182
Job completed
BFGS: 15 02:06:46 -7505.467489 0.058289
Job completed
BFGS: 16 02:06:54 -7505.467928 0.044353
Job completed
BFGS: 17 02:07:03 -7505.468194 0.040088 BFGS: 18 02:07:10 -7505.468340 0.023008
Job completed
Final geometry
from ase.build import graphene
from ase.visualize import view
from ase.calculators.siesta import Siesta
import os
from ase.optimize.bfgs import BFGS
from ase.units import Ry
from ase import Atom
from ase.io import write
from ase.io import read
atoms = read('gr1_siesta.traj')
mg_index = [i for i, a in enumerate(atoms) if a.symbol == 'Ni'][0]
dx = 1.5
dy = 2
atoms[mg_index].position += [dx, dy, 0.0]
os.mkdir("gr2_siesta_opt")
os.chdir("gr2_siesta_opt")
atoms.calc = Siesta(label='Graphene_Ni',
xc='PBE',
pseudo_path='/Users/staykov/pseudodojo',
pseudo_qualifier='gga',
symlink_pseudos=True,
mesh_cutoff=230 * Ry,
energy_shift=0.01 * Ry,
basis_set='DZP',
spin='non-polarized',
kpts=(3, 3, 1),
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 200},
)
energy = atoms.get_total_energy()
opt = BFGS(atoms, trajectory='atoms.traj')
opt.run(fmax=0.03)
os.chdir("../")
write('gr2_siesta.traj', atoms)
view(atoms, viewer='x3d')
Job completed
Step Time Energy fmax BFGS: 0 02:08:42 -7505.101137 2.653569
Job completed
BFGS: 1 02:08:55 -7505.291189 1.264890
Job completed
BFGS: 2 02:09:06 -7505.341076 0.840949
Job completed
BFGS: 3 02:09:16 -7505.363142 0.781935
Job completed
BFGS: 4 02:09:27 -7505.401438 0.780404
Job completed
BFGS: 5 02:09:40 -7505.429630 0.648716
Job completed
BFGS: 6 02:09:53 -7505.444627 0.404577
Job completed
BFGS: 7 02:10:05 -7505.454556 0.226003
Job completed
BFGS: 8 02:10:18 -7505.459737 0.140047
Job completed
BFGS: 9 02:10:29 -7505.461282 0.118983
Job completed
BFGS: 10 02:10:40 -7505.462648 0.098937
Job completed
BFGS: 11 02:10:50 -7505.464240 0.072891
Job completed
BFGS: 12 02:11:03 -7505.465232 0.067889
Job completed
BFGS: 13 02:11:14 -7505.465929 0.075571
Job completed
BFGS: 14 02:11:25 -7505.466648 0.083473
Job completed
BFGS: 15 02:11:34 -7505.467440 0.078246
Job completed
BFGS: 16 02:11:44 -7505.468072 0.039563 BFGS: 17 02:11:54 -7505.468292 0.023693
Job completed
The NEB calculation proceeds by first reading the initial and final geometries and constructing a series of intermediate images using a linear interpolation algorithm. These images define the discrete points along the reaction pathway. Each image is then assigned its own SIESTA calculator so that forces and energies are evaluated consistently and independently. With the band fully defined, the NEB optimizer adjusts the images by removing the parallel component of the true force and adding spring forces along the path, gradually pushing the band toward the minimum‑energy pathway. Once the optimization converges, the code extracts the energies of all images, plots the NEB energy profile, and identifies the transition state as the image with the highest energy. Finally, the geometry of this transition‑state image is written out for visualization and further analysis, providing a complete picture of the reaction barrier and the atomic configuration at the saddle point.
import numpy as np
import matplotlib.pyplot as plt
import os
from ase.io import read, write
from ase.mep import NEB
from ase.optimize import BFGS
from ase.calculators.siesta import Siesta
from ase.units import Ry
from ase.visualize import view
# ------------------------------------------------------------
# 1. Read initial and final geometries
# ------------------------------------------------------------
initial = read('gr1_siesta.traj') # or initial.xyz / initial.traj
final = read('gr2_siesta.traj')
os.mkdir("neb_siesta")
os.chdir("neb_siesta")
# ------------------------------------------------------------
# 2. Create NEB images (initial + 3 intermediates + final)
# ------------------------------------------------------------
n_images = 5
images = [initial]
for i in range(n_images - 2):
images.append(initial.copy())
images.append(final)
neb = NEB(images)
neb.interpolate() # linear interpolation
# ------------------------------------------------------------
# 3. Attach a *separate* SIESTA calculator to each image
# ------------------------------------------------------------
def make_siesta_calc(label):
return Siesta(label='neb',
xc='PBE',
pseudo_path='/Users/staykov/pseudodojo',
pseudo_qualifier='gga',
symlink_pseudos=True,
mesh_cutoff=230 * Ry,
energy_shift=0.01 * Ry,
basis_set='DZP',
spin='non-polarized',
kpts=(3, 3, 1),
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 200},
)
for i, img in enumerate(images):
img.calc = make_siesta_calc(f'neb_img_{i}')
# ------------------------------------------------------------
# 4. Optimize the NEB path
# ------------------------------------------------------------
opt = BFGS(neb, logfile='neb.log')
opt.run(fmax=0.05) # convergence criterion for forces on images
# ------------------------------------------------------------
# 5. Collect energies and compute relative profile
# ------------------------------------------------------------
energies = [img.get_potential_energy() for img in images]
E0 = energies[0]
rel_energies = [E - E0 for E in energies]
# Identify TS (highest energy image)
ts_index = int(np.argmax(energies))
ts_image = images[ts_index]
write('TS.traj', ts_image)
os.chdir("../")
# ------------------------------------------------------------
# 6. Plot NEB energy profile
# ------------------------------------------------------------
plt.figure()
plt.plot(range(n_images), rel_energies, '-o')
plt.xlabel('Image index')
plt.ylabel('Relative energy (eV)')
plt.title('NEB Energy Profile')
plt.grid(True)
plt.tight_layout()
plt.savefig('neb_profile.png', dpi=200)
print("Absolute energies (eV):", energies)
print("TS image index:", ts_index)
print("TS geometry saved as TS.traj")
view(ts_image, viewer='x3d')
/opt/miniconda3/lib/python3.13/site-packages/ase/mep/neb.py:329: UserWarning: The default method has changed from 'aseneb' to 'improvedtangent'. The 'aseneb' method is an unpublished, custom implementation that is not recommended as it frequently results in very poor bands. Please explicitly set method='improvedtangent' to silence this warning, or set method='aseneb' if you strictly require the old behavior (results may vary). See: https://gitlab.com/ase/ase/-/merge_requests/3952 warnings.warn( Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed
Absolute energies (eV): [-7505.454216, -7505.082417, -7504.759844, -7505.081602, -7505.454246] TS image index: 2 TS geometry saved as TS.traj
Job completed
Basis Set Superposition Error (BSSE) in SIESTA arises from the use of localized basis functions, such as Numerical Atomic Orbitals (NAOs), where the spatial overlap of basis functions from neighboring atoms provides an artificial improvement in the description of an individual atom's electronic wavefunction. When two atoms are close together, each atom "borrows" the basis functions of its neighbor to better represent its electron density, leading to a lower total energy; as they move apart, this overlapping basis is lost, making the separated atoms appear higher in energy. This results in a systematic overestimation of binding energies or interaction energies because the energy gain is partially due to the mathematical improvement of the basis set rather than purely physical bonding.
Consequently, the overestimation of the Ni binding energy at the graphene resting sites due to BSSE may lead to an artificially strong interaction. This, in turn, can result in inflated activation barriers for the migration of Ni atoms across the surface, potentially leading to an inaccurate representation of the diffusion kinetics.
In systems utilizing localized basis sets like SIESTA, the ghost atom method is the standard correction for BSSE. Ghost atoms provide the 'missing' basis functions at neighboring sites during the calculation of isolated atoms. By ensuring that the orbital space available to the Ni atom is consistent whether it is alone or sitting on the graphene surface, the ghost atoms prevent the artificial lowering of energy caused by overlapping basis functions. This results in a stabilized calculation of the adsorption energy, ensuring that the predicted migration barriers are a result of the potential energy surface rather than numerical artifacts.
In the subsequent calculation, we incorporate carbon ghost atoms at the positions of the nearest graphene sites. By providing these ghost orbitals, the Ni atom's basis set remains consistent whether it is isolated or adsorbed. This correction mitigates the BSSE, resulting in a more accurate—and typically lower—adsorption energy by removing the artificial energy gains from basis function overlap.
The second graphene layer shown in the figure below is a 'ghost' layer. It consists entirely of basis functions (orbitals) without associated nuclei or electrons. As such, it remains fixed during the geometry optimization and serves exclusively to provide the necessary orbital space to mitigate basis set superposition error.
from ase.build import graphene
from ase.visualize import view
from ase.calculators.siesta import Siesta
import os
from ase.optimize.bfgs import BFGS
from ase.units import Ry
from ase import Atom
from ase.io import write
from ase.calculators.siesta.parameters import Species
from ase.constraints import FixAtoms
# Build graphene with total 8 Å vacuum along z
atoms = graphene(a=2.46, vacuum=8.0)
# Optional: center the sheet exactly in the cell along z
atoms.center(axis=2)
atoms = atoms.repeat((3,3,1))
ghost_layer = atoms.copy()
ghost_layer.positions[:, 2] += 2
ghost_layer.set_tags([1] * len(ghost_layer))
atoms = atoms + ghost_layer
# Append Ni adatom
atoms.append(Atom('Ni', (2.5, 1.5, 10.0)))
# Find ghost atom indices (tag == 1) and fix them so optimizer won't move them
ghost_indices = [i for i, t in enumerate(atoms.get_tags()) if t == 1]
constraint = FixAtoms(indices=ghost_indices)
atoms.set_constraint(constraint)
species = [
Species(
symbol='C',
basis_set='DZP'
),
Species(
symbol='Ni',
basis_set='DZP'
),
Species(
symbol='C',
basis_set='DZ',
tag=1,
ghost=True
)
]
# Create output directory safely
os.makedirs("gr1_siesta_ghost_opt", exist_ok=True)
os.chdir("gr1_siesta_ghost_opt")
atoms.calc = Siesta(label='Graphene_Ni_ghost',
xc='PBE',
species=species,
pseudo_path='/Users/staykov/pseudodojo',
pseudo_qualifier='gga',
symlink_pseudos=True,
mesh_cutoff=230 * Ry,
energy_shift=0.01 * Ry,
basis_set='DZP',
spin='non-polarized',
kpts=(3, 3, 1),
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 200},
)
energy = atoms.get_total_energy()
opt = BFGS(atoms, trajectory='atoms.traj')
opt.run(fmax=0.03)
os.chdir("../")
write('gr1_siesta_ghost_opt.traj', atoms)
view(atoms, viewer='x3d')
Job completed
Step Time Energy fmax BFGS: 0 02:34:02 -7506.179218 2.366846
Job completed
BFGS: 1 02:34:27 -7506.281305 2.344015
Job completed
BFGS: 2 02:35:07 -7506.724608 1.592271
Job completed
BFGS: 3 02:36:08 -7506.796014 0.986888
Job completed
BFGS: 4 02:36:53 -7506.867876 0.469328
Job completed
BFGS: 5 02:37:23 -7506.883072 0.337218
Job completed
BFGS: 6 02:37:58 -7506.893434 0.220944
Job completed
BFGS: 7 02:38:22 -7506.897540 0.222113
Job completed
BFGS: 8 02:38:53 -7506.912079 0.331010
Job completed
BFGS: 9 02:39:23 -7506.922800 0.246036
Job completed
BFGS: 10 02:39:50 -7506.928428 0.097307
Job completed
BFGS: 11 02:40:13 -7506.930046 0.100639
Job completed
BFGS: 12 02:40:32 -7506.931608 0.112107
Job completed
BFGS: 13 02:40:50 -7506.932708 0.070960
Job completed
BFGS: 14 02:41:09 -7506.933378 0.079281
Job completed
BFGS: 15 02:41:24 -7506.933898 0.070072
Job completed
BFGS: 16 02:41:45 -7506.934525 0.068020
Job completed
BFGS: 17 02:42:07 -7506.935260 0.076913
Job completed
BFGS: 18 02:42:30 -7506.936216 0.102184
Job completed
BFGS: 19 02:42:55 -7506.937450 0.104209
Job completed
BFGS: 20 02:43:20 -7506.938887 0.096719
Job completed
BFGS: 21 02:43:44 -7506.940509 0.100503
Job completed
BFGS: 22 02:44:05 -7506.941813 0.113295
Job completed
BFGS: 23 02:44:24 -7506.943363 0.132146
Job completed
BFGS: 24 02:44:45 -7506.946508 0.167140
Job completed
BFGS: 25 02:45:09 -7506.954672 0.364167
Job completed
BFGS: 26 02:45:46 -7506.975087 0.925262
Job completed
BFGS: 27 02:46:16 -7506.998383 1.004333
Job completed
BFGS: 28 02:46:51 -7507.082669 0.560871
Job completed
BFGS: 29 02:47:18 -7507.102233 0.270001
Job completed
BFGS: 30 02:47:39 -7507.113473 0.194098
Job completed
BFGS: 31 02:48:06 -7507.118171 0.165110
Job completed
BFGS: 32 02:48:34 -7507.121043 0.099777 BFGS: 33 02:49:00 -7507.121489 0.020628
Job completed
from ase.build import graphene
from ase.visualize import view
from ase.calculators.siesta import Siesta
import os
from ase.optimize.bfgs import BFGS
from ase.units import Ry
from ase import Atom
from ase.io import write
from ase.io import read
from ase.calculators.siesta.parameters import Species
from ase.constraints import FixAtoms
atoms = read('gr1_siesta_ghost_opt.traj')
mg_index = [i for i, a in enumerate(atoms) if a.symbol == 'Ni'][0]
dx = 1.5
dy = 2
atoms[mg_index].position += [dx, dy, 0.0]
# Find ghost atom indices (tag == 1) and fix them so optimizer won't move them
ghost_indices = [i for i, t in enumerate(atoms.get_tags()) if t == 1]
constraint = FixAtoms(indices=ghost_indices)
atoms.set_constraint(constraint)
species = [
Species(
symbol='C',
basis_set='DZP'
),
Species(
symbol='Ni',
basis_set='DZP'
),
Species(
symbol='C',
basis_set='DZ',
tag=1,
ghost=True
)
]
# Create output directory safely
os.makedirs("gr2_siesta_ghost_opt", exist_ok=True)
os.chdir("gr2_siesta_ghost_opt")
atoms.calc = Siesta(label='Graphene_Ni_ghost',
xc='PBE',
species=species,
pseudo_path='/Users/staykov/pseudodojo',
pseudo_qualifier='gga',
symlink_pseudos=True,
mesh_cutoff=230 * Ry,
energy_shift=0.01 * Ry,
basis_set='DZP',
spin='non-polarized',
kpts=(3, 3, 1),
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 200},
)
energy = atoms.get_total_energy()
opt = BFGS(atoms, trajectory='atoms.traj')
opt.run(fmax=0.03)
os.chdir("../")
write('gr2_siesta_ghost_opt.traj', atoms)
view(atoms, viewer='x3d')
Job completed
Step Time Energy fmax BFGS: 0 02:52:18 -7506.871804 1.692035
Job completed
BFGS: 1 02:52:42 -7506.986702 1.003210
Job completed
BFGS: 2 02:53:05 -7507.018208 0.610953
Job completed
BFGS: 3 02:53:31 -7507.029922 0.523387
Job completed
BFGS: 4 02:54:08 -7507.048125 0.562125
Job completed
BFGS: 5 02:54:46 -7507.062208 0.643801
Job completed
BFGS: 6 02:55:31 -7507.082783 0.660669
Job completed
BFGS: 7 02:56:11 -7507.098035 0.503809
Job completed
BFGS: 8 02:56:40 -7507.109672 0.319639
Job completed
BFGS: 9 02:57:14 -7507.117732 0.095287
Job completed
BFGS: 10 02:57:37 -7507.118219 0.082081
Job completed
BFGS: 11 02:58:02 -7507.119056 0.066823
Job completed
BFGS: 12 02:58:22 -7507.119770 0.058475
Job completed
BFGS: 13 02:58:44 -7507.120353 0.073538
Job completed
BFGS: 14 02:59:05 -7507.120935 0.075341
Job completed
BFGS: 15 02:59:26 -7507.121366 0.052643 BFGS: 16 02:59:43 -7507.121600 0.028841
Job completed
import numpy as np
import matplotlib.pyplot as plt
import os
from ase.io import read, write
from ase.mep import NEB
from ase.optimize import BFGS
from ase.calculators.siesta import Siesta
from ase.units import Ry
from ase.visualize import view
from ase.constraints import FixAtoms
from ase.calculators.siesta.parameters import Species
# ------------------------------------------------------------
# 1. Read initial and final geometries
# ------------------------------------------------------------
initial = read('gr1_siesta_ghost_opt.traj') # or initial.xyz / initial.traj
final = read('gr2_siesta_ghost_opt.traj')
os.mkdir("neb_siesta_ghost")
os.chdir("neb_siesta_ghost")
# ------------------------------------------------------------
# 2. Create NEB images (initial + 3 intermediates + final)
# ------------------------------------------------------------
n_images = 5
images = [initial]
for i in range(n_images - 2):
images.append(initial.copy())
images.append(final)
neb = NEB(images)
neb.interpolate() # linear interpolation
# ------------------------------------------------------------
# 3. Attach a *separate* SIESTA calculator to each image
# ------------------------------------------------------------
def make_siesta_calc(label):
return Siesta(label='c2',
xc='PBE',
species=species,
pseudo_path='/Users/staykov/pseudodojo',
pseudo_qualifier='gga',
symlink_pseudos=True,
mesh_cutoff=230 * Ry,
energy_shift=0.01 * Ry,
basis_set='DZP',
spin='non-polarized',
kpts=(3, 3, 1),
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 200},
)
for i, img in enumerate(images):
ghost_indices = [i for i, t in enumerate(img.get_tags()) if t == 1]
constraint = FixAtoms(indices=ghost_indices)
img.set_constraint(constraint)
species = [
Species(
symbol='C',
basis_set='DZP'
),
Species(
symbol='Ni',
basis_set='DZP'
),
Species(
symbol='C',
basis_set='DZ',
tag=1,
ghost=True
)
]
img.calc = make_siesta_calc(f'neb_img_{i}')
# ------------------------------------------------------------
# 4. Optimize the NEB path
# ------------------------------------------------------------
opt = BFGS(neb, logfile='neb.log')
opt.run(fmax=0.05) # convergence criterion for forces on images
# ------------------------------------------------------------
# 5. Collect energies and compute relative profile
# ------------------------------------------------------------
energies = [img.get_potential_energy() for img in images]
E0 = energies[0]
rel_energies = [E - E0 for E in energies]
# Identify TS (highest energy image)
ts_index = int(np.argmax(energies))
ts_image = images[ts_index]
write('TS.traj', ts_image)
os.chdir("../")
# ------------------------------------------------------------
# 6. Plot NEB energy profile
# ------------------------------------------------------------
plt.figure()
plt.plot(range(n_images), rel_energies, '-o')
plt.xlabel('Image index')
plt.ylabel('Relative energy (eV)')
plt.title('NEB Energy Profile')
plt.grid(True)
plt.tight_layout()
plt.savefig('neb_profile.png', dpi=200)
print("Absolute energies (eV):", energies)
print("TS image index:", ts_index)
print("TS geometry saved as TS.traj")
view(ts_image, viewer='x3d')
Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed
Absolute energies (eV): [-7507.111451, -7506.905358, -7506.818043, -7506.907304, -7507.111617] TS image index: 2 TS geometry saved as TS.traj
Job completed
The inclusion of the graphene ghost layer resulted in a 50% reduction in the activation energy for Ni migration on the graphene surface, bringing the results into excellent agreement with values obtained via plane-wave calculations.