Basis sets - Optimization¶

Using the default basis sets generated by SIESTA might be enough for some applications, but doing some degree of manual optimization of the basis sets may help to achieve better results with similar computational costs. This is especially advisable when dealing with extremely large systems, where going for high basis set cardinality (triple-zeta, quadruple-zeta) is not really an option.

Note that the concepts presented in this tutorial apply to both molecules and crystals. Don’t hesitate to use them for your own systems, whatever they are!

Take me to the exercises!

background General Concepts¶

Cardinality: Multiple-zeta and polarization¶

SIESTA uses atomic-centered basis sets, which are the product of a radial function (R(r)) and a spherical harmonic (Y(m,l)). The basis set cardinality is what determines the total amount of basis functions that we have per atom.

In a minimal basis set, one would have one basis function per valence orbital. So for a Carbon atom, we would get one 2s function and three 2p functions. This is called a single-zeta (SZ) basis set, which contains one function per orbital.

We can, of course, increase the quality of this basis set by increasing the amount of radial functions that describe the same orbital. So in a double-zeta basis set (DZ), we have two functions per orbital. This means that for C, we would now have eight functions instead of the minimal four. This can go up indefinitely, to triple-zeta (TZ), quadruple-zeta (QZ), and so on.

The last part that contributes to cardinality are the polarization orbitals. These increase the “flexibility” of a basis set, providing and extra set of orbitals of a higher l. So in the case of a double-zeta polarized basis set (DZP), our C atom will now add 5 d-type functions to the mix (one l higher than the p-orbitals), thus increasing the total to 13 functions.

Save for a few exceptions, a DZP basis is more than enough to get a good description. In general, increasing the cardinality of a basis set greatly increases the cost, since the most time-consuming part of SIESTA has a cost proportional to the cube of the amount of basis functions.

Cut-off radii¶

Remember that a basis function in SIESTA are the product of a radial function and a spherical harmonic. The spherical harmonic is determined by the quantum numbers of the orbital; as such, it is not something we can really modify. However, we can do work on the radial part.

The cut-off radius of a basis function indicates the point after which the radial function becomes zero. For multiple-zeta basis sets, the matching radius indicates the point after which the n-zeta orbitals coincide with the first-zeta.

In general, the larger the cut-off radii, the better described will be our system. However, this also comes at a cost; in particular, an increase in computational costs. So how do we choose a cut-off radius that is large enough to ensure good quality results, but also does not increase our simulation times unnecessarily?

Basis Enthalpy¶

The basis enthalpy is a fictitious magnitude that adds an extra term to the energy, which depends on the volume of the orbitals; essentially, this can tell us when the basis orbitals are getting needlessly large. This extra term is essentially P*V, where V is the orbital volume and P is a fictitious “pressure”. This pressure is an input we will use in SIESTA.

For bulk materials such as metals (Cu, Fe) or semiconductors (Si, Ge), the default value of 0.2 GPa for the basis pressure produces reasonable results. However, for small atoms in isolated molecules (C, O, H, etc.) this pressure results in extremely short cut-off radii. One typically reduces the basis pressure to a tenth of the default value, i.e. 0.02 GPa.

How do you know whether your orbitals are large enough? A good rule of thumb is to consider that orbitals with a cut-off radius of less than 5.0 Bohr are very small, while a cut-off radius larger than 10.0 Bohr is unnecessarily large. Of course, there are a few exceptions: 1s and 2s orbitals of the first two rows of the periodic table (specially H) can be slightly smaller than 5.0 Bohr, but going below 4.0 Bohr should be avoided. In the oposite case, when dealing with surfaces or 2D materials one might need larger, diffuse orbitals in order to better describe the area around the exposed atoms: these diffuse orbitals might easily reach 10 to 12 Bohr. More on this topic is covered in Basis sets - Special cases.

Optimization and transferability¶

When optimizing basis sets, the concept of transferability is one of the most important ideas to keep in mind. You might get the ultimate basis set that gives the perfect results for a given set of atoms in a given set of coordinates; however, if even a slight re-arrangement of those atoms produces terrible results, that means your basis set is not transferable. Ideally, one should aim for highly transferable basis sets that also produce reasonable results.

This ties in with a key aspect of basis set optimization: defining your optimization function, i.e., the magnitude that you will be using to monitor how well your basis is performing. One good magnitude to look at is the total energy of a given system, but even better is to use the concept of basis enthalpy.

While it might be tempting to optimize basis sets using bond lengths, interaction energies, or even band structures, these typically result in poorly transferable basis sets. The use of these magnitudes as optimization functions is not encouraged; however, these can be used after the optimization is done in order to benchmark the behaviour of the newly optimized basis set.

Having defined the optimization function, one needs to define which are the optimization variables, i.e. the things that we will modify manually to improve our basis sets. While SIESTA has a ton of possibilities to modify your basis set, the most important targets for optimization are the orbital radii: the cut-off radii for first-zeta orbitals, and the matching radii for multiple-zeta. Polarization orbitals can also be optimized separately, but the way these are generated in SIESTA is generally good and they usually do not require further optimization.

Note

To summarize:

We aim for good-quality, highly-transferable basis sets.
Typically we will be modifying the cut-off radii of first-zeta and the matching radii of multiple-zeta orbitals.
We will take the basis enthalpy (or the total energy) of a system as our main way of optimizing the basis sets.
We need to have a way to test if our basis set is transferable enough.

basic Rough optimization of the basis set¶

Note

This tutorial shows just one possible method to optimize your basis sets. Different people tend to use different techniques, there is no “one true way”. This method does, however, work reliably across different kinds of systems.

Note

We will stick to double-zeta polarized (DZP) basis sets for this tutorial, which should be enough for most calculations with SIESTA.

Fig. 1 For SIESTA, the water molecule is three points in space: one representing oxygen (red) and two representing hydrogen (white).¶

Molecular dynamics of water (or aqueous systems in general) is still one of the most common types of simulations. In this tutorial, we will optimize the basis set for a single water molecule, and we will use the energy of a water dimer as a way to benchmark the quality of our basis functions.

You can find the input files for this first section here:

Pseudos: O.psml, H.psml
Input files: water.fdf

It is also important to note that we have set up a relatively big simulation box. Since water is an isolated molecule and not a bulk metal, we want to avoid orbital overlaps between one molecule and its periodic images:

LatticeConstant 15.0 Ang
%block LatticeVectors
 1.0 0.0 0.0
 0.0 1.0 0.0
 0.0 0.0 1.0
%endblock LatticeVectors

Take a look at the following lines in the water.fdf file:

PAO.EnergyShift     0.01 Ry
PAO.SplitNorm       0.15

These two options allow us a rough control of the quality of the basis set. PAO.EnergyShift allows us to control the cut-off radii of the first-zeta orbitals: the lower the energy shift, the larger the cut-off radii. PAO.SplitNorm allows us to control the matching radii of the multiple-zeta orbitals: the lower the value, the larger the matching radii (i.e., more similar to the first-zeta orbitals). The values presented in the input file are the defaults.

We want orbitals that are good quality but not needlessly large, so we will be relying on the Basis Enthalpy (as defined above) to guide our optimization. Since we are using small atoms, the variable BasisPressure will be set to 0.02 GPa, as you can see in the following line:

BasisPressure  0.02 GPa

You can run this fdf as usual:

$ siesta < water.fdf > water.out

Now that we have all set, lets start by running a single-point calculation for different values of PAO.EnergyShift. We will compare the results for the basis enthalpy in each case. You can get the basis enthalpy from the general output under “(Free)E+ p_basis*V_orbitals”, or from the first line of the <systemlabel>.BASIS_ENTHALPY file:

$ grep "E + p_basis" water.out

$ cat h2o.BASIS_ENTHALPY

Test the following values for the energy shift 0.01 Ry, 0.001 Ry, and 0.0001 Ry. Compare the basis enthalpy for all three values. Which of the three performs better?.

Now we will test different values for PAO.SplitNorm. Set the energy shift to the one we found above, 0.001 Ry. Test the values 0.45, 0.55 and 0.65 for the Split Norm. Which of the three gives the better basis enthalpy?

Note

In real-life applications, you would typically optimize both the energy shift and the split norm iteratively, as they might not be totally independent of each other. Be sure to explore a larger array of values for both parameters, although differences below 1 meV in the basis enthalpy are usually not that relevant.

basic Calculating the binding energy of a water dimer¶

Fig. 2 A water dimer, i.e. two water molecules sharing a hydrogen bond.¶

Now is the time to put our optimized basis set to test. We will use as an example the binding energy of a water dimer, i.e.:

\(E_{\rm binding} = E_{\rm dimer} - 2E_{\rm monomer}\)

The input files are the same as before, with monomer.fdf being the single water molecule we used before:

Pseudos: O.psml, H.psml
Input files: dimer.fdf and monomer.fdf

You will notice that both files have the default values for the energy shift and split norm:

PAO.EnergyShift     0.01 Ry
PAO.SplitNorm       0.15

Run both files as usual:

$ siesta < monomer.fdf > monomer.out
$ siesta < dimer.fdf > dimer.out

Take note of the total energy for both systems (remember to look for Total in the output file). Then, change the PAO.SplitNorm and PAO.EnergyShift values in the monomer.fdf and dimer.fdf to the values obtained in the previous section (if you didn’t do it, pick 0.001 Ry for the energy shift and 0.55 for the split norm). Run the calculations again. How different is our binding energy for the dimer?

You can see now that the binding energy with the lightly optimized basis setis significantly different. Does this mean that the default basis sets are awful? Not really, because we are not yet accounting the Basis Set Superposition Error (BSSE), which will be covered in a later tutorial: Basis sets - Special cases. However, it highlights the impact of even a slight optimization of the basis set.

background Manual optimization of the cut-off radii¶

Just relying on PAO.EnergyShift and PAO.SplitNorm is surely convenient when we have to work on a large amount of very different systems. But if our focus is to get the best quality results (for example when working with very subtle effects in energy or fine differences in band structure), then we need to go a step further and optimize the cut-off and matching radii manually.

In the following sections, we will start with a basis set generated with a value of 0.001 Ry for PAO.EnergyShift and the default value for the PAO.SplitNorm. However, this is just the starting point, we will be modifying these values the values for the orbital radii manually. By the end of this tutorial, the energy shift will become redundant, and can be actually removed from the input file.

But first, we need to have a look at the PAO.Basis block. Now, instead of defining the cardinality just as PAO.BasisSize, we will make use of the PAO.Basis block which expands the information as follows:

%block PAO.Basis
  H 1
    n=1 0 2 P 1     # n, l (=0, s-orbital), number of zetas, include 1 set of polarization orbitals
        0.0 0.0     # rcut for first zeta, rmatch for second zeta
  O 2               # Number of different l
    n=2 0 2         # n, l (=0, s-orbital), number of zetas
        0.0 0.0     # rcut for first zeta, rmatch for second zeta
    n=2 1 2 P 1     # n, l (=1, p-orbital), number of zetas, include 1 set of polarization orbitals
        0.0 0.0     # rcut for first zeta, rmatch for second zeta
%endblock PAO.Basis

Notice how the H atom has a single 1s level (n=1/l=0) with 2 functions (zetas) and a single set of polarization orbitals (P 1). Meanwhile, Oxygen has two levels, 2s and 2p (n=2/l=0 and n=2/l=1), with an additional set polarization orbitals on the 2p level (P 1).

In this example, the input is the same as specifying the flag PAO.BasisSize DZP, since we have a double-zeta basis set for all valence orbitals, and we are including a single polarization orbital set for each atom.

The values for rcut and rmatch are all set to 0.0; this tells SIESTA that it should automatically find appropriate values for each of those radii using the energy shift. By the end of this tutorial, we will have set each of these radii manually. It is worth noting that all values for radii here should be in Bohr.

Note

As with the first part of the tutorial, in real applications you should converge all rcut values either iteratively or simultaneously.

basic Optimizing the first-zeta cut-off radii¶

You will now need the following input files:

Pseudos: O.psml, H.psml
Input files: water.fdf

We will first start by optimizing the cut-off radii for the first zeta orbitals. This means we will be changing the first 0.0 we find for each of the orbitals (the * * *are for emphasis and should not be included in the actual input file):

%block PAO.Basis
H 1
  n=1 0 2 P 1     # n, l (=0, s-orbital), number of zetas, include 1 set of polarization orbitals\\
    *0.0* 0.0     # rcut for first zeta, rmatch for second zeta
O 2               # Number of different l
  n=2 0 2         # n, l (=0, s-orbital), number of zetas
    *0.0* 0.0     # rcut for first zeta, rmatch for second zeta
  n=2 1 2 P 1     # n, l (=1, p-orbital), number of zetas, include 1 set of polarization orbitals
    *0.0* 0.0     # rcut for first zeta, rmatch for second zeta
%endblock PAO.Basis

Modify the value for rcut in the 1s orbital of Hydrogen, going from 4.0 Bohr to 6.0 Bohr by 0.5 Bohr increases. Take note of the basis enthalpy of the system for each value of rcut.

Now we choose the cut-off radii for which the basis enthalpy is minimized. Assuming the value is 4.5 Bohr, your basis block will look somewhat similar to this:

%block PAO.Basis
  H 1
    n=1 0 2 P 1     # n, l (=0, s-orbital), number of zetas, include 1 set of polarization orbitals
        4.5 0.0     # rcut for first zeta, rmatch for second zeta
  O 2               # Number of different l
    n=2 0 2         # n, l (=0, s-orbital), number of zetas
        0.0 0.0     # rcut for first zeta, rmatch for second zeta
    n=2 1 2 P 1     # n, l (=1, p-orbital), number of zetas, include 1 set of polarization orbitals
        0.0 0.0     # rcut for first zeta, rmatch for second zeta
%endblock PAO.Basis

After choosing a value for the Hydrogen 1s orbital, do the same procedure for the 2s orbital of Oxygen, and then for the 2p orbital. Explore values from 4.5 Bohr to 7.5 Bohr for both orbitals.

Do the values for the optimized cut-off radii coincide with what you would expect from chemical intuition? Which is the largest orbital and which one is the shortest?

For the purposes of this tutorial, we will stop here the optimization of rcut. In a more general case, it would be better to re-optimize all the orbitals iteratively to guarantee that the rcut values are truly converged.

basic Optimizing the second-zeta matching radii¶

The input files are essentially the same as in the previous section:

Pseudos: O.psml, H.psml
Input files: water.fdf

We will now expand the previous block to the one we had at the beginning for the DZP basis set, but instead of zeroes in the rcut, we will place our optimized values for the single zeta basis:

%block PAO.Basis
  H 1
    n=1 0 2 P 1      # n, l (=0, s-orbital), number of zetas, include 1 set of polarization orbitals
        4.5 *0.0*    # rcut for first zeta, rmatch for second zeta
  O 2                # Number of different l
    n=2 0 2          # n, l (=0, s-orbital), number of zetas
        5.5 *0.0*    # rcut for first zeta, rmatch for second zeta
    n=2 1 2 P 1      # n, l (=1, p-orbital), number of zetas, include 1 set of polarization orbitals
        7.0 *0.0*    # rcut for first zeta, rmatch for second zeta

%endblock PAO.Basis

We will do the same as before, but now changing the value for the matching radii (which are still 0.0 in the block). Using a similar procedure as before, find the value of rmatch for the 1s orbital of Hydrogen that minimizes the basis enthalpy, evaluating it from 1.5 Bohr to 3.0 Bohr. Do the same for Oxygen 2s and then 2p, changing the values from 2.5 Bohr to 4.0 Bohr for each of them.*

Congratulations! You optimized your basis set for a water molecule. In the same way as in the previous section, we will stop here. However, remember that in actual production runs you should verify that all values for the cut-off radii and matching radii are properly converged.

basic Checking the result with the water dimer¶

It’s the water dimer again! We will use the same input files as before, but now we will be using the optimized basis set with the PAO.Basis block. The input files are:

Pseudos: O.psml, H.psml
Input files: dimer.fdf and monomer.fdf

Now is the time to put our optimized basis set to test. Remember the binding energy of a water dimer, i.e.:

\(E_{\rm binding} = E_{\rm dimer} - 2E_{\rm monomer}\)

Calculate the binding energy of the dimer again. How different is our result when compared to the previous ones? Consider both the total energies for each system, and the total binding energy.

advanced Optional: Optimizing the polarization orbitals¶

Note

This step is tipically not necessary, and it may help only in very specific cases. SIESTA already generates good polarization orbitals.

One last step in the orbital optimization could be to optimize the polarization orbitals. This is conceptually a bit different from what we have been doing up to now, even if the procedure is somewhat similar. Here, we will not be changing the cut-off radius of the orbital as above; instead, we will be altering the overall shape of the polarization orbital.

For the input files:

Pseudos: O.psml, H.psml
Input files: water.fdf

For that purpose, we will expand the basis block to add one more orbital set per atom:

%block PAO.Basis
  H *2*
    n=1 0 2              # n, l (=0, s-orbital), number of zetas
        *RCUT1 RMATCH1*  # rcut for first zeta, rmatch for second zeta
    n=2 1 1 *Q 0.0*      # n, l (=1, p-orbital), number of zetas, charge value for confinement
        *RCUT1*          # rcut for first zeta, rmatch for second zeta
  O *3*
    n=2 0 2              # n, l (=0, s-orbital), number of zetas
        *RCUT2 RMATCH2*  # rcut for first zeta, rmatch for second zeta
    n=2 1 2              # n, l (=1, p-orbital), number of zetas
        *RCUT3 RMATCH3*  # rcut for first zeta, rmatch for second zeta
    n=3 2 1 *Q 0.0*      # n, l (=2, d-orbital), number of zetas, charge value for confinement
        *RCUT3*          # rcut for first zeta
%endblock PAO.Basis

Note that, in both cases, the polarization orbital belongs to the following principal number (n=3 for Oxygen, n=2 for Hydrogen). In addition, all polarization orbitals have the same cut-off radius as the orbital they are polarizing. So, using the results from the previous sections of this tutorial, our starting PAO.Basis block looks like this:

%block PAO.Basis
H 2
 n=1 0 2
  4.5 2.0
 n=2 1 1 Q 0.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 0.0
  7.0
%endblock PAO.Basis

One of the recommended methods to generate appropriate polarization orbitals is the charge confinement method (more details in the manual), but suffice to say that this method relies on seeing how the atomic wavefunctions change in the presence of a charge confinement potential, \({ V(r;Q) = Q / r^{2}}\). In order to optimize our polarization orbitals, we just need to change the value of the charge Q.

Modify the charge used in the Hydrogen 2p polarization orbital (the number following Q) from 3.5 to 6.0, and choose a value that minimizes the basis enthalpy. Do the same for the Oxygen 3d polarization orbital, modifying the charge from 5.0 to 8.0.

Repeat the calculations for the water dimer using your newly optimized basis set. How do results change? Use the same input as in the previous section and modify the PAO.Basis block accordingly.