Frozen-Core FCI Calculations using CIPSI
HOME
Table of Contents
In this part of the tutorial, you will learn how to perform excited state calculations using the Configuration Interaction using a Perturbative Selection made Iteratively (CIPSI) algorithm in Quantum Package (QP). As an illustration, we will compute the energies of the two lowest excited states of the formaldehyde molecule in the aug-cc-pVDZ basis set.
The following table displays, in eV, the theoretical best estimates (TBE) of the excitation energies of formaldehyde. \(V\) denotes valence states, and \(R\) denotes Rydberg states. The data are taken from https://arxiv.org/abs/1807.02045 .
| Transition | TBE (eV) |
|---|---|
| \(1\:^{1}\mathrm{A_1} \rightarrow 1\:^{1}\mathrm{A_2} (\mathrm{V};n \rightarrow \pi^\star)\) | 3.97 |
| \(1\:^{1}\mathrm{A_1} \rightarrow 1\:^{1}\mathrm{B_2} (\mathrm{R};n \rightarrow 3s)\) | 7.3 |
| \(1\:^{1}\mathrm{A_1} \rightarrow 2\:^{1}\mathrm{B_2} (\mathrm{R};n \rightarrow 3p)\) | 8.14 |
| \(1\:^{1}\mathrm{A_1} \rightarrow 2\:^{1}\mathrm{A_1} (\mathrm{R};n \rightarrow 3p)\) | 8.27 |
| \(1\:^{1}\mathrm{A_1} \rightarrow 2\:^{1}\mathrm{A_2} (\mathrm{R};n \rightarrow 3p)\) | 8.5 |
| \(1\:^{1}\mathrm{A_1} \rightarrow 1\:^{1}\mathrm{B_1} (\mathrm{V})\) | 9.21 |
| \(1\:^{1}\mathrm{A_1} \rightarrow 3\:^{1}\mathrm{A_1} (\mathrm{V};\pi \rightarrow \pi^\star)\) | 9.26 |
| \(1\:^{1}\mathrm{A_1} \rightarrow 1\:^{3}\mathrm{A_2} (\mathrm{V};n \rightarrow \pi^\star)\) | 3.58 |
| \(1\:^{1}\mathrm{A_1} \rightarrow 1\:^{3}\mathrm{A_1} (\mathrm{V};\pi \rightarrow \pi^\star)\) | 6.07 |
| \(1\:^{1}\mathrm{A_1} \rightarrow 1\:^{3}\mathrm{B_2} (\mathrm{R};n \rightarrow 3s)\) | 7.14 |
| \(1\:^{1}\mathrm{A_1} \rightarrow 2\:^{3}\mathrm{B_2} (\mathrm{R};n \rightarrow 3p)\) | 7.96 |
| \(1\:^{1}\mathrm{A_1} \rightarrow 2\:^{3}\mathrm{A_1} (\mathrm{R};n \rightarrow 3p)\) | 8.15 |
| \(1\:^{1}\mathrm{A_1} \rightarrow 1\:^{3}\mathrm{B_1} (\mathrm{R})\) | 8.42 |
| \(1\:^{1}\mathrm{A^\prime} \rightarrow 1\:^{1}\mathrm{A^{\prime\prime}} [\mathrm{F}] (\mathrm{V};n \rightarrow \pi^\star)\) | 2.8 |
1 Preparing the EZFIO
First, let's create an xyz file for formaldehyde, in Angstroms:
4 Formaldehyde C 0.00000000 0.00000000 -0.60298508 O 0.00000000 0.00000000 0.60539399 H 0.00000000 0.93467313 -1.18217476 H 0.00000000 -0.93467313 -1.18217476
Next, create an EZFIO database containing the geometry and the basis set parameters for the aug-cc-pVDZ basis set, and run the Hartree-Fock calculation.
qp create_ezfio -b aug-cc-pvdz formaldehyde.xyz qp run scf | tee formaldehyde.scf.out qp set_frozen_core
The expected Hartree-Fock energy is -113.8850441 au.
Excited states generally have diffuse character. To accurately describe excited states, it is necessary to use a basis set that includes diffuse atomic basis functions. Diffuse functions are orbitals that are spatially extended beyond the region where the valence electrons are found.
When excited states are computed using a basis set that lacks diffuse functions, the computed wavefunctions may have insufficiently diffuse character and be of lower quality, resulting in overestimated excited-state energies and inaccurate properties.
2 Multi-state CIPSI calculation with Hartree-Fock orbitals
2.1 Blind test
To enable the calculation of multiple states, you can use qp
set determinants n_states. In this section, we will compute the
ground state and the two lowest states of formaldehyde, namely three
states:
qp set determinants n_states 3 qp set determinants n_det_max 50000 qp run fci | tee formaldehyde.fci1.out
In multi-state calculations, QP selects a common set of determinants for all the states. The consistency between the states is ensured by selecting determinants such that the PT2 correction is comparable within the different states.
Summary at N_det = 51242
-----------------------------------
# ============ ============================= ============================= =============================
State 1 State 2 State 3
# ============ ============================= ============================= =============================
# E -114.15561119 -113.83545337 -113.76434931
# Excit. (au) 0.00000000 0.32015782 0.39126188
# Excit. (eV) 0.00000000 8.71194149 10.64678213
# PT2 -0.08061958 0.00014194 -0.08570764 0.00015779 -0.09765686 0.00018134
# rPT2 -0.07907595 0.00013922 -0.08402890 0.00015470 -0.09537572 0.00017710
#
# E+PT2 -114.23623077 0.00014194 -113.92116101 0.00015779 -113.86200617 0.00018134
# E+rPT2 -114.23468714 0.00013922 -113.91948227 0.00015470 -113.85972503 0.00017710
# Excit. (au) 0.00000000 0.00020073 0.31506977 0.00021223 0.37422460 0.00023028
# Excit. (eV) 0.00000000 0.00546209 8.57348838 0.00577519 10.18317401 0.00626629
# ============ ============================= ============================= =============================
In this calculation, the PT2-corrected excitation energies are around 8.5 and 10.1 eV. This result is far from the expected values of 3.97 and 7.3 eV.
Let us look at the wave functions with qp edit to understand what
has gone wrong. Go to the end of the file, and search backwards for
= to get to the Determinants section:
Determinants ::
-0.956190388929 0.00808349705018 0.0442038650847
++++++++--------------------------------------------------------
++++++++--------------------------------------------------------
-0.00466646692438 -0.656599820416 -0.0219363805532
++++++++--------------------------------------------------------
+++++++--+------------------------------------------------------
-0.00466647325716 -0.656599812797 -0.021934926707
+++++++--+------------------------------------------------------
++++++++--------------------------------------------------------
-0.0160794523102 0.0200969800694 -0.614572492065
++++++++--------------------------------------------------------
++++++-+--+-----------------------------------------------------
-0.0160794624442 0.020096909905 -0.614572284411
++++++-+--+-----------------------------------------------------
++++++++--------------------------------------------------------
0.0216298014759 -0.00369205179465 0.221702915328
++++++++--------------------------------------------------------
++++++-+----+---------------------------------------------------
The determinants are shown by rotating 90 degrees clockwise the traditional spin-orbital diagrams. The first determinant is the Hartree-Fock determinant, and for each determinant, three numbers are given which correspond to their CI coefficients in states 1, 2, and 3.
The first state has a dominant Hartree-Fock character, as expected. The second state is a 8-10 singly excited state. Looking at the MO coefficients earlier in the file, we can identify this excitation as \(\pi \rightarrow \pi^*\), which might be the \(3 ^1 A_1\) state. The third state is a 7-11 single excitation, which is also a \(\pi \rightarrow \pi^*\).
Here, QP has not been able to obtain any of the two lowest states, due to the symmetry of the system. Formaldehyde has \(C_{2v}\) symmetry, and the lowest excited states, are in the irreducible representations \(A_2\) and \(B_2\), whereas the Hartree-Fock determinant and the ground state are in \(A_1\). As the \(C_{2v}\) symmetry group is Abelian, \(\langle I | H | J \rangle = 0\) if \(|I\rangle\) and \(|J\rangle\) are in different irreducible representations. Hence, only determinants in \(A_1\) can be selected starting from the Hartree-Fock determinant, and it is impossible to describe any of the \(A_2\), \(B_1\) and \(B_2\) states.
2.2 A correct run
A simple solution is to start from a set of determinants containing
at least one determinant of each irreducible representation.
The cis program performs the CI in the determinant space
containing the Hartree-Fock determinant and all possible singly
excited determinants.
qp run cis | tee formaldehyde.cis1.out
The CIS program obtains a first excited state at 4.55 eV, and a second one at 8.57 eV:
Excitation energies (au) (eV)
2 0.16732247703939152 4.5529677944444256
3 0.31507863549308013 8.5735216541156269
Looking at the wave function with qp edit, we can confirm we have
obtained different states:
Determinants ::
0.999999999917 7.74912835069e-11 -9.90449932734e-10
++++++++--------------------------------------------------------
++++++++--------------------------------------------------------
5.38429711827e-10 -1.9819140786e-10 0.618374351559
++++++++--------------------------------------------------------
+++++++-+-------------------------------------------------------
5.38460906745e-10 -1.98192930704e-10 0.618374351559
+++++++-+-------------------------------------------------------
++++++++--------------------------------------------------------
-3.62264652628e-11 0.494296392547 7.50508298599e-10
++++++++--------------------------------------------------------
+++++++---+-----------------------------------------------------
-3.62142767305e-11 0.494296392547 7.52683353491e-10
+++++++---+-----------------------------------------------------
++++++++--------------------------------------------------------
The second state is a 8-11 excitation, and the third state is a 8-13 excitation. Both are \(n \rightarrow \pi^*\) excitations. This looks correct!
A last point to address is that the orbitals are biased in favor of the ground state. We can make the orbitals more neutral by taking the state-average natural orbitals of the CIS calculation, and running the CIS again:
qp run save_natorb
The occupation number is now close to 4/3 in orbital 8, and close to 1/3 in orbitals 9 and 10:
======== ================ ================
MO Eigenvalue Cumulative
======== ================ ================
1 2.0000000000 2.0000000000
2 2.0000000000 4.0000000000
3 1.9999818795 5.9999818795
4 1.9999753489 7.9999572284
5 1.9999558595 9.9999130878
6 1.9992870769 11.9992001647
7 1.9928747860 13.9920749507
8 1.3412583826 15.3333333333
9 0.3319424230 15.6652757563
10 0.3305709648 15.9958467211
11 0.0020803062 15.9979270272
12 0.0012731375 15.9992001647
13 0.0006298587 15.9998300235
14 0.0000968613 15.9999268848
15 0.0000526128 15.9999794976
16 0.0000069726 15.9999864703
17 0.0000065101 15.9999929804
18 0.0000035978 15.9999965782
19 0.0000028594 15.9999994375
20 0.0000005625 16.0000000000
Running the CIS again gives more compact wave functions:
qp run cis | tee formaldehyde.cis2.out
Determinants ::
0.999999999983 7.38937386846e-11 -1.01495203432e-10
++++++++--------------------------------------------------------
++++++++--------------------------------------------------------
4.38247190443e-11 -0.70362182149 4.4184513023e-10
++++++++--------------------------------------------------------
+++++++-+-------------------------------------------------------
4.37452813195e-11 -0.70362182149 4.42073794179e-10
+++++++-+-------------------------------------------------------
++++++++--------------------------------------------------------
-1.27668931035e-11 -4.50577175529e-10 -0.702139343308
++++++++--------------------------------------------------------
+++++++--+------------------------------------------------------
-1.27533052582e-11 -4.50583098529e-10 -0.702139343307
+++++++--+------------------------------------------------------
++++++++--------------------------------------------------------
We are now ready to run the CIPSI calculation, using the CIS as a starting point:
qp set determinants read_wf true qp run fci | tee formaldehyde.fci2.out
The obtained energies are
Summary at N_det = 61167
-----------------------------------
# ============ ============================= ============================= =============================
State 1 State 2 State 3
# ============ ============================= ============================= =============================
# E -114.15552016 -114.00871579 -113.89649664
# Excit. (au) 0.00000000 0.14680438 0.25902352
# Excit. (eV) 0.00000000 3.99475213 7.04839183
# PT2 -0.08757543 0.00017014 -0.08655526 0.00016992 -0.08227528 0.00012616
# rPT2 -0.08574253 0.00016658 -0.08478030 0.00016643 -0.08072771 0.00012379
#
# E+PT2 -114.24309559 0.00017014 -114.09527104 0.00016992 -113.97877192 0.00012616
# E+rPT2 -114.24126270 0.00016658 -114.09349609 0.00016643 -113.97722435 0.00012379
# Excit. (au) 0.00000000 0.00024062 0.14782455 0.00024046 0.26432367 0.00021181
# Excit. (eV) 0.00000000 0.00654755 4.02251249 0.00654319 7.19261633 0.00576373
# ============ ============================= ============================= =============================
The excitation energies seem reasonable, but the run is too short to get accurate enough extrapolated energies:
Extrapolated energies
------------------------
State 1
=========== ===================
minimum PT2 Extrapolated energy
=========== ===================
-0.1279 -114.24302968
-0.1706 -114.24100537
-0.2101 -114.23929233
-0.2473 -114.23696140
-0.2967 -114.23323620
-0.3874 -114.22613671
=========== ===================
State 2
=========== =================== =================== ===================
minimum PT2 Extrapolated energy Excitation (a.u) Excitation (eV)
=========== =================== =================== ===================
-0.1256 -114.09427941 0.14875027 4.04770262
-0.1664 -114.09306651 0.14793886 4.02562301
-0.2077 -114.09216726 0.14712508 4.00347883
-0.2502 -114.09091014 0.14605127 3.97425895
-0.3019 -114.08733065 0.14590555 3.97029379
-0.3892 -114.08668158 0.13945513 3.79476873
=========== =================== =================== ===================
State 3
=========== =================== =================== ===================
minimum PT2 Extrapolated energy Excitation (a.u) Excitation (eV)
=========== =================== =================== ===================
-0.1146 -113.98182916 0.26120052 7.10763096
-0.1507 -113.97980692 0.26119845 7.10757457
-0.1930 -113.97838501 0.26090733 7.09965283
-0.2362 -113.97701405 0.25994736 7.07353069
-0.2733 -113.97578459 0.25745160 7.00561764
-0.4078 -113.97487597 0.25126074 6.83715554
=========== =================== =================== ===================
2.3 An improved run
To improve the results, we can take the state-average natural orbitals of the previous CIPSI calculation, and run another CIPSI calculation after a CIS. Also, it is worth noting that 50000 determinants is a very small calculation, and a few million determinants should be considered.
qp run save_natorb qp run cis | tee formaldehyde.cis3.out
The CIS energies are already greatly improved:
Excitation energies (au) (eV)
2 0.150674637886709 4.09996783382701
3 0.278687823609005 7.58330086935594
qp run fci | tee formaldehyde.fci3.out
We now obtain:
Summary at N_det = 102890
-----------------------------------
# ============ ============================= ============================= =============================
State 1 State 2 State 3
# ============ ============================= ============================= =============================
# E -114.19125322 -114.04674536 -113.93267218
# Excit. (au) 0.00000000 0.14450786 0.25858104
# Excit. (eV) 0.00000000 3.93226072 7.03635129
# PT2 -0.05717516 0.00011029 -0.05359972 0.00010551 -0.05141479 0.00010125
# rPT2 -0.05644379 0.00010888 -0.05296138 0.00010425 -0.05083531 0.00010011
#
# E+PT2 -114.24842838 0.00011029 -114.10034508 0.00010551 -113.98408697 0.00010125
# E+rPT2 -114.24769701 0.00010888 -114.09970674 0.00010425 -113.98350749 0.00010011
# Excit. (au) 0.00000000 0.00015597 0.14808331 0.00015263 0.26434141 0.00014972
# Excit. (eV) 0.00000000 0.00424423 4.02955356 0.00415329 7.19309903 0.00407396
# ============ ============================= ============================= =============================
and the extrapolated energies are:
State 1 =========== =================== minimum PT2 Extrapolated energy =========== =================== -0.0813 -114.24479277 -0.1173 -114.24404595 -0.1711 -114.24426568 -0.2445 -114.24194770 -0.3985 -114.23699775 =========== =================== State 2 =========== =================== =================== =================== minimum PT2 Extrapolated energy Excitation (a.u) Excitation (eV) =========== =================== =================== =================== -0.0755 -114.09662130 0.14817146 4.03195238 -0.1113 -114.09703285 0.14701310 4.00043179 -0.1659 -114.09685928 0.14740640 4.01113395 -0.2404 -114.09557286 0.14637485 3.98306404 -0.3913 -114.09220186 0.14479589 3.94009845 =========== =================== =================== =================== State 3 =========== =================== =================== =================== minimum PT2 Extrapolated energy Excitation (a.u) Excitation (eV) =========== =================== =================== =================== -0.0727 -113.98287958 0.26191319 7.12702357 -0.1036 -113.98180071 0.26224524 7.13605928 -0.1534 -113.98137173 0.26289394 7.15371140 -0.2233 -113.98010171 0.26184599 7.12519512 -0.4013 -113.97583603 0.26116172 7.10657503 =========== =================== =================== ===================
Leaving the calculation run much longer by considering a larger n_det_max,
one should obtain 3.99 eV for the first state, and 7.11 eV for the second state.