Fix MCSCF Final CI Energy Mismatch

[lcyyork]

Description

Fix the mismatch issue on Github (#2258) and on the forum http://forum.psicode.org/t/slight-mismatching-of-casscf-between-total-mcscf-energy-and-mcscf-root-0-energy/1693/4.
A more efficient way to obtain the CI coefficients after a redundant orbital rotation would be a direct transformation without solving the CI again, see https://doi.org/10.1063/1.479573.

There might be two related issue:

• H0block_ in DETCI keeps filling up during the MCSCF macro iteration. This will lead to repeated entries during the MCSCF iterations and sometimes lead to uncertain numerical behavior.
• The "DFILE" CI guess for the CI after orbital canonicalization is directly related to this issue. Using this guess leads to very slow convergence for the CI.

Checklist

• Tests added for any new features
• All or relevant fraction of full tests run

Status

• Ready for review
• Ready for merge

[fevangelista]

The changes look reasonable to me. It would probably be good for @CDSherrill to take a look at these changes.

[CDSherrill]

I put these same comments on the forum.

believe @lcyyork is correct in identifying a problem with an ill-defined H0 (small block of H in the determinant basis, formed out of the ~400 “most important” determinants). The problem appears to be that this H0 block was computed in MCSCF macroiteration 0 but never updated… it should be updated each MCSCF macroiteration. I’m not sure why this didn’t cause even more serious MCSCF problems before now (not just for CAS/MCSCF canonical/semicanonical orbitals). Anyway, York’s PR recomputes H0 each macroiteration, which is a good improvement. This is especially true because some of the CI preconditioner options in DETCI use H0 explicitly.

York also modified the code to discard the previous best-guess for the CI vector when the new semicanonical orbitals are obtained. I don’t see why this was necessary. But maybe it is, for some subtle reason. At worst it increases the time to solution a little. If York has time to see what happens NOT including line 371 (ciwfn.set_ci_guess("H0_BLOCK")), I would be curious to see if this still works. (Possibly there are problems because the correction vector is basically zero?). But otherwise I think this is ok as-is.

[lcyyork]

@CDSherrill Yes, I did try that. Removing line 371 yields the original problem (using the input in the forum, CI space: 792 determinants):

   ==> Starting DF-MCSCF iterations <==

           Iter         Total Energy       Delta E   Orb RMS    CI RMS  NCI NORB
   @DF-MCSCF  1:   -112.452176153325   -3.9110e-01  3.88e-02  1.11e-11    1    1  Initial CI
      Warning! Maxstep = 0.10, scaling to 0.10
   @DF-MCSCF  2:   -112.515943250091   -6.3767e-02  1.81e-02  0.00e+00   12    1  TS
   @DF-MCSCF  3:   -112.522488962067   -6.5457e-03  6.14e-03  0.00e+00   12    1  TS
   @DF-MCSCF  4:   -112.523535480248   -1.0465e-03  2.40e-03  0.00e+00   12    1  TS
   @DF-MCSCF  5:   -112.523757306710   -2.2183e-04  1.06e-03  0.00e+00   12    1  TS
   @DF-MCSCF  6:   -112.523809657459   -5.2351e-05  4.44e-04  0.00e+00   12    1  TS
   @DF-MCSCF  7:   -112.523823448504   -1.3791e-05  3.18e-04  0.00e+00   12    1  TS
   @DF-MCSCF  8:   -112.523827480498   -4.0320e-06  1.33e-04  0.00e+00   12    1  TS
   @DF-MCSCF  9:   -112.523828813577   -1.3331e-06  1.79e-04  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 10:   -112.523829321491   -5.0791e-07  6.79e-05  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 11:   -112.523829546801   -2.2531e-07  1.38e-04  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 12:   -112.523829698540   -1.5174e-07  9.78e-05  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 13:   -112.523829801775   -1.0324e-07  1.30e-05  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 14:   -112.523829823173   -2.1397e-08  1.93e-05  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 15:   -112.523829813293    9.8798e-09  5.67e-05  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 16:   -112.523829832174   -1.8881e-08  2.90e-05  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 17:   -112.523829839177   -7.0027e-09  7.35e-06  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 18:   -112.523829840344   -1.1669e-09  3.02e-06  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 19:   -112.523829840592   -2.4842e-10  3.34e-07  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 20:   -112.523829840596   -3.6948e-12  2.10e-07  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 21:   -112.523829840597   -8.2423e-13  1.56e-07  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 22:   -112.523829840597   -6.2528e-13  7.02e-08  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 23:   -112.523829840598   -1.5632e-13  6.16e-08  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 24:   -112.523829840598   -2.8422e-14  4.59e-08  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 25:   -112.523829840598   -1.1369e-13  2.09e-08  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 26:   -112.523829840598    0.0000e+00  3.15e-09  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 27:   -112.523829840598    2.8422e-14  1.21e-09  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 28:   -112.523829840598   -4.2633e-14  2.51e-10  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 29:   -112.523829840598   -4.2633e-14  2.29e-10  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 30:   -112.523829840588    4.2633e-14  1.05e-10  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 31:   -112.523829840588   -2.8422e-14  2.75e-11  0.00e+00   12    1  TS, DIIS

          @DF-MCSCF has converged!

   @DF-MCSCF Final Energy: -112.523829840588462

   Computing CI Semicanonical Orbitals

   ==> Starting CI iterations <==

    H0 Block Eigenvalue = -112.52382984

    Simultaneous Expansion Method (Block Davidson Method)
    Attempting to use 1 previous converged vectors

     Iter   Root       Total Energy       Delta E      C RMS

   @CI  0:     0   -111.417799029788   -2.4659E+01   1.1817E+00  
   @CI  1:     0   -112.189379709492   -7.7158E-01   3.5706E-01  
   @CI  2:     0   -112.250278204452   -6.0898E-02   3.1809E-01  
   @CI  3:     0   -112.324883713784   -7.4606E-02   2.2146E-01  
   @CI  4:     0   -112.352310238893   -2.7427E-02   1.2530E-01  
   @CI  5:     0   -112.358631917780   -6.3217E-03   9.7660E-02  
   @CI  6:     0   -112.365095716692   -6.4638E-03   7.3202E-02  
   @CI  7:     0   -112.366948614338   -1.8529E-03   5.3732E-02  
   @CI  8:     0   -112.372824198774   -5.8756E-03   1.1022E-01  
   @CI  9:     0   -112.387582513115   -1.4758E-02   1.8337E-01  
   @CI 10:     0   -112.410692911567   -2.3110E-02   2.3866E-01  
   @CI 11:     0   -112.444565012608   -3.3872E-02   2.3335E-01  
   @CI 12:     0   -112.483195864884   -3.8631E-02   1.9437E-01  

    Warning! CI diagonalization did not fully converge!


   ==> Energetics <==

    SCF energy =         -112.061077881343110
    Total MCSCF energy = -112.523829840588462

   ==> MCSCF root 0 information <==

    MCSCF Root 0 energy =  -112.483195864884152

However, with line 371, I get:

   @DF-MCSCF 30:   -112.523829840588    4.2633e-14  1.05e-10  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 31:   -112.523829840588   -2.8422e-14  2.75e-11  0.00e+00   12    1  TS, DIIS

          @DF-MCSCF has converged!

   @DF-MCSCF Final Energy: -112.523829840588462

   Computing CI Semicanonical Orbitals

   ==> Starting CI iterations <==

    H0 Block Eigenvalue = -112.52382984

    Simultaneous Expansion Method (Block Davidson Method)
    Using 1 initial trial vectors

     Iter   Root       Total Energy       Delta E      C RMS

   @CI  0:     0   -112.523829840638   -2.5765E+01   4.3475E-11  
    Warning: Norm of correction (root 0) is < 1.0E-13
   @CI  1:     0   -112.523829840638   0.0000E+00   2.2840E-10 c

   ==> Energetics <==

    SCF energy =         -112.061077881343110
    Total MCSCF energy = -112.523829840588462

   ==> MCSCF root 0 information <==

    MCSCF Root 0 energy =  -112.523829840637774

I am not sure why reading the old guess leads to very slow convergence. It seems the old vector is not a good guess anymore: the energy from the first sigma build is way off. I actually wonder how the old vector remains to be a good guess during the MCSCF orbital optimization procedure.

There is a subtle issue using H0BLOCK for the final CASCI. The MCSCF sets CI max iteration to 12 by default, which is probably not enough to converge a large CI problem. This could be overwritten by the keyword CI_MAXITER. However, doing so will also increase the CI iterations during the MCSCF orbital optimization.