CamCASP/Numerical

CamCASP => Numerical Issues

Introduction

This page contains information related to the numerical accuracy of CamCASP.

Integral Switch

Integrals can be calculated using density-fitting (the default). But for a few kinds of 2-index integrals (nuclear, overlap, etc.) there is the possibility of calculating them without density-fitting. This can often be advantageous, in fact, for nuclear and overlap integrals, this is probably the thing to do whenever you can.

The energy modules all have the optional command

 INTEGRAL Switch = <switch>

that can be used to set the method used to calculate the integrals. In general, Switch = 0 means use density-fitting, and Switch = 1 means use the exact evaluation (no density-fitting) if possible. Not all integrals can be evaluated without density-fitting. See df_integrals.F90 for details.

The 4-index Coulomb integrals are robust (see Robust Integrals for details.). So these can be obtained quite accurately using density-fitting. In any case, there is no other way of calculating these integrals, so the Switch command has no effect on these.

On the other hand, the nuclear and overlap integrals are not robust when density-fitting is used (the error made in these integrals is linear in the error made in the density. So these should be evaluated using Switch = 1. This is not the default in version 5.4.00 of the code, so it needs to be put in manually.

Here are some examples:

• Helium dimer: aTZ/aTZ MC+ PBE0/AC R = 5.6 Bohr. All energies in Kelvin.

First the SAPT(KS) energies. I.e. second-order energies calculated using the un-coupled approximation. The reference energies don't use density-fitting. Auxiliary basis: aug-cc-pVTZ.

            E1elst    E1exch    E2ind(UC)  E2exind(UC)  E2disp(UC)  E2exdisp(UC)
===============================================================================
Ref.        -1.82018  12.79049  -0.27725   0.22934    -23.03610     0.65406
Switch = 0   0.32482  12.79195  -0.37673   0.23211    -23.18343     0.67636
Switch = 1  -1.81505  12.79118  -0.27725   0.22929    -23.18343     0.67590
===============================================================================

The improvement in E1elst and E2ind(UC) and E2exind(UC) is phenomenal! The other energy components are largely unchanged. To improve these you need to increase the size of your auxiliary basis:

              E1elst    E1exch    E2ind(UC) E2exind(UC)  E2disp(UC)  E2exdisp(UC)
===============================================================================
aQZ/Switch=1 -1.81947  12.79028  -0.27725  0.22928    -23.094059    0.65693
===============================================================================

Using the aug-cc-pVQZ auxiliary basis with switch = 1 gives us energies with the worst-case error of 0.3%. Further improvements may be possible by using a better mid-bond set for the auxiliary basis.

And now the SAPT(DFT) energies. There is no reference for the second-order energies so I've used our early SAPT(DFT) calculations here.

             E1elst    E1exch    E2ind   E2exind  E2disp   E2exdisp  E2int     %Diff
============================================================================
Ref.         -1.82018  12.79049 -0.38953 0.32221 -22.72746 0.64530 -11.17916     0.0
aTZ/switch=0  0.32482  12.79195 -0.38953 0.23999 -22.72746 0.66306  -9.09717    18.62
aTZ/switch=1 -1.81505  12.79118 -0.38953 0.32214 -22.72746 0.66260 -11.15612     0.21
aQZ/switch=1 -1.81947  12.79028 -0.31768 0.26272 -22.63060 0.64375 -11.07100     0.97
============================================================================

The larger percentage difference in the last case (aQZ/Switch = 1) is not a measure of any inaccuracy in this calculation because the reference calculation also has numerical errors associated with density-fitting. What is great is the excellent agreement with the reference and the aTZ/swtich=1 result. These two should agree as both use the aug-cc-pVTZ auxiliary basis.

Auxiliary Basis Sets

Mid-bond set

So far, we have used the same mid-bond set for the main and auxiliary bases. This is not right as the mid-bond set for the auxiliary basis should be larger. Exactly how large is the question, so let's look at some energies.

• Helium dimer: aTZ/aTZ MC+ PBE0/AC R = 5.6 Bohr. All energies in Kelvin.

Since we have a reference for the SAPT(KS) energies (i.e., the un-coupled approximation) we will look at these energies first. All calculations use switch = 1.

• Mb = Same mid-bond set as used in the main basis.
• Mb1 = Mid-bond set created from Mb. Approximately spans the direct product space Mb x Mb.

            E1elst    E1exch    E2ind(UC)  E2exind(UC)  E2disp(UC)  E2exdisp(UC)
===============================================================================
Ref.        -1.82018  12.79049  -0.27725   0.22934    -23.03610     0.65406
aTZ/  -1.81505  12.79118  -0.27725   0.22929    -23.18343     0.67590
===============================================================================

MC or DC?

What sort of basis set