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)
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: 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.0378433408259
- 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)
Show 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.33104335562553233
Electrons Excitation
0.391 : 3 -> 8
0.574 : 3 -> 7
CI weight: 0.616 Orb weight: 0.384
* State 2
- Exc. energy: 0.40904333562086964
Electrons Excitation
0.275 : 3 -> 9
0.664 : 3 -> 6
CI weight: 0.700 Orb weight: 0.300
* State 3
- Exc. energy: 0.4311622010138668
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.5113196797605557
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.5783859960264506
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.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)
Show 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.3305287922814588
Electrons Excitation
0.579 : 3 -> 7
0.393 : 3 -> 8
CI weight: 0.000 Orb weight: 1.000
* State 2
- Exc. energy: 0.40846093347948503
Electrons Excitation
0.670 : 3 -> 6
0.275 : 3 -> 9
CI weight: 0.000 Orb weight: 1.000
* State 3
- Exc. energy: 0.430746842356351
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.5110821086026898
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
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]