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 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()

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

# 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
                                                                                                                          
        Final results
        -------------
                                                                                                                          
                                                                                                                          
* State 1
- Energy: -76.03784333970783
- 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)                              
                                                                                                                          
                                                                                                                          
                                                Dipole moment for state 1                                                 
                                               ---------------------------                                                
                                                                                                                          
                                   X   :         0.000000 a.u.         0.000000 Debye                                     
                                   Y   :         0.000000 a.u.         0.000000 Debye                                     
                                   Z   :        -0.981805 a.u.        -2.495500 Debye                                     
                                 Total :         0.981805 a.u.         2.495500 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
                                                                                                                          
Warning, significant negative value in orbital S(2): -0.3720628025672088
          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.663245798      2.9e-01      0/  5     5 +   5      0:00:00    
        2       0.530193941      1.1e-01      0/  5     5 +  15      0:00:00    
        3       0.477684538      2.8e-02      0/  5    11 +  19      0:00:00    
        4       0.463500859      4.1e-02      0/  5    15 +  24      0:00:00    
        5       0.455006499      1.3e-02      0/  5    19 +  30      0:00:00    
        6       0.452743653      9.0e-04      0/  5    27 +  32      0:00:00    
        7       0.452274946      4.2e-04      2/  5    29 +  40      0:00:00    
        8       0.452221617      2.2e-04      3/  5    35 +  40      0:00:00    
        9       0.452193120      3.0e-05      4/  5    37 +  42      0:00:00    
       10       0.452191485      6.4e-06      4/  5    37 +  44      0:00:00    
       11       0.452190487      3.9e-07      4/  5    39 +  44      0:00:00    
       12       0.452190455      4.2e-08      4/  5    39 +  46      0:00:00    
       13       0.452190450      1.1e-08      4/  5    41 +  46      0:00:00    
       14       0.452190449      1.1e-09      5/  5    41 +  48      0:00:00    
                                                                                                                          
** Convergence reached in 14 iterations
                                                                                                                          
                                                                                                                          
                                                                                                                          
                                                 Oscillator strength (au)
 Transition   Energy (Ha)    Energy (eV)    Velocity gauge     Length gauge 
  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
          1     0.3310413      9.00809          0.0389375        0.0122857
          2     0.4090370     11.13046          0.0000000        0.0000000
          3     0.4311667     11.73264          0.1736421        0.1284502
          4     0.5113203     13.91373          0.1444200        0.1449619
          5     0.5783868     15.73871          0.3471772        0.4072682
                                                                                                                          
                                                                                                                          
                                                 Rotatory strength (cgs)
 Transition   Energy (Ha)    Energy (eV)    Velocity gauge     Length gauge 
 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
          1     0.3310413      9.00809         -0.0000000       -0.0000000
          2     0.4090370     11.13046          0.0000000        0.0000000
          3     0.4311667     11.73264         -0.0000000       -0.0000000
          4     0.5113203     13.91373         -0.0000000        0.0000000
          5     0.5783868     15.73871         -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.00809375 11.13046476 11.73264413 13.91373449]

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.3310413425168086
 Electrons     Excitation
 0.391    :    3  ->    8
 0.574    :    3  ->    7
 CI weight: 0.616  Orb weight: 0.384
                                                                                                                          
* State 2
- Exc. energy: 0.40903703907882366
 Electrons     Excitation
 0.275    :    3  ->    9
 0.664    :    3  ->    6
 CI weight: 0.700  Orb weight: 0.300
                                                                                                                          
* State 3
- Exc. energy: 0.431166723077617
 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.5113203161722619
 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.5783868256049282
 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.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.03727403103234
- 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)                                                        
                                                                                                                          
                                                                                                                          
                                                Dipole moment for state 1                                                 
                                               ---------------------------                                                
                                                                                                                          
                                   X   :        -0.000000 a.u.        -0.000000 Debye                                     
                                   Y   :        -0.000000 a.u.        -0.000000 Debye                                     
                                   Z   :        -0.983192 a.u.        -2.499025 Debye                                     
                                 Total :         0.983192 a.u.         2.499025 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.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
                                                                                                                          
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.00809375 11.13046476 11.73264413 13.91373449]
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.3305460361595499
 Electrons     Excitation
 0.578    :    3  ->    7
 0.393    :    3  ->    8
 CI weight: 0.000  Orb weight: 1.000
                                                                                                                          
* State 2
- Exc. energy: 0.4084732614620102
 Electrons     Excitation
 0.670    :    3  ->    6
 0.275    :    3  ->    9
 CI weight: 0.000  Orb weight: 1.000
                                                                                                                          
* State 3
- Exc. energy: 0.4307438863711155
 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.5110799773765211
 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
                                                                                                                          
          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.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
                                                                                                                          
Total MC-RPA 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)
[ 8.99461578 11.11512359 11.72113816 13.90719454]