Skip to content

Latest commit

 

History

History
471 lines (376 loc) · 21.1 KB

03-Formaldehyde.org

File metadata and controls

471 lines (376 loc) · 21.1 KB

Frozen-Core FCI Calculations using CIPSI

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.

./images/formaldehyde.png

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} → 1\:1\mathrm{A_2} (\mathrm{V};n → π^*)$ 3.97
$1\:1\mathrm{A_1} → 1\:1\mathrm{B_2} (\mathrm{R};n → 3s)$ 7.3
$1\:1\mathrm{A_1} → 2\:1\mathrm{B_2} (\mathrm{R};n → 3p)$ 8.14
$1\:1\mathrm{A_1} → 2\:1\mathrm{A_1} (\mathrm{R};n → 3p)$ 8.27
$1\:1\mathrm{A_1} → 2\:1\mathrm{A_2} (\mathrm{R};n → 3p)$ 8.5
$1\:1\mathrm{A_1} → 1\:1\mathrm{B_1} (\mathrm{V})$ 9.21
$1\:1\mathrm{A_1} → 3\:1\mathrm{A_1} (\mathrm{V};π → π^*)$ 9.26
$1\:1\mathrm{A_1} → 1\:3\mathrm{A_2} (\mathrm{V};n → π^*)$ 3.58
$1\:1\mathrm{A_1} → 1\:3\mathrm{A_1} (\mathrm{V};π → π^*)$ 6.07
$1\:1\mathrm{A_1} → 1\:3\mathrm{B_2} (\mathrm{R};n → 3s)$ 7.14
$1\:1\mathrm{A_1} → 2\:3\mathrm{B_2} (\mathrm{R};n → 3p)$ 7.96
$1\:1\mathrm{A_1} → 2\:3\mathrm{A_1} (\mathrm{R};n → 3p)$ 8.15
$1\:1\mathrm{A_1} → 1\:3\mathrm{B_1} (\mathrm{R})$ 8.42
$1\:1\mathrm{A^′} → 1\:1\mathrm{A′\prime} [\mathrm{F}] (\mathrm{V};n → π^*)$ 2.8

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.

Multi-state CIPSI calculation with Hartree-Fock orbitals

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 $π → π^*$, which might be the $3 ^1 A_1$ state. The third state is a 7-11 single excitation, which is also a $π → π^*$.

Here, QP has not been able to obtain any of the two lowest states, due to the symmetry of the system. Formaldehyde has $C2v$ 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 $C2v$ 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.

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 → π^*$ 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
 =========== =================== =================== ===================

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.