Slight mismatching of CASSCF between "Total MCSCF energy" and "MCSCF root 0 energy"

1

Hello,

I made a strange observation that there is a small discrepancy between the printed energies and psi4’s CURRENT ENERGY in CASSCF. The following is the input I use.

molecule pbenzyne{
  0 1
   C   0.0000000000  -2.5451795941   0.0000000000
   C   0.0000000000   2.5451795941   0.0000000000
   C  -2.2828001669  -1.3508352528   0.0000000000
   C   2.2828001669  -1.3508352528   0.0000000000
   C   2.2828001669   1.3508352528   0.0000000000
   C  -2.2828001669   1.3508352528   0.0000000000
   H  -4.0782187459  -2.3208602146   0.0000000000
   H   4.0782187459  -2.3208602146   0.0000000000
   H   4.0782187459   2.3208602146   0.0000000000
   H  -4.0782187459   2.3208602146   0.0000000000

  units bohr
}

set globals{
  basis                  cc-pvdz
  df_basis_scf           cc-pvdz-jkfit
  reference              rhf
  scf_type               df
  d_convergence          8
  e_convergence          12
  mcscf_type             df
  mcscf_e_convergence    12
  mcscf_r_convergence    10
  mcscf_diis_start       8
  mcscf_max_rot          0.1
  docc                   [5,3,1,1,0,1,5,4]
  restricted_docc        [5,3,0,0,0,0,4,4]
  active                 [1,0,1,2,1,2,1,0]
}
Ecas, wfn = energy('casscf', return_wfn=True)

Here is the part of the output that bothers me.

...
   @DF-MCSCF 16:   -229.495450424074   -1.3379e-09  7.13e-07  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 17:   -229.495450424145   -7.0315e-11  8.76e-10  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 18:   -229.495450424145    5.6843e-14  3.93e-10  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 19:   -229.495450424145   -5.6843e-14  6.15e-11  0.00e+00   12    1  TS, DIIS

          @DF-MCSCF has converged!

   @DF-MCSCF Final Energy: -229.495450424144593

   Computing CI Semicanonical Orbitals

   ==> Energetics <==

    SCF energy =         -229.283874146841754
    Total MCSCF energy = -229.495450424144593

   ==> MCSCF root 0 information <==

    MCSCF Root 0 energy =  -229.495450379691391

The Total MCSCF energy differs from MCSCF Root 0 energy by 4 * 10^-8 and the CURRENT ENERGY uses MCSCF Root 0 energy. It seems to me the problem is the integral transformation step to semicanonical orbitals and re-diagonalizing the Hamiltonian (the latter is probably not necessary), but I am not sure. I tried other molecules (e.g., NH3, N2) but I have not found similar behaviors. Any ideas why this is happening?

Edits: I am using Psi4 1.4a2.dev313.

2

@francesco anything that you can think of with the semicanonical orbitals here?

3

Thank you for bringing this up @york0822: I am not sure where the discrepancy comes from. Let me investigate the code that saves the energy to “MCSCF Root 0 energy.”

4

I found the same problem, but this time the difference is as large as 40 mHa.

Here’s the input file:

molecule mol {
0 1
C  0.0 0.0 0.0
O  0.0 0.0 2.7
}

set {
  ints_tolerance 0.0
  guess sad
  basis cc-pvqz #cc-pvdz
  scf_type df
  maxiter 500
  reference rhf
  restricted_docc [ 2, 0, 0, 0 ]
  active          [ 4, 0, 2, 2 ]
  num_roots 1
  mcscf_maxiter 100
  d_convergence          8
  e_convergence          12
  mcscf_type             df
  mcscf_e_convergence    12
  mcscf_r_convergence    10
  mcscf_diis_start       8
  mcscf_max_rot          0.1
}
en=energy('casscf')

Here’s the output:

   @DF-MCSCF 29:   -112.523829840602    2.8422e-14  2.29e-10  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 30:   -112.523829840602    4.2633e-14  1.05e-10  0.00e+00   12    1  TS, DIIS
   @DF-MCSCF 31:   -112.523829840602   -5.6843e-14  2.78e-11  0.00e+00   12    1  TS, DIIS

          @DF-MCSCF has converged!

   @DF-MCSCF Final Energy: -112.523829840601877

   Computing CI Semicanonical Orbitals

   ==> Energetics <==

    SCF energy =         -112.061077881342072
    Total MCSCF energy = -112.523829840601877

   ==> MCSCF root 0 information <==

    MCSCF Root 0 energy =  -112.483195864901745

Many nearby geometries (r = 2.5 - 3.0) have similar issue, r = 2.7 is the worst.
My psi4 version should be 1.4a3.dev18, (commit a610e3b7bee8b1937f175843d7d4591e2866e983)

5

Thanks for the report. I’ve filed an issue. I’m currently working on some other issues, but MCSCF is important for the group I’m in, so I’ll be sure to give it a look when my schedule clears up.

6

The error comes from the last CI diagonalization step, which is not converged due to an ill-defined H0 for Davidson-Liu. I am not sure how to fix it yet but a workaround for the calculation is to do

Ecas, wfn = energy('casscf', return_wfn=True)
en = energy('detci', ref_wfn=wfn)
7

If you are in a hurry, few lines can be added/modified to your local branch of Psi4 according to this PR https://github.com/psi4/psi4/pull/2259. I tested locally that your example returns the correct energy.

8

I believe York 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.