A First Encounter with SIESTA¶

This first tutorial provides a first acquaintance with SIESTA by studying one simple molecule, CH4. It covers most of the basic features and concepts, without worrying too much about issues of convergence or fine-tuning.

Take me to the exercises!

background Background Information¶

The inputs of SIESTA¶

The essential inputs of SIESTA are the FDF (Flexible Data Format) file and the pseudopotential file.

FDF file¶

It is the core input, since it provides information about the system and sets the values of several parameters regarding the quality of the calculation, as well as the pseudopotential information.

The FDF file (input.fdf) is a plain-text input file. It is structured in name-value pairs—that is, a name characterizing the data and its value. Units can be specified after the value. These name-value pairs can be part of more complex structures called blocks (%block label [...] %endblock). The FDF file is case-insensitive, and characters -, _, and . in a data label are ignored.

The list of available descriptors is extensive; however, most of them have default values. To describe the system for a calculation, the minimum information needed is provided by the following descriptors:

General system information: SystemLabel, NumberOfSpecies, NumberOfAtoms
Structural information: LatticeConstant, LatticeVectors, AtomicCoordinatesFormat, AtomicCoordinatesAndAtomicSpecies.

As a part of the FDF file we can find the information regarding the basis set used. SIESTA offers many options for the specification of the basis set. Two of the key concepts are the cardinality of the set and range of the orbitals; both can be specified, in most cases, simply with a couple of lines.

The cardinality of the set, or its size, can be set by PAO.BasisSize. It follows the usual convention of amount of functions per valence orbital, the possible options being SZ, SZP, DZ, DZP… For example, for C, a SZ basis set has 4 functions (s+p) and a DZP basis set has 13 functions (2*s, 2*p, and a d polarization shell).
The energy shift is a physical concept determining the increase in energy brought about by the confinement of the orbital; the lower the energy shift, the larger the cut-off radii of basis orbitals. It is set using PAO.EnergyShift.

As in any DFT code, the exchange-correlation functional must be specified; the default option is using the local density approximation (LDA) for our calculations. However, it is also possible to use other functionals, such as those of GGA type. The information of the exchange-correlation functional should be provided like this:

XC.functional GGA
XC.authors    PBE

Notice that both XC.Authors and XC.Functional are needed, and they must be consistent, check the manual. For XC.Functional GGA, the most used variants are BLYP, PBE and PBEsol. For XC.Functional VDW, frequently used functionals are DRSLL, LMKLL and VV. See the manual for more options on XC functionals.

Pseudopotentials¶

The other essential input of SIESTA is the pseudopotential information. SIESTA uses pseudopotentials to represent the core electron-ion interaction, as do most plane-wave codes (in contrast to so-called “all-electron” programs). In particular, the pseudopotentials are of the “norm-conserving” kind. The choice of a pseudopotential significantly impacts the quality of the calculation.

SIESTA reads pseudopotentials from different files based on the species information in %block ChemicalSpeciesLabel. By default, files are searched for in the current directory. SIESTA accepts different formats, but PSML is the one recommended. Curated databases of high-quality PSML files are available. In particular, the Pseudo-Dojo project offers PSML files for almost the entire periodic table, along with reports detailing the tests performed during their generation. The database includes variants for scalar and fully relativistic pseudopotentials using LDA, PBE, and PBEsol functionals.

Generating a new pseudopotential from scratch is beyond the scope of this tutorial, but more details can be found in Pseudopotentials. Note that ATOM generates psf files, which are also supported by SIESTA.

The pseudopotential should be tested. The method for testing depends on the type of system being simulated, so reviewing existing literature is strongly recommended.

The output¶

SIESTA writes a log of its workings to standard output, which is recommended to be redirected to an output file. In the standard output, the entire thinking process of SIESTA can be followed: starting from reading the input parameters, defining the system, and proceeding to the self-consistent-field convergence.

System properties can be found in it; however, not all of them are written by default. More detailed information will be given in the exercises: see below.

basic Look at the input files¶

You will find an input file named ch4.fdf, along with the files C.psml and H.psml which contain the information about the pseudopotentials.

Pseudos: C.psml, H.psml
Input file: ch4.fdf

We find first the descriptors that specify the system

SystemName          CH4 molecule
SystemLabel         ch4
NumberOfAtoms       5
NumberOfSpecies     2

%block ChemicalSpeciesLabel
 1  6 C   # Species index, atomic number, species label
 2  1 H   # Species index, atomic number, species label
%endblock ChemicalSpeciesLabel

Note

You will notice several lines starting with #. These are comments on the file and will not be read by SIESTA.

The SystemName is only used for tagging purposes, while the SystemLabel will determine the root name for all output files (i.e., they will all begin by ch4.* in our case). NumberOfAtoms indicates the total number of atoms in our simulation box, and NumberOfSpecies indicates the different kinds of atoms present in the system.

Pay special attention to the %block ChemicalSpeciesLabel. In this block, you assign an index, an atomic number, and a label to each atomic species. The atomic number is, as usual, the nuclear charge of the atom. The label will allow SIESTA to recognize the files containing the information about the pseudopotential and the basis set (when provided); this means that the name of the label must coincide with the name of the pseudopotential file. The amount of species in this block must be the one specified by the NumberOfSpecies input.

Note

Options in the fdf file are case-insensitive. This means that NumberOfAtoms is the same as numberofatoms or NUMBEROFATOMS.

A few lines below, we find the coordinates for the atoms in the system

AtomicCoordinatesFormat  Ang

%block AtomicCoordinatesAndAtomicSpecies
  0.000      0.000      0.000   1
  1.219     -0.284     -0.377   2
 -0.284      1.219     -0.377   2
 -0.140     -0.140      1.219   2
 -0.833     -0.833     -0.503   2
%endblock AtomicCoordinatesAndAtomicSpecies

In this block, each of the lines corresponds to an atom. The first three numbers are the X, Y and Z coordinates for each atom, and the last (integer) number indicates the species index in the %block ChemicalSpeciesLabel. In this case, there is a Carbon atom (species index = 1) in the position (0.0,0.0,0.0), and the rest are the four Hydrogens (species index = 2).

The AtomicCoordinatesFormat option can also be set to Bohr (which is the default value), or to Fractional if one desires to set the coordinates as a fraction of the lattice vectors (see below).

You might have wondered about the existence of this block in the input file:

#Unit cell for the calculation
LatticeConstant 15 Ang
%block LatticeVectors
  1.000 0.000 0.000
  0.000 1.000 0.000
  0.000 0.000 1.000
%endblock LatticeVectors

It does not seem to make sense for a ‘molecule’ calculation. In fact, SIESTA always uses periodic boundary conditions (PBC); this means that we are doing a calculation for an infinite collection of regularly spaced molecules. If we want to simulate an isolated molecule, it is important to have enough distance between the molecule and its neighboring images. At the very minimum, there should not be overlap between the orbitals on different image molecules; however, this size should be bigger for molecules with a dipole, which will have a long-range interaction with their images in PBC.

The input option LatticeConstant is mandatory and multiplies the lattice vectors by that constant. In this case, we will end up with a cubic simulation box of 15 x 15 x 15 Angstrom.

In the rest of the file you will find options regarding the quality of the calculation and other performance related variables:

# Basis set definition
PAO.EnergyShift 250 meV
PAO.SplitNorm   0.15
PAO.BasisSize   SZ

#Real space grid
MeshCutoff 125.0 Ry

# Convergence of SCF
MaxSCFIterations 50
DM.MixingWeight 0.4
DM.NumberPulay  2

# Type of solution
SolutionMethod diagon

The specifics of these will be explored in other tutorials; but to provide a brief overview:

The first three variables are a rough definition the quality of the basis set. The number of orbitals per atom is defined by the parameter PAO.BasisSize, and we have set it to select a minimal basis (SZ) for a quick, cheap calculation disregarding quality.
The parameter MeshCutoff controls the fineness of the real-space grid used to compute the integrals for the matrix elements of the Hamiltonian. A higher value corresponds to a finer grid.
The next three options control the performance of the self-consistent field (SCF) cycle: MaxSCFIterations, DM.MixingWeight, DM.NumberPulay.
The kind of solver used is set by SolutionMethod.

While the parameters specifying the system are mandatory, all other parameters have some default values and, in principle, it is not necessary to explicitly include them in the input file. However, it is important to note that the default values do not always guarantee a fast or a good quality calculation.

basic Running SIESTA¶

To run SIESTA via command line, you just need to do

$ siesta < ch4.fdf > ch4.out

or simply

$ siesta ch4.fdf > ch4.out

You will see that a lot of output files are generated in the simulation directory. We will begin by having a look at the main output file, ch4.out (or any name that you gave to the standard output).

basic A look at the output¶

Take a look at ch4.out once the program has finished. The file contains a human-readable account of the doings of SIESTA for this calculation:

A header with the version and compilation options, and a copy (in the first mode above) or a mention (in the second mode) of the input file.

A summary of existing species and pseudopotential information.

Details on the basis set generation.

A summary of the values used in the calculation for the most important parameters.

Information of the geometry and matrix sparsity for a given geometry.

Followed by a full breakdown of the initial, non-converged energy.

A log of the self-consistent-field (SCF) cycle.

A full breakdown of the total energy decomposition and forces, and stress.

The last part of the output file includes information on the cell pressure and the total electric dipole of the system.

And to finish, literature regarding the simulation done.

As you can see, SIESTA generates a ton of output files besides the main log/out file. Each of these might interest you depending on what you intend to do with your simulation, but some of the most commonly relevant files are:

*.DM files (ch4.DM) contain a restart of the Density Matrix. Useful to restart already converged or close-to-convergence calculations without doing so from scratch.
*.XV files (ch4.XV) contain a restart of the current geometry and cell size. This is particularly useful to restart or continue molecular dynamics runs.
*.FA files (ch4.FA) contain total forces over all atoms. Each line corresponds to the cartesian force (x, y and z) over each atom (in eV/Ang).
*.EIG files (ch4.EIG) contain the eigenvalues of the Hamiltonian for all k-points involved. k-points are explored in more detail in Sampling of the BZ with k-points.
*.times and TIMES files (ch4.times and TIMES) contain timing information, which is invaluable when optimizing parameters for fast performance.

Now that you had a brief overview of how running SIESTA works, let’s play a bit with your CH4 inputs. See what happens if…:

You change the fineness of the real-space grid. Try an unreasonably low value for the the parameter MeshCutoff (around 10 Ry) and check the resulting total energy and forces (you can also find the forces in the file ch4.FA). Then experiment with a much higher one (300 Ry). More on this topic will be covered in this tutorial.
Answer: Change of MeshCutoff.
Everytime that you run SIESTA in a folder, all the output will be overwritten, thus to be able to compare the results of the diferent runs, a specific directory should be created for each of them.

The information or the total energy can be extracted in a more automatic way doing a grep of Total = of the standard output file and to extract the total force look for the line that contains siesta: Tot; for example if you have structured the calculation in different folders:
$ grep -i  "siesta:    Tot" */output.dat
10ry/output.dat:siesta:    Tot    1.798364    1.798364   -3.227789
300ry/output.dat:siesta:    Tot   -0.000058   -0.000058   -0.000291
$ grep -i  "siesta:    Total =" */output.dat
10ry/output.dat:siesta:         Total =    -230.840670
300ry/output.dat:siesta:         Total =    -222.538059
You can play with the size of the lattice parameters to go from ‘interacting’ molecules to effectively isolated ones. Look at the variation in the total energy as a function of the cell size, to see how the interaction between molecules decreases with increasing distance between images. For this non-polar molecule, the interaction should be very small.
Answer: Change of lattice parameter.
Everytime that you run SIESTA in a folder, all the output will be overwritten, thus to be able to compare the results of the diferent runs, a specific directory should be created for each of them.

The information or the total energy can be extracted in a more automatic way doing a grep of Total = of the standard output file and to extract the total force look for the line that contains siesta: Tot. For example, if you have structured the calculation in different folders:
$ grep -i  "Total =" */output.dat
15AA/output.dat:siesta:         Total =    -222.538221
35AA/output.dat:siesta:         Total =    -222.538164
7_5AA/output.dat:siesta:         Total =    -222.538303

basic Simple basis-sets options¶

Now let’s change some options related to basis-sets.

You can find the input files here:

Pseudos psml: C.psml and H.psml
Input file: ch4.fdf

Change basis-set size or cardinality of the basis. Set PAO.BasisSize to any of: DZ, SZP, DZP, TZP. See what happens with the total energy.

Let’s look at the effect of changing the basis set in some parameters with physical meaning, like the bond length. In this tutorial you will be doing a geometry optimization for CH4.

The input file ch4.fdf is already set up to optimize the geometry of a single methane molecule, but until now it was comemented. Let’s uncomment the corresponding lines by removing # (line 39 and 40).

# Variables for Geometry optimization
# MD.TypeOfRun   CG
# MD.NumCGsteps  50

At the bottom of the file, you will find two options used to define the basis set:

PAO.BasisSize   SZ
PAO.EnergyShift 250 meV

The idea of this exercise is to see what happens with the total energy and the bond lengths of the system when modifying these options. The bond lengths can be found in the ch4.BONDS_FINAL file.

The SystemLabel.BONDS_FINAL contains the information regarding the atomic position of each atom of the unit cell and their distance with respect to each of the other ones.

See what happens when you increase the size of the basis set to SZP and DZP.

For DZP, try changing the energy shift to 100 meV, 50 meV and 10 meV.

Consider that the experimental bond-length for gas methane is 1.087 Å.

Which has the biggest effect in the energy, increasing the cardinality/size of the basis set, or changing the energy shift?

Answer: Effect in the energy.

From the first experiment, we can see that the cardinality of the basis set induces a change in the range of a few eV, which is much larger than that induced by the energy shift. Be aware that this happens when the basis set is relatively small; for basis sets larger than DZP (such as TZP) you don’t usually get a huge improvement. However, going to SZ to SZP, or from SZP to DZP, the change in flexibility of the basis set is a determining factor in the quality of the results.
In which case would you consider the basis set is “good enough”?

Answer: How do we select the good one?

We definitely need the DZP basis set to start getting reasonable results. But with regards to the energy shift, we can see that 50 meV is already reasonable enough to get good results. The total energy seems to be already converged to the meV level, and the bond lengths are also reasonable. That said, increasing the energy shift is not detrimental in this case, it will just increase the time of the calculation.
Is it ok to compare the results to an experimental bond length? Or would it be better to check against an already converged basis set?

Answer: Is the comparison with the experimental possible?

The question is tricky and depends on what you are going after. If you are testing convergence of parameters, then by all means, compare to an already converged case, even if it is with another type of basis set (such as gaussians or plane-waves). That said, if you are already sure that your input parameters are good enough, then comparing to experimental values is a good idea: it will give you a hint on how your level of theory (i.e., DFT with a given functional) is able to model the real world data. In most cases, the best approach is to compare trends between experiments and theory, rather than absolute values.

basic DFT functional¶

You can find the input files here:

Pseudos LDA: C.psml and H.psml
Pseudos GGA: C.gga.psml and H.gga.psml
Input file: ch4.fdf

Up to now we have been implicitly using the local density approximation (LDA) for our calculations, which is the default option in SIESTA. However, it is also possible to use other functionals, such as those of GGA type. Edit the ch4.fdf file to include this block:

XC.functional GGA
XC.authors    PBE

Notice that both XC.Authors and XC.Functional are needed, and they must be consistent. For XC.Functional GGA, the most used variants are BLYP, PBE and PBEsol. For XC.Functional VDW, frequently used functionals are DRSLL, LMKLL and VV. See the manual for more options on XC functionals.

Run SIESTA and look for possible lines with WARNING or ERROR in them.

Modify the input file again to use these files by simply changing the species strings:

%block ChemicalSpeciesLabel
1  6 C.gga   # Species index, atomic number, species label
2  1 H.gga   # Species index, atomic number, species label
%endblock ChemicalSpeciesLabel

Remember, the species name can be anything you want, but it must be consistent with the name of the pseudopotential files. Run the program again and check whether the warning disappears from the output.

It is worth noting that, as with pseudopotentials, DFT functionals are not easy to compare between each other systematically. Again, it is worth expending some time searching the literature for the type of simulations you are trying to run.