CamCASP/Numerical

From CUC3
Revision as of 15:33, 14 May 2009 by import>Am592 (→‎Water Dimer)
Jump to navigation Jump to search

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. This is a 3s3p3d2f1g set.
         Limit   E1elst    E1exch    E2ind(UC)  E2exind(UC)  E2disp(UC)  E2exdisp(UC)
===============================================================================
Ref/Mb     D   -1.82018  12.79049  -0.27725   0.22934    -23.03610     0.65406
-------------------------------------------------------------------------------
aTZ/Mb     D   -1.81505  12.79118  -0.27725   0.22929    -23.18343     0.67590
aTZ/Mb1    G   -1.83438  12.79007  -0.27724   0.22929    -23.03462     0.65994 <---*
           F   -1.83665  12.79021  -0.27724   0.22929    -23.04413     0.64328
           D   -1.81445  12.79021  -0.27725   0.22928    -23.09893     0.69828
           P   -1.82227  12.79127  -0.27725   0.22930    -23.38051     0.61662
-------------------------------------------------------------------------------
aQZ/Mb     D   -1.81947  12.79028  -0.27725   0.22928    -23.09406     0.65693
aQZ/Mb1    G   -1.82034  12.79005  -0.27724   0.22928    -23.03755     0.65513 <---*
           F   -1.82037  12.79007  -0.27724   0.22928    -23.03577     0.65197
           D   -1.82026  12.79007  -0.27724   0.22928    -23.04761     0.66652
           P   -1.82066  12.79031  -0.27724   0.22928    -23.06353     0.63258
===============================================================================

The best results in each basis are marked. The aQZ/Mb1 basis results in a 0.006% error in E2disp(UC) and 0.008% error in E1elst. These are negligible. More realistically, we would use the aTZ/Mb1(Limit F) basis: this results in errors of 0.9% in E1elst and 0.03% in E2disp(UC).

The mid-bond set Mb1 has not been optimized, so we should expect even better energies if we do so - possibly with a smaller basis! This should be explored at some point.

Now for SAPT(DFT) energies. Keep in mind that we have no SAPT(DFT) benchmark here. The Korona et al. benchmark is Eint(2) = -11.06 K. Given the numerical accuracies in the SAPT(KS) energies, we should take the SAPT(DFT) energies calculated in the aQZ/Mb1 basis as our benchmark. This is the first row.All calculations use Switch = 1. 

             Limit   E1elst    E1exch   E2ind     E2exind    E2disp    E2exdisp  E2int     %Diff
========================================================================================
Ref: aQZ/Mb1   G    -1.82034  -0.28699 -22.56953  12.79005   0.23734   0.64182   -11.00765  0.0
----------------------------------------------------------------------------------------
aTZ/Mb         D    -1.81505  -0.38953 -22.72746  12.79118   0.32214   0.66260   -11.15612  1.35
aTZ/Mb1        G    -1.83438  -0.31149 -22.56627  12.79007   0.25761   0.64652   -11.01794  0.09
               F    -1.83665  -0.38420 -22.57699  12.79021   0.31775   0.63024   -11.05964  0.47
               D    -1.81445  -0.34006 -22.63486  12.79021   0.28123   0.68426   -11.03368  0.24
               P    -1.82227  -0.76152 -22.94349  12.79127   0.62980   0.60509   -11.50111  4.48
----------------------------------------------------------------------------------------
aQZ/Mb         D    -1.81946  -0.31768 -22.63060  12.79028   0.26272   0.64375   -11.07100  0.58
aQZ/Mb1        G    -1.82034  -0.28699 -22.56953  12.79005   0.23734   0.64182   -11.00765  0.0
               F    -1.82037  -0.29048 -22.56772  12.79007   0.24023   0.63872   -11.00955  0.02
               D    -1.82026  -0.29644 -22.58014  12.79007   0.24516   0.65301   -11.00861  0.009
               P    -1.82066  -0.29462 -22.59850  12.79031   0.24365   0.61983   -11.05999  0.48
========================================================================================

MC or DC?

We use main basis sets of different types: MC, MC+, DC and DC+. What sort of basis type should we use for the corresponding auxiliary basis set?

So far, I've been using the same type as used for the main basis. But there is numerical evidence that this may not be the best idea.

  • MC+/DC/DC+: No problem here. I always used the DC or DC+ types for the auxiliary basis set. From the Helium dimer example above, you can see that the DC+ basis type works beautifully.
  • MC: Here's where there is a problem. So let's have a look at some calculations.

Helium dimer

aTZ/MC R=5.6 Bohr All energies in Kelvin.

As usual, we begin with the SAPT(KS) energies. The reference has been obtained without density-fitting. The auxiliary basis description is given in the first column. In all cases we have used Switch = 0 and Cartesian GTOs in the auxiliary basis. 

          E1elst  E1exch  E2ind(UC)  E2exind(UC) E2disp(UC) E2exdisp(UC)
===============================================================================
Ref     -1.43891   8.21752  -0.00302   0.00647 -19.16999   0.14762
-------------------------------------------------------------------------------
aTZ/MC   3.88603   8.30448  -0.06989   0.00339 -18.96535   0.11470
    DC  -1.43919   8.47950  -0.00303   0.00658 -19.16981   0.11523
aQZ/MC -11.40375   8.30964  -0.42265  -0.00253 -19.78456   0.13666  <---??????
    DC  -1.43899   8.31032  -0.00303   0.00652 -19.16594   0.13385
-------------------------------------------------------------------------------
aDZ/MC -4293.69336  8.84783 -126.16097   0.56257 -15.87909   0.03960  Just for fun!
    DC  -24.50340   9.09293   -0.02501   0.01275 -18.94748   0.05356
===============================================================================
  • I can't believe the aQZ/MC result. E1elst is nonsense in this basis. But I do not know why this is so.I included the aDZ/MC/DC calculations just for fun. It looks like E1elst is very sensitive to the auxiliary basis used if we use the MC basis type. You do have to use an auxiliary basis consistent with the main basis. Of course, this could be because of Helium. Will look at water later.
  • The DC basis types are clearly much better.
  • E1exch seems hard to get correct to more than 1-3%. I wonder why! I think we are seeing the difference in the full E1exch expression and the S2 approximation used by CamCASP.
  • E2disp(UC) is very accurate in the DC basis types.

So, I would conclude that we need the DC basis types even if we use the MC basis type for the main basis. Why is this so?

Water Dimer

  • Minimum energy dimer geometry from Mas and Szalewicz, J. Chem. Phys. 104 (1006)
  • Main Basis: Sadlej pVTZ : MC+ and MC types
  • Aux basis: aTZ Basis type as indicated. New mid-bond set used: 3s3p3d2f1g
  • Reference calculations without DF are in /home/am592/DistProp/systems/water/sapt-dft/outputs/sadlej/ Look at geometry (3).

(A) MC+ type for main basis

First, SAPT(KS) energies. Note that E1exch is calculated in the S2 approximation in CamCASP but not in the reference calculations from SAPT. 

                   E1elst   E1exch    E2ind(UC)  E2exind(UC) E2disp(UC) E2exdisp(UC)
=======================================================================================
Ref.              -29.0115   23.8384  -10.2594    5.3546  -12.1981    2.0515
---------------------------------------------------------------------------------------
aTZ/MC+ Switch=0  -29.06974  23.72704 -10.32765   5.38184 -12.32278   2.08945
        Switch=1  -28.94095  23.73961 -10.30060   5.38430 -12.32278   2.08879
=======================================================================================
  • Why the somewhat large error in E2disp(UC)? I think it is because the CamCASP code uses the Weighted algorithm for the position of the mid-bond set. So the mid-bonds are located in a different position: (0.0000, 0.0000, 2.8345) in the reference calculation versus (0.0038, 0.0000, 3.6479) in CamCASP. This difference could account for all the small differences in energies.

So let's say that there are no surprises here and that we have rather good energies.

Here are the SAPT(DFT) energies: 

                   E1elst   E1exch    E2ind      E2exind   E2disp    E2exdisp   E2int
=======================================================================================
Ref.              -29.0115   23.8384  -10.615     5.540    -9.630     1.619     -18.258
---------------------------------------------------------------------------------------
aTZ/MC+ Switch=0  -29.0697   23.7270  -10.5667    5.4881   -9.7279    1.6495    -18.4998
        Switch=1  -28.9410   23.7396  -10.5667    5.5068   -9.7279    1.6489    -18.3402
=======================================================================================

(B) MC type for main basis

Now let's look at the MC type. This should be more interesting.

Retrieved from "http://wikis.ch.cam.ac.uk/cuc3/wiki/index.php?title=CamCASP/Numerical&oldid=3828"