A First Encounter - Part 2: Choosing your level of theory

Have you set up the local environment?

If not, do that now before proceeding.

Choosing a Pseudopotential

An important thing that determines the quality of the calculation is the choice of the pseudopotential. As mentioned in A First Encounter - Part 1: Running SIESTA, pseudopotentials in SIESTA are provided as external input files, usually in psf or psml formats.

Generating a new pseudopotential from scratch is outside the scope of this tutorial, but more information can be seen in Pseudopotentials. For most cases, you can instead download pre-generated pseudopotentials at Pseudo-Dojo, which covers variants for scalar and fully relativistic pseudos for LDA, PBE and PBEsol functionals.

The pseudopotential should be tested. How to test it will depend on the type of system you are running, so diving in to the existing literature is strongly recommended.

DFT functional

Hint

Enter the directory CH4-XC-Functional

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. Even though the program ran successfully, you will find the following warning in the output:

xc_check: WARNING: Pseudopotential generated with LDA CA functional

You are using a GGA functional with a pseudopotential generated using LDA, as this is not consistent!. Fortunately, we have produced also the pseudopotentials using GGA for you. They are in the files C.gga.psf and H.gga.psf. You can 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

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 in a systematic way. Again, it is worth expending some time searching the literature for the type of simulations you are trying to run.

Simple basis-sets options

Note

As a general background for this topic, you can watch the 2021 talks by Emilio Artacho:

Note

This tutorial will only cover basic, general options for basis set generation. We encourage you to jump into the dedicated tutorial for basis set optimization afterwards.

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, which can be specified, in most cases, simply with a couple of lines. For example:

PAO.BasisSize   DZP
PAO.EnergyShift 0.01 Ry

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.

Other interesting options are:

PAO.SplitNorm       0.15
PAO.SoftDefault     T
PAO.OldStylePolOrbs F

PAO.SplitNorm esentially controls the matching radii for higher-zeta orbitals. PAO.SoftDefault enables the use of soft-confinement potentials for basis set generation (ensuring the derivatives of the radial functions become zero at the orbital cut-off radii), and PAO.OldStylePolOrbs can be used as a fallback to an older, different way of generating polarization orbitals.

Note

The values for PAO.SoftDefault and PAO.OldStylePolOrbs are set to true and false by default starting in SIESTA 5. These are the recommended values; however, previous versions of SIESTA used other defaults.

Playing with CH4

Note

In this tutorial you will be doing a geometry optimization for CH4. This method is covered in Structural optimization using forces and stresses.

Hint

Enter the directory CH4-Basis

The input file ch4.fdf is already set up to optimize the geometry of a single methane molecule. 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 tutorials is to see what happens with the total energy and the bond lengths of the system when modifying these options. The total energy is found in the standard output, look for siesta:         Total = near the bottom of the file, under siesta: Final energy (eV):. The bond lengths can be found in the ch4.BONDS_FINAL file.

  • 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?

  • In which case would you consider the basis set is “good enough”?

  • Is it ok to compare the results to an experimental bond lenght? Or would it be better to check against an already converged basis set?

Optional

Here are a few extra things for you to explore:

  • For DZP at a given energy shift, try different values for PAO.SplitNorm (0.10, 0.15 and 0.20). See how that affects your results.

  • See what happens when you disable PAO.SoftDefault. In particular, you could have a look at the cut-off radii for basis orbitals, which are present in the standard output (search for rc =). Try again different values for the energy shift and see if your previous conclusions change.