Basis sets - Special cases¶

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 tutorial on basis set optimization. This optimization usually comes before considering what we will cover here.

Take me to the exercises!

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

basic Fixing the Basis Set Superposition Error¶

(A water molecule and a ghost water molecule.)

Fig. 3 A water molecule (right) with another ghost water molecule (left) to represent what we have in the dimer.¶

If you ran the tutorial on 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₁+ (Ghost)Water₂

(Ghost)Water₁+ Water₂

We’ll now have a look at the input files:

Pseudos: O.psml, H.psml, O_ghost.psml, H_ghost.psml
Input files: monomer.fdf, dimer.fdf, ghost1.fdf, 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:

%block ChemicalSpeciesLabel
 1  8  O
 2  1  H
 3 -8  O_ghost
 4 -1  H_ghost
%endblock ChemicalSpeciesLabel

AtomicCoordinatesFormat  Ang
%block AtomicCoordinatesAndAtomicSpecies
 0.000000    0.000000    0.000000 4
 0.000000    0.000000    1.522047 4
 0.436628    0.395162    0.761024 3
 2.341468    0.000000    0.760784 2
 3.795635    0.459564    0.760954 2
 3.251445   -0.331791    0.760767 1
%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?

Second, if we use the PAO.Basis block (which, we will), we also need to generate a new section for this new species!

%block PAO.Basis
H 1
 n=1 0 2 P 1
  0.0 0.0
O 2
 n=2 0 2
  0.0 0.0
 n=2 1 2 P 1
  0.0 0.0
H_ghost 1
 n=1 0 2 P 1
  0.0 0.0
O_ghost 2
 n=2 0 2
  0.0 0.0
 n=2 1 2 P 1
  0.0 0.0
%endblock PAO.Basis

Now that you’ve seen the important parts of it, What is the difference between ghost1.fdf and ghost2.fdf ?

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:

\(E_{\rm binding} = E_{\rm dimer} - 2E_{\rm 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:

\(E_{\rm binding} = E_{\rm dimer} - E_{\rm ghost1} - E_{\rm ghost2}\)

Now we can run our siesta calculation as usual:

$ 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:

$ grep "Total =" *.out

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?

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 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 PAO.Basis block for ghost1.fdf and ghost2.fdf. They should look something similar to the following block:

basic Basis sets for 2D materials: improving vacuum description¶

Fig. 4 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: C.psml, C_ghost.psml
Input files: graphene.fdf, graphene.diffuse.fdf, 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 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 K-point sampling and Analysis tools tutorials).

Just to characterize the problem a bit, run the graphene system varying the energy shift (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:

$ siesta < graphene.fdf > graphene.out

And that you can easily grep the total energy and the basis enthalpy from the output:

$ grep "E + p_basis" graphene.out
$ grep "Total =" graphene.out

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. 5 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:

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

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.

The file graphene.diffuse.fdf includes 3s and 3p orbitals with very long tails. These are enough to lower the energy of the \(\sigma^*\) states below 5 eV. You can see these orbitals in the PAO.Basis block section of the input:

%block PAO.Basis
C        4            # We include 4 shells
 n=2   0   2          # 2s (double zeta)
   0.000      0.000
   1.000      1.000
 n=2   1   2 P   1    # 2p (double zeta) + polarization (single zeta)
   0.000      0.000
   1.000      1.000
 n=3   0   1          # 3s diffuse orbital (single zeta)
  10.000
   1.000
 n=3   1   1          # 3p diffuse orbital (single zeta)
  10.000
   1.000
%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 \(\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?

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:

%block ChemicalSpeciesLabel
    1   6  C
    2  -6  C_ghost  # Ghost species!
%endblock ChemicalSpeciesLabel

%block AtomicCoordinatesAndAtomicSpecies
33333333  0.66666667  0.5        1     # C
66666667  0.33333333  0.5        1     # C
33333333  0.66666667  0.55       2     # Ghost-C above
66666667  0.33333333  0.55       2     # Ghost-C above
33333333  0.66666667  0.45       2     # Ghost-C below
66666667  0.33333333  0.45       2     # Ghost-C below
%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?

Diffuse orbitals versus Ghosts¶

You have a summary of the comparisons in the following image:

Fig. 6 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?