Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Configuration interaction

Full-CI response

As we saw in the previous section, for an exact state, we have the formula

⟨ ⁣⟨Ω^;V^ω⟩ ⁣⟩=−1ℏ∑n>0[⟨0∣Ω^∣n⟩⟨n∣V^ω∣0⟩ωn0−ω+⟨0∣V^ω∣n⟩⟨n∣Ω^∣0⟩ωn0+ω]\langle \! \langle \hat{\Omega}; \hat{V}^{\omega} \rangle \! \rangle = - \frac{1}{\hbar} \sum_{n>0} \left[ \frac{\langle 0 | \hat{\Omega} | n \rangle \langle n | \hat{V}^{\omega} | 0 \rangle } {\omega_{n0}-\omega} + \frac{\langle 0 | \hat{V}^{\omega} | n \rangle \langle n | \hat{\Omega} | 0 \rangle } {\omega_{n0}+\omega} \right]
(1)

For example, in the case of a molecule irradiated by light, we are interested in computing the polarizability ααβ(−ω;ω)\alpha_{\alpha\beta}(-\omega; \omega), corresponding to the response with Ω^=μ^α\hat{\Omega} = \hat{\mu}_\alpha and V^ω=−μ^β\hat{V}^\omega = - \hat{\mu}_\beta, where α,β∈{x,y,z}\alpha, \beta \in \{x,y,z\}. Looking more precisely at the isotropic component, αˉ(−ω;ω)=13∑αααα(−ω;ω)\bar\alpha (-\omega; \omega) = \frac{1}{3} \sum_\alpha \alpha_{\alpha\alpha} (-\omega; \omega), we obtain

αˉ(−ω;ω)=23∑n>0∑αωn0∣⟨0∣μ^α∣n⟩∣2ωn02−ω2\bar\alpha (-\omega; \omega) = \frac{2}{3} \sum_{n>0} \sum_\alpha \frac{ \omega_{n0} | \langle 0 | \hat{\mu}_\alpha | n \rangle |^2 }{ \omega_{n0}^2 - \omega^2 }
(2)

We now define the oscillator strength fnf_n for a transition 0→n0 \rightarrow n as

fn=23∑αωn0∣⟨0∣μ^α∣n⟩∣2f_n = \frac{2}{3} \sum_\alpha \omega_{n0} | \langle 0 | \hat{\mu}_\alpha | n \rangle |^2
(3)

such that

αˉ(−ω;ω)=∑n>0fnωn02−ω2\bar\alpha (-\omega; \omega) = \sum_{n>0} \frac{ f_n }{ \omega_{n0}^2 - \omega^2 }
(4)

Those formulas can be applied directly for a configuration interaction wave function, where it is easy to get the excited state energies simply by obtaining additional roots in the diagonalization of the Hamiltonian. We then only need to compute the matrix element ⟨0∣μ^α∣n⟩\langle 0 | \hat{\mu}_\alpha | n \rangle. In second quantization, the α\alpha component of the dipole operator can be written as

μ^α=∑pqμα,pqa^p†a^q\hat{\mu}_{\alpha} = \sum_{pq} \mu_{\alpha, pq} \hat{a}_p^\dagger \hat{a}_q
(5)

We thus need to compute the matrix γ\boldsymbol{\gamma}, with elements γpq=⟨0∣a^p†a^q∣n⟩ \gamma_{pq} =\langle 0 |\hat{a}_p^\dagger \hat{a}_q | n \rangle, which is called the transition density matrix, and contract it with the dipole moment integrals.

Let’s apply these for full CI (FCI) of a water molecule.

import veloxchem as vlx
import multipsi as mtp
import numpy as np
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)
Loading...
# Compute the 5 lowest states of water using full CI

space=mtp.OrbSpace(molecule,scf_drv.mol_orbs)
# Full CI with frozen 1s orbital
space.fci(n_frozen=1)

nstates=5
cidrv=mtp.CIDriver()
ci_results = cidrv.compute(molecule,basis,space,nstates)
Output
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Cell In[4], line 9
      5 space.fci(n_frozen=1)
      6 
      7 nstates=5
      8 cidrv=mtp.CIDriver()
----> 9 ci_results = cidrv.compute(molecule,basis,space,nstates)

File ~/miniconda3/envs/echem/lib/python3.12/site-packages/multipsi/ci_driver.py:173, in CIDriver.compute(self, molecule, basis, orbital_space, n_states)
    171 else:
    172     mcham = MC_Hamiltonian()
--> 173     mcham.init_integrals(molecule, basis, old_space)
    174     energy, e_inactive, Ftu, tuvw = mcham.ci_gradients(None, None, 1.0e-12)
    175     ham_operator.twoel_operator(e_inactive, Ftu, tuvw)

File ~/miniconda3/envs/echem/lib/python3.12/site-packages/multipsi/mc_hamiltonian.py:74, in MC_Hamiltonian.init_integrals(self, molecule, basis, orbital_space)
     71     do_mcscf = self.functional.do_mcscf
     73 # Nuclear term
---> 74 self.nuc_energy = molecule.nuclear_repulsion_energy()
     76 # 1el integrals
     77 ## Do nuclear in parallel
     78 if molecule.number_of_atoms() >= self.MPI_size and self.MPI_size > 1:

File ~/miniconda3/envs/echem/lib/python3.12/site-packages/veloxchem/molecule.py:548, in _Molecule_nuclear_repulsion_energy(self, basis)
    538 def _Molecule_nuclear_repulsion_energy(self, basis=None):
    539     """
    540     Deprecated compatibility shim for effective nuclear repulsion energy.
    541 
   (...)    545         This function always raises.
    546     """
--> 548     assert_msg_critical(
    549         False,
    550         'Molecule.nuclear_repulsion_energy is deprecated; use Molecule.effective_nuclear_repulsion_energy(basis)'
    551     )

File ~/miniconda3/envs/echem/lib/python3.12/site-packages/veloxchem/errorhandler.py:48, in assert_msg_critical(condition, msg)
     38 """
     39 Asserts that the condition is true. Otherwise terminates the program with
     40 an error message.
   (...)     45     The error message.
     46 """
     47 if __debug__ and MPI.COMM_WORLD.Get_size() == 1:
---> 48     assert condition, msg
     49 else:
     50     if not condition:

AssertionError: Molecule.nuclear_repulsion_energy is deprecated; use Molecule.effective_nuclear_repulsion_energy(basis)
au2ev = 27.211386

# Get the molecular orbitals within the active space
C=space._get_active_mos()

# Compute dipole integrals in AO basis
dipole_drv = vlx.ElectricDipoleIntegralsDriver()
dipole_mats = dipole_drv.compute(molecule, basis)
dipole=[dipole_mats.x_to_numpy(),dipole_mats.y_to_numpy(),dipole_mats.z_to_numpy()]

# Transform to MO basis
dipole_mo=[]
for icomp in range(0,3):
    dipole_mo.append(np.einsum('pq,pt,qu->tu', dipole[icomp], C, C))
    

# Initialize the CIOperator class to compute the transition densities.
expansion=mtp.CIExpansion(space)
DenDriver=mtp.CIOperator(expansion, cidrv.comm)
energies = ci_results["energies"]
ci_vectors = ci_results["ci_vectors"]
        
# Compute all 0->n transitions
for n in range(1,nstates):
    dE=energies[n]-energies[0]
    tden=DenDriver.get_1dm(ci_vectors, 0, ci_vectors, n) #Transition density matrix
    dx=np.tensordot(tden, dipole_mo[0]) #<0|x|n>
    dy=np.tensordot(tden, dipole_mo[1]) #<0|y|n>
    dz=np.tensordot(tden, dipole_mo[2]) #<0|z|n>
    F=2/3*dE*(dx*dx+dy*dy+dz*dz)
    print(f"Excitation 0->{n} energy: {dE*au2ev:5.3f} eV oscillator strength: {F:.5f}")

We thus obtain 3 visible transitions and one dark state (near-zero dipole oscillator strength). We obtain the exact same result by using the built-in class of multipsi:

SI = mtp.StateInteraction()
si_results = SI.compute(molecule, basis, ci_results)

Truncated CI

In truncated CI, the energy does not just depend on the CI coefficients but also on the molecular orbitals, which makes the truncated CI response more complicated. However, it is interesting to look at the results we obtain by simply applying the above equations. This gives a hierarchy of method from CIS to full CI. However, the excitation energies do not necessarily improve from order to order. Let us demonstrate this on the water example.

nstates=5
# nstates-1 transitions for CIS, CISD, CISDT and CISDTQ and FCI
Energies=np.empty((5,nstates-1))
Energies[4,:]=au2ev*np.array(si_results['energies']) # Save the FCI result

#CIS to CISDTQ
space=mtp.OrbSpace(molecule,scf_drv.mol_orbs)
cidrv=mtp.CIDriver()
SI=mtp.StateInteraction()

for exc in range(1,5):
    space.ci(exc,n_frozen=1)
    ci_results = cidrv.compute(molecule,basis,space, nstates)
    si_results = SI.compute(molecule,basis,ci_results)
    Energies[exc-1,:]=au2ev*np.array(si_results['energies'])
Output
                                                                                                                          
List of oscillator strengths greather than 1e-10
                                                                                                                          
  From     to       Energy (eV)    Oscillator strength (length and velocity)
     1       2       10.29159         1.548866e-02    3.658376e-02
     1       4       12.78501         1.247451e-01    1.355482e-01
     1       5       15.12866         1.136810e-01    9.274229e-02
                                                                                                                          
List of rotatory strengths greater than 1e-10
                                                                                                                          
  From     to       Energy (eV)    Rot. strength (a.u. and 10^-40 cgs)
                          Configuration Interaction Driver
                         ==================================
                                                                                                                          

          Active space definition:
          ------------------------
Number of inactive (occupied) orbitals: 1
Number of active orbitals:              12
Number of virtual orbitals:             0

    This is a GASSCF wavefunction

         Cumulated   Min cumulated    Max cumulated 
 Space    orbitals      occupation       occupation
     1           4               5                8
     2          12               8                8


          CI expansion:
          -------------
Number of determinants:      6329


        ╭────────────────────────────────────╮
        │          Driver settings           │
        ╰────────────────────────────────────╯
          Max. iterations         : 40
          Initial diagonalization : 200
          Max subspace size       : 50
          Convergence thresholds:
            - Energy change       : 1e-08
            - Residual square norm: 1e-08
          Olsen step
                                                                                                                          
                                                                                                                          
        CI Iterations
        -------------
                                                                                                                          
     Iter. | Average Energy | E. Change | Grad. Norm |   Time
     ----------------------------------------------------------
        1     -75.708185001     0.0e+00      3.5e-01   0:00:04
        2     -75.787071882    -9.3e-02      1.4e-02   0:00:04
        3     -75.791912542    -5.3e-03      1.2e-03   0:00:04
        4     -75.792184181    -3.9e-04      8.9e-05   0:00:04
        5     -75.792206026    -4.2e-05      1.1e-05   0:00:04
        6     -75.792207797    -4.0e-06      8.3e-07   0:00:04
        7     -75.792207936    -3.5e-07      7.9e-08   0:00:04
        8     -75.792207948    -2.7e-08      7.6e-09   0:00:05
        9     -75.792207950    -7.4e-09      6.3e-09   0:00:04
                                                                                                                          
** Convergence reached in 9 iterations
                                                                                                                          
        Final results
        -------------
                                                                                                                          
                                                                                                                          
* State 1
- Energy: -76.11438963659869
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.98960 1.97168 1.97493 1.98327 0.00047 0.02478 0.00289 0.01562 0.02338 0.01071 0.00060 0.00208
                                                                                                                          
* State 2
- Energy: -75.80003060665913
- S^2   : -0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99159 1.97002 1.98187 0.99537 0.99592 0.03029 0.00318 0.00554 0.01330 0.00820 0.00318 0.00237
                                                                                                                          
* State 3
- Energy: -75.71893849552846
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99111 1.97998 1.97319 0.99537 0.03033 0.99427 0.00339 0.00552 0.03033 0.00953 0.01223 0.00311
                                                                                                                          
* State 4
- Energy: -75.70731804012323
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.98804 1.96044 1.12926 1.98430 0.86217 0.03726 0.00263 0.01467 0.01016 0.00542 0.00348 0.00217
                                                                                                                          
* State 5
- Energy: -75.62036297238913
- S^2   : -0.00  (multiplicity = 1.0 )
- Natural orbitals
1.98665 1.95751 1.01386 1.98540 0.04016 0.97604 0.00341 0.01331 0.00418 0.00547 0.01214 0.00418
                                                                                                                          
List of oscillator strengths greather than 1e-10
                                                                                                                          
  From     to       Energy (eV)    Oscillator strength (length and velocity)
     1       2        8.55414         1.328951e-02    4.580239e-02
     1       4       11.07698         1.154300e-01    1.705756e-01
     1       5       13.44315         1.119484e-01    1.180816e-01
                                                                                                                          
List of rotatory strengths greater than 1e-10
                                                                                                                          
  From     to       Energy (eV)    Rot. strength (a.u. and 10^-40 cgs)
                          Configuration Interaction Driver
                         ==================================
                                                                                                                          

          Active space definition:
          ------------------------
Number of inactive (occupied) orbitals: 1
Number of active orbitals:              12
Number of virtual orbitals:             0

    This is a GASSCF wavefunction

         Cumulated   Min cumulated    Max cumulated 
 Space    orbitals      occupation       occupation
     1           4               4                8
     2          12               8                8


          CI expansion:
          -------------
Number of determinants:      27763


        ╭────────────────────────────────────╮
        │          Driver settings           │
        ╰────────────────────────────────────╯
          Max. iterations         : 40
          Initial diagonalization : 200
          Max subspace size       : 50
          Convergence thresholds:
            - Energy change       : 1e-08
            - Residual square norm: 1e-08
          Olsen step
                                                                                                                          
                                                                                                                          
        CI Iterations
        -------------
                                                                                                                          
     Iter. | Average Energy | E. Change | Grad. Norm |   Time
     ----------------------------------------------------------
        1     -75.708185001     0.0e+00      3.6e-01   0:00:07
        2     -75.793118310    -9.5e-02      3.6e-02   0:00:07
        3     -75.798908867    -7.2e-03      2.0e-03   0:00:08
        4     -75.799401560    -5.8e-04      2.2e-04   0:00:24
        5     -75.799440899    -6.4e-05      2.9e-05   0:00:24

        6     -75.799446125    -9.5e-06      3.3e-06   0:00:08
        7     -75.799446623    -9.0e-07      5.5e-07   0:00:08
        8     -75.799446711    -2.5e-07      2.1e-07   0:00:08
        9     -75.799446734    -9.2e-08      6.1e-08   0:00:07
       10     -75.799446738    -2.0e-08      7.4e-09   0:00:07
       11     -75.799446739    -2.1e-09      8.5e-10   0:00:10
                                                                                                                          
** Convergence reached in 11 iterations
                                                                                                                          
        Final results
        -------------
                                                                                                                          
                                                                                                                          
* State 1
- Energy: -76.12003281964515
- S^2   : -0.02  (multiplicity = 1.0 )
- Natural orbitals
1.98832 1.96820 1.97163 1.98078 0.00049 0.02805 0.00308 0.01801 0.02650 0.01210 0.00062 0.00221
                                                                                                                          
* State 2
- Energy: -75.80828932732831
- S^2   : -0.03  (multiplicity = 0.9 )
- Natural orbitals
1.99041 1.96630 1.97922 0.99448 0.99598 0.03407 0.00351 0.00646 0.01523 0.00935 0.00243 0.00257
                                                                                                                          
* State 3
- Energy: -75.72622990160053
- S^2   : -0.03  (multiplicity = 0.9 )
- Natural orbitals
1.98982 1.97785 1.96971 0.99453 0.03391 0.99421 0.00372 0.00639 0.00215 0.01098 0.01350 0.00324
                                                                                                                          
* State 4
- Energy: -75.71537749268845
- S^2   : -0.03  (multiplicity = 0.9 )
- Natural orbitals
1.98672 1.95763 1.14165 1.98172 0.84888 0.04031 0.00280 0.01710 0.01128 0.00608 0.00347 0.00235
                                                                                                                          
* State 5
- Energy: -75.6273041535162
- S^2   : -0.03  (multiplicity = 0.9 )
- Natural orbitals
1.98504 1.95735 1.01140 1.98291 0.04058 0.97774 0.00364 0.01556 0.00427 0.00615 0.01336 0.00427
                                                                                                                          
List of oscillator strengths greather than 1e-10
                                                                                                                          
  From     to       Energy (eV)    Oscillator strength (length and velocity)
     1       2        8.48297         1.323123e-02    4.607713e-02
     1       4       11.01123         1.154562e-01    1.737591e-01
     1       5       13.40783         1.150580e-01    1.218071e-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

plt.figure(figsize=(6,4))
x = np.array(range(1,6))
plt.plot(x,Energies[:,0], label='0->1')
plt.plot(x,Energies[:,1], label='0->2')
plt.plot(x,Energies[:,2], label='0->3')
plt.plot(x,Energies[:,3], label='0->4')

plt.title('Convergence of excitation energies')
plt.xticks([1, 2, 3, 4, 5], ['CIS', 'CISD', 'CISDT', 'CISDTQ', 'FCI'])
plt.ylabel("Energies (eV)")
plt.legend()
plt.tight_layout(); plt.show()
Energies[0,:]

As we can see, we do not have a monotonous convergence towards the full CI result. In particular, CISD is actually worse than CIS. The reason for this is that CIS has a more even-handed treatment of the ground and excited states than CISD. In CIS, due to the Brillouin theorem, the ground state does not improve with the inclusion of single excitations, but the single excitations are needed to generate the dominant configuration for the (singly) excited states. Both ground and excited states are thus treated at a mean-field (HF-like) level.

However, in CISD, the ground state contains double excitation, and thus dynamical correlation, while the excited states would need up to triple excitations to both generate the singly-excited dominant configuration and correlating double excitations on top. There is an unbalance, the ground state being described more accurately than the excited states, leading to too high excitation energies. Adding the triple excitations bring correlation also for excited states, which brings the result closer to the final FCI result, and in this simple case is essentially converged.

For this reason, while CIS is sometimes used to compute excitation energies (and as we can see in the next section, is actually related to HF response), the other truncated CIs are not.

Spectra and properties
Quasi-energy determination of parameters
Spectra and properties
Algebraic-diagrammatic construction