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 SIESTA_COMMAND must hold the command to invoke the siesta calculation. The variable must be a python format string with exactly two string fields for the input and output files. Examples: siesta < %s > %s, mpirun -np 4 /bin/siesta3.2 < %s > %s.

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 SIESTA_COMMAND="siesta < ./%s > ./%s"
$ export SIESTA_PP_PATH=$HOME/mypps

Alternatively, the path to the pseudopotentials can be given in the calculator initialization.

keyword	type	default value	description
`pseudo_path`	`str`	`None`	Directory for pseudopotentials to use None means using $SIESTA_PP_PATH

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`	`float`	`COLLINEAR`	The spin approximation used, must be either `UNPOLARIZED`, `COLLINEAR` or `FULL`
`species`	`list`	`[]`	A method for specifying a specific description for some atoms.
`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, i.e. `H.lda.psf`

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.set_calculator(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 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.set_calculator(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.

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 considered relative to the default 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.

TDDFT Calculations¶

It is possible to run Time Dependent Density Functional Theory (TDDFT) using the PYSCF-NAO code together with the SIESTA code. This code allows to run TDDFT up to thousand atoms with small computational ressources. Visit the github webpage for further informations about PYSCF-NAO.

Example of code to calculate polarizability of Na8 cluster,:

from ase.units import Ry, eV, Ha
from ase.calculators.siesta import Siesta
from ase import Atoms
import numpy as np
import matplotlib.pyplot as plt

# Define the systems
Na8 = Atoms('Na8',
             positions=[[-1.90503810, 1.56107288, 0.00000000],
                        [1.90503810, 1.56107288, 0.00000000],
                        [1.90503810, -1.56107288, 0.00000000],
                        [-1.90503810, -1.56107288, 0.00000000],
                        [0.00000000, 0.00000000, 2.08495836],
                        [0.00000000, 0.00000000, -2.08495836],
                        [0.00000000, 3.22798122, 2.08495836],
                        [0.00000000, 3.22798122, -2.08495836]],
             cell=[20, 20, 20])

# Siesta input
siesta = Siesta(
            mesh_cutoff=150 * Ry,
            basis_set='DZP',
            pseudo_qualifier='',
            energy_shift=(10 * 10**-3) * eV,
            fdf_arguments={
                'SCFMustConverge': False,
                'COOP.Write': True,
                'WriteDenchar': True,
                'PAO.BasisType': 'split',
                'DM.Tolerance': 1e-4,
                'DM.MixingWeight': 0.01,
                'MaxSCFIterations': 300,
                'DM.NumberPulay': 4,
                'XML.Write': True})

Na8.set_calculator(siesta)
e = Na8.get_potential_energy()
freq, pol = siesta.get_polarizability_pyscf_inter(label="siesta",
                                                  jcutoff=7,
                                                  iter_broadening=0.15/Ha,
                                                  xc_code='LDA,PZ',
                                                  tol_loc=1e-6,
                                                  tol_biloc=1e-7,
                                                  freq = np.arange(0.0, 5.0, 0.05))
# plot polarizability
plt.plot(freq, pol[:, 0, 0].imag)
plt.show()

Remark:¶

The PYSCF-NAO code is still under active development and to have access to it with ASE you will need to use this PYSCF fork and use the branch nao. To summarize:

git clone https://github.com/cfm-mpc/pyscf
git fetch
git checkout nao

Then you can follow the instruction of the README. The installation is relatively easy, go to the lib directory:

cd pyscf/pyscf/lib
cp cmake_arch_config/cmake.arch.inc-your-config cmake.arch.inc
mkdir build
cd build
cmake ..
make

Then you need to add the pyscf directory to your PYTHONPATH

export PYTHONPATH=/PATH-TO-PYSCF/pyscf:$PYTHONPATH

Raman Calculations with SIESTA and PYSCF-NAO¶

It is possible to calulate the Raman spectra with SIESTA, PYSCF-NAO anf the vibration module from ASE. Example with CO2,:

from ase.units import Ry, eV, Ha
from ase.calculators.siesta import Siesta
from ase.calculators.siesta.siesta_raman import SiestaRaman
from ase import Atoms
import numpy as np

# Define the systems
# example of Raman calculation for CO2 molecule,
# comparison with QE calculation can be done from
# https://github.com/maxhutch/quantum-espresso/blob/master/PHonon/examples/example15/README

CO2 = Atoms('CO2',
            positions=[[-0.009026, -0.020241, 0.026760],
                       [1.167544, 0.012723, 0.071808],
                       [-1.185592, -0.053316, -0.017945]],
            cell=[20, 20, 20])

# enter siesta input
# To perform good vibrational calculations it is strongly advised
# to relax correctly the molecule geometry before to actually run the
# calculations. Then to use a large mesh_cutoff and to have the option
# PAO.SoftDefault turned on
siesta = Siesta(
    mesh_cutoff=450 * Ry,
    basis_set='DZP',
    xc="GGA",
    pseudo_qualifier='gga',
    energy_shift=(10 * 10**-3) * eV,
    fdf_arguments={
        'SCFMustConverge': False,
        'COOP.Write': True,
        'WriteDenchar': True,
        'PAO.BasisType': 'split',
        "PAO.SoftDefault": True,
        'DM.Tolerance': 1e-4,
        'DM.MixingWeight': 0.01,
        'MaxSCFIterations': 300,
        'DM.NumberPulay': 4,
        'XML.Write': True,
        'DM.UseSaveDM': True})

CO2.set_calculator(siesta)

ram = SiestaRaman(CO2, siesta, nfree=4, label="siesta", jcutoff=7, iter_broadening=0.15/Ha,
        xc_code='LDA,PZ', tol_loc=1e-6, tol_biloc=1e-7, freq = np.arange(0.0, 5.0, 0.05))

ram.run()
ram.summary(intensity_unit_ram='A^4 amu^-1')
ram.write_spectra(start=200, intensity_unit_ram='A^4 amu^-1')

Further Examples¶

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

Siesta Calculator Class¶

class ase.calculators.siesta.base_siesta.BaseSiesta(**kwargs)[source]¶

Calculator interface to the SIESTA code.

ASE interface to the SIESTA code.

Parameters:

label : The base head of all created files.
mesh_cutoff : Energy in eV.

The mesh cutoff energy for determining number of grid points.
energy_shift : Energy in eVV

The confining energy of the basis sets.
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”, strings which specify the

type of functions basis set.
spin : “UNPOLARIZED”|”COLLINEAR”|”FULL”. 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.
atoms : The Atoms object.
restart : str. Prefix for restart file.

May contain a directory. Default is None, don’t restart.
siesta_default: Use siesta default parameter if the parameter

is not explicitly set.
ignore_bad_restart_file: bool.

Ignore broken or missing restart file. By default, it is an error if the restart file is missing or broken.
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.