02 – First simulation with SIESTA

Example 01 (Water)

Download the pseudopotential files H.gga.psf and O.gga.psf in https://departments.icmab.es/leem/siesta/Databases/Pseudopotentials/periodictable-gga-abinit.html:

$ cd $HOME/siesta/pseudos
$ wget https://departments.icmab.es/leem/siesta/Databases/Pseudopotentials/Pseudos_GGA_Abinit/H_html/H.psf
$ wget https://departments.icmab.es/leem/siesta/Databases/Pseudopotentials/Pseudos_GGA_Abinit/O_html/O.psf

Copy the input.fdf file with the git clone command:

$ cd $HOME/siesta/IO
$ git clone https://github.com/QuantumMackGraphe/siesta-tutorials.git
$ cd siesta-tutorials/example01

Now, type the following command and cross your fingers.

$ siesta-omp input.fdf > job.out &

if all goes well, you will see a message “Job completed” in the last line of your output file (job.out). The “siesta input.fdf” command execute the SIESTA program reading the input file (input.fdf). The following character “>” indicates to write the output data to “job.out”. The last character “&” is used to leave the job in the background until the end of execution.

Understanding your input file

The input file (input.fdf) is a plain text file with many keys and values. To simplify our input file, we will group some keys and values in “blocks”. These blocks will be separated by comments. Comments are texts in the input file following the hash character (#). These texts will be ignored by the SIESTA program, and will be introduced only to facilitate understanding of the input file. For example, in our input.fdf file in example01, we use keys and values grouped in two blocks: Self-consistent field and Material.

The first block is given by:

# in file: input.fdf

         # -- SELF-CONSISTENT FIELD --

MeshCutoff                   300 Ry
PAO.EnergyShift              30 meV
XC.Functional                GGA
XC.Authors                   PBE
MaxSCFIterations             1000
  • The MeshCutoff key is used to control the real-space grid for charge density. The energy value is used to define the plane wave cutoff for the grid. From the version 4.1, the default value is 300 Ry. Increase this value to make the grid finer, but keep in mind that this increases the computational cost considerably. We will optimized the value in the following tutorial section;
  • The PAO.EnergyShift key is used to control the basis function range. The SIESTA program is based on linear combination of atomic orbitals (LCAO). However, the atomic orbitals are actually numerical orbitals confined in a real-space range. This confinement shift the energy levels as an infinite spherical potential well. Smaller PAO.EnergyShift indicate far-reaching, and hence better basis functions. The value will also be optimized in the following tutorial section;
  • The XC.Functional and XC.Authors are two keys to choose the exchange-correlation (XC) functional in the Kohn-Sham equations. The XC functionals more popular by now is the Generalized Gradient Approximation (GGA) parameterized by Perdew, Burke and Ernzerhof (PBE) [Phys. Rev. Lett. 77, 3865 (2016)];
  • The MaxSCFIterations key is just the maximum number of iterations in Kohn-Sham self-consistent field to converge the density matrix and hamiltonian. The default value is 50, but this may be too low. We will increase to a high value as 1000 in our simulations.

The second block in our input file is given by:

# in file: input.fdf
        # -- MATERIAL --

SystemLabel                  H2O
WriteCoorStep                true
WriteCoorXmol                true

%block ChemicalSpeciesLabel
 1   8  O.gga
 2   1  H.gga
%endblock ChemicalSpeciesLabel

AtomicCoordinatesFormat      Ang

LatticeConstant              1.000 Ang

%block LatticeVectors
       20.0000000000         0.0000000000         0.0000000000
        0.0000000000        20.0000000000         0.0000000000
        0.0000000000         0.0000000000        20.0000000000
%endblock LatticeVectors

%block AtomicCoordinatesAndAtomicSpecies
   10.05053122   10.05053123    9.99999999   1
   11.01902700    9.93038215   10.00000000   2
    9.93038215   11.01902699   10.00000000   2
%endblock AtomicCoordinatesAndAtomicSpecies
  • SystemLabel is a string key to define the filename of some output files. Use something simple with no spaces and special characters;
  • WriteCoorStep and WriteCoorXmol are two keys to write the atomic coordinates in the main output (WriteCoorStep) and in a xyz file (WriteCoorXmol);
  • ChemicalSpeciesLabel is key-value in a block format. Each line of this block receives an index for a chemical species, the atomic number, and the label for the pseudopotential file (label.psf);
  • AtomicCoordinatesFormat is the key to define the format of the atomic coordinates. We usually use Ang or ScaledByLatticeVectors. For molecules, the former is more intuitive;
  • LatticeConstant is the key that multiplies the lattice vectors;
  • LatticeVectors are keys to define a unit cell that repeats periodically in space. For molecules, these lattice vectors are just large vectors that define a large box;
  • AtomicCoordinatesAndAtomicSpecies is another key-value in a block format that defines the positions and chemical species of each atom. Each row in this block is an atom for the unit cell.

More details of the keys and values that can be used in the input file can be found in the SIESTA manual.

Checking your output file

To show the contents of the output file, use the cat command:

$ cat job.out

The output file usually are very long, with a lot of information about the simulation. For now, let’s focus on self-consistent field (scf) part. For each iteration, the Kohn-Sham energy, Free energy, and Fermi level are showed. The self-consistent convergence is controlled by the difference density matrix (dDmax) and difference of Hamiltonian (dHmax).

        iscf     Eharris(eV)        E_KS(eV)     FreeEng(eV)     dDmax    Ef(eV) dHmax(eV)
   scf:    1     -474.831923     -476.892759     -476.892759  0.902382 -2.734301 11.509096
timer: Routine,Calls,Time,% = IterSCF        1       3.753  53.65
   scf:    2     -480.167452     -478.826398     -478.826398  0.176264 -7.267793  4.413552
   scf:    3     -479.128249     -479.090327     -479.090327  0.075113 -5.951575  0.356110
   scf:    4     -479.100644     -479.096260     -479.096260  0.014051 -5.916251  0.292365
   scf:    5     -479.097704     -479.097035     -479.097035  0.004317 -5.896236  0.268409
   scf:    6     -479.098074     -479.097594     -479.097594  0.006155 -5.839471  0.213176
   scf:    7     -479.097937     -479.097845     -479.097845  0.017635 -5.612172  0.050964
   scf:    8     -479.098079     -479.097983     -479.097983  0.003525 -5.638938  0.018375
   scf:    9     -479.097998     -479.097994     -479.097994  0.001679 -5.655289  0.002709
   scf:   10     -479.097994     -479.097994     -479.097994  0.000099 -5.655716  0.000802

Note that this scf is considered converged if dDmax < 1.0E-4 and dHmax < 1.0E-3 eV. In our example, the scf converged after 10 iterations.

Another important part, is the final energy decomposition, as shown below:

siesta: Final energy (eV):
siesta:  Band Struct. =    -109.265444
siesta:       Kinetic =     329.160653
siesta:       Hartree =     512.246332
siesta:       Eldau   =       0.000000
siesta:       Eso     =       0.000000
siesta:    Ext. field =       0.000000
siesta:   Exch.-corr. =    -128.943101
siesta:  Ion-electron =   -1312.809437
siesta:       Ion-ion =     121.247559
siesta:       Ekinion =       0.000000
siesta:         Total =    -479.097994
siesta:         Fermi =      -5.655716

This part include the final total energy, and this value can be used to calculate various other physical properties. Here. the total energy is -479.097994 eV. However, we must remember that this value is not absolute because we are using a pseudopotential approach. Only differences in total energies will be compared with experimental data.

A quick way to see the total energy in the output file is by using the grep command:

$ grep "Total =" job.out

This command filters a text file so that only the lines containing a specific string appear. In our case, filter the job.out file with lines containing the string “Total =”.

Previous: 01 – Installing SIESTA

Next: 03 – Geometry optimization