:sequential_nav: next .. _tutorial-basic-first-encounter-theorylevel: A First Encounter - Part 2: Choosing your level of theory ********************************************************* .. sidebar:: **Have you set up the local environment?** If not, :ref:`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 :ref:`tutorial-basic-first-encounter`, 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 :ref:`tutorial-basic-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: * Part I: https://www.youtube.com/watch?v=puweAzW9ZNs * Part II: https://www.youtube.com/watch?v=AlY2SDb1ta0 .. note:: This tutorial will only cover basic, general options for basis set generation. We encourage you to jump into the dedicated tutorial for :ref:`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 :ref:`tutorial-basic-structure-optimization`. .. 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? .. Changing basis: SZ : Total: -215.244568 eV, Bonds: 1.1809/1.1809/1.1816/1.1820 SZP: Total: -215.797752 eV, Bonds: 1.1506/1.1509/1.1509/1.1512 DZP: Total: -218.068475 eV, Bonds: 1.1040/1.1044/1.1044/1.1058 Changing eshift for DZP: 100: Total: -218.162319 eV, Bonds: 1.1126/1.1126/1.1128/1.1130 50 : Total: -218.114205 eV, Bonds: 1.1161/1.1161/1.1163/1.1164 10 : Total: -218.006615 eV, Bonds: 1.1191/1.1191/1.1193/1.1194 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.