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)
Show code cell output
Multi-Configurational Self-Consistent Field Driver
====================================================
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: 100
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
Number of states : 1
State-specific calculation
- State of interest : 1
Max. iterations : 50
BFGS window : 5
Convergence thresholds:
- Energy change : 1e-08
- Gradient sq. norm : 1e-08
MCSCF Iterations
-------------------
Iter. | Average Energy | E. Change | Grad. Norm | CI Iter. | Trust rad. | Time
---------------------------------------------------------------------------------
1 -75.991312026 0.0e+00 5.7e-04 0 0.40 0:00:00
2 -76.006166338 -1.5e-02 1.9e-03 0 0.40 0:00:00
3 -76.032227066 -2.6e-02 6.1e-04 0 0.56 0:00:00
4 -76.031421785 8.1e-04 1.3e-03 0 0.78 0:00:00
5 -76.037415898 -6.0e-03 1.2e-04 0 0.39 0:00:00
6 -76.037698943 -2.8e-04 2.2e-05 0 0.47 0:00:00
7 -76.037780505 -8.2e-05 7.7e-06 0 0.47 0:00:00
8 -76.037840886 -6.0e-05 3.2e-06 0 0.66 0:00:00
9 -76.037842682 -1.8e-06 7.6e-07 0 0.66 0:00:00
10 -76.037843309 -6.3e-07 4.4e-08 0 0.79 0:00:00
11 -76.037843332 -2.3e-08 1.2e-08 0 0.79 0:00:00
12 -76.037843341 -8.6e-09 4.5e-10 0 0.80 0:00:00
** Convergence reached in 12 iterations
13 -76.037843341 -3.5e-10 7.9e-11 0 0.80 0:00:00
Final results
-------------
* State 1
- S^2 : 0.00 (multiplicity = 1.0 )
- Energy : -76.03784334095118
- 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.07302 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.50216 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.75291 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.73045 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.68035 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.78476 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)
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)
Show code cell output
Multi-Configurational Linear response Eigenvalue Solver
=========================================================
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
Number of states : 5
Full MC-RPA calculation
Max. iterations : 50
Residual sq. norm threshold : 1e-08
Trial vector settings:
- Lin. dep. threshold : 1e-10
- Max. trial vectors : 20
- Max. trial vectors : 60
Number of parameters:
- CI : 99
- Orbital : 52
# 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)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[6], line 9
6 au2ev = 27.211386
8 # Multiconfigurational linear response for CAS(6,5)
----> 9 MCRPA_6_5_excitation_energies = au2ev * rsp_dict["eigenvalues"][:4]
11 print("Full CI Energies")
12 print(FCI_excitation_energies)
TypeError: 'NoneType' object is not subscriptable
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.3310413425168189
Electrons Excitation
0.391 : 3 -> 8
0.574 : 3 -> 7
CI weight: 0.616 Orb weight: 0.384
* State 2
- Exc. energy: 0.40903703907882255
Electrons Excitation
0.275 : 3 -> 9
0.664 : 3 -> 6
CI weight: 0.700 Orb weight: 0.300
* State 3
- Exc. energy: 0.43116672307762466
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.5113203161722654
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.5783868256049328
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)
Show 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)
Total MCSCF time: 00:00:00
mcrpa = mtp.Mclr_EigenSolver()
rsp_dict = mcrpa.compute(molecule, basis, mcscf_drv, n_states=4)
Show code cell output
Initialization time: 00:00:00
Linear eigensolver
------------------
CI parameters: 20
Orbital parameters: 54
MC-RPA Iterations
-----------------
Iter. | Average Energy | Grad. Norm | Converged | CI+Orb vec. | Time
-----------------------------------------------------------------------
1 0.710081794 2.4e-01 0/ 4 4 + 4 0:00:00
2 0.504438122 7.0e-02 0/ 4 4 + 12 0:00:00
3 0.460463894 4.4e-02 1/ 4 4 + 20 0:00:00
4 0.438810895 6.4e-03 1/ 4 6 + 24 0:00:00
5 0.436944673 2.4e-03 2/ 4 8 + 27 0:00:00
6 0.423462081 1.1e-02 2/ 4 8 + 31 0:00:00
7 0.421356042 3.8e-03 2/ 4 10 + 33 0:00:00
8 0.420278497 4.0e-04 2/ 4 12 + 35 0:00:00
9 0.420223667 3.4e-05 2/ 4 12 + 39 0:00:00
10 0.420210977 2.3e-06 2/ 4 16 + 39 0:00:00
11 0.420210806 8.3e-08 2/ 4 16 + 43 0:00:00
12 0.420210790 6.9e-09 4/ 4 18 + 45 0:00:00
** Convergence reached in 12 iterations
Oscillator strength (au)
Transition Energy (Ha) Energy (eV) Velocity gauge Length gauge
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
1 0.3305460 8.99462 0.0389589 0.0121782
2 0.4084733 11.11512 0.0000000 0.0000000
3 0.4307439 11.72114 0.1738522 0.1284474
4 0.5110800 13.90719 0.1450456 0.1453738
Rotatory strength (cgs)
Transition Energy (Ha) Energy (eV) Velocity gauge Length gauge
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
1 0.3305460 8.99462 0.0000000 0.0000000
2 0.4084733 11.11512 0.0000000 0.0000000
3 0.4307439 11.72114 -0.0000000 -0.0000000
4 0.5110800 13.90719 -0.0000000 0.0000000
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.3305460361595478
Electrons Excitation
0.578 : 3 -> 7
0.393 : 3 -> 8
CI weight: 0.000 Orb weight: 1.000
* State 2
- Exc. energy: 0.40847326146200613
Electrons Excitation
0.670 : 3 -> 6
0.275 : 3 -> 9
CI weight: 0.000 Orb weight: 1.000
* State 3
- Exc. energy: 0.4307438863711089
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.5110799773765147
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)
Show code cell output
Initialization time: 00:00:00
Linear eigensolver
------------------
CI parameters: 20
Orbital parameters: 54
MC-RPA Iterations
-----------------
Iter. | Average Energy | Grad. Norm | Converged | CI+Orb vec. | Time
-----------------------------------------------------------------------
1 0.710081794 2.4e-01 0/ 4 4 + 4 0:00:00
2 0.504438122 7.0e-02 0/ 4 4 + 12 0:00:00
3 0.460463894 4.4e-02 1/ 4 4 + 20 0:00:00
4 0.438810895 6.4e-03 1/ 4 6 + 24 0:00:00
5 0.436944673 2.4e-03 2/ 4 8 + 27 0:00:00
6 0.423462081 1.1e-02 2/ 4 8 + 31 0:00:00
7 0.421356042 3.8e-03 2/ 4 10 + 33 0:00:00
8 0.420278497 4.0e-04 2/ 4 12 + 35 0:00:00
9 0.420223667 3.4e-05 2/ 4 12 + 39 0:00:00
10 0.420210977 2.3e-06 2/ 4 16 + 39 0:00:00
11 0.420210806 8.3e-08 2/ 4 16 + 43 0:00:00
12 0.420210790 6.9e-09 4/ 4 18 + 45 0:00:00
** Convergence reached in 12 iterations
Oscillator strength (au)
Transition Energy (Ha) Energy (eV) Velocity gauge Length gauge
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
1 0.3305460 8.99462 0.0389589 0.0121782
2 0.4084733 11.11512 0.0000000 0.0000000
3 0.4307439 11.72114 0.1738522 0.1284474
4 0.5110800 13.90719 0.1450456 0.1453738
Rotatory strength (cgs)
Transition Energy (Ha) Energy (eV) Velocity gauge Length gauge
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
1 0.3305460 8.99462 0.0000000 0.0000000
2 0.4084733 11.11512 0.0000000 0.0000000
3 0.4307439 11.72114 -0.0000000 -0.0000000
4 0.5110800 13.90719 -0.0000000 0.0000000
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]