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 =”.