CamCASP/Bugs/9

CamCASP => Bugs => Errors in Second-order Energies

Background

Way back in 2010, while developing a potential for molecule XVI (part of the 2010 Blind test), we had encountered large errors in the second-order energies for dimers of this molecule. Sally had alerted us to this (she was developing the potential with my help). She had realised there was a problem as there were significant outliers in the total energies and dispersion energies when compared with our fits. Closer investigation revealed that the contributions to the dispersion energy at imaginary frequency were erratic, and in some cases, even changed sign (they should all be of one sign and decay as <math> 1/(w^2+w_0^2)</math>), where <math>w_0</math> is a characteristic frequency). At that point we brushed aside the problem as there wasn't time to investigate it. But we knew that it was related to the new propagator module that allows us to calculate the FDDS entirely within CamCASP using the ALDAX+CHF model.

Recently, I have encountered the problem yet again, this time with the pyridine dimer. As before, the calculations involve the new propagator module with the ALDAX+CHF model. Unfortunately, this seems to be a large-system problem. Calculations on smaller systems or large systems with smaller (less-diffuse) basis sets have not exhibited the problem. Nor has there been a problem with calculations using the ALDAX propagator.

Here I present a summary of my investigations with a few pyridine dimer systems.

Details

Unless stated, the calculations use the parameters defined here.

CamCASP:

 trunk  rev 22218
 

DALTON

 Version 2 with patches from Rafal. 
 Also includes changes suggested by Anthony to get it to compile with ifort. 
 See CamCASP/Compilation  for details.

Default run parameters:

 Integral switch = 1
 PBE0
 ALDAX+CHF
 Standard integration grids

Machines used:

 Tati:
    ifort Version 11.1
    BLAS/LAPACK: Goto2 
 mmp12:
    ifort Version 12.0.5
    BLAS/LAPACK: ATLAS 3.8.4 with LAPACK additions from Netlib
       Compiled with Gfortran 4.5 or 4.6 (probably 4.5)
    Gfortran: Versions 4.5.4 and 4.6.1

Example 1: Pyridine

Click here for a generic CLUSTER file

Sadlej/MC+/AC

Second-order energies

mmp12:
     ALDA+CHF :  -1628.32   -2316.07
     ALDAX+CHF:  -1546.89   -2324.97
     ALDAX    :  -1460.61   -2065.68
    ---------------------------------
tati*:
     ALDAX+CHF:  -311.77    -1138.50
    ---------------------------------
    * Original run: PBE0/AC Integral switch = 0 ifort 11.1
      CamCASP: Compiled >= Aug 05 2011  Rev. 22175

• The calculations on mmp12 are sensible.
  • E2disp: Very good agreement between ALDA+CHF and ALDAX+CHF (0.4%). Errors large for ALDAX, but this is to be expected as we have used PBE0 for MOs, so we must have 25% of CHF.
  • E2ind: This time ALDAX+CHF is not a good approximation. It turns out that the total induction energy from ALDA+CHF and ALDAX+CHF is not too bad: -549.63 cm-1 from ALDA+CHF and -522.14 cm-1 from ALDAX+CHF. A -5% error with ALDAX+CHF. I had overlooked this in my tests of ALDAX+CHF.
• The results on Tati are clearly nonsense. Here are some details of that run:

DATADIR: CamCASP/tests/bugs/kernel-bugs/bug-1-aldx-chf/pyr2/pyr2-sadlej-mcP-original/

File: $DATADIR/OUT/pot1721-Min13-s0.95-aldax-chf.out

  INDUCTION CALCULATION

  Parameters same. Skipping DF for pyr_1
OV DF finished for pyr_1
 E^{2}_{ind}(UC) =     -1634.523269 CM-1
     E^{2}_{ind}(UC)(A) =      -817.260620 CM-1
     E^{2}_{ind}(UC)(B) =      -817.262648 CM-1
 E^{2}_{ind} =      -311.774657 CM-1
     E^{2}_{ind}(A) =      -194.440126 CM-1
     E^{2}_{ind}(B) =      -117.334532 CM-1
  Density-fitted propagator was used.
  Propagator type : cks
 ---------------------------------

 E^{2}_{ind,exch} =       207.525965 CM-1
     E^{2}_{ind,exch}(A) =       129.424747 CM-1
     E^{2}_{ind,exch}(B) =        78.101219 CM-1
 E^{2}_{ind,exch}(UC) =      1087.984745 CM-1
     E^{2}_{ind,exch}(UC)(A) =       543.991362 CM-1
     E^{2}_{ind,exch}(UC)(B) =       543.993383 CM-1
 ---------------------------------

 ---------------------------------
 Quad points and XI at each point:
 i     OMEGA             XI
 ---------------------------------
    1      0.006610   -0.903720E-02
    2      0.036175    0.100213E+00
    3      0.095447   -0.197757E-01
    4      0.197644   -0.154010E-01
    5      0.370417   -0.103565E-01
    6      0.674915   -0.495012E-02
    7      1.264899   -0.114417E-02
    8      2.619245   -0.200620E-03
    9      6.910886   -0.698310E-05
   10     37.823762   -0.109340E-07
 ---------------------------------
 E^{2}_{disp} =     -1138.500450 CM-1
 E^{2}_{disp}(UC) =     -3550.241639 CM-1
     Density-fitted propagator was used.
     Propagator type : cks
 ---------------------------------

 E^{2}_{disp,exch} =       180.486265 CM-1
 E^{2}_{disp,exch}(UC) =       562.819148 CM-1
 ---------------------------------

It is worth looking at the details of this calculation:

• XI(w) should be monotonic with w. It is not.
• The induction energy in the coupled and un-coupled approximations are typically similar for systems like pyridine with a moderately sized HOMO-LUMO gap, but they are dramatically off.
• Symmetry has broken. This is a <math>D_{2h}</math> symmetry dimer. E2ind(A) and E2ind(B) should be the same. The uncoupled energies are, but not the coupled.

PUZZLE 1: The last observation is puzzling. Both the coupled and uncoupled expressions for the induction involve virtual orbitals. why has only the coupled calculation broken symmetry?

First-order energies

Here are first-order energies from runs on the two machines. They should be identical. Same inputs/settings. Only compilers and libraries differ.

          ==========================
                   E1elst    E1exch
          --------------------------
          mmp12   -2792.27   4812.59
          Tati    -2798.89   4947.78
          ==========================

How do we understand this? Look at the DALTON DFT energies...

DALTON energies

                  MonA      MonB
          --------------------------
          mmp12  -248.0367  -248.0367 (yes, this is a symmetric config)
          tati   -249.3680  -249.3680 (D2h symm so monomers are identical)
          ==========================

This is amazing. The DFT solutions on these two machines are really very different. And these results are reproducible. Possible causes:

• Basis too diffuse/large. Uses mid-bonds. Perhaps this leads to linear-dependencies and oddities in the MOs. But how does this explain PUZZLE 1?
• DALTON parameters (in the *.dal files) could be the cause of the problem. We should allow orbital-space truncation by reducing .CMOMAX and increasing .AO DELETE.
• Compile on Tati with the new version of ifort and the ATLAS library.

Cross-runs: MOs from Tati with CamCASP on mmp12

As a cross-check I ran CamCASP on mmp12 using MOs from Tati. If there was a problem with the MOs from Tati I should have obtained the same CamCASP results on mmp12 as on Tati. This is the case:

Results in
           /bug-1-aldx-chf/pyr2/pyr2-sadlej-mcP-original/
          Two runs: 2GB and 4GB memory.
                E2ind    E2disp
          ----------------------------
          2GB   -307.57  -1050.56
          4GB   -307.57  -1050.56
          ----------------------------

I ran this job with two memory settings to ensure that this was not a large-job-small-memory problem as we have encountered on some occasions.