SIESTA

Introduction

SIESTA is a density-functional theory code for very large systems based on atomic orbital (LCAO) basis sets.

Environment variables

The environment variable ASE_SIESTA_COMMAND must hold the command to invoke the siesta calculation. The variable must be a string where PREFIX.fdf/PREFIX.out are the placeholders for the input/output files. This variable allows you to specify serial or parallel execution of SIESTA. Examples: siesta < PREFIX.fdf > PREFIX.out and mpirun -np 4 /bin/siesta4.0 < PREFIX.fdf > PREFIX.out.

A default directory holding pseudopotential files .vps/.psf can be defined to avoid defining this every time the calculator is used. This directory can be set by the environment variable SIESTA_PP_PATH.

Set both environment variables in your shell configuration file:

$ export ASE_SIESTA_COMMAND="siesta < PREFIX.fdf > PREFIX.out"
$ export SIESTA_PP_PATH=$HOME/mypps

Alternatively, the path to the pseudopotentials can be given in the calculator initialization. Listed below all parameters related to pseudopotential control.

keyword

type

default value

description

pseudo_path

str

None

Directory for pseudopotentials to use None means using $SIESTA_PP_PATH

pseudo_qualifier

str

None

String for picking out specific type type of pseudopotentials. Giving example means that H.example.psf or H.example.vps will be used. None means that the XC.functional keyword is used, e.g. H.lda.psf

symlink_pseudos

bool

---

Whether psedos will be sym-linked into the execution directory. If False they will be copied in stead. Default is True on Unix and False on Windows.

SIESTA Calculator

These parameters are set explicitly and overrides the native values if different.

keyword

type

default value

description

label

str

'siesta'

Name of the output file

mesh_cutoff

float

200*Ry

Mesh cut-off energy in eV

xc

str

'LDA'

Exchange-correlation functional. Corresponds to either XC.functional or XC.authors keyword in SIESTA

energy_shift

float

100 meV

Energy shift for determining cutoff radii

kpts

list

[1,1,1]

Monkhorst-Pack k-point sampling

basis_set

str

DZP

Type of basis set (‘SZ’, ‘DZ’, ‘SZP’, ‘DZP’)

spin

str

'non-polarized'

The spin approximation used, must be either 'non-polarized', 'collinear', 'non-collinear' or 'spin-orbit'.

species

list

[]

A method for specifying the basis set for some atoms.

Most other parameters are set to the default values of the native interface.

Extra FDF parameters

The SIESTA code reads the input parameters for any calculation from a .fdf file. This means that you can set parameters by manually setting entries in this input .fdf file. This is done by the argument:

>>> Siesta(fdf_arguments={'variable_name': value, 'other_name': other_value})

For example, the DM.MixingWeight can be set using

>>> Siesta(fdf_arguments={'DM.MixingWeight': 0.01})

The explicit fdf arguments will always override those given by other keywords, even if it breaks calculator functionality. The complete list of the FDF entries can be found in the official SIESTA manual.

Example

Here is an example of how to calculate the total energy for bulk Silicon, using a double-zeta basis generated by specifying a given energy-shift:

>>> from ase import Atoms
>>> from ase.calculators.siesta import Siesta
>>> from ase.units import Ry
>>>
>>> a0 = 5.43
>>> bulk = Atoms('Si2', [(0, 0, 0),
...                      (0.25, 0.25, 0.25)],
...              pbc=True)
>>> b = a0 / 2
>>> bulk.set_cell([(0, b, b),
...                (b, 0, b),
...                (b, b, 0)], scale_atoms=True)
>>>
>>> calc = Siesta(label='Si',
...               xc='PBE',
...               mesh_cutoff=200 * Ry,
...               energy_shift=0.01 * Ry,
...               basis_set='DZ',
...               kpts=[10, 10, 10],
...               fdf_arguments={'DM.MixingWeight': 0.1,
...                              'MaxSCFIterations': 100},
...               )
>>> bulk.calc = calc
>>> e = bulk.get_potential_energy()

Here, the only input information on the basis set is, that it should be double-zeta (basis='DZP') and that the confinement potential should result in an energy shift of 0.01 Rydberg (the energy_shift=0.01 * Ry keyword). Sometimes it can be necessary to specify more information on the basis set.

Defining Custom Species

Standard basis sets can be set by the keyword basis_set directly, but for anything more complex than one standard basis size for all elements, a list of species must be defined. Each specie is identified by atomic element and the tag set on the atom.

For instance if we wish to investigate a H2 molecule and put a ghost atom (the basis set corresponding to an atom but without the actual atom) in the middle with a special type of basis set you would write:

>>> from ase.calculators.siesta.parameters import Specie, PAOBasisBlock
>>> from ase import Atoms
>>> from ase.calculators.siesta import Siesta
>>> atoms = Atoms(
...     '3H',
...     [(0.0, 0.0, 0.0),
...      (0.0, 0.0, 0.5),
...      (0.0, 0.0, 1.0)],
...     cell=[10, 10, 10])
>>> atoms.set_tags([0, 1, 0])
>>>
>>> basis_set = PAOBasisBlock(
... """1
... 0  2 S 0.2
... 0.0 0.0""")
>>>
>>> siesta = Siesta(
...     species=[
...         Specie(symbol='H', tag=None, basis_set='SZ'),
...         Specie(symbol='H', tag=1, basis_set=basis_set, ghost=True)])
>>>
>>> atoms.calc = siesta

When more species are defined, species defined with a tag has the highest priority. General species with tag=None has a lower priority. Finally, if no species apply to an atom, the general calculator keywords are used.

Pseudopotentials

Pseudopotential files in the .psf or .vps formats are needed. Pseudopotentials generated from the ABINIT code and converted to the SIESTA format are available in the SIESTA website. A database of user contributed pseudopotentials is also available there.

Optimized GGA–PBE pseudos and DZP basis sets for some common elements are also available from the SIMUNE website.

You can also find an on-line pseudopotential generator from the OCTOPUS code.

Species can also be used to specify pseudopotentials:

>>> specie = Specie(symbol='H', tag=1, pseudopotential='H.example.psf')

When specifying the pseudopotential in this manner, both absolute and relative paths can be given. Relative paths are interpreted as relative to the set pseudopotential path.

Restarting from an old Calculation

If you want to rerun an old SIESTA calculation, whether made using the ASE interface or not, you can set the keyword restart to the siesta .XV file. The keyword ignore_bad_restart (True/False) will decide whether a broken file will result in an error(False) or the whether the calculator will simply continue without the restart file.

Choosing the coordinate format

If you are mainly using ASE to generate SIESTA files for relaxation with native SIESTA relaxation, you may want to write the coordinates in the Z-matrix format which will best allow you to use the advanced constraints present in SIESTA.

keyword

type

default value

description

atomic_coord_format

str

'xyz'

Choose between 'xyz' and 'zmatrix' for the format that coordinates will be written in.

Siesta Calculator Class

class ase.calculators.siesta.siesta.Siesta(command=None, profile=None, directory='.', **kwargs)[source]

Calculator interface to the SIESTA code.

ASE interface to the SIESTA code.

Parameters:
  • label (-) – The basename of all files created during calculation.

  • mesh_cutoff (-) – Energy in eV. The mesh cutoff energy for determining number of grid points in the matrix-element calculation.

  • energy_shift (-) – Energy in eV The confining energy of the basis set generation.

  • kpts (-) – Tuple of 3 integers, the k-points in different directions.

  • xc (-) – The exchange-correlation potential. Can be set to any allowed value for either the Siesta XC.funtional or XC.authors keyword. Default “LDA”

  • basis_set (-) – “SZ”|”SZP”|”DZ”|”DZP”|”TZP”, strings which specify the type of functions basis set.

  • spin (-) – “non-polarized”|”collinear”| “non-collinear|spin-orbit”. The level of spin description to be used.

  • species (-) – None|list of Species objects. The species objects can be used to to specify the basis set, pseudopotential and whether the species is ghost. The tag on the atoms object and the element is used together to identify the species.

  • pseudo_path (-) – None|path. This path is where pseudopotentials are taken from. If None is given, then then the path given in $SIESTA_PP_PATH will be used.

  • pseudo_qualifier (-) – None|string. This string will be added to the pseudopotential path that will be retrieved. For hydrogen with qualifier “abc” the pseudopotential “H.abc.psf” will be retrieved.

  • symlink_pseudos (-) – None|bool If true, symlink pseudopotentials into the calculation directory, else copy them. Defaults to true on Unix and false on Windows.

  • atoms (-) – The Atoms object.

  • restart (-) – str. Prefix for restart file. May contain a directory. Default is None, don’t restart.

  • fdf_arguments (-) – Explicitly given fdf arguments. Dictonary using Siesta keywords as given in the manual. List values are written as fdf blocks with each element on a separate line, while tuples will write each element in a single line. ASE units are assumed in the input.

  • atomic_coord_format (-) – “xyz”|”zmatrix”, strings to switch between the default way of entering the system’s geometry (via the block AtomicCoordinatesAndAtomicSpecies) and a recent method via the block Zmatrix. The block Zmatrix allows to specify basic geometry constrains such as realized through the ASE classes FixAtom, FixedLine and FixedPlane.

Excited states calculations

The PyNAO code can be used to access excited state properties after having obtained the ground state properties with SIESTA. PyNAO allows to perform

  • Time Dependent Density Functional Theory (TDDFT) calculations

  • GW approximation calculations

  • Bethe-Salpeter equation (BSE) calculations.

Example of code to calculate polarizability of CH4 molecule:

from ase.calculators.siesta.siesta_lrtddft import SiestaLRTDDFT
from ase.build import molecule
import numpy as np

# Define the systems
ch4 = molecule('CH4')
lrtddft = SiestaLRTDDFT(label="siesta", xc_code='LDA,PZ')

# run DFT with siesta
lrtddft.get_ground_state(ch4)

# Run TDDFT with PyNAO
freq = np.arange(0.0, 25.0, 0.5)
pmat = lrtddft.get_polarizability(freq)

import matplotlib.pyplot as plt

plt.plot(freq, pmat[0, 0, :].imag)
plt.show()

Raman Calculations with SIESTA and PyNAO

It is possible to calculate the Raman spectra with SIESTA and PyNAO using the Raman function of the vibration module:

from ase.calculators.siesta.siesta_lrtddft import RamanCalculatorInterface
from ase.calculators.siesta import Siesta
from ase.vibrations.raman import StaticRamanCalculator
from ase.vibrations.placzek import PlaczekStatic
from ase.build import molecule

n2 = molecule('N2')

# enter siesta input
n2.calc = Siesta(
    basis_set='DZP',
    fdf_arguments={
        'COOP.Write': True,
        'WriteDenchar': True,
        'XML.Write': True})

name = 'n2'
pynao_args = dict(label="siesta", jcutoff=7, iter_broadening=0.15,
                  xc_code='LDA,PZ', tol_loc=1e-6, tol_biloc=1e-7)
rm = StaticRamanCalculator(n2, RamanCalculatorInterface, name=name,
                           delta=0.011, exkwargs=pynao_args)
rm.run()

Pz = PlaczekStatic(n2, name=name)
e_vib = Pz.get_energies()
Pz.summary()

Further Examples

See also ase/test/calculators/siesta/lrtddft for further examples on how the calculator can be used.

Siesta lrtddft Class

class ase.calculators.siesta.siesta_lrtddft.SiestaLRTDDFT(initialize=False, **kw)[source]

Interface for linear response TDDFT for Siesta via PyNAO

When using PyNAO please cite the papers indicated in the documentation

Parameters:
  • initialize (bool) – To initialize the tddft calculations before calculating the polarizability Can be useful to calculate multiple frequency range without the need to recalculate the kernel

  • kw (dictionary) – keywords for the tddft_iter function from PyNAO

Siesta RamanCalculatorInterface Calculator Class

class ase.calculators.siesta.siesta_lrtddft.RamanCalculatorInterface(omega=0.0, **kw)[source]

Raman interface for Siesta calculator. When using the Raman calculator, please cite

M. Walter and M. Moseler, Ab Initio Wavelength-Dependent Raman Spectra: Placzek Approximation and Beyond, J. Chem. Theory Comput. 2020, 16, 1, 576–586

Parameters:
  • omega (float) – frequency at which the Raman intensity should be computed, in eV

  • kw (dictionary) – The parameter for the siesta_lrtddft object