CamCASP/Numerical

From CUC3
Revision as of 10:41, 10 November 2010 by import>Am592 (→‎Water Dimer)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
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. All calculations use Switch = 1. 

             Limit   E1elst    E2ind    E2disp    E1exch     E2exind   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 : DC+ 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).One problem: The main basis in these calculations used Cartesian GTOs! Will cause differences in all energies.

(A) DC+ 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. All energies in kJ/mol. 

                   E1elst   E1exch    E2ind(UC)  E2exind(UC) E2disp(UC) E2exdisp(UC)
=======================================================================================
Ref.              -29.0160   23.9420  -10.3709    5.4518  -12.3344    2.0928
---------------------------------------------------------------------------------------
aTZ/DC+ 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 E1exch? Error due to S2 approximation in CamCASP?
  • 0.7% error in E2ind(UC). Also keep in mind that the reference calculation used Cartesian GTOs in the main basis.
  • Some differences should arise because CamCASP 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: (kJ/mol) 

                   E1elst   E1exch    E2ind      E2exind   E2disp    E2exdisp   E2int
=======================================================================================
Ref.              -29.0160   23.9420  -10.7739    5.6634   -9.7518    1.6547    -18.2815
---------------------------------------------------------------------------------------
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
=======================================================================================
  • The difference in Switch=0 and 1 energies is almost entirely due to differences in E1elst. Which one is more accurate? They are both about 0.2% from the reference.
  • The reference induction energies are, as explained above, not really reference energies as there are errors due to the DF in this energies. In particular, the nuclear integrals were obtained with density-fitting.

So I'd say that the overall energy is rather good.

(B) MC type for main basis

Now let's look at the MC type. This should be more interesting. We do have a problem because the reference calculation had used Cartesian GTOs in the main basis. What I'd like to emphasise here is the stability of the different basis types. As we shall see, the stability of our results will also depend on the integral switch used.

First, SAPT(KS) energies. (kJ/mol).

  • Aux basis: MC type : Integral switch = 0
           E1elst     E1exch  E2ind(UC)  E2exind(UC) E2disp(UC) E2exdisp(UC)
=======================================================================================
Ref.       -28.8435   24.1824   -5.0944    2.0529   -9.2025    0.9758
---------------------------------------------------------------------------------------
aDZ        -26.6585   24.2225   -5.0145    2.0169   -9.1100    0.9630
aTZ        -28.4528   24.1657   -5.0624    2.0840   -9.1697    0.9824
aQZ        -28.6441   24.1766   -5.1033    2.1152   -9.1846    0.9864
JK-tzvpp   -28.9813   24.1948   -4.0452    1.3108   -8.5681    0.8010
=======================================================================================
  • Aux basis: MC type : Integral switch = 1
           E1elst     E1exch  E2ind(UC)  E2exind(UC) E2disp(UC) E2exdisp(UC)
=======================================================================================
Ref.       -28.8435   24.1824   -5.0944    2.0529   -9.2025    0.9758
---------------------------------------------------------------------------------------
aDZ        -36.8057   24.1126  -10.1633    2.6087   -9.1100    0.9598
aTZ        -29.1432   24.1849   -6.7141    2.4606   -9.1697    0.9813
aQZ        -29.5866   24.1755   -5.4282    2.2136   -9.1846    0.9859
JK-tzvpp   -34.7426   24.1369  -54.8832    9.1899   -8.5681    0.7985
=======================================================================================
  • Aux basis: DC type : Integral switch = 0
           E1elst     E1exch  E2ind(UC)  E2exind(UC) E2disp(UC) E2exdisp(UC)
=======================================================================================
Ref.       -28.8435   24.1824   -5.0944    2.0529   -9.2025    0.9758
---------------------------------------------------------------------------------------
aDZ        -29.0056   24.2002   -5.1613    2.1379   -9.1873    0.9860
aTZ        -29.7837   24.1658   -5.0887    2.1104   -9.1865    0.9881
aQZ        -29.1700   24.1750   -5.1621    2.1354   -9.1865    0.9887
JK-tzvpp   -29.0325   24.1952   -5.1972    2.1510   -9.1833    0.9834
=======================================================================================
  • Aux basis: DC type : Integral switch = 1
           E1elst     E1exch  E2ind(UC)  E2exind(UC) E2disp(UC) E2exdisp(UC)
=======================================================================================
Ref.       -28.8435   24.1824   -5.0944    2.0529   -9.2025    0.9758
---------------------------------------------------------------------------------------
aDZ        -29.0921   24.1837   -5.1580    2.1343   -9.1873    0.9864
aTZ        -29.0949   24.1769   -5.1574    2.1345   -9.1865    0.9883
aQZ        -29.0921   24.1837   -5.1580    2.1343   -9.1873    0.9864
JK-tzvpp   -29.0948   24.1844   -5.1573    2.1345   -9.1833    0.9842
=======================================================================================


And now SAPT(DFT) energies: (kJ/mol)

  • Aux basis: MC type : Integral switch = 0
           E1elst     E2ind      E2disp   E1exch    E2exind    E2exdisp    E2int
=======================================================================================
Ref.       -28.8435   -4.8671   -6.8557   24.1824    1.9613    0.7270     -13.6956 
---------------------------------------------------------------------------------------
aDZ        -26.6585   -4.8811   -6.7738   24.2225    1.9640    0.7160     -11.4109
aTZ        -28.4528   -4.9476   -6.8365   24.1657    2.0360    0.7324     -13.3027
aQZ        -28.6441   -4.9948   -6.8482   24.1766    2.0695    0.7355     -13.5056
JK-tzvpp   -28.9813   -3.8415   -6.3363   24.1948    1.2562    0.5923     -13.1158
=======================================================================================
  • Aux basis: MC type : Integral switch = 1
           E1elst     E2ind      E2disp   E1exch    E2exind    E2exdisp   E2int
=======================================================================================
Ref.       -28.8435   -4.8671   -6.8557   24.1824    1.9613    0.7270     -13.6956
---------------------------------------------------------------------------------------
aDZ        -36.8057   -4.8811   -6.7738   24.1126    1.2568    0.7136     -22.3775
aTZ        -29.1432   -4.9476   -6.8365   24.1849    1.8300    0.7316     -14.1808
aQZ        -29.5866   -4.9948   -6.8482   24.1755    2.0365    0.7351     -14.4825
JK-tzvpp   -34.7426   -3.8415   -6.3363   24.1369    0.6432    0.5905     -19.5499
=======================================================================================
  • Aux basis: DC type : Integral switch = 0
           E1elst     E2ind      E2disp   E1exch    E2exind    E2exdisp   E2int
=======================================================================================
Ref.       -28.8435   -4.8671   -6.8557   24.1824    1.9613    0.7270    -13.6956
---------------------------------------------------------------------------------------
aDZ        -29.0056   -5.0598   -6.8513   24.2002    2.0953    0.7353    -13.8860
aTZ        -29.7837   -4.9782   -6.8502   24.1658    2.0644    0.7368    -14.6452
aQZ        -29.1700   -5.0598   -6.8502   24.1750    2.0924    0.7372    -14.0754 
JK-tzvpp   -29.0325   -5.0986   -6.8474   24.1952    2.1094    0.7332    -13.9408
=======================================================================================
  • Aux basis: DC type : Integral switch = 1
           E1elst     E2ind      E2disp   E1exch    E2exind    E2exdisp   E2int
=======================================================================================
Ref.       -28.8435   -4.8671   -6.8557   24.1824    1.9613    0.7270    -13.6956
---------------------------------------------------------------------------------------
aDZ        -29.0921   -5.0598   -6.8513   24.1837    2.0930    0.7356    -13.9909
aTZ        -29.0949   -4.9782   -6.8502   24.1769    2.0598    0.7370    -13.9497
  *10Nov10 -29.0948   -5.0648   -6.8502   24.1769    2.0955    0.7370    -14.0006 <-- SVD instead of LU in DF
aQZ        -29.0948   -5.0598   -6.8502   24.1761    2.0934    0.7372    -13.9982
JK-tzvpp   -29.0948   -5.0986   -6.8474   24.1844    2.1095    0.7339    -14.0130
=======================================================================================
  • This is a lot of data, but what (I believe) is important is the very small variation in the energies when we use the DC/switch=1 combination. The interaction energy varies by less than 0.5% for the four auxiliary basis sets when we use this combination.
  • Compare this with the 15% variation if the MC/switch=0 combination is used. Even if we skip the aug-cc-pVDZ auxiliary basis (as it is small), the variation in E2int is still about 3% for this case.
  • But this still large error would be even larger if not for a cancellation in errors from the e-e and e-n integrals. There is no such cancellation in the DC/switch=1 case: here all integrals are accurate.

Conclusions?

  • So can we conclude that we have good evidence that the best way to calculate interaction energies is to use the DC auxiliary basis with Integral Switch = 1?
  • I think we have good reason to believe that this is the case, but we need to use an energy scan with multiple orientations to make sure.
  • What are the consequences of this observation for our calculations of density overlaps? The density-fitted density must introduce a good bit of noise into the overlap model. And I cannot advocate calculating the first-order energies using the present ENERGY-SCAN module. We probably need to use the DC auxiliary basis and repeat DF at each geometry. In which case, CamCASP needs to know how to handle DC auxiliary bases efficiently.
Retrieved from "http://wikis.ch.cam.ac.uk/cuc3/wiki/index.php?title=CamCASP/Numerical&oldid=3840"