Linear response

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 parameterization

\[ | \bar{\psi}(t) \rangle = e^{-i\hat{\kappa}(t)} e^{-i\hat{R}(t)} |0\rangle \]

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

\[ \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] \]

but also the CI state-transfer operator

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

For CASSCF wave functions, 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
mol_str = """3

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_xyz_string(mol_str)
basis = vlx.MolecularBasis.read(molecule, "6-31g")
scf_drv = vlx.ScfRestrictedDriver()
scf_results = scf_drv.compute(molecule, basis)
import py3Dmol as p3d
viewer = p3d.view(width=300, height=200)
viewer.addModel(mol_str, 'xyz')
viewer.setViewStyle({"style": "outline", "width": 0.05})
viewer.setStyle({"stick":{},"sphere": {"scale":0.25}})
viewer.rotate(-90,'x')
viewer.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

# CAS(6,5) calculation of water
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_results = mcscf_drv.compute(molecule, basis, space)
Hide code cell output
                                                                                                                          
                          Multi-Configurational Self-Consistent Field Driver
                         ====================================================
                                                                                                                          
        ╭────────────────────────────────────╮
        │          Driver settings           │
        ╰────────────────────────────────────╯
          State-specific calculation
          Max. iterations         : 50
          BFGS window             : 5
          Convergence thresholds:
            - Energy change       : 1e-08
            - Gradient norm       : 0.0001
                                                                                                                          
          Integrals in memory
                                                                                                                          
          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
       13     -76.037843341    -1.1e-09      1.6e-05          1   0:00:00
                                                                                                                          
        Final results
        -------------
                                                                                                                          
                                                                                                                          
* State 1
- Energy: -76.03784334082589
- 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.83591 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.75289 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.73042 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.68038 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.78478 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)                              
                                                                                                                          
                                                                                                                          
                                                Dipole moment for state 1                                                 
                                               ---------------------------                                                
                                                                                                                          
                                   X   :         0.000000 a.u.         0.000000 Debye                                     
                                   Y   :        -0.000000 a.u.        -0.000000 Debye                                     
                                   Z   :        -0.981834 a.u.        -2.495574 Debye                                     
                                 Total :         0.981834 a.u.         2.495574 Debye                                     
                                                                                                                          
                                                                                                                          
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, n_states=5)
Hide code cell output
                                                                                                                          
                          Multi-Configurational Linear response Eigenvalue Solver
                         =========================================================
                                                                                                                          
        ╭────────────────────────────────────╮
        │          Driver settings           │
        ╰────────────────────────────────────╯
          Number of states        : 5
                                                                                                                          
          Max. iterations         : 50
          Residual norm threshold : 1e-08
                                                                                                                          
          Integrals in memory
          Trial vector settings:
            - Linear dependence threshold :1e-10
            - Min. trial vectors          :20
            - Max. trial vectors          :60
                                                                                                                          
          CI parameters:      54
          Orbital parameters: 52
                                                                                                                          
                                                                                                                          
Initialization time: 00:00:00
                                                                                                                          
                                                                                                                          
        MC-RPA Iterations
        -----------------
                                                                                                                          
     Iter. | Average Energy | Grad. Norm | Converged | CI+Orb vec. | Time
     -----------------------------------------------------------------------
        1       0.647019008      2.9e-01      0/  5     5 +   5      0:00:00    
        2       0.513232268      1.1e-01      0/  5     5 +  15      0:00:00    
        3       0.473594122      4.4e-02      0/  5    13 +  17      0:00:00    
        4       0.458084076      1.9e-02      0/  5    15 +  24      0:00:00    
        5       0.453452539      3.1e-03      0/  5    23 +  26      0:00:00    
        6       0.452548808      1.4e-03      0/  5    25 +  34      0:00:00    
        7       0.452313010      2.9e-04      1/  5    31 +  38      0:00:00    
        8       0.452201613      7.0e-05      3/  5    35 +  41      0:00:00    
        9       0.452196083      2.6e-05      4/  5    37 +  43      0:00:00    
       10       0.452191087      9.9e-07      4/  5    39 +  43      0:00:00    
       11       0.452190945      2.0e-07      4/  5    39 +  45      0:00:00    
       12       0.452190916      5.1e-08      4/  5    41 +  45      0:00:00    
       13       0.452190914      4.0e-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.3310434      9.00815          0.0389358        0.0122841
          2     0.4090433     11.13064          0.0000000        0.0000000
          3     0.4311622     11.73252          0.1736409        0.1284413
          4     0.5113197     13.91372          0.1444330        0.1449627
          5     0.5783860     15.73868          0.3471652        0.4072428
                                                                                                                          
                                                                                                                          
                                                 Rotatory strength (cgs)
 Transition   Energy (Ha)    Energy (eV)    Velocity gauge     Length gauge 
 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
          1     0.3310434      9.00815          0.0000000        0.0000000
          2     0.4090433     11.13064          0.0000000        0.0000000
          3     0.4311622     11.73252         -0.0000000       -0.0000000
          4     0.5113197     13.91372          0.0000000       -0.0000000
          5     0.5783860     15.73868          0.0000000        0.0000000
                                                                                                                          
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.00814853 11.1306361  11.73252108 13.91371718]

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.33104335562554305
 Electrons     Excitation
 0.391    :    3  ->    8
 0.574    :    3  ->    7
 CI weight: 0.616  Orb weight: 0.384
                                                                                                                          
* State 2
- Exc. energy: 0.40904333562089923
 Electrons     Excitation
 0.275    :    3  ->    9
 0.664    :    3  ->    6
 CI weight: 0.700  Orb weight: 0.300
                                                                                                                          
* State 3
- Exc. energy: 0.43116220101387687
 Electrons     Excitation
 0.279    :    4  ->    8
 0.207    :    2  ->    7
 0.158    :    2  ->    8
 0.332    :    4  ->    7
 CI weight: 0.335  Orb weight: 0.665
                                                                                                                          
* State 4
- Exc. energy: 0.5113196797605605
 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.5783859960264695
 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_results = mcscf_drv.compute(molecule, basis, space)
Hide code cell output
                                                                                                                          
                          Multi-Configurational Self-Consistent Field Driver
                         ====================================================
                                                                                                                          
        ╭────────────────────────────────────╮
        │          Driver settings           │
        ╰────────────────────────────────────╯
          State-specific calculation
          Max. iterations         : 50
          BFGS window             : 5
          Convergence thresholds:
            - Energy change       : 1e-08
            - Gradient norm       : 0.0001
                                                                                                                          
          Integrals in memory
                                                                                                                          
          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.037271082     0.0e+00      2.6e-03          1   0:00:00
        2     -76.037273353    -2.3e-06      6.4e-04          1   0:00:00
        3     -76.037273800    -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.6e-09      6.2e-05          1   0:00:00
** Convergence reached in 6 iterations
        7     -76.037274032    -9.4e-10      3.1e-05          1   0:00:00
                                                                                                                          
        Final results
        -------------
                                                                                                                          
                                                                                                                          
* State 1
- Energy: -76.03727403201913
- 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.83330 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.75294 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.73232 a.u.                                                                  
               (   1 O   2s  :     0.49) (   1 O   3s  :     0.18) (   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.67760 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.78416 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)                                                        
                                                                                                                          
                                                                                                                          
                                                Dipole moment for state 1                                                 
                                               ---------------------------                                                
                                                                                                                          
                                   X   :        -0.000000 a.u.        -0.000000 Debye                                     
                                   Y   :         0.000000 a.u.         0.000000 Debye                                     
                                   Z   :        -0.983159 a.u.        -2.498941 Debye                                     
                                 Total :         0.983159 a.u.         2.498941 Debye                                     
                                                                                                                          
                                                                                                                          
Total MCSCF time: 00:00:00
mcrpa = mtp.Mclr_EigenSolver()
rsp_dict = mcrpa.compute(molecule, basis, mcscf_drv, n_states=4)
Hide code cell output
                                                                                                                          
                          Multi-Configurational Linear response Eigenvalue Solver
                         =========================================================
                                                                                                                          
        ╭────────────────────────────────────╮
        │          Driver settings           │
        ╰────────────────────────────────────╯
          Number of states        : 4
                                                                                                                          
          Max. iterations         : 50
          Residual norm threshold : 1e-08
                                                                                                                          
          Integrals in memory
          Trial vector settings:
            - Linear dependence threshold :1e-10
            - Min. trial vectors          :16
            - Max. trial vectors          :48
                                                                                                                          
          CI parameters:      20
          Orbital parameters: 54
                                                                                                                          
                                                                                                                          
Initialization time: 00:00:00
                                                                                                                          
                                                                                                                          
        MC-RPA Iterations
        -----------------
                                                                                                                          
     Iter. | Average Energy | Grad. Norm | Converged | CI+Orb vec. | Time
     -----------------------------------------------------------------------
        1       0.710067027      2.4e-01      0/  4     4 +   4      0:00:00    
        2       0.504433348      7.0e-02      0/  4     4 +  12      0:00:00    
        3       0.460454267      4.4e-02      1/  4     4 +  20      0:00:00    
        4       0.438804358      6.4e-03      1/  4     6 +  24      0:00:00    
        5       0.436938061      2.4e-03      2/  4     8 +  27      0:00:00    
        6       0.423456314      1.1e-02      2/  4     8 +  31      0:00:00    
        7       0.421352332      3.8e-03      2/  4    10 +  33      0:00:00    
        8       0.420272535      4.0e-04      2/  4    12 +  35      0:00:00    
        9       0.420217560      3.4e-05      2/  4    12 +  39      0:00:00    
       10       0.420204856      2.3e-06      2/  4    16 +  39      0:00:00    
       11       0.420204685      8.3e-08      2/  4    16 +  43      0:00:00    
       12       0.420204669      6.8e-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.3305288      8.99415          0.0389580        0.0121777
          2     0.4084609     11.11479          0.0000000        0.0000000
          3     0.4307468     11.72122          0.1738542        0.1284536
          4     0.5110821     13.90725          0.1450310        0.1453630
                                                                                                                          
                                                                                                                          
                                                 Rotatory strength (cgs)
 Transition   Energy (Ha)    Energy (eV)    Velocity gauge     Length gauge 
 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
          1     0.3305288      8.99415         -0.0000000       -0.0000000
          2     0.4084609     11.11479         -0.0000000       -0.0000000
          3     0.4307468     11.72122          0.0000000        0.0000000
          4     0.5110821     13.90725          0.0000000       -0.0000000
                                                                                                                          
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.00814853 11.1306361  11.73252108 13.91371718]
MC-RPA Energies CAS(4,4)
[ 8.99414655 11.11478813 11.7212186  13.90725253]

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.33052879228146137
 Electrons     Excitation
 0.579    :    3  ->    7
 0.393    :    3  ->    8
 CI weight: -0.000  Orb weight: 1.000
                                                                                                                          
* State 2
- Exc. energy: 0.4084609334794797
 Electrons     Excitation
 0.670    :    3  ->    6
 0.275    :    3  ->    9
 CI weight: 0.000  Orb weight: 1.000
                                                                                                                          
* State 3
- Exc. energy: 0.43074684235635774
 Electrons     Excitation
 0.281    :    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.5110821086026888
 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, n_states=4)
Hide code cell output
                                                                                                                          
                          Multi-Configurational Linear response Eigenvalue Solver
                         =========================================================
                                                                                                                          
        ╭────────────────────────────────────╮
        │          Driver settings           │
        ╰────────────────────────────────────╯
          Number of states        : 4
          Tamm-Damcoff approximation
                                                                                                                          
          Max. iterations         : 50
          Residual norm threshold : 1e-08
                                                                                                                          
          Integrals in memory
          Trial vector settings:
            - Linear dependence threshold :1e-10
            - Min. trial vectors          :8
            - Max. trial vectors          :24
                                                                                                                          
          CI parameters:      20
          Orbital parameters: 54
                                                                                                                          
                                                                                                                          
Initialization time: 00:00:00
                                                                                                                          
                                                                                                                          
        MC-TDA Iterations
        -----------------
                                                                                                                          
     Iter. | Average Energy | Grad. Norm | Converged | CI+Orb vec. | Time
     -----------------------------------------------------------------------
        1       0.710192193      9.4e-01      0/  4     4 +   4      0:00:00    
        2       0.505688090      2.7e-01      0/  4     4 +   8      0:00:00    
        3       0.467941118      1.6e-01      0/  4     4 +  12      0:00:00    
        4       0.447283455      1.3e-02      1/  4     5 +  15      0:00:00    
        5       0.445599132      6.4e-03      1/  4     6 +  17      0:00:00    
        6       0.434442055      3.1e-02      2/  4     6 +  20      0:00:00    
        7       0.432814242      1.1e-02      2/  4     7 +  21      0:00:00    
        8       0.431128477      9.6e-04      2/  4     8 +  22      0:00:00    
        9       0.430925260      1.3e-04      2/  4     9 +  23      0:00:00    
       10       0.430911781      2.4e-05      2/  4    10 +  24      0:00:00    
       11       0.430909612      3.1e-06      2/  4    10 +  10      0:00:00    Orb subspace collapsed 
       12       0.430909517      8.1e-07      2/  4    11 +  11      0:00:00    
       13       0.430909448      1.4e-07      3/  4    12 +  12      0:00:00    
       14       0.430909446      6.0e-09      4/  4    12 +  13      0:00:00    
                                                                                                                          
** Convergence reached in 14 iterations
                                                                                                                          
                                                                                                                          
                                                                                                                          
                                                 Oscillator strength (au)
 Transition   Energy (Ha)    Energy (eV)    Velocity gauge     Length gauge 
  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
          1     0.3459460      9.41367          0.0403367        0.0142157
          2     0.4209258     11.45397          0.0000000        0.0000000
          3     0.4376907     11.91017          0.1841986        0.1222064
          4     0.5190753     14.12476          0.1187446        0.1440273
                                                                                                                          
                                                                                                                          
                                                 Rotatory strength (cgs)
 Transition   Energy (Ha)    Energy (eV)    Velocity gauge     Length gauge 
 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
          1     0.3459460      9.41367         -0.0000000       -0.0000000
          2     0.4209258     11.45397         -0.0000000       -0.0000000
          3     0.4376907     11.91017          0.0000000        0.0000000
          4     0.5190753     14.12476         -0.0000000        0.0000000
                                                                                                                          
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.99414655 11.11478813 11.7212186  13.90725253]
MC-TDA Energies CAS(4,4)
[ 9.41366984 11.45397333 11.9101715  14.12475838]