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

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:

- Calculate the total energy with both reactants, including all basis functions and nuclei. This total energy is referred to as \(E_{12}\).
- 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}\).
- 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:

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:

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