:sequential_nav: next .. _tutorial-basic-basis-extra: Basis sets - Special cases ========================== .. dropdown:: **Do you want to get all tutorial materials on your local machine?** :color: primary You can get all files and set up your local environment. :ref:`Do that now ` before proceeding. In this tutorial we will cover a few additional considerations on basis sets when using SIESTA, particularly in the context of the *Basis Set Superposition Error* and properties that require a better description of the vacuum region. If you have not done so already, we recommend that you first read the :ref:`tutorial on basis set optimization `. This optimization usually comes before considering what we will cover here. .. dropdown:: Additional information :color: info 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 .. button-ref:: basis-sets-ghosts :color: success :outline: :expand: :shadow: Take me to the exercises! :bdg-info:`background` Staring at the void: the vacuum region ------------------------------------------------------------- Contrary to plane-wave descriptions, atomic-centered basis sets do not fill the entire space of your simulation. This is key to their efficiency, but it also means that the description of the vacuum region might be lacking in some cases. The two tricks we will be covering here (diffuse orbitals and ghost atoms) are meant to solve this issue with little effort. Diffuse orbitals ^^^^^^^^^^^^^^^^ The term *diffuse orbitals* refers to basis functions that have a long cut-off radius. These orbitals are **added** after the basis set has been optimized, and we will use an arbitrarily large cut-off radius. So suppose we have a PAO.Basis block that looks like this for a carbon atom:: %block PAO.Basis C 2 n=2 0 2 0.000 0.000 1.000 1.000 n=2 1 2 P 1 0.000 0.000 1.000 1.000 %endblock PAO.Basis We have here two shells with n=2: a double-zeta s-type shell and a double-zeta p-type shell with a single polarization orbital. We can then add a diffuse orbital by adding a new shell with a higher n, for example:: %block PAO.Basis C 3 ! Note the increase 2->3 n=2 0 2 0.000 0.000 1.000 1.000 n=2 1 2 P 1 0.000 0.000 1.000 1.000 n=3 0 1 10.000 1.000 %endblock PAO.Basis Note that we used an arbitrary value for the cut-off radius, of 10 Bohr. Usually, values chosen range for 10 to 15 Bohr, depending on the system. If 10 is enough to get a given property converged, then you're good to go! A thing worth noting is that you can always have any amount of extra diffuse orbitals. Usually, you add diffuse orbitals until the property you want to calculate is converged enough, but a little chemical or physical intuition also helps. Basis Set Superposition Error ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ A concept that is key for atomic-centered basis in general is the *Basis Set Superposition Error* or *BSSE*. All kinds of DFT codes are subject to this error, so long as they use atom-centered basis sets. The general idea is that, given an atom A, the basis functions coming from neighbouring atoms B also help improve the description of the electronic density around A. If we remove those neighbours, we effectively end up with a poorer description. These are some examples in which BSSE is surely present: * Dimer energy calculations: we will cover this example below. Essentially, each single monomer molecule has a better numerical description in the dimer than alone. * Adsorption energy: similar to the dimer above. * Vacancy formation energy: by removing an entire atom from the system, we have effectively removed a place in space where functions from other atoms might not reach. The better the basis set, the less important this issue becomes. However, one way to completely remove this issue in some cases is by using *ghost atoms*. Ghost atoms ^^^^^^^^^^^ SIESTA can place basis orbitals anywhere, not just in the position of atoms. The main way to do this is by using *ghost atoms*. Ghost atoms are not real atoms, in the sense that they do not have a nuclear charge, and they do not add electrons to the system. They are just a convenient way to place basis functions in places where there are no atoms. In order to do this in SIESTA, we just add an atom and set their atomic number to a negative value: -8, for example, will place a ghost oxygen, while -6 will place a ghost carbon (this will be detailed later). For this, you will need to generate an extra atom species in the :ref:`ChemicalSpeciesLabel` block. Ghost atoms can be useful in a variety of cases, for example: * When fixing the BSSE in dimer, binding or vacancy calculations. * To improve the treatment for the electronic density in slabs, increasing the quality of the description in the near-vacuum region. Both of these cases will be considered below. .. _basis-sets-ghosts: :bdg-success:`basic` Fixing the Basis Set Superposition Error ------------------------------------------------------------- .. dropdown:: **Are you under a local environment?** :color: primary Enter the directory *01-WaterDimer-BSSE* .. figure:: images/dimerghost.png :width: 300px :align: center :alt: (A water molecule and a ghost water molecule.) A water molecule (right) with another ghost water molecule (left) to represent what we have in the dimer. If you ran the tutorial on :ref:`basis set optimization `, you might have noticed that the binding energy can be twice as large as the planewaves result; at least, when choosing the default settings for SIESTA. Even for our optimized basis set, we might be up to 50 meV away from it. Is this really the best result we can get? It just so happens that in our case, the dimer, each of the water molecules is slightly better described if we take in to account the basis functions of the other water molecule. Thus, it is as if our monomer had a poorer basis set than the dimer, and we are in a case in which the BSSE is important. Luckily enough, there are a few ways so solve this issue, the most common one involving ghost atoms. As mentioned before, ghost atoms are just points in space in which we add basis functions that would belong to a given atom if it were there. In our case, we will turn each of the water molecules in a ghost ending up with the following systems: * Water\ :sub:`1`\ + (Ghost)Water\ :sub:`2`\ * (Ghost)Water\ :sub:`1`\ + Water\ :sub:`2`\ We'll now have a look at the input files: * Pseudos: :download:`O.psml`, :download:`H.psml`, :download:`O_ghost.psml`, :download:`H_ghost.psml` * Input files: :download:`monomer.fdf`, :download:`dimer.fdf`, :download:`ghost1.fdf`, :download:`ghost2.fdf` We will use four different input files: *monomer.fdf* has a single water molecule, while *dimer.fdf* has the water dimer. *ghost1.fdf* and *ghost2.fdf* have the geometries of the dimer, but with one of the water molecules as a set of ghost atoms. Take a look, for example, at the *ghost1.fdf* file, you will notice a few distinct things: .. literalinclude:: work-files/01-WaterDimer-BSSE/ghost1.fdf :linenos: :lineno-match: :start-at: %block ChemicalSpeciesLabel :end-at: %endblock AtomicCoordinatesAndAtomicSpecies Aside from the normal oxygen and hydrogen, do you see the **two additional species**? Note O_ghost and H_ghost, which have negative atomic numbers. This indicates that those species will be ghosts! If you look below in the coordinates, you will see that the first water molecule is the ghost (the atoms have species indexes 3 and 4). There are two more things we have to take into account now: first, the pseudopotential files. Since these ghosts are techincally a new species, we need to have pseudos for them. So what do we do? We just copy and rename the pseudopotential files for the normal atoms. See how *H.psml* and *H_ghost.psml* are the same file? .. dropdown:: Additional information: Reusing the same pseudopotential file. :color: info Instead of copying and pasting the pseudopotential file with a new name, you can re-use the same file by adding an extra term to the :ref:`ChemicalSpeciesLabel` block:: %block ChemicalSpeciesLabel 1 8 O 2 1 H 3 -8 O_ghost O.psml 4 -1 H_ghost H.psml %endblock ChemicalSpeciesLabel In this case, O_ghost will use the *O.psml* file, and H_ghost will use the *H.psml* file. Second, if we use the :ref:`PAO.Basis` block (which, we will), we also need to generate a new section for this new species! .. literalinclude:: work-files/01-WaterDimer-BSSE/ghost1.fdf :linenos: :lineno-match: :start-at: %block PAO.Basis :end-at: %endblock Now that you've seen the important parts of it, *What is the difference between* ghost1.fdf *and* ghost2.fdf *?* .. dropdown:: Answer: Difference between ghost1.fdf and ghost2.fdf. :color: success The only difference is in the block :ref:`AtomicCoordinatesAndAtomicSpecies`. For *ghost2.fdf*, what changes is the order of the species: .. literalinclude:: work-files/01-WaterDimer-BSSE/ghost1.fdf :linenos: :lineno-match: :start-at: %block AtomicCoordinatesAndAtomicSpecies :end-at: %endblock AtomicCoordinatesAndAtomicSpecies This means that the other water molecule is the ghost. Now for the last part of the puzzle, how do we calculate the binding energy when having ghost atoms? Remember that the original definition is: :math:`E_{binding} = E_{dimer} - 2*E_{monomer}` But now we have two possible escenarios for a monomer with ghosts: water 1 as a ghost with water 2 as the monomer (ghost1); or water 1 as the monomer, with water 2 as the ghost (ghost2). So now, the definition of the binding energy becomes: :math:`E_{binding} = E_{dimer} - E_{ghost1} - E_{ghost2}` Now we can run our siesta calculation as usual: .. code-block:: console $ siesta < monomer.fdf > monomer.out $ siesta < ghost1.fdf > ghost1.out $ siesta < ghost2.fdf > ghost2.out $ siesta < dimer.fdf > dimer.out And we can get the total energy for each of the systems: .. code-block:: console $ grep "Total =" *.out .. dropdown:: Output .. code-block:: dimer.out:siesta: Total = -965.509762 ghost1.out:siesta: Total = -482.645640 ghost2.out:siesta: Total = -482.646198 monomer.out:siesta: Total = -482.594903 For this case, with 0.001 Ry, we would get a binding energy of **-320.0 meV** if we do not consider the BSSE, but we get **-217.9 meV** if we take it into account with the ghost atoms. That's more than a 100 meV difference! This error (102.1 meV) accounts for as much as the **32% of the binding energy** if we do not take it into account. If we look in the literature for reported binding energy values for the dimer, we find a value of **228.6 meV** using planewaves, which is fairly close to the value we're finding roughly if we take BSSE into account. (`Corsetti et al., 2013 `__). Now, try different values for the energy shift, namely 0.01 Ry and 0.0001 Ry. Compare the binding energies with and without the BSSE correction. *How much of the original binding energy can be adjudicated to BSSE?* .. dropdown:: Answer: Basis enthalpy for different values of energy shift. :color: success We have now the energies for monomer, dimer, and ghost combinations for different energy shifts:: Energy Shift Monomer Dimer Ghost1 Ghost2 0.01 Ry -482.136281 -964.717861 -482.174049 -482.329897 0.001 Ry -482.594903 -965.509762 -482.645640 -482.646198 0.0001 Ry -482.632604 -965.578226 -482.684753 -482.671886 If we calculate the binding energies with and without BSSE correction, we should get something like this:: Energy Shift No BSSE With BSSE BSSE BSSE% 0.01 Ry -471.9 -213.9 -258.0 55 0.001 Ry -320.0 -217.9 -102.1 32 0.0001 Ry -313.0 -221.5 -91.5 29 Note how, with a lower energy shift (i.e., larger orbitals), we get a much lower contribution with regards to BSSE. In the highest value for energy shift, the BSSE is actually larger than the binding energy! In general, seeing how much of BSSE is present in things like binding energy is a good indicator of how transferrable your basis set is. If you did the :ref:`basis set optimization ` tutorial, you might want to check how each of your basis sets perform when checking for BSSE! Remember that you will have to duplicate the sections in the :ref:`PAO.Basis` block for *ghost1.fdf* and *ghost2.fdf*. They should look something similar to the following block: .. dropdown:: Optimized PAO.Basis with ghosts. .. code-block:: %block PAO.Basis H 1 n=1 0 2 P 1 4.5 2.0 O 2 n=2 0 2 5.5 3.0 n=2 1 2 P 1 7.0 3.0 H_ghost 1 n=1 0 2 P 1 4.5 2.0 O_ghost 2 n=2 0 2 5.5 3.0 n=2 1 2 P 1 7.0 3.0 %endblock PAO.Basis .. dropdown:: Answer: Basis enthalpy for optimized basis sets. :color: success We have picked here different stages from the basis set optimization tutorial. .. code-block:: Basis Set Monomer Dimer Ghost1 Ghost2 Default -482.136281 -964.717861 -482.174049 -482.329897 Rough -482.884135 -966.016972 -482.888456 -482.906653 OptimRcut -482.936001 -966.119097 -482.941834 -482.948588 OptimPolar -482.980804 -966.205087 -482.987083 -482.986743 The *Default* refers to the siesta default basis set, while *Rough* refers to the optimization done only with :ref:`PAO.EnergyShift` and :ref:`PAO.SplitNorm`, using 0.001 Ry and 0.55 respectively. *OptimRcut* refers to the optimization done on cut-off and matching radii, which results in the following block: .. dropdown:: Optimized PAO.Basis radii .. code-block:: %block PAO.Basis H 1 n=1 0 2 P 1 4.5 2.0 O 2 n=2 0 2 5.5 3.0 n=2 1 2 P 1 7.0 3.0 %endblock PAO.Basis Finally, *OptimPolar* refers to an additional optimization performed on the Polarization orbitals, resulting in this last block: .. dropdown:: Optimized PAO.Basis radii .. code-block:: %block PAO.Basis H 2 n=1 0 2 4.5 2.0 n=2 1 1 Q 4.0 4.5 O 3 n=2 0 2 5.5 3.0 n=2 1 2 7.0 3.0 n=3 2 1 Q 6.5 7.0 %endblock PAO.Basis Translating these into binding energies, we get:: Basis Set No BSSE With BSSE BSSE BSSE% Default -445.3 -213.9 -231.4 51.96 Rough -248.7 -221.9 -26.8 10.79 OptimRcut -247.1 -228.7 -18.4 7.45 OptimPolar -243.5 -231.3 -12.2 5.02 Note how the fully optimized basis sets not only reproduce the binding energy more accurately, but also the contribution of BSSE is greatly reduced. Note, also, that there is very little gain when optimizing the polarization orbitals over just optimizing the cut-off radii. .. _basis-sets-2D-materials: :bdg-success:`basic` Basis sets for 2D materials: improving vacuum description ------------------------------------------------------------------------------ .. dropdown:: **Are you under a local environment?** :color: primary Enter the directory *02-Graphene* .. figure:: images/graphene.png :width: 300px :align: center :alt: (Graphene) Graphene, the quintessential 2D system. In this example, we will illustrate the tuning of the basis set orbitals for a slab system. 2D systems in general represent a scenario in which where the vacuum must be properly described in one direction (perpendicular to the plane), specially if you want to capture specific surface-related properties. We will use graphene (the prototypical 2D system), but similar effects are important for other layered materials and surfaces (see for example `García-Gil et al., 2009 `__). First, let's have a look at the input files: * Pseudos: :download:`C.psml`, :download:`C_ghost.psml` * Input files: :download:`graphene.fdf`, :download:`graphene.diffuse.fdf`, :download:`graphene.ghosts.fdf` .. note:: In this section of the tutorial, we will be using *gnubands* and *gnuplot* to obtain plots of the band structure. *gnubands* is a utility within SIESTA, so you do not need to install anything else. Let's first have a look at *graphene.fdf*. As you can see, it is a conventional fdf file, and it doesn't even have a :ref:`PAO.Basis` block (yet). The only parts that might not sound familiar are the options regarding k-point sampling and band structure output, which are covered in other tutorials (see the :ref:`K-point sampling ` and :ref:`Analysis tools ` tutorials). .. dropdown:: A peek at the fdf file .. literalinclude:: work-files/02-Graphene/graphene.fdf :linenos: Just to characterize the problem a bit, run the graphene system varying the energy shift (:ref:`PAO.EnergyShift`) across 300 meV, 100 meV, 50 meV and 10 meV. *Take note of the total energy and the basis enthalpy*. Remember that you can run this system as usual: .. code-block:: console $ siesta < graphene.fdf > graphene.out And that you can easily grep the total energy and the basis enthalpy from the output: .. code-block:: console $ grep "E + p_basis" graphene.out $ grep "Total =" graphene.out .. dropdown:: Answer: Energy and Basis Enthalpy for different Energy Shifts :color: success You should get something similar to these results:: Eshift BasisEnthalpy TotalEnergy 10 meV -325.843477 eV -326.027291 eV 50 meV -325.992000 eV -326.102045 eV 100 meV -326.025687 eV -326.111045 eV 300 meV -325.912154 eV -325.965841 eV At first look, 100 meV looks like a reasonable value for the energy shift, by virtue of being a minimum of both total energy and basis enthalpy. If we are following the same criteria as in the basis set optimization tutorial, it would seem that 100 meV is a good energy shift to get a reasonable electronic structure. So now let's have a look at the bands. We have here the bands obtained with a planewaves code, which we asume is a reasonable description for this kind of system: .. _fig-sigma-star: .. figure:: images/Band-structure-of-graphene-VASP.png :width: 400px :align: center :alt: Graphene band structure with VASP. Graphene band structure with VASP. Now, let's try to see if SIESTA produces a structure similar to the one in the picture above. We will rely on gnubands and gnuplot for this purpose: .. code-block:: console $ gnubands -G -F -o graph -E 10 -e -20 graphene.bands $ gnuplot --persist -e "set grid" graph.gplot The output of gnubands is a datafile with the band structure and a *graph.plot* file that can be used by gnuplot. How does your band structure compare to the one obtained by planewaves? .. dropdown:: Answer: Bands for the current basis set :color: success .. figure:: images/bands-siesta-default.png :width: 400px :align: center :alt: Graphene band structure SIESTA We can see that the bands that ar below the Fermi level (0 in the Y axis) are very well descripted; however, in planewaves we can see a very different structure for those bands 4-5 eV above the Fermi level (:math:`\sigma^*`). You can see that the description is not as good as we hoped for. Moreover, the patterns we see will remain for even lower values for energy shifts (i.e., longer cutoff radii). We will now see two different ways in which this problem is tackled. Using diffuse orbitals ^^^^^^^^^^^^^^^^^^^^^^ We can now have a look at the *graphene.diffuse.fdf* file. .. dropdown:: A peek at the fdf file .. literalinclude:: work-files/02-Graphene/graphene.diffuse.fdf :linenos: The file *graphene.diffuse.fdf* includes 3s and 3p orbitals with very long tails. These are enough to lower the energy of the :math:`\sigma^*` states below 5 eV. You can see these orbitals in the :ref:`PAO.Basis` block section of the input: .. literalinclude:: work-files/02-Graphene/graphene.diffuse.fdf :linenos: :lineno-match: :start-at: %block PAO.Basis :end-at: %endblock PAO.Basis Note that the carbon atom now has *four* shells (``C 4``). We add both s and p orbitals as diffuse since we expect both to participate in a :math:`\sigma^*` band! Run the *graphene.diffuse.fdf* file and compute the bands (you will now have to use gnubands with *graphene.diffuse.bands*). Repeat the calculation, but increasing the cut-off radii of the diffuse orbitals to 12 Bohr. *How do the bands behave?* .. dropdown:: Answer: Bands with diffuse orbitals. :color: success We have here the bands with 10 Bohr-long diffuse orbitals. .. figure:: images/bands-siesta-diffuse10.png :width: 400px :align: center :alt: Graphene band structure SIESTA with diffuse orbitals (10 Bohr). And here, the bands with 12 Bohr-long diffuse orbitals: .. figure:: images/bands-siesta-diffuse12.png :width: 400px :align: center :alt: Graphene band structure SIESTA with diffuse orbitals (12 Bohr). .. dropdown:: Additional information: Using only S diffuse orbitals. :color: info In case you are wondering what happens if you only add a single diffuse S orbital, you can see the bands here: .. figure:: images/bands-siesta-diffuse-onlyS.png :width: 400px :align: center :alt: Graphene band structure SIESTA with a single S diffuse orbital. Only the first of the :math:`\sigma^*` is affected. To fully capture the correct band structure, we need the p orbitals indeed. We can see in this case that adding diffuse orbitals solves our problem. But we still have another way to fix this... using ghost atoms again! Using ghost atoms ^^^^^^^^^^^^^^^^^ How do we use ghost atoms to solve the same problem? Let's have a look at the *graphene.ghosts.fdf* file. In particular, the lines on species and geometry: .. literalinclude:: work-files/02-Graphene/graphene.ghosts.fdf :linenos: :lineno-match: :start-at: %block ChemicalSpeciesLabel :end-at: %endblock ChemicalSpeciesLabel .. literalinclude:: work-files/02-Graphene/graphene.ghosts.fdf :linenos: :lineno-match: :start-at: %block AtomicCoordinatesAndAtomicSpecies :end-at: %endblock AtomicCoordinatesAndAtomicSpecies You can see that it includes a layer of ghost orbitals (C ghosts) slightly above and below the graphene plane. By doing this we are able to describe electronic states that extend farther away from the atomic layer into the vacuum (at the cost of including a significant number of orbitals in the basis). .. note:: We use here the same carbon atoms as ghosts for convenience, but we are not aiming to correct BSSE here! Your ghost atoms could be of any other type, you could have for example ghost hydrogens, or sulfur, or whatever atom you want. Run the *graphene.ghosts.fdf* file and compute the bands (you will now have to use gnubands with *graphene.ghosts.bands*). *How do the bands compare to the cases before?* .. dropdown:: Answer: Bands with ghost atoms. :color: success We have here the bands with ghosts atoms on top of the surface. .. figure:: images/bands-siesta-ghosts.png :width: 400px :align: center :alt: Graphene band structure SIESTA with ghost atoms. Diffuse orbitals versus Ghosts ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ You have a summary of the comparisons in the following image: .. _fig-sigma-star-multi: .. figure:: images/bands.png :width: 600px :align: center :alt: Graphene bands with Siesta with various basis sets. Graphene bands with SIESTA with various basis sets. The "DZP long" part refers to a normal DZP basis set but with extremely long cut-off radii. You can see that both the ghost atoms and the diffuse orbitals are good solutions, but which one would you pick? For your own systems, consider the following: * Which of the approaches adds more basis functions to the system? How does this reflect on the total computational time? * Which of the approaches makes more sense for the technique you are using? Using ghost atoms is perfectly fine for RT-TDDFT, but it's not ideal for molecular dynamics, for example. * If you are using ghost atoms, where do you place them? And do you use the same atoms as in the surface? Or a different atom for the ghosts? * If you use diffuse orbitals, do you apply the diffuse orbital to all atoms of the same species? Or do you create a specific atomic species only for those atoms on the surface?