Skip to content

The Counterpoise Correction

The convergence of observables with the number basis functions has been demonstrated to be very slow for electronic structure methods that depend also on the unoccupied states (cf. the chapter for GW). This has an immediate consequence when studying bond breaking or forming between two or more reactants: The basis functions localized on one reactant can act as additional functions for the other reactant, and vise versa. This error is called basis set superposition error (BSSE). For finite basis-set sizes, this can introduce a significant error for the bonding curves – the smaller the basis set the larger the BSSE error.

The brute force way to account for the BSSE is to add more and more basis functions to the system. However in practice, this is not always possible as the computational cost of such calculations can explode rapidly. Alternatively, it is possible to calculate the counterpoise correction, which approximates the bias to the observables of each reactant due to the additional basis functions in the system with both reactants included.

Here, we will study the BSSE and the corresponding counterpoise correction for the sodium dimer for different basis sets. We will use the "RPA + Exact Exchange" total energy method (that is, the correlation energy in the random-phase approximation plus the exact exchange energy) to demonstrate this effect. The physical observable of interest in this part is the binding energy of the sodium dimer:

\[ E_\text{bind} = E_\text{Na-Na} - 2*E_\text{Na} \]

where \(E_\text{Na-Na}\) refers to the total energy of the Na dimer and \(E_\text{Na}\) to the total energy of the single Na atom.

How to perform counterpoise correction

We will explain the counterpoise correction for the case of two reactants. The calculation of the counterpoise correction requires three steps:

  1. Calculate the total energy with both reactants, including all basis functions and nuclei. This total energy is referred to as \(E_{12}\).
  2. Perform a calculation for each reactant using the same geometry as in the complex for each of the parts. This will give the energies for \(E_{1}\) and \(E_{2}\).
  3. Repeat the calculations of step 2. However, this time including the basis functions of the other reactant at their original positions from step 1, but without including the charge of the nuclei. This will give the energies \(E_{1}^\ast\) and \(E_{2}^\ast\). In FHI-aims, step 3 can be realized by using the empty keyword instead of the atom keyword in the geometry.in file.

The counterpoise correction is the difference of total energies of the single reactants with and without the additional basis function included:

\[ \Delta E_\text{cpc} = (E_{1}^\ast - E_{1}) + (E_{2}^\ast - E_{2}) \]

The final counterpoise correction for the dimer reduces to (since both reactants are identical): $$ \Delta E_\text{cpc} = 2 * (E_\text{atom}^\ast - E_\text{atom}) $$ where \(E_\text{atom}\) is the total energy for the single atom.

Let us get more specific. We want to discuss the effect of the counterpoise correction for the Na dimer. In case of a dimer, the counterpoise corrected binding energy is:

\[ E_\text{bind, cpc} = E_\text{Na-Na} - 2*E_\text{Na} - \Delta E_\text{cpc} = E_\text{Na-Na} - 2*E_\text{Na}^\ast \]

The geometry.in files

To calculate the counterpoise correction we need to setup three different geometries:

For step 1, we need the following geometry.in file:

atom 0.0 0.0 0.0 Na
atom 0.0 0.0 <z> Na

The <z> will be bonding distance of the Na atoms in the dimer, which will be varied for a certain range.

For step 2:

atom 0.0 0.0 0.0 Na
    initial_moment         1.000000

The Na atom naturally has one unpaired electron, so we initialize its spin state with an initial moment of 1.

And finally for step 3:

empty 0.0 0.0 0.0 Na
atom  0.0 0.0 <z> Na
    initial_moment         1.000000

Here, we replace atom with the empty keyword for the first Na atom. This implies that no nuclei charge will be considered, but all of the basis functions will be included at this coordinate.

The control.in files

The corresponding control.in file for step 1 is:

xc                                 pbe
relativistic                       atomic_zora scalar
total_energy_method                rpa

This triggers a non-spinpolarized PBE calculation (xc pbe) with a subsequent RPA calculation (total_energy_method rpa). Beforehand, we tested different spin-configurations for the Na dimer (up-up, up-down), but it turned out that for various distances the total spin moment results to zero.

For the steps 2 and 3 (spin collinear calculations) the control.in file is:

xc                                 pbe
spin                               collinear
relativistic                       atomic_zora scalar
total_energy_method                rpa

Very large basis sets and nearly singular overlap matrix

For very large basis sets it might be needed to include the keyword override_illconditioning .true. in the control.in file, which allows to override a built-in stop in FHI-aims and run with a nearly singular overlap matrix. If included, you should take special care that the results appear reasonable (e.g. by cross-checking for smaller basis sets).

The binding curve of the Na dimer

To demonstrate the drastic impact of the BSSE, we will perform single-point calculations for the following distances of the Na dimer (in Å):

2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6

For each distance, you should create a separate folder. For all distances, you will need to perform two calculations: raw and counterpoise, where raw refers to step 1 (the dimer geometry with all basis functions and the nuclei charge included) and counterpoise refers to step 3 (the dimer geometry, but one atom with empty site, i.e. no nucleus charge).

Overall, you should have folder structure that looks as follows:

light
├── atom
├── counterpoise
│   ├── 2.0
│   ├── 2.2
│   ├── 2.4
│   ├── 2.6
│   ├── 2.8
│   ├── 3.0
│   ├── 3.2
│   ├── 3.4
│   └── 3.6
└── raw
    ├── 2.0
    ├── 2.2
    ├── 2.4
    ├── 2.6
    ├── 2.8
    ├── 3.0
    ├── 3.2
    ├── 3.4
    └── 3.6

As a first test we will use the light species defaults. Attach the light species defaults to all control.in files, and for all folders create the corresponding geometry.in files (with the correct dimer distance) and control.in file.

Since we want to calculate the counterpoise correction for different basis sets, it is advisable to create a script for the creation of the folder structure. In total, you need to run 19 calculations for a single basis set.

Repeat the above 19 calculations, for the following species defaults:

  • light: defaults_2020/light
  • tight: defaults_2020/tight
  • cc-pV5Z (Dunning correlation-consistent Gaussian basis sets): non-standard/gaussian_tight_770/cc-pV5Z
  • NAO-VCC-5Z (Dunning-like valence-correlation-consistent NAO basis sets): NAO-VCC-nZ/NAO-VCC-5Z

Warning

For the basis sets cc-pV5Z and NAO-VCC-5Z you have to set override_illconditioning .true. in all of the corresponding control.in files.

Plot the raw (uncorrected) and counterpoise-corrected binding energies as a function of the dimer distance. You should be able to reproduce a plot similar to the one shown below.

For the discussion of this plot, let us first focus on the light results (blue curves). The dashed lines show the raw (uncorrected) binding curve. Comparing to the counterpoise-corrected binding, we find that for the the minimum of the binding curve, the binding energy increases by ~0.2 eV and the dimer distance moves ~0.2Å to the right due to the counterpoise-correction.

Now look at the tight result. The effect of the BSSE is enormous: The minimum is off by 0.7Å and 1.1eV too low compared to the corrected curve. At a first glance, this may appear inconsistent to the light results. We promised at the beginning that the BSSE will decrease by increasing the basis set. However, the FHI-aims NAO species defaults not only affect the number of basis function, it also affects the real-space integration grid, which in case of tight settings is also tightened (the cut-off radius and the number of grid points are increased). Overall, the NAO FHI-aims basis sets are optimized for DFT calculations and do not guarantee systematic convergence for beyond-DFT methods.

Lets look at the non-standard NAO and Gaussian basis sets NAO-VCC-5Z and cc-pV5Z. We find that despite a difference in the uncorrected data, the results for the counterpoise-corrected binding energies agree very well. In these cases, although we used the largest available basis set but the BSSE is still significant.