Multiconfigurational linear response#

Concept#

So-far we have seen CI response and how pure CASCI response can fail dramatically due to the neglect of orbital relaxation. We fixed it in a simple way using state-averaging, but this also introduced some difficulties, such as the dependency of the result on the number of states and the need to include all orbitals involved in the excitation in the active space.

These difficulties can be avoided altogether by deriving a proper linear response for MCSCF. For this, we introduce our parametrization:

\[\begin{equation*} \label{mclr-param} | \bar{\psi}(t) \rangle = e^{-i\hat{\kappa}(t)} e^{-i\hat{R}(t)} |0\rangle \end{equation*}\]

which contains both the orbital rotation operator similar to HF or DFT response:

\[\begin{equation*} \hat{\kappa}(t) = \sum_a^\mathrm{act,unocc} \sum_i^\mathrm{occ,act} \left[ \kappa_{ai}(t) \hat{a}^\dagger_a \hat{a}_i + \kappa_{ai}^*(t) \hat{a}^\dagger_i \hat{a}_a \right] , \end{equation*}\]

but also the CI state-transfer operator:

\[\begin{equation*} \hat{R}(t) = \sum_n R_{n0}(t) | n \rangle \langle 0 | + R^*_{0n}(t) | 0 \rangle \langle n | . \end{equation*}\]

For CASSCF wavefunctions, there is a redundancy between the active-active orbital response parameters and the CI response parameters, so these orbital response terms are dropped. The orbital rotations are thus restricted to inactive-active, inactive-virtual and active-virtual.

This theory can thus be seen as CASCI response but with added orbital excitations to and from non-active orbitals. Hence, it fixes the lack of orbital relaxation of CASCI and also lifts the need to include all orbitals involved in the excitations within our active space. The only difference between having the orbitals inside or outside the active space is the description of the reference state (typically the ground state) and the fact that the CI state-transfer includes all order of excitations while the orbital operators are restricted to single excitations.

Illustrative calculations#

Let us go back to our water example.

import multipsi as mtp
import veloxchem as vlx
* Warning * Environment variable OMP_NUM_THREADS not set.
* Warning * Setting OMP_NUM_THREADS to 4.
# CAS(6,5) calculation of water

mol_str = """
O       0.0000000000     0.0000000000     0.1178336003
H      -0.7595754146    -0.0000000000    -0.4713344012
H       0.7595754146     0.0000000000    -0.4713344012
"""

molecule = vlx.Molecule.read_str(mol_str, units="angstrom")
basis = vlx.MolecularBasis.read(molecule, "6-31g")

scf_drv = vlx.ScfRestrictedDriver()
scf_drv.compute(molecule, basis)

space = mtp.OrbSpace(molecule, scf_drv.mol_orbs)
space.CAS(6, 5)  # 3 O_2p and 2 H_1s
# CASSCF calculation
mcscf_drv = mtp.McscfDriver()
mcscf_drv.compute(molecule, basis, space)
* Info * Reading basis set from file: /home/thomas/Notebook/anaconda/envs/echem/lib/python3.9/site-packages/veloxchem/basis/6-31G
                                                                                                                          
                                              Molecular Basis (Atomic Basis)                                              
                                             ================================                                             
                                                                                                                          
                                  Basis: 6-31G                                                                            
                                                                                                                          
                                  Atom Contracted GTOs          Primitive GTOs                                            
                                                                                                                          
                                   O   (3S,2P)                  (10S,4P)                                                  
                                   H   (2S)                     (4S)                                                      
                                                                                                                          
                                  Contracted Basis Functions : 13                                                         
                                  Primitive Basis Functions  : 30                                                         
                                                                                                                          
                                                                                                                          
                                            Self Consistent Field Driver Setup                                            
                                           ====================================                                           
                                                                                                                          
                   Wave Function Model             : Spin-Restricted Hartree-Fock                                         
                   Initial Guess Model             : Superposition of Atomic Densities                                    
                   Convergence Accelerator         : Two Level Direct Inversion of Iterative Subspace                     
                   Max. Number of Iterations       : 50                                                                   
                   Max. Number of Error Vectors    : 10                                                                   
                   Convergence Threshold           : 1.0e-06                                                              
                   ERI Screening Scheme            : Cauchy Schwarz + Density                                             
                   ERI Screening Mode              : Dynamic                                                              
                   ERI Screening Threshold         : 1.0e-12                                                              
                   Linear Dependence Threshold     : 1.0e-06                                                              
                                                                                                                          
* Info * Nuclear repulsion energy: 9.1561447194 a.u.                                                                      
                                                                                                                          
* Info * Overlap matrix computed in 0.00 sec.                                                                             
                                                                                                                          
* Info * Kinetic energy matrix computed in 0.00 sec.                                                                      
                                                                                                                          
* Info * Nuclear potential matrix computed in 0.00 sec.                                                                   
                                                                                                                          
* Info * Orthogonalization matrix computed in 0.00 sec.                                                                   
                                                                                                                          
* Info * SAD initial guess computed in 0.00 sec.                                                                          
                                                                                                                          
* Info * Starting Reduced Basis SCF calculation...                                                                        
* Info * ...done. SCF energy in reduced basis set: -75.983870205311 a.u. Time: 0.05 sec.                                  
                                                                                                                          
* Info * Overlap matrix computed in 0.00 sec.                                                                             
                                                                                                                          
* Info * Kinetic energy matrix computed in 0.00 sec.                                                                      
                                                                                                                          
* Info * Nuclear potential matrix computed in 0.00 sec.                                                                   
                                                                                                                          
* Info * Orthogonalization matrix computed in 0.00 sec.                                                                   
                                                                                                                          
                                                                                                                          
               Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change               
               --------------------------------------------------------------------------------------------               
                  1       -75.983870373952    0.0000000000      0.00006826      0.00001638      0.00000000                
                  2       -75.983870375702   -0.0000000017      0.00002636      0.00000475      0.00006304                
                  3       -75.983870375765   -0.0000000001      0.00000396      0.00000061      0.00000524                
                  4       -75.983870375769   -0.0000000000      0.00000029      0.00000008      0.00000296                
                                                                                                                          
               *** SCF converged in 4 iterations. Time: 0.04 sec.                                                         
                                                                                                                          
               Spin-Restricted Hartree-Fock:                                                                              
               -----------------------------                                                                              
               Total Energy                       :      -75.9838703758 a.u.                                              
               Electronic Energy                  :      -85.1400150952 a.u.                                              
               Nuclear Repulsion Energy           :        9.1561447194 a.u.                                              
               ------------------------------------                                                                       
               Gradient Norm                      :        0.0000002874 a.u.                                              
                                                                                                                          
                                                                                                                          
               Ground State Information                                                                                   
               ------------------------                                                                                   
               Charge of Molecule            :  0.0                                                                       
               Multiplicity (2S+1)           :  1.0                                                                       
               Magnetic Quantum Number (M_S) :  0.0                                                                       
                                                                                                                          
                                                                                                                          
                                                 Spin Restricted Orbitals                                                 
                                                 ------------------------                                                 
                                                                                                                          
               Molecular Orbital No.   1:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:  -20.56128 a.u.                                                                  
               (   1 O   1s  :    -1.00)                                                                                  
                                                                                                                          
               Molecular Orbital No.   2:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -1.35434 a.u.                                                                  
               (   1 O   1s  :     0.21) (   1 O   2s  :    -0.47) (   1 O   3s  :    -0.48)                              
                                                                                                                          
               Molecular Orbital No.   3:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -0.70774 a.u.                                                                  
               (   1 O   1p+1:    -0.51) (   1 O   2p+1:    -0.27) (   2 H   1s  :     0.26)                              
               (   3 H   1s  :    -0.26)                                                                                  
                                                                                                                          
               Molecular Orbital No.   4:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -0.56024 a.u.                                                                  
               (   1 O   2s  :    -0.18) (   1 O   3s  :    -0.31) (   1 O   1p0 :    -0.55)                              
               (   1 O   2p0 :    -0.40)                                                                                  
                                                                                                                          
               Molecular Orbital No.   5:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -0.50122 a.u.                                                                  
               (   1 O   1p-1:     0.64) (   1 O   2p-1:     0.51)                                                        
                                                                                                                          
               Molecular Orbital No.   6:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.20280 a.u.                                                                  
               (   1 O   3s  :     1.17) (   1 O   1p0 :    -0.23) (   1 O   2p0 :    -0.48)                              
               (   2 H   2s  :    -0.99) (   3 H   2s  :    -0.99)                                                        
                                                                                                                          
               Molecular Orbital No.   7:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.29874 a.u.                                                                  
               (   1 O   1p+1:    -0.34) (   1 O   2p+1:    -0.82) (   2 H   2s  :    -1.38)                              
               (   3 H   2s  :     1.38)                                                                                  
                                                                                                                          
               Molecular Orbital No.   8:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    1.05548 a.u.                                                                  
               (   1 O   2p+1:    -0.74) (   2 H   1s  :    -0.97) (   2 H   2s  :     0.50)                              
               (   3 H   1s  :     0.97) (   3 H   2s  :    -0.50)                                                        
                                                                                                                          
               Molecular Orbital No.   9:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    1.16442 a.u.                                                                  
               (   1 O   1p-1:    -0.96) (   1 O   2p-1:     1.04)                                                        
                                                                                                                          
               Molecular Orbital No.  10:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    1.18237 a.u.                                                                  
               (   1 O   2s  :     0.26) (   1 O   3s  :    -0.18) (   1 O   1p0 :    -0.67)                              
               (   1 O   2p0 :     0.22) (   2 H   1s  :    -0.84) (   2 H   2s  :     0.55)                              
               (   3 H   1s  :    -0.84) (   3 H   2s  :     0.55)                                                        
                                                                                                                          
          Active space definition:
          ------------------------
Number of inactive (occupied) orbitals: 2
Number of active orbitals:              5
Number of virtual orbitals:             6

    This is a CASSCF wavefunction: CAS(6,5)

          CI expansion:
          -------------
Number of determinants:      55


                                                                                                                          
        MCSCF Iterations
        ----------------
                                                                                                                          
     Iter. | Average Energy | E. Change | Grad. Norm | CI Iter. |   Time
     ---------------------------------------------------------------------
        1     -75.991312026     0.0e+00      3.6e-02          1   0:00:00
        2     -76.006166338    -1.5e-02      6.3e-02          1   0:00:00
        3     -76.032227092    -2.6e-02      3.2e-02          1   0:00:00
        4     -76.034907536    -2.7e-03      3.7e-02          1   0:00:00
        5     -76.037509034    -2.6e-03      1.5e-02          1   0:00:00
        6     -76.037708221    -2.0e-04      7.3e-03          1   0:00:00
        7     -76.037802925    -9.5e-05      4.9e-03          1   0:00:00
        8     -76.037839963    -3.7e-05      3.1e-03          1   0:00:00
        9     -76.037842633    -2.7e-06      1.4e-03          1   0:00:00
       10     -76.037843292    -6.6e-07      3.9e-04          1   0:00:00
       11     -76.037843330    -3.8e-08      2.0e-04          1   0:00:00
       12     -76.037843340    -9.4e-09      5.7e-05          1   0:00:00
                                                                                                                          
** Convergence reached in 12 iterations
                                                                                                                          
        Final results
        -------------
                                                                                                                          
                                                                                                                          
* State 1
- Energy: -76.03784333970785
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99931 1.97806 1.97495 0.02400 0.02368
                                                                                                                          
                                                 Spin Restricted Orbitals                                                 
                                                 ------------------------                                                 
                                                                                                                          
               Molecular Orbital No.   1:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:  -20.55485 a.u.                                                                  
               (   1 O   1s  :    -1.00)                                                                                  
                                                                                                                          
               Molecular Orbital No.   2:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -1.07299 a.u.                                                                  
               (   1 O   1s  :     0.22) (   1 O   2s  :    -0.48) (   1 O   3s  :    -0.57)                              
               (   1 O   1p0 :    -0.24) (   1 O   2p0 :    -0.19)                                                        
                                                                                                                          
               Molecular Orbital No.   3:                                                                                 
               --------------------------                                                                                 
               Occupation: 1.999 Energy:   -0.50217 a.u.                                                                  
               (   1 O   1p-1:    -0.64) (   1 O   2p-1:    -0.51)                                                        
                                                                                                                          
               Molecular Orbital No.   4:                                                                                 
               --------------------------                                                                                 
               Occupation: 1.978 Energy:   -0.83593 a.u.                                                                  
               (   1 O   1p0 :     0.51) (   1 O   2p0 :     0.35) (   2 H   1s  :    -0.21)                              
               (   3 H   1s  :    -0.21)                                                                                  
                                                                                                                          
               Molecular Orbital No.   5:                                                                                 
               --------------------------                                                                                 
               Occupation: 1.975 Energy:   -0.70221 a.u.                                                                  
               (   1 O   1p+1:     0.51) (   1 O   2p+1:     0.25) (   2 H   1s  :    -0.27)                              
               (   3 H   1s  :     0.27)                                                                                  
                                                                                                                          
               Molecular Orbital No.   6:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.024 Energy:    0.75288 a.u.                                                                  
               (   1 O   1p+1:    -0.87) (   2 H   1s  :    -0.46) (   2 H   2s  :    -0.30)                              
               (   3 H   1s  :     0.46) (   3 H   2s  :     0.30)                                                        
                                                                                                                          
               Molecular Orbital No.   7:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.024 Energy:    0.73040 a.u.                                                                  
               (   1 O   2s  :     0.50) (   1 O   3s  :     0.17) (   1 O   1p0 :    -0.77)                              
               (   1 O   2p0 :     0.18) (   2 H   1s  :    -0.43) (   2 H   2s  :    -0.20)                              
               (   3 H   1s  :    -0.43) (   3 H   2s  :    -0.20)                                                        
                                                                                                                          
               Molecular Orbital No.   8:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.68041 a.u.                                                                  
               (   1 O   2s  :    -0.19) (   1 O   3s  :     1.21) (   1 O   1p0 :     0.41)                              
               (   1 O   2p0 :    -0.76) (   2 H   1s  :     0.40) (   2 H   2s  :    -1.13)                              
               (   3 H   1s  :     0.40) (   3 H   2s  :    -1.13)                                                        
                                                                                                                          
               Molecular Orbital No.   9:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.78480 a.u.                                                                  
               (   1 O   1p+1:    -0.28) (   1 O   2p+1:     0.73) (   2 H   1s  :    -0.62)                              
               (   2 H   2s  :     1.60) (   3 H   1s  :     0.62) (   3 H   2s  :    -1.60)                              
                                                                                                                          
               Molecular Orbital No.  10:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    1.16356 a.u.                                                                  
               (   1 O   1p-1:    -0.96) (   1 O   2p-1:     1.04)                                                        
                                                                                                                          
                                                                                                                          
 MTP execution time table
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
   Active integrals
           │     CI step
           │             │     Density matrices
           │             │             │     Active Fock
           │             │             │             │     Quasi-Newton
It.        │             │             │             │             │             │     Wrap-up
           │             │             │             │             │             │  
   Σ      0.0s (15%)    0.1s (33%)    0.0s ( 6%)    0.0s ( 6%)    0.0s (23%)    0.0s ( 8%)
                                                                                                                          
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
 Total time: 0.15 s
                                                                                                                          
Total MCSCF time: 00:00:00
mcrpa = mtp.Mclr_EigenSolver()
# We only need 4 states but use 5 states to help a bit with convergence
rsp_dict = mcrpa.compute(molecule, basis, mcscf_drv, nstates=5)
                                                                                                                          
Initialization time: 00:00:00
                                                                                                                          
          Linear eigensolver
          ------------------
                                                                                                                          
CI parameters:      54
Orbital parameters: 52
                                                                                                                          
                                                                                                                          
                                                                                                                          
        MC-RPA Iterations
        -----------------
                                                                                                                          
     Iter. | Average Energy | Grad. Norm | Converged | CI+Orb vec. | Time
     -----------------------------------------------------------------------
        1       0.647025237      2.9e-01      0/  5     5 +   5      0:00:00    
        2       0.513233576      1.1e-01      0/  5     5 +  15      0:00:00    
        3       0.473596395      4.4e-02      0/  5    13 +  17      0:00:00    
        4       0.458091234      1.9e-02      0/  5    15 +  24      0:00:00    
        5       0.453458989      3.1e-03      0/  5    23 +  26      0:00:00    
        6       0.452555176      1.4e-03      0/  5    25 +  34      0:00:00    
        7       0.452312403      2.9e-04      1/  5    31 +  38      0:00:00    
        8       0.452201158      6.9e-05      3/  5    35 +  41      0:00:00    
        9       0.452195626      2.6e-05      4/  5    37 +  43      0:00:00    
       10       0.452190637      9.9e-07      4/  5    39 +  43      0:00:00    
       11       0.452190496      2.0e-07      4/  5    39 +  45      0:00:00    
       12       0.452190468      5.5e-08      4/  5    41 +  45      0:00:00    
       13       0.452190465      5.3e-09      5/  5    41 +  47      0:00:00    
                                                                                                                          
** Convergence reached in 13 iterations
                                                                                                                          
                                                                                                                          
                                                                                                                          
                                                 Oscillator strength (au)
 Transition   Energy (Ha)    Energy (eV)    Velocity gauge     Length gauge 
  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
          1     0.3310416      9.00810          0.0389377        0.0122857
          2     0.4090370     11.13046          0.0000000        0.0000000
          3     0.4311669     11.73265          0.1736441        0.1284482
          4     0.5113192     13.91371          0.1444119        0.1449526
          5     0.5783876     15.73873          0.3471868        0.4072747
                                                                                                                          
                                                                                                                          
                                                 Rotatory strength (cgs)
 Transition   Energy (Ha)    Energy (eV)    Velocity gauge     Length gauge 
 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
          1     0.3310416      9.00810          0.0000000        0.0000000
          2     0.4090370     11.13046         -0.0000000       -0.0000000
          3     0.4311669     11.73265          0.0000000        0.0000000
          4     0.5113192     13.91371         -0.0000000        0.0000000
          5     0.5783876     15.73873         -0.0000000       -0.0000000
                                                                                                                          
 MTP execution time table
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
   Initialization
           │     Sigma
           │             │     Reduced equations
           │             │             │     New trial vectors
It.        │             │             │             │             │     Properties
           │             │             │             │             │  
   Σ      0.0s ( 1%)    0.3s (31%)    0.4s (40%)    0.2s (24%)    0.0s ( 1%)
                                                                                                                          
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
 Total time: 0.97 s
                                                                                                                          
Total MC-RPA time: 00:00:00
# Results computed in CI response and SA-MCSCF
FCI_excitation_energies = [8.44992, 10.68451, 10.97654, 13.37692]
CIS_excitation_energies = [9.38059, 11.30739, 11.83172, 13.90237]
CASCI_excitation_energies = [18.61974291, 18.73087038, 26.61579239, 28.51654574]
SA_CASSCF_excitation_energies = [7.55787583, 9.83566209, 10.32868402, 12.56323891]
au2ev = 27.211386
# Multiconfigurational linear response for CAS(6,5)
MCRPA_6_5_excitation_energies = au2ev * rsp_dict["eigenvalues"][:4]

print("Full CI Energies")
print(FCI_excitation_energies)
print("CIS Energies")
print(CIS_excitation_energies)
print("SA-CASSCF Energies")
print(SA_CASSCF_excitation_energies)
print("CASCI Energies")
print(CASCI_excitation_energies)
print("MC-RPA Energies")
print(MCRPA_6_5_excitation_energies)
Full CI Energies
[8.44992, 10.68451, 10.97654, 13.37692]
CIS Energies
[9.38059, 11.30739, 11.83172, 13.90237]
SA-CASSCF Energies
[7.55787583, 9.83566209, 10.32868402, 12.56323891]
CASCI Energies
[18.61974291, 18.73087038, 26.61579239, 28.51654574]
MC-RPA Energies
[ 9.0080996  11.13046476 11.73264782 13.91370526]

We can see that this multiconfigurational response results is significantly better than CASCI, and is actually even better than SA-CASSCF and CIS.

It is also useful to look at the output of the calculation to understand what is happening. In particular, the code prints the number of trial vectors that are used to converge the calculation. We can see that MC-RPA uses both CI vectors and orbital vectors, corresponding to the 2 parts of our parametrization. The final response vector is a linear combination of both. We can also print the excitation character analysis to see this:

mcrpa.get_excitation_character()
                                                                                                                          
    Excitation characters
    ---------------------
                                                                                                                          
* State 1
- Exc. energy: 0.33104155724950823
 Electrons     Excitation
 0.391    :    3  ->    8
 0.574    :    3  ->    7
 CI weight: 0.616  Orb weight: 0.384
                                                                                                                          
* State 2
- Exc. energy: 0.4090370390848041
 Electrons     Excitation
 0.275    :    3  ->    9
 0.664    :    3  ->    6
 CI weight: 0.700  Orb weight: 0.300
                                                                                                                          
* State 3
- Exc. energy: 0.431166858545419
 Electrons     Excitation
 0.279    :    4  ->    8
 0.208    :    2  ->    7
 0.158    :    2  ->    8
 0.332    :    4  ->    7
 CI weight: 0.335  Orb weight: 0.665
                                                                                                                          
* State 4
- Exc. energy: 0.5113192419996073
 Electrons     Excitation
 0.223    :    2  ->    6
 0.215    :    4  ->    9
 0.106    :    2  ->    9
 0.386    :    4  ->    6
 CI weight: 0.415  Orb weight: 0.585
                                                                                                                          
* State 5
- Exc. energy: 0.5783876272044596
 Electrons     Excitation
 0.437    :    5  ->    8
 0.519    :    5  ->    7
 CI weight: 0.541  Orb weight: 0.459
                                                                                                                          

The orbital analysis shows the weights of the CI and orbital contribution to each response vector. Even if all of those states correspond broadly to excitations between valence orbitals which are all in the active space, all contain a sizeable contribution to the orbital response, which explains why they were so poorly described by CASCI.

We can go even further and drop the oxygen lone pair from our active space. As seen in the occupation number of the CASSCF calculation, the first orbital (the lone pair) is almost perfectly doubly occupied, and was only added in our active space because it is the HOMO of our system and is thus excited in the lowest excited states. Since we now also include excitations outside of the active space, we can safely remove this orbital and focus on the 2 pairs of O-H \(\sigma\) and \(\sigma^*\) (4 electrons in 4 orbitals):

space.CAS(4, 4)  # Remove the first orbital
# CASSCF calculation
mcscf_drv = mtp.McscfDriver()
mcscf_drv.compute(molecule, basis, space)
          Active space definition:
          ------------------------
Number of inactive (occupied) orbitals: 3
Number of active orbitals:              4
Number of virtual orbitals:             6

    This is a CASSCF wavefunction: CAS(4,4)

          CI expansion:
          -------------
Number of determinants:      21


                                                                                                                          
        MCSCF Iterations
        ----------------
                                                                                                                          
     Iter. | Average Energy | E. Change | Grad. Norm | CI Iter. |   Time
     ---------------------------------------------------------------------
        1     -76.037271059     0.0e+00      2.6e-03          1   0:00:00
        2     -76.037273348    -2.3e-06      6.5e-04          1   0:00:00
        3     -76.037273801    -4.5e-07      3.2e-04          1   0:00:00
        4     -76.037273988    -1.9e-07      2.6e-04          1   0:00:00
        5     -76.037274022    -3.4e-08      1.6e-04          1   0:00:00
        6     -76.037274031    -8.7e-09      6.4e-05          1   0:00:00
                                                                                                                          
** Convergence reached in 6 iterations
                                                                                                                          
        Final results
        -------------
                                                                                                                          
                                                                                                                          
* State 1
- Energy: -76.03727403103228
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.97794 1.97490 0.02386 0.02330
                                                                                                                          
                                                 Spin Restricted Orbitals                                                 
                                                 ------------------------                                                 
                                                                                                                          
               Molecular Orbital No.   1:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:  -20.55475 a.u.                                                                  
               (   1 O   1s  :    -1.00)                                                                                  
                                                                                                                          
               Molecular Orbital No.   2:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -1.07554 a.u.                                                                  
               (   1 O   1s  :     0.22) (   1 O   2s  :    -0.49) (   1 O   3s  :    -0.58)                              
               (   1 O   1p0 :    -0.24) (   1 O   2p0 :    -0.19)                                                        
                                                                                                                          
               Molecular Orbital No.   3:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -0.50212 a.u.                                                                  
               (   1 O   1p-1:    -0.64) (   1 O   2p-1:    -0.51)                                                        
                                                                                                                          
               Molecular Orbital No.   4:                                                                                 
               --------------------------                                                                                 
               Occupation: 1.978 Energy:   -0.83329 a.u.                                                                  
               (   1 O   1p0 :     0.51) (   1 O   2p0 :     0.35) (   2 H   1s  :    -0.21)                              
               (   3 H   1s  :    -0.21)                                                                                  
                                                                                                                          
               Molecular Orbital No.   5:                                                                                 
               --------------------------                                                                                 
               Occupation: 1.975 Energy:   -0.70218 a.u.                                                                  
               (   1 O   1p+1:    -0.51) (   1 O   2p+1:    -0.25) (   2 H   1s  :     0.27)                              
               (   3 H   1s  :    -0.27)                                                                                  
                                                                                                                          
               Molecular Orbital No.   6:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.024 Energy:    0.75297 a.u.                                                                  
               (   1 O   1p+1:    -0.87) (   2 H   1s  :    -0.46) (   2 H   2s  :    -0.30)                              
               (   3 H   1s  :     0.46) (   3 H   2s  :     0.30)                                                        
                                                                                                                          
               Molecular Orbital No.   7:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.023 Energy:    0.73237 a.u.                                                                  
               (   1 O   2s  :    -0.49) (   1 O   3s  :    -0.17) (   1 O   1p0 :     0.78)                              
               (   1 O   2p0 :    -0.18) (   2 H   1s  :     0.43) (   2 H   2s  :     0.20)                              
               (   3 H   1s  :     0.43) (   3 H   2s  :     0.20)                                                        
                                                                                                                          
               Molecular Orbital No.   8:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.67755 a.u.                                                                  
               (   1 O   2s  :    -0.18) (   1 O   3s  :     1.21) (   1 O   1p0 :     0.41)                              
               (   1 O   2p0 :    -0.76) (   2 H   1s  :     0.40) (   2 H   2s  :    -1.13)                              
               (   3 H   1s  :     0.40) (   3 H   2s  :    -1.13)                                                        
                                                                                                                          
               Molecular Orbital No.   9:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.78413 a.u.                                                                  
               (   1 O   1p+1:    -0.27) (   1 O   2p+1:     0.72) (   2 H   1s  :    -0.62)                              
               (   2 H   2s  :     1.60) (   3 H   1s  :     0.62) (   3 H   2s  :    -1.60)                              
                                                                                                                          
               Molecular Orbital No.  10:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    1.16377 a.u.                                                                  
               (   1 O   1p-1:    -0.96) (   1 O   2p-1:     1.04)                                                        
                                                                                                                          
                                                                                                                          
 MTP execution time table
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
   Active integrals
           │     CI step
           │             │     Density matrices
           │             │             │     Active Fock
           │             │             │             │     Quasi-Newton
It.        │             │             │             │             │             │     Wrap-up
           │             │             │             │             │             │  
   Σ      0.0s (12%)    0.0s (34%)    0.0s ( 2%)    0.0s ( 3%)    0.0s (18%)    0.0s (15%)
                                                                                                                          
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
 Total time: 0.08 s
                                                                                                                          
Total MCSCF time: 00:00:00
mcrpa = mtp.Mclr_EigenSolver()
rsp_dict = mcrpa.compute(molecule, basis, mcscf_drv, nstates=4)
                                                                                                                          
Initialization time: 00:00:00
                                                                                                                          
          Linear eigensolver
          ------------------
                                                                                                                          
CI parameters:      20
Orbital parameters: 54
                                                                                                                          
                                                                                                                          
                                                                                                                          
        MC-RPA Iterations
        -----------------
                                                                                                                          
     Iter. | Average Energy | Grad. Norm | Converged | CI+Orb vec. | Time
     -----------------------------------------------------------------------
        1       0.710081794      2.4e-01      0/  4     4 +   4      0:00:00    
        2       0.504438122      7.0e-02      0/  4     4 +  12      0:00:00    
        3       0.460463894      4.4e-02      1/  4     4 +  20      0:00:00    
        4       0.438810895      6.4e-03      1/  4     6 +  24      0:00:00    
        5       0.436944673      2.4e-03      2/  4     8 +  27      0:00:00    
        6       0.423462081      1.1e-02      2/  4     8 +  31      0:00:00    
        7       0.421356042      3.8e-03      2/  4    10 +  33      0:00:00    
        8       0.420278497      4.0e-04      2/  4    12 +  35      0:00:00    
        9       0.420223667      3.4e-05      2/  4    12 +  39      0:00:00    
       10       0.420210977      2.3e-06      2/  4    16 +  39      0:00:00    
       11       0.420210806      8.3e-08      2/  4    16 +  43      0:00:00    
       12       0.420210790      6.9e-09      4/  4    18 +  45      0:00:00    
                                                                                                                          
** Convergence reached in 12 iterations
                                                                                                                          
                                                                                                                          
                                                                                                                          
                                                 Oscillator strength (au)
 Transition   Energy (Ha)    Energy (eV)    Velocity gauge     Length gauge 
  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
          1     0.3305460      8.99462          0.0389589        0.0121782
          2     0.4084733     11.11512          0.0000000        0.0000000
          3     0.4307439     11.72114          0.1738522        0.1284474
          4     0.5110800     13.90719          0.1450456        0.1453738
                                                                                                                          
                                                                                                                          
                                                 Rotatory strength (cgs)
 Transition   Energy (Ha)    Energy (eV)    Velocity gauge     Length gauge 
 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
          1     0.3305460      8.99462          0.0000000        0.0000000
          2     0.4084733     11.11512         -0.0000000       -0.0000000
          3     0.4307439     11.72114          0.0000000        0.0000000
          4     0.5110800     13.90719          0.0000000        0.0000000
                                                                                                                          
 MTP execution time table
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
   Initialization
           │     Sigma
           │             │     Reduced equations
           │             │             │     New trial vectors
It.        │             │             │             │             │     Properties
           │             │             │             │             │  
   Σ      0.0s ( 3%)    0.2s (42%)    0.2s (34%)    0.1s (16%)    0.0s ( 1%)
                                                                                                                          
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
 Total time: 0.55 s
                                                                                                                          
Total MC-RPA time: 00:00:00
MCRPA_4_4_excitation_energies = au2ev * rsp_dict["eigenvalues"]

print("Full CI Energies")
print(FCI_excitation_energies)
print("CIS Energies")
print(CIS_excitation_energies)
print("SA-CASSCF Energies")
print(SA_CASSCF_excitation_energies)
print("MC-RPA Energies CAS(6,5)")
print(MCRPA_6_5_excitation_energies)
print("MC-RPA Energies CAS(4,4)")
print(MCRPA_4_4_excitation_energies)
Full CI Energies
[8.44992, 10.68451, 10.97654, 13.37692]
CIS Energies
[9.38059, 11.30739, 11.83172, 13.90237]
SA-CASSCF Energies
[7.55787583, 9.83566209, 10.32868402, 12.56323891]
MC-RPA Energies CAS(6,5)
[ 9.0080996  11.13046476 11.73264782 13.91370526]
MC-RPA Energies CAS(4,4)
[ 8.99461578 11.11512359 11.72113816 13.90719454]

As expected, since this orbital was almost perfectly doubly occupied, removing it from the active space made almost no difference in the excitation energies (this would not be true in SA-CASSCF).

As a result, the MC-RPA calculation used fewer CI trial vectors, and we can also see it in the excitation character analysis where the first 2 states are now purely described by the orbital response parameters:

mcrpa.get_excitation_character()
                                                                                                                          
    Excitation characters
    ---------------------
                                                                                                                          
* State 1
- Exc. energy: 0.330546036159546
 Electrons     Excitation
 0.578    :    3  ->    7
 0.393    :    3  ->    8
 CI weight: 0.000  Orb weight: 1.000
                                                                                                                          
* State 2
- Exc. energy: 0.40847326146200813
 Electrons     Excitation
 0.670    :    3  ->    6
 0.275    :    3  ->    9
 CI weight: 0.000  Orb weight: 1.000
                                                                                                                          
* State 3
- Exc. energy: 0.43074388637109207
 Electrons     Excitation
 0.282    :    4  ->    8
 0.205    :    2  ->    7
 0.157    :    2  ->    8
 0.332    :    4  ->    7
 CI weight: 0.335  Orb weight: 0.665
                                                                                                                          
* State 4
- Exc. energy: 0.5110799773765028
 Electrons     Excitation
 0.221    :    2  ->    6
 0.216    :    4  ->    9
 0.105    :    2  ->    9
 0.388    :    4  ->    6
 CI weight: 0.417  Orb weight: 0.583
                                                                                                                          

Finally, as for HF or DFT response, we can invoke the Tamm-Dancoff approximation (TDA), which typically gives rather moderate error compared to the full RPA response although leads to only modest gains in computational time:

mcrpa = mtp.Mclr_EigenSolver()
mcrpa.TDA = True
rsp_dict = mcrpa.compute(molecule, basis, mcscf_drv, nstates=4)
                                                                                                                          
Initialization time: 00:00:00
                                                                                                                          
          Linear eigensolver
          ------------------
                                                                                                                          
Tamm-Dancoff approximation
CI parameters:      20
Orbital parameters: 54
                                                                                                                          
                                                                                                                          
                                                                                                                          
        MC-TDA Iterations
        -----------------
                                                                                                                          
     Iter. | Average Energy | Grad. Norm | Converged | CI+Orb vec. | Time
     -----------------------------------------------------------------------
        1       0.710206962      9.4e-01      0/  4     4 +   4      0:00:00    
        2       0.505692806      2.7e-01      0/  4     4 +   8      0:00:00    
        3       0.467944503      1.6e-01      0/  4     4 +  12      0:00:00    
        4       0.447285216      1.3e-02      1/  4     5 +  15      0:00:00    
        5       0.445600907      6.4e-03      1/  4     6 +  17      0:00:00    
        6       0.434444757      3.1e-02      2/  4     6 +  20      0:00:00    
        7       0.432815917      1.1e-02      2/  4     7 +  21      0:00:00    
        8       0.431130278      9.6e-04      2/  4     8 +  22      0:00:00    
        9       0.430927141      1.3e-04      2/  4     9 +  23      0:00:00    
       10       0.430913653      2.4e-05      2/  4    10 +  24      0:00:00    
       11       0.430909925      3.1e-06      2/  4    10 +  10      0:00:00    Orb subspace collapsed 
       12       0.430909829      8.2e-07      3/  4    11 +  11      0:00:00    
       13       0.430909761      1.5e-07      3/  4    12 +  11      0:00:00    
       14       0.430909758      9.2e-09      4/  4    12 +  12      0:00:00    
                                                                                                                          
** Convergence reached in 14 iterations
                                                                                                                          
                                                                                                                          
                                                                                                                          
                                                 Oscillator strength (au)
 Transition   Energy (Ha)    Energy (eV)    Velocity gauge     Length gauge 
  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
          1     0.3459476      9.41371          0.0403390        0.0142156
          2     0.4209270     11.45401          0.0000000        0.0000000
          3     0.4376897     11.91014          0.1841958        0.1221934
          4     0.5190748     14.12474          0.1187558        0.1440216
                                                                                                                          
                                                                                                                          
                                                 Rotatory strength (cgs)
 Transition   Energy (Ha)    Energy (eV)    Velocity gauge     Length gauge 
 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
          1     0.3459476      9.41371         -0.0000000       -0.0000000
          2     0.4209270     11.45401          0.0000000        0.0000000
          3     0.4376897     11.91014          0.0000000        0.0000000
          4     0.5190748     14.12474          0.0000000        0.0000000
                                                                                                                          
 MTP execution time table
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
   Initialization
           │     Sigma
           │             │     Reduced equations
           │             │             │     New trial vectors
It.        │             │             │             │             │     Properties
           │             │             │             │             │  
   Σ      0.0s ( 3%)    0.2s (52%)    0.0s (11%)    0.1s (23%)    0.0s ( 2%)
                                                                                                                          
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
 Total time: 0.44 s
                                                                                                                          
Total MC-TDA time: 00:00:00
MCTDA_4_4_excitation_energies = au2ev * rsp_dict["eigenvalues"]

print("Full CI Energies")
print(FCI_excitation_energies)
print("MC-RPA Energies CAS(4,4)")
print(MCRPA_4_4_excitation_energies)
print("MC-TDA Energies CAS(4,4)")
print(MCTDA_4_4_excitation_energies)
Full CI Energies
[8.44992, 10.68451, 10.97654, 13.37692]
MC-RPA Energies CAS(4,4)
[ 8.99461578 11.11512359 11.72113816 13.90719454]
MC-TDA Energies CAS(4,4)
[ 9.41371372 11.45400689 11.91014268 14.12474375]