Siesta

In this section:

Description
Available Versions
Loading the Module
Serial Job Submission
MPI Job Submission
Example Input File (FDF)
Geometry Optimization
Band Structure Calculation
Density of States
Job Array for Parameter Scan
Results Analysis
References

Description 

According to the page of Siesta, SIESTA is a computational implementation method for performing efficient electronic structure calculations and ab initio molecular dynamics simulations of molecules and solids, using strictly localized bases.

Available Versions 

siesta/4.0.2
siesta/4.1-b3 (default)
siesta/4.1.5.mpi
siesta/4.1.5.serial
siesta/5.2.0.mpi
siesta/5.2.0.serial

Loading the Module 

# Serial version
module load siesta/5.2.0.serial

# MPI version
module load siesta/5.2.0.mpi

Serial Job Submission 

submit_siesta_serial.sh

#!/bin/bash
#SBATCH -J siesta_serial
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -t 24:00:00
#SBATCH --mem=8G

export INPUT="h2o.fdf"
export OUTPUT="h2o.out"

module load siesta/5.2.0.serial

job-nanny siesta < h2o.fdf > h2o.out

MPI Job Submission 

submit_siesta_mpi.sh

#!/bin/bash
#SBATCH -J siesta_mpi
#SBATCH -N 2
#SBATCH --ntasks-per-node=28
#SBATCH -t 48:00:00
#SBATCH --mem-per-cpu=2G

export INPUT="si.fdf"
export OUTPUT="si.out"

module load siesta/5.2.0.mpi

job-nanny mpirun -np $SLURM_NTASKS siesta < si.fdf > si.out

Example Input File (FDF)

h2o.fdf

# Siesta input for H2O molecule

SystemName          Water molecule
SystemLabel         h2o

# Lattice parameters (large box for molecule)
LatticeConstant     20.0 Bohr
%block LatticeVectors
  1.0  0.0  0.0
  0.0  1.0  0.0
  0.0  0.0  1.0
%endblock LatticeVectors

# Atomic coordinates
AtomicCoordinatesFormat  Ang
%block AtomicCoordinatesAndAtomicSpecies
  0.000   0.000   0.000  1  # O
  0.757   0.586   0.000  2  # H
 -0.757   0.586   0.000  2  # H
%endblock AtomicCoordinatesAndAtomicSpecies

# Atomic species
NumberOfSpecies        2
%block ChemicalSpeciesLabel
  1  8  O
  2  1  H
%endblock ChemicalSpeciesLabel

# Basis set
PAO.BasisSize      DZP
PAO.EnergyShift    0.02 Ry

# Exchange-correlation
XC.functional      GGA
XC.authors         PBE

# Mesh
MeshCutoff         200 Ry

# SCF options
MaxSCFIterations    200
DM.MixingWeight     0.1
DM.Tolerance        1.d-6

# Output
WriteForces         T
WriteCoorStep       T
WriteMDhistory      F

# Geometry optimization
MD.TypeOfRun        CG
MD.NumCGsteps       100
MD.MaxForceTol      0.01 eV/Ang

Geometry Optimization 

submit_siesta_opt.sh

#!/bin/bash
#SBATCH -J siesta_opt
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -t 36:00:00
#SBATCH --mem=8G

export INPUT="optimize.fdf"
export OUTPUT="optimize.out"

module load siesta/5.2.0.serial

job-nanny siesta < optimize.fdf > optimize.out

optimize.fdf

# Siesta geometry optimization

SystemName          Si bulk optimization
SystemLabel         si_opt

# Lattice parameters for Si
LatticeConstant     10.20 Bohr
%block LatticeVectors
  0.0  0.5  0.5
  0.5  0.0  0.5
  0.5  0.5  0.0
%endblock LatticeVectors

# Atomic positions
AtomicCoordinatesFormat  Fractional
NumberOfAtoms        2
%block AtomicCoordinatesAndAtomicSpecies
  0.00   0.00   0.00  1
  0.25   0.25   0.25  1
%endblock AtomicCoordinatesAndAtomicSpecies

NumberOfSpecies      1
%block ChemicalSpeciesLabel
  1  14  Si
%endblock ChemicalSpeciesLabel

# Basis
PAO.BasisSize        DZP
PAO.EnergyShift      0.02 Ry

# k-points
%block kgrid_Monkhorst_Pack
  4 0 0 0.0
  0 4 0 0.0
  0 0 4 0.0
%endblock kgrid_Monkhorst_Pack

# Optimization
MD.TypeOfRun         CG
MD.NumCGsteps        100
MD.MaxForceTol       0.01 eV/Ang

Band Structure Calculation 

submit_siesta_bands.sh

#!/bin/bash
#SBATCH -J siesta_bands
#SBATCH -N 2
#SBATCH --ntasks-per-node=28
#SBATCH -t 48:00:00
#SBATCH --mem-per-cpu=2G

export INPUT="bands.fdf"
export OUTPUT="bands.out"

module load siesta/5.2.0.mpi

job-nanny mpirun -np $SLURM_NTASKS siesta < bands.fdf > bands.out

bands.fdf

# Siesta band structure calculation

SystemName          Si bands
SystemLabel         si_bands

# Use optimized geometry
MD.TypeOfRun        BS   # Band structure calculation

# Path in k-space
%block BandLines
  1  0.0 0.0 0.0  \Gamma
 20  0.5 0.0 0.0  X
 15  0.5 0.5 0.0  W
 10  0.0 0.0 0.0  \Gamma
 15  0.5 0.5 0.5  L
%endblock BandLines

# Include eigenvalues in output
WriteEigenvalues    T

# Rest of input (similar to previous)

Density of States 

submit_siesta_dos.sh

#!/bin/bash
#SBATCH -J siesta_dos
#SBATCH -N 2
#SBATCH --ntasks-per-node=28
#SBATCH -t 36:00:00
#SBATCH --mem-per-cpu=2G

export INPUT="dos.fdf"
export OUTPUT="dos.out"

module load siesta/5.2.0.mpi

job-nanny mpirun -np $SLURM_NTASKS siesta < dos.fdf > dos.out

dos.fdf

# Siesta DOS calculation

SystemName          Si DOS
SystemLabel         si_dos

# Projected DOS
ProjectedDensity    T
PDOS.Kgrid          10 10 10
PDOS.Emin           -20.0 eV
PDOS.Emax           10.0 eV
PDOS.Broaden        0.2 eV

# Rest of input (similar to previous)

Job Array for Parameter Scan 

submit_siesta_array.sh

#!/bin/bash
#SBATCH -J siesta_array
#SBATCH --array=1-5
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -t 24:00:00
#SBATCH --mem=8G

MESH_VALUES=(100 150 200 250 300)
MESH=${MESH_VALUES[$SLURM_ARRAY_TASK_ID-1]}

export INPUT="template.fdf"
export OUTPUT="mesh_${MESH}/"

module load siesta/5.2.0.serial

mkdir -p mesh_${MESH}
cd mesh_${MESH}

# Modify MeshCutoff
sed "s/MeshCutoff.*/MeshCutoff $MESH Ry/" ../template.fdf > input.fdf

job-nanny siesta < input.fdf > siesta.out

# Extract total energy
grep "siesta: Total =" siesta.out | tail -1 >> ../energies.txt

Results Analysis 

analyze_siesta.sh

#!/bin/bash
#SBATCH -J analyze_siesta
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -t 01:00:00
#SBATCH --mem=2G

# Extract total energy
grep "siesta: Total =" *.out > energies.txt

# Extract final forces
grep "siesta: Atomic forces" -A 5 *.out > forces.txt

# Extract optimized geometry
for out in *.out; do
    echo "Final geometry for $out:" >> final_geometry.txt
    awk '/outcoor: Atomic coordinates/,/outcoor: End/' $out >> final_geometry.txt
    echo "" >> final_geometry.txt
done

# Extract gaps (if any)
grep "gap" *.out > gaps.txt

# Plot with Python
cat > plot_energy.py << 'EOF'
import numpy as np
import matplotlib.pyplot as plt

data = np.loadtxt('energies.txt', usecols=(0,1))
mesh = data[:,0]
energy = data[:,1]

plt.figure()
plt.plot(mesh, energy, 'bo-')
plt.xlabel('MeshCutoff (Ry)')
plt.ylabel('Total Energy (eV)')
plt.title('Energy convergence')
plt.grid(True)
plt.savefig('convergence.png')
EOF

python plot_energy.py

References 

Documentation: https://siesta-project.org/siesta/Documentation/
Tutorial: https://docs.siesta-project.org/projects/siesta/en/stable/tutorials/index.html
Pseudopotentials: https://docs.siesta-project.org/projects/siesta/en/stable/reference/pseudopotentials/