State-averaged MCSCF#
CASCI response#
As for truncated CI, the exact state response formulas also apply to CASCI and other active space CIs. The method can thus create all excited states that arise from excitations within the active space. This can also be used as an approximation for a MCSCF response, neglecting the orbital part of the response.
Let us try it on our usual water case. We’ll use an active space including the 2 pairs of O-H \(\sigma\) and \(\sigma^*\) (4 electrons in 4 orbitals). However, since the main excitation also involves the oxygen lone pair, we add this orbital too, for a total CAS(6,5).
First, we try the intuitive concept of first optimizing the CASSCF wavefunction for the ground state, and then applying CASCI response with the 5 states.
import veloxchem as vlx
import multipsi as mtp
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)
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)
nstates=5
# CASCI excitation energies
cidrv=mtp.CIDriver()
ci_results = cidrv.compute(molecule,basis,space,nstates)
# Transition properties
SI=mtp.InterState()
DipOsc=SI.compute(molecule, basis, ci_results)
# Results computed in CI response
FCI_excitation_energies=[8.44992, 10.68451, 10.97654, 13.37692]
CIS_excitation_energies=[9.38059, 11.30739, 11.83172, 13.90237]
au2ev = 27.211386
print("Full CI Energies")
print(FCI_excitation_energies)
print("CIS Energies")
print(CIS_excitation_energies)
print("CASCI Energies")
print(au2ev*DipOsc['energies'])
Full CI Energies
[8.44992, 10.68451, 10.97654, 13.37692]
CIS Energies
[9.38059, 11.30739, 11.83172, 13.90237]
CASCI Energies
[18.61974291 18.73087038 26.61579239 28.51654574]
As we can see, the result is very disappointing, very far from the full CI result, and even from the CIS result. The reason for that is the orbitals. The orbitals were optimized for the ground state, and not for the excited states. Like for truncated CI, we thus have an imbalance between the description of the ground and excited states leading to (significantly) too high excitation energies. In particular, in this case, the \(\sigma^*\) orbitals are optimal to describe correlation in the O-H bonds, but not for the electronic excitation.
State-specific optimization#
A way to improve over the previous situation is to optimize every state separately with their own orbitals. This is called state-specific optimization. To ensure that we are doing an excited state, we still need to compute the ground states (or in general all states below the targeted state), but tell the program to optimize the orbitals specifically for the n_th state. In MultiPsi, this is done by providing “weights” (this concept will become more clear in the next section), a weight of 1 for the state of interest and 0 for all the other, meaning they will not influence the orbital optimization.
mcscf_results = mcscf_drv.compute(molecule,basis,space)
e0 = mcscf_drv.get_energy()
mcscf_results = mcscf_drv.compute(molecule,basis,space.copy(), n_states = 2 , weights = [0,1])
e1 = mcscf_drv.get_energy()
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 -76.037843340 0.0e+00 5.7e-05 1 0:00:00
** Convergence reached in 1 iterations
Final results
-------------
* State 1
- Energy: -76.03784333970778
- 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
Multi-Configurational Self-Consistent Field Driver
====================================================
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
State-specific calculation
Number of states : 2
State of interest : 2
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.353580402 0.0e+00 4.4e-01 1 0:00:00
2 -75.593582599 -2.4e-01 2.8e-01 1 0:00:00
3 -75.632467002 -3.9e-02 1.3e-01 1 0:00:00
4 -75.642345183 -9.9e-03 1.7e-02 1 0:00:00
5 -75.642439224 -9.4e-05 5.8e-03 1 0:00:00
6 -75.642461614 -2.2e-05 2.4e-03 1 0:00:00
7 -75.642463965 -2.4e-06 1.3e-03 1 0:00:00
8 -75.642464740 -7.7e-07 4.8e-04 1 0:00:00
9 -75.642464800 -6.0e-08 2.2e-04 1 0:00:00
10 -75.642464810 -1.1e-08 5.7e-05 1 0:00:00
11 -75.642464811 -9.0e-10 2.1e-05 1 0:00:00
** Convergence reached in 11 iterations
Final results
-------------
* State 1
- Energy: -75.99867639594106
- S^2 : -0.00 (multiplicity = 1.0 )
- Natural orbitals
1.98064 1.98062 0.01721 1.99899 0.02255
* State 2
- Energy: -75.64246481134326
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.99396 1.97991 1.00091 1.00000 0.02521
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 1:
--------------------------
Occupation: 2.000 Energy: -20.71709 a.u.
( 1 O 1s : 1.00)
Molecular Orbital No. 2:
--------------------------
Occupation: 2.000 Energy: -1.18832 a.u.
( 1 O 1s : -0.21) ( 1 O 2s : 0.50) ( 1 O 3s : 0.54)
( 1 O 1p0 : 0.27) ( 1 O 2p0 : 0.19)
Molecular Orbital No. 3:
--------------------------
Occupation: 1.994 Energy: -0.84981 a.u.
( 1 O 1p+1: -0.59) ( 1 O 2p+1: -0.31) ( 2 H 1s : 0.22)
( 3 H 1s : -0.22)
Molecular Orbital No. 4:
--------------------------
Occupation: 1.980 Energy: -0.98577 a.u.
( 1 O 2s : 0.17) ( 1 O 1p0 : -0.54) ( 1 O 2p0 : -0.34)
( 2 H 1s : 0.18) ( 3 H 1s : 0.18)
Molecular Orbital No. 5:
--------------------------
Occupation: 1.001 Energy: 0.12396 a.u.
( 1 O 1p+1: -0.37) ( 1 O 2p+1: -0.71) ( 2 H 1s : -0.16)
( 2 H 2s : -1.25) ( 3 H 1s : 0.16) ( 3 H 2s : 1.25)
Molecular Orbital No. 6:
--------------------------
Occupation: 1.000 Energy: -0.35135 a.u.
( 1 O 1p-1: 0.70) ( 1 O 2p-1: 0.44)
Molecular Orbital No. 7:
--------------------------
Occupation: 0.025 Energy: 0.52545 a.u.
( 1 O 2s : -0.42) ( 1 O 3s : -0.41) ( 1 O 1p0 : 0.67)
( 2 H 1s : 0.47) ( 2 H 2s : 0.33) ( 3 H 1s : 0.47)
( 3 H 2s : 0.33)
Molecular Orbital No. 8:
--------------------------
Occupation: 0.000 Energy: 0.73664 a.u.
( 1 O 3s : -1.03) ( 1 O 1p0 : -0.49) ( 1 O 2p0 : 0.75)
( 2 H 1s : -0.45) ( 2 H 2s : 1.06) ( 3 H 1s : -0.45)
( 3 H 2s : 1.06)
Molecular Orbital No. 9:
--------------------------
Occupation: 0.000 Energy: 0.99710 a.u.
( 1 O 2p+1: 0.65) ( 2 H 1s : 0.97) ( 2 H 2s : -0.62)
( 3 H 1s : -0.97) ( 3 H 2s : 0.62)
Dipole moment for state 1
---------------------------
X : 0.000000 a.u. 0.000000 Debye
Y : 0.000000 a.u. 0.000000 Debye
Z : -0.878307 a.u. -2.232433 Debye
Total : 0.878307 a.u. 2.232433 Debye
Dipole moment for state 2
---------------------------
X : -0.000000 a.u. -0.000000 Debye
Y : -0.000000 a.u. -0.000000 Debye
Z : -0.001750 a.u. -0.004447 Debye
Total : 0.001750 a.u. 0.004447 Debye
Total MCSCF time: 00:00:00
print(f"First state-specific MCSCF excitation energy = {(e1-e0)*au2ev:.4f}")
print(f"Corresponding (2nd) FCI excitation energy = {FCI_excitation_energies[1]:.4f}")
First state-specific MCSCF excitation energy = 10.7588
Corresponding (2nd) FCI excitation energy = 10.6845
As hoped, the first excitation energy is now a much better match to the FCI result. Actually, if we analyzed the orbitals, we would see that the state we converged is actually the second excited state, and the result is thus very close to the FCI result. We can also note that the second calculation still prints the ground state energy, but on the orbitals optimized for the excited state, and because of this, the energy of this state is higher than in the first calculation.
While in principle the most accurate description of the state that can be obtained at the MCSCF level, state-specific optimization has severe shortcomings hindering broader applicability. First, they way the calculation is done, the excited state is not necessarily orthogonal to the ground state, and it may therefore be artificially lower than it should. Second, if there is a small energy gap between states, the optimization may have trouble converging. We can see it if we try to optimize a second excited state:
mcscf_results = mcscf_drv.compute(molecule,basis,space.copy(), n_states = 3 , weights = [0,0,1])
Multi-Configurational Self-Consistent Field Driver
====================================================
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
State-specific calculation
Number of states : 3
State of interest : 3
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.349496543 0.0e+00 5.0e-01 1 0:00:00
2 -75.357366238 -7.9e-03 4.5e-01 1 0:00:00
3 -75.498311381 -1.4e-01 3.6e-01 1 0:00:00
4 -75.612877454 -1.1e-01 1.4e-01 1 0:00:00
5 -75.623966474 -1.1e-02 4.3e-02 1 0:00:00
6 -75.624835914 -8.7e-04 2.6e-02 1 0:00:00
7 -75.612608663 1.2e-02 4.3e-01 1 0:00:00
8 -75.621434275 -8.8e-03 3.7e-01 1 0:00:00
9 -75.624429758 -3.0e-03 1.7e-01 1 0:00:00
10 -75.629065196 -4.6e-03 1.1e-01 1 0:00:00
11 -75.630933085 -1.9e-03 4.1e-02 1 0:00:00
12 -75.631870791 -9.4e-04 3.2e-02 1 0:00:00
13 -75.635001940 -3.1e-03 1.0e-01 1 0:00:00
14 -75.474764698 1.6e-01 5.0e-01 1 0:00:00
15 -75.559698693 -8.5e-02 4.2e-01 1 0:00:00
16 -75.636132197 -7.6e-02 1.4e-01 1 0:00:00
17 -75.632633824 3.5e-03 3.3e-01 1 0:00:00
18 -75.628219538 4.4e-03 2.9e-01 1 0:00:00
19 -75.634712298 -6.5e-03 1.3e-01 1 0:00:00
20 -75.633640246 1.1e-03 3.0e-01 1 0:00:00
21 -75.632310724 1.3e-03 2.6e-01 1 0:00:00
22 -75.635207566 -2.9e-03 2.2e-01 1 0:00:00
23 -75.636894185 -1.7e-03 1.1e-01 1 0:00:00
24 -75.637407852 -5.1e-04 2.8e-01 1 0:00:00
25 -75.637461038 -5.3e-05 1.1e-01 1 0:00:00
26 -75.637599340 -1.4e-04 2.3e-01 1 0:00:00
27 -75.637647869 -4.9e-05 8.8e-02 1 0:00:00
28 -75.638008599 -3.6e-04 5.0e-02 1 0:00:00
29 -75.633848849 4.2e-03 3.1e-01 1 0:00:00
30 -75.635054284 -1.2e-03 3.1e-01 1 0:00:00
31 -75.637214078 -2.2e-03 3.1e-01 1 0:00:00
32 -75.637444925 -2.3e-04 1.1e-01 1 0:00:00
33 -75.637664463 -2.2e-04 7.1e-02 1 0:00:00
34 -75.637864064 -2.0e-04 3.7e-02 1 0:00:00
35 -75.637908727 -4.5e-05 2.1e-02 1 0:00:00
36 -75.637953607 -4.5e-05 3.4e-02 1 0:00:00
37 -75.636609098 1.3e-03 2.8e-01 1 0:00:00
38 -75.637232557 -6.2e-04 2.8e-01 1 0:00:00
39 -75.638062077 -8.3e-04 6.1e-02 1 0:00:00
40 -75.637591655 4.7e-04 1.7e-01 1 0:00:00
41 -75.637843060 -2.5e-04 1.5e-01 1 0:00:00
42 -75.637634075 2.1e-04 3.0e-01 1 0:00:00
43 -75.637841526 -2.1e-04 2.7e-01 1 0:00:00
44 -75.637931356 -9.0e-05 1.6e-01 1 0:00:00
45 -75.637956249 -2.5e-05 1.4e-01 1 0:00:00
46 -75.638010069 -5.4e-05 5.1e-02 1 0:00:00
47 -75.638030298 -2.0e-05 4.4e-02 1 0:00:00
48 -75.638040439 -1.0e-05 6.2e-02 1 0:00:00
49 -75.637775518 2.6e-04 2.9e-01 1 0:00:00
50 -75.637814890 -3.9e-05 2.8e-01 1 0:00:00
** Convergence not reached after 50 iterations
Final results
-------------
* State 1
- Energy: -75.9919107570913
- S^2 : -0.00 (multiplicity = 1.0 )
- Natural orbitals
1.98030 1.97986 0.02309 1.99882 0.01794
* State 2
- Energy: -75.63815587448187
- S^2 : -0.00 (multiplicity = 1.0 )
- Natural orbitals
1.98477 1.99709 0.01787 1.00000 1.00026
* State 3
- Energy: -75.63781488990989
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.99766 1.98523 1.00136 1.00000 0.01575
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 1:
--------------------------
Occupation: 2.000 Energy: -20.65202 a.u.
( 1 O 1s : -1.00)
Molecular Orbital No. 2:
--------------------------
Occupation: 2.000 Energy: -1.19753 a.u.
( 1 O 1s : 0.22) ( 1 O 2s : -0.51) ( 1 O 3s : -0.54)
( 1 O 1p0 : -0.25) ( 1 O 2p0 : -0.18)
Molecular Orbital No. 3:
--------------------------
Occupation: 1.998 Energy: -0.88016 a.u.
( 1 O 2s : -0.16) ( 1 O 1p+1: -0.30) ( 1 O 1p0 : 0.52)
( 1 O 2p+1: -0.17) ( 1 O 2p0 : 0.32) ( 3 H 1s : -0.24)
Molecular Orbital No. 4:
--------------------------
Occupation: 1.985 Energy: -0.83932 a.u.
( 1 O 1p+1: 0.49) ( 1 O 1p0 : 0.24) ( 1 O 2p+1: 0.22)
( 1 O 2p0 : 0.16) ( 2 H 1s : -0.31) ( 2 H 2s : -0.16)
Molecular Orbital No. 5:
--------------------------
Occupation: 1.001 Energy: 0.15875 a.u.
( 1 O 2s : 0.21) ( 1 O 3s : 0.73) ( 1 O 1p+1: 0.17)
( 1 O 1p0 : -0.36) ( 1 O 2p+1: 0.37) ( 1 O 2p0 : -0.19)
( 2 H 1s : -0.18) ( 3 H 1s : -0.38) ( 3 H 2s : -1.26)
Molecular Orbital No. 6:
--------------------------
Occupation: 1.000 Energy: -0.31746 a.u.
( 1 O 1p-1: 0.70) ( 1 O 2p-1: 0.44)
Molecular Orbital No. 7:
--------------------------
Occupation: 0.016 Energy: 0.24804 a.u.
( 1 O 3s : 0.46) ( 1 O 1p+1: -0.36) ( 1 O 1p0 : -0.26)
( 1 O 2p+1: -0.63) ( 2 H 1s : -0.28) ( 2 H 2s : -1.43)
( 3 H 2s : 0.69)
Molecular Orbital No. 8:
--------------------------
Occupation: 0.000 Energy: 0.96888 a.u.
( 1 O 3s : 0.44) ( 1 O 1p0 : 0.54) ( 1 O 2p+1: 0.43)
( 1 O 2p0 : -0.53) ( 2 H 1s : 1.08) ( 2 H 2s : -1.02)
( 3 H 2s : -0.27)
Molecular Orbital No. 9:
--------------------------
Occupation: 0.000 Energy: 1.04914 a.u.
( 1 O 3s : -0.28) ( 1 O 1p0 : -0.59) ( 1 O 2p+1: 0.51)
( 1 O 2p0 : 0.78) ( 2 H 1s : 0.62) ( 3 H 1s : -0.80)
( 3 H 2s : 0.82)
Dipole moment for state 1
---------------------------
X : -0.000218 a.u. -0.000553 Debye
Y : -0.000000 a.u. -0.000001 Debye
Z : -0.823779 a.u. -2.093838 Debye
Total : 0.823779 a.u. 2.093838 Debye
Dipole moment for state 2
---------------------------
X : 1.254525 a.u. 3.188685 Debye
Y : 0.000001 a.u. 0.000002 Debye
Z : 0.022982 a.u. 0.058414 Debye
Total : 1.254736 a.u. 3.189220 Debye
Dipole moment for state 3
---------------------------
X : -1.250024 a.u. -3.177245 Debye
Y : 0.000003 a.u. 0.000006 Debye
Z : 0.026023 a.u. 0.066143 Debye
Total : 1.250295 a.u. 3.177933 Debye
Total MCSCF time: 00:00:00
The optimization oscillates as optimizing the second excited state brings its energy very close to the first, and thus the order of the states keep switching. While there exists “root following” techniques to overcome this, state-specific optimizations are usually fairly cumbersome.
For all these reasons, state-specific optimization is usually replaced by a so-called state-averaged optimization.
State-averaging#
Ideally, we would want each state to have its own orbitals while maintaining orthogonality betweem states. In the next chapter, we will see how this can be achieved within the framework of full MCSCF response. However, there is a way to improve the result without including the orbital response. The concept is to use only one orbital set for all states, but to make these orbitals a “compromise” between the needs of the ground state and that of the excited states. The description of each state would be worse but the hope is that this error would be consistent across all states and mostly cancel out in the excitation energies.
The way to do this is to construct a “state-average” MCSCF (SA-MCSCF), that is optimize the orbitals to minimize the average energy of all states. Doing this requires only a trivial modification of the code, namely the CI part of the MCSCF includes all states (here 5) and then the MCSCF energy is replaced by the average energy of all states and the CI densities are replaced by the average density of all states. Note that this is simply an extension of the state-specific optimization mentioned above, but using equal weights for all states.
Let’s see it in practice:
# CASSCF calculation
# 5 states included already in the MCSCF.
mcscf_results = mcscf_drv.compute(molecule,basis,space,nstates)
# Transition properties
DipOsc=SI.compute(molecule, basis, mcscf_results)
print("Full CI Energies")
print(FCI_excitation_energies)
print("CIS Energies")
print(CIS_excitation_energies)
print("SA-CASSCF Energies")
print(au2ev*DipOsc['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]
The result is now significantly closer to the converged full CI result. In this particular case, the result is not better than CIS, but of course, one the main motivation to use MCSCF is for systems with strong correlation, which was not the case here.
This method is of course not without issues. In particular, unlike standard responses calculation, the introduction of state-averaging means that the excitation energies depend on the total number of states, since the state-averaged orbitals depend on all states. In particular, we can expect the results to deteriorate with increasing number of states as the orbitals begin to compromise many states of widely different natures:
first_exc_energy=[]
for nstates in range(2,6):
# CASSCF calculation
# 5 states included already in the MCSCF
mcscf_results = mcscf_drv.compute(molecule,basis,space,nstates)
# Transition properties
DipOsc=SI.compute(molecule, basis, mcscf_results)
first_exc_energy.append(au2ev*DipOsc['energies'][0])
Show code cell output
Multi-Configurational Self-Consistent Field Driver
====================================================
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
State-averaged calculation
Number of states : 2
Equal-weights
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.843270715 0.0e+00 1.1e-01 1 0:00:00
2 -75.858829274 -1.6e-02 6.9e-02 1 0:00:00
3 -75.864455907 -5.6e-03 2.8e-02 1 0:00:00
4 -75.864858794 -4.0e-04 9.1e-03 1 0:00:00
5 -75.864894174 -3.5e-05 3.0e-03 1 0:00:00
6 -75.864897105 -2.9e-06 7.1e-04 1 0:00:00
7 -75.864897238 -1.3e-07 1.2e-04 1 0:00:00
8 -75.864897243 -4.9e-09 2.3e-05 1 0:00:00
** Convergence reached in 8 iterations
Final results
-------------
* State 1
- Energy: -76.01351196497133
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.98571 1.97737 1.99925 0.01513 0.02254
* State 2
- Energy: -75.71628252119801
- S^2 : -0.00 (multiplicity = 1.0 )
- Natural orbitals
1.99623 1.97627 1.00000 1.00146 0.02604
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 1:
--------------------------
Occupation: 2.000 Energy: -20.64857 a.u.
( 1 O 1s : 1.00)
Molecular Orbital No. 2:
--------------------------
Occupation: 2.000 Energy: -1.19270 a.u.
( 1 O 1s : -0.22) ( 1 O 2s : 0.51) ( 1 O 3s : 0.56)
( 1 O 1p0 : 0.22) ( 1 O 2p0 : 0.16)
Molecular Orbital No. 3:
--------------------------
Occupation: 1.985 Energy: -0.86833 a.u.
( 1 O 1p0 : 0.55) ( 1 O 2p0 : 0.35) ( 2 H 1s : -0.19)
( 3 H 1s : -0.19)
Molecular Orbital No. 4:
--------------------------
Occupation: 1.973 Energy: -0.76522 a.u.
( 1 O 1p+1: 0.54) ( 1 O 2p+1: 0.27) ( 2 H 1s : -0.25)
( 3 H 1s : 0.25)
Molecular Orbital No. 5:
--------------------------
Occupation: 1.500 Energy: -0.43591 a.u.
( 1 O 1p-1: 0.66) ( 1 O 2p-1: 0.49)
Molecular Orbital No. 6:
--------------------------
Occupation: 0.515 Energy: 0.12362 a.u.
( 1 O 2s : -0.18) ( 1 O 3s : -1.04) ( 1 O 1p0 : 0.33)
( 1 O 2p0 : 0.38) ( 2 H 2s : 0.89) ( 3 H 2s : 0.89)
Molecular Orbital No. 7:
--------------------------
Occupation: 0.028 Energy: 0.60893 a.u.
( 1 O 1p+1: 0.82) ( 2 H 1s : 0.43) ( 2 H 2s : 0.44)
( 3 H 1s : -0.43) ( 3 H 2s : -0.44)
Molecular Orbital No. 8:
--------------------------
Occupation: 0.000 Energy: 0.84499 a.u.
( 1 O 1p+1: 0.29) ( 1 O 2p+1: -0.56) ( 2 H 1s : 0.70)
( 2 H 2s : -1.51) ( 3 H 1s : -0.70) ( 3 H 2s : 1.51)
Molecular Orbital No. 9:
--------------------------
Occupation: 0.000 Energy: 1.10548 a.u.
( 1 O 1p-1: -0.95) ( 1 O 2p-1: 1.05)
Dipole moment for state 1
---------------------------
X : -0.000000 a.u. -0.000000 Debye
Y : -0.000000 a.u. -0.000000 Debye
Z : -0.914416 a.u. -2.324213 Debye
Total : 0.914416 a.u. 2.324213 Debye
Dipole moment for state 2
---------------------------
X : -0.000000 a.u. -0.000000 Debye
Y : 0.000000 a.u. 0.000000 Debye
Z : 0.066682 a.u. 0.169488 Debye
Total : 0.066682 a.u. 0.169488 Debye
Total MCSCF time: 00:00:00
List of oscillator strengths greather than 1e-10
From to Energy (eV) Oscillator strength (length and velocity)
1 2 8.08803 8.153140e-03 3.791438e-02
List of rotatory strengths greater than 1e-10
From to Energy (eV) Rot. strength (a.u. and 10^-40 cgs)
Multi-Configurational Self-Consistent Field Driver
====================================================
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
State-averaged calculation
Number of states : 3
Equal-weights
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.724068112 0.0e+00 4.6e-02 1 0:00:00
2 -75.725200061 -1.1e-03 1.2e-02 1 0:00:00
3 -75.725331152 -1.3e-04 2.9e-03 1 0:00:00
4 -75.725334435 -3.3e-06 8.8e-04 1 0:00:00
5 -75.725334998 -5.6e-07 3.2e-04 1 0:00:00
6 -75.725335073 -7.5e-08 1.4e-04 1 0:00:00
7 -75.725335090 -1.7e-08 2.4e-05 1 0:00:00
8 -75.725335090 -5.3e-10 9.8e-06 1 0:00:00
** Convergence reached in 8 iterations
Final results
-------------
* State 1
- Energy: -76.00986345021587
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.98623 1.99931 1.97762 0.01412 0.02272
* State 2
- Energy: -75.717605128665
- S^2 : -0.00 (multiplicity = 1.0 )
- Natural orbitals
1.99620 1.00000 1.97731 1.00100 0.02548
* State 3
- Energy: -75.4485366920499
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.99442 1.99771 1.00266 0.99771 0.00751
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 1:
--------------------------
Occupation: 2.000 Energy: -20.66135 a.u.
( 1 O 1s : 1.00)
Molecular Orbital No. 2:
--------------------------
Occupation: 2.000 Energy: -1.25118 a.u.
( 1 O 1s : -0.22) ( 1 O 2s : 0.52) ( 1 O 3s : 0.56)
( 1 O 1p0 : 0.18)
Molecular Orbital No. 3:
--------------------------
Occupation: 1.988 Energy: -0.84539 a.u.
( 1 O 1p0 : -0.56) ( 1 O 2p0 : -0.36) ( 2 H 1s : 0.19)
( 3 H 1s : 0.19)
Molecular Orbital No. 4:
--------------------------
Occupation: 1.666 Energy: -0.50141 a.u.
( 1 O 1p-1: 0.66) ( 1 O 2p-1: 0.48)
Molecular Orbital No. 5:
--------------------------
Occupation: 1.648 Energy: -0.70435 a.u.
( 1 O 1p+1: -0.55) ( 1 O 2p+1: -0.28) ( 2 H 1s : 0.25)
( 3 H 1s : -0.25)
Molecular Orbital No. 6:
--------------------------
Occupation: 0.675 Energy: 0.08480 a.u.
( 1 O 2s : 0.16) ( 1 O 3s : 1.07) ( 1 O 1p0 : -0.31)
( 1 O 2p0 : -0.39) ( 2 H 2s : -0.92) ( 3 H 2s : -0.92)
Molecular Orbital No. 7:
--------------------------
Occupation: 0.023 Energy: 0.55592 a.u.
( 1 O 1p+1: 0.77) ( 1 O 2p+1: 0.18) ( 2 H 1s : 0.44)
( 2 H 2s : 0.51) ( 3 H 1s : -0.44) ( 3 H 2s : -0.51)
Molecular Orbital No. 8:
--------------------------
Occupation: 0.000 Energy: 0.84549 a.u.
( 1 O 1p+1: 0.26) ( 1 O 2p+1: -0.40) ( 2 H 1s : 0.75)
( 2 H 2s : -1.44) ( 3 H 1s : -0.75) ( 3 H 2s : 1.44)
Molecular Orbital No. 9:
--------------------------
Occupation: 0.000 Energy: 1.08882 a.u.
( 1 O 1p-1: -0.95) ( 1 O 2p-1: 1.05)
Dipole moment for state 1
---------------------------
X : 0.000000 a.u. 0.000000 Debye
Y : 0.000000 a.u. 0.000000 Debye
Z : -0.902882 a.u. -2.294897 Debye
Total : 0.902882 a.u. 2.294897 Debye
Dipole moment for state 2
---------------------------
X : -0.000000 a.u. -0.000000 Debye
Y : -0.000000 a.u. -0.000000 Debye
Z : 0.071922 a.u. 0.182808 Debye
Total : 0.071922 a.u. 0.182808 Debye
Dipole moment for state 3
---------------------------
X : -0.000000 a.u. -0.000000 Debye
Y : 0.000000 a.u. 0.000000 Debye
Z : -0.220124 a.u. -0.559500 Debye
Total : 0.220124 a.u. 0.559500 Debye
Total MCSCF time: 00:00:00
List of oscillator strengths greather than 1e-10
From to Energy (eV) Oscillator strength (length and velocity)
1 2 7.95275 8.899154e-03 4.003592e-02
1 3 15.27448 3.502949e-01 3.169114e-01
List of rotatory strengths greater than 1e-10
From to Energy (eV) Rot. strength (a.u. and 10^-40 cgs)
Multi-Configurational Self-Consistent Field Driver
====================================================
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
State-averaged calculation
Number of states : 4
Equal-weights
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.650640385 0.0e+00 1.1e-01 1 0:00:00
2 -75.692220054 -4.2e-02 5.0e-02 1 0:00:00
3 -75.694350434 -2.1e-03 1.3e-02 1 0:00:00
4 -75.694439713 -8.9e-05 3.7e-03 1 0:00:00
5 -75.694445644 -5.9e-06 1.2e-03 1 0:00:00
6 -75.694446718 -1.1e-06 4.2e-04 1 0:00:00
7 -75.694446789 -7.2e-08 1.9e-04 1 0:00:00
8 -75.694446804 -1.5e-08 2.5e-05 1 0:00:00
9 -75.694446804 -2.8e-10 5.8e-06 1 0:00:00
** Convergence reached in 9 iterations
Final results
-------------
* State 1
- Energy: -75.99698106363174
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.98702 1.98210 1.99942 0.01415 0.01731
* State 2
- Energy: -75.70831737893658
- S^2 : -0.00 (multiplicity = 1.0 )
- Natural orbitals
1.99704 1.98460 1.00000 1.00024 0.01812
* State 3
- Energy: -75.62405350917159
- S^2 : -0.00 (multiplicity = 1.0 )
- Natural orbitals
1.98743 1.99379 1.00000 0.02056 0.99822
* State 4
- Energy: -75.44843526583908
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.99702 1.00111 1.99835 0.99825 0.00526
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 1:
--------------------------
Occupation: 2.000 Energy: -20.67802 a.u.
( 1 O 1s : -1.00)
Molecular Orbital No. 2:
--------------------------
Occupation: 2.000 Energy: -1.24525 a.u.
( 1 O 1s : 0.22) ( 1 O 2s : -0.52) ( 1 O 3s : -0.55)
( 1 O 1p0 : -0.20)
Molecular Orbital No. 3:
--------------------------
Occupation: 1.988 Energy: -0.87394 a.u.
( 1 O 1p0 : 0.56) ( 1 O 2p0 : 0.36) ( 2 H 1s : -0.18)
( 3 H 1s : -0.18)
Molecular Orbital No. 4:
--------------------------
Occupation: 1.734 Energy: -0.74018 a.u.
( 1 O 1p+1: 0.56) ( 1 O 2p+1: 0.28) ( 2 H 1s : -0.24)
( 3 H 1s : 0.24)
Molecular Orbital No. 5:
--------------------------
Occupation: 1.499 Energy: -0.46605 a.u.
( 1 O 1p-1: -0.67) ( 1 O 2p-1: -0.48)
Molecular Orbital No. 6:
--------------------------
Occupation: 0.512 Energy: 0.10271 a.u.
( 1 O 2s : 0.16) ( 1 O 3s : 1.07) ( 1 O 1p0 : -0.30)
( 1 O 2p0 : -0.39) ( 2 H 2s : -0.92) ( 3 H 2s : -0.92)
Molecular Orbital No. 7:
--------------------------
Occupation: 0.267 Energy: 0.23802 a.u.
( 1 O 1p+1: 0.45) ( 1 O 2p+1: 0.67) ( 2 H 1s : 0.18)
( 2 H 2s : 1.19) ( 3 H 1s : -0.18) ( 3 H 2s : -1.19)
Molecular Orbital No. 8:
--------------------------
Occupation: 0.000 Energy: 0.99722 a.u.
( 1 O 2p+1: 0.62) ( 2 H 1s : 0.97) ( 2 H 2s : -0.66)
( 3 H 1s : -0.97) ( 3 H 2s : 0.66)
Molecular Orbital No. 9:
--------------------------
Occupation: 0.000 Energy: 1.08042 a.u.
( 1 O 1p-1: -0.94) ( 1 O 2p-1: 1.05)
Dipole moment for state 1
---------------------------
X : 0.000000 a.u. 0.000000 Debye
Y : 0.000000 a.u. 0.000000 Debye
Z : -0.865078 a.u. -2.198810 Debye
Total : 0.865078 a.u. 2.198810 Debye
Dipole moment for state 2
---------------------------
X : -0.000000 a.u. -0.000000 Debye
Y : -0.000000 a.u. -0.000000 Debye
Z : 0.085986 a.u. 0.218555 Debye
Total : 0.085986 a.u. 0.218555 Debye
Dipole moment for state 3
---------------------------
X : 0.000000 a.u. 0.000000 Debye
Y : -0.000000 a.u. -0.000000 Debye
Z : -0.042740 a.u. -0.108634 Debye
Total : 0.042740 a.u. 0.108634 Debye
Dipole moment for state 4
---------------------------
X : -0.000000 a.u. -0.000000 Debye
Y : 0.000000 a.u. 0.000000 Debye
Z : -0.237138 a.u. -0.602745 Debye
Total : 0.237138 a.u. 0.602745 Debye
Total MCSCF time: 00:00:00
List of oscillator strengths greather than 1e-10
From to Energy (eV) Oscillator strength (length and velocity)
1 2 7.85494 7.812110e-03 4.149883e-02
1 4 14.92669 4.583994e-01 3.096527e-01
List of rotatory strengths greater than 1e-10
From to Energy (eV) Rot. strength (a.u. and 10^-40 cgs)
Multi-Configurational Self-Consistent Field Driver
====================================================
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
State-averaged calculation
Number of states : 5
Equal-weights
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.630250938 0.0e+00 7.0e-02 1 0:00:00
2 -75.672532201 -4.2e-02 9.5e-02 1 0:00:00
3 -75.683260759 -1.1e-02 4.7e-02 1 0:00:00
4 -75.686039130 -2.8e-03 5.0e-03 1 0:00:00
5 -75.686051087 -1.2e-05 1.1e-03 1 0:00:00
6 -75.686051424 -3.4e-07 2.2e-04 1 0:00:00
7 -75.686051445 -2.1e-08 2.9e-05 1 0:00:00
8 -75.686051445 -2.8e-10 8.7e-06 1 0:00:00
** Convergence reached in 8 iterations
Final results
-------------
* State 1
- Energy: -75.98214464848321
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.98461 1.99145 1.99958 0.01560 0.00875
* State 2
- Energy: -75.70439712311851
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.98783 1.99722 1.00000 0.01582 0.99913
* State 3
- Energy: -75.62069009030203
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.99499 1.99411 1.00000 0.99764 0.01326
* State 4
- Energy: -75.60257199547218
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.96066 1.22366 1.99968 0.04213 0.77386
* State 5
- Energy: -75.52045336778457
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.97487 1.01895 1.99956 0.97610 0.03052
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 1:
--------------------------
Occupation: 2.000 Energy: -20.68941 a.u.
( 1 O 1s : -1.00)
Molecular Orbital No. 2:
--------------------------
Occupation: 2.000 Energy: -1.44678 a.u.
( 1 O 1s : 0.22) ( 1 O 2s : -0.49) ( 1 O 3s : -0.48)
Molecular Orbital No. 3:
--------------------------
Occupation: 1.976 Energy: -0.80534 a.u.
( 1 O 1p+1: -0.55) ( 1 O 2p+1: -0.28) ( 2 H 1s : 0.23)
( 3 H 1s : -0.23)
Molecular Orbital No. 4:
--------------------------
Occupation: 1.604 Energy: -0.58138 a.u.
( 1 O 2s : -0.18) ( 1 O 3s : -0.24) ( 1 O 1p0 : -0.60)
( 1 O 2p0 : -0.40)
Molecular Orbital No. 5:
--------------------------
Occupation: 1.600 Energy: -0.50270 a.u.
( 1 O 1p-1: -0.67) ( 1 O 2p-1: -0.48)
Molecular Orbital No. 6:
--------------------------
Occupation: 0.414 Energy: 0.21172 a.u.
( 1 O 1p+1: -0.43) ( 1 O 2p+1: -0.69) ( 2 H 2s : -1.24)
( 3 H 2s : 1.24)
Molecular Orbital No. 7:
--------------------------
Occupation: 0.406 Energy: 0.11950 a.u.
( 1 O 2s : -0.17) ( 1 O 3s : -1.09) ( 1 O 1p0 : 0.25)
( 1 O 2p0 : 0.41) ( 2 H 2s : 0.93) ( 3 H 2s : 0.93)
Molecular Orbital No. 8:
--------------------------
Occupation: 0.000 Energy: 1.00726 a.u.
( 1 O 2p+1: 0.69) ( 2 H 1s : 0.97) ( 2 H 2s : -0.59)
( 3 H 1s : -0.97) ( 3 H 2s : 0.59)
Molecular Orbital No. 9:
--------------------------
Occupation: 0.000 Energy: 1.07263 a.u.
( 1 O 1p-1: 0.94) ( 1 O 2p-1: -1.05)
Dipole moment for state 1
---------------------------
X : -0.000000 a.u. -0.000000 Debye
Y : -0.000000 a.u. -0.000000 Debye
Z : -0.847572 a.u. -2.154313 Debye
Total : 0.847572 a.u. 2.154313 Debye
Dipole moment for state 2
---------------------------
X : 0.000000 a.u. 0.000000 Debye
Y : -0.000000 a.u. -0.000000 Debye
Z : 0.071706 a.u. 0.182259 Debye
Total : 0.071706 a.u. 0.182259 Debye
Dipole moment for state 3
---------------------------
X : 0.000000 a.u. 0.000000 Debye
Y : 0.000000 a.u. 0.000000 Debye
Z : -0.023143 a.u. -0.058824 Debye
Total : 0.023143 a.u. 0.058824 Debye
Dipole moment for state 4
---------------------------
X : -0.000000 a.u. -0.000000 Debye
Y : -0.000000 a.u. -0.000000 Debye
Z : 0.110859 a.u. 0.281774 Debye
Total : 0.110859 a.u. 0.281774 Debye
Dipole moment for state 5
---------------------------
X : 0.000000 a.u. 0.000000 Debye
Y : -0.000000 a.u. -0.000000 Debye
Z : 0.135445 a.u. 0.344267 Debye
Total : 0.135445 a.u. 0.344267 Debye
Total MCSCF time: 00:00:00
List of oscillator strengths greather than 1e-10
From to Energy (eV) Oscillator strength (length and velocity)
1 2 7.55790 1.182597e-02 2.629915e-02
1 4 10.32870 1.101138e-01 1.302110e-01
1 5 12.56326 1.491866e-01 1.439092e-01
List of rotatory strengths greater than 1e-10
From to Energy (eV) Rot. strength (a.u. and 10^-40 cgs)
import matplotlib.pyplot as plt
import numpy as np
plt.figure(figsize=(6,4))
x = np.array(range(2,6))
plt.plot(x,np.array(first_exc_energy), label='First SA-CASSCF excitation energy')
plt.title('Deteriotation of state-averaged CASSCF excitation energies')
plt.ylabel("Energies (eV)")
plt.xlabel("Number of states")
plt.axhline(y=FCI_excitation_energies[0],label ='FCI',color='red')
plt.xticks(x)
plt.legend()
plt.tight_layout(); plt.show()
While for 2 states, the first excitation energy was relatively close to the fullCI, it degrades significantly as more states are added, with a difference of more than 0.5 eV going from 2 to 5 states. Adding dynamical correlation on-top of the CASSCF typically reduces this error, but it is still an undesirable behaviour.
Still, because of its simplicity and its reasonable accuracy, SA-MCSCF is one of the most frequently used technique to compute properties at the MCSCF level.