Configuration interaction#
Multipsi provides a very general configuration code, allowing a wide variety of CI expansions, including FullCI and truncated CI (for example CI with single and double excitations CISD). Note however that the code is mostly designed around Full CI (for use in particular in complete active space CAS methods, see following chapters) and is thus not as efficient as a dedicated CIS or CISD code.
Full CI#
Let’s start with a simple FullCI of water.
MultiPsi relies on VeloxChem to define the molecule and basis set. For more informations about how to run Veloxchem, please go to the Veloxchem webpage.
Here we will start from a simple Hartree-Fock calculation in Veloxchem:
import veloxchem as vlx
water_str='''
O 0.0 0.0 0.0
H 0.0 1.4 1.1
H 0.0 -1.4 1.1
'''
water=vlx.Molecule.read_str(water_str)
basis = vlx.MolecularBasis.read(water,"6-31G")
scfdrv = vlx.ScfRestrictedDriver()
hf_results = scfdrv.compute(water, basis)
Self Consistent Field Driver Setup
====================================
Wave Function Model : Spin-Restricted Hartree-Fock
Initial Guess Model : Superposition of Atomic Densities
Convergence Accelerator : Two Level Direct Inversion of Iterative Subspace
Max. Number of Iterations : 50
Max. Number of Error Vectors : 10
Convergence Threshold : 1.0e-06
ERI Screening Threshold : 1.0e-12
Linear Dependence Threshold : 1.0e-06
* Info * Nuclear repulsion energy: 4.9444403801 a.u.
* Info * Overlap matrix computed in 0.00 sec.
* Info * Kinetic energy matrix computed in 0.00 sec.
* Info * Nuclear potential matrix computed in 0.00 sec.
* Info * Orthogonalization matrix computed in 0.00 sec.
* Info * Starting Reduced Basis SCF calculation...
* Info * ...done. SCF energy in reduced basis set: -75.641519996755 a.u. Time: 0.02 sec.
* Info * Overlap matrix computed in 0.00 sec.
* Info * Kinetic energy matrix computed in 0.00 sec.
* Info * Nuclear potential matrix computed in 0.00 sec.
* Info * Orthogonalization matrix computed in 0.00 sec.
Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change
--------------------------------------------------------------------------------------------
1 -75.641523271211 0.0000000000 0.00087825 0.00019340 0.00000000
2 -75.641525364813 -0.0000020936 0.00084045 0.00014743 0.00353953
3 -75.641525451665 -0.0000000869 0.00004595 0.00000631 0.00023152
4 -75.641525453425 -0.0000000018 0.00000898 0.00000138 0.00009522
5 -75.641525453439 -0.0000000000 0.00000040 0.00000012 0.00000630
*** SCF converged in 5 iterations. Time: 0.03 sec.
Spin-Restricted Hartree-Fock:
-----------------------------
Total Energy : -75.6415254534 a.u.
Electronic Energy : -80.5859658336 a.u.
Nuclear Repulsion Energy : 4.9444403801 a.u.
------------------------------------
Gradient Norm : 0.0000003971 a.u.
Ground State Information
------------------------
Charge of Molecule : 0.0
Multiplicity (2S+1) : 1
Magnetic Quantum Number (M_S) : 0.0
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 1:
--------------------------
Occupation: 2.000 Energy: -20.66858 a.u.
( 1 O 1s : -1.00)
Molecular Orbital No. 2:
--------------------------
Occupation: 2.000 Energy: -1.20252 a.u.
( 1 O 1s : 0.23) ( 1 O 2s : -0.53) ( 1 O 3s : -0.54)
Molecular Orbital No. 3:
--------------------------
Occupation: 2.000 Energy: -0.50810 a.u.
( 1 O 1p+1: -0.68) ( 1 O 2p+1: -0.47)
Molecular Orbital No. 4:
--------------------------
Occupation: 2.000 Energy: -0.45244 a.u.
( 1 O 1p-1: 0.45) ( 1 O 2p-1: 0.36) ( 2 H 1s : 0.16)
( 2 H 2s : 0.28) ( 3 H 1s : -0.16) ( 3 H 2s : -0.28)
Molecular Orbital No. 5:
--------------------------
Occupation: 2.000 Energy: -0.43893 a.u.
( 1 O 1p0 : 0.47) ( 1 O 2p0 : 0.38) ( 2 H 1s : 0.15)
( 2 H 2s : 0.27) ( 3 H 1s : 0.15) ( 3 H 2s : 0.27)
Molecular Orbital No. 6:
--------------------------
Occupation: 0.000 Energy: -0.00093 a.u.
( 1 O 3s : 0.23) ( 1 O 1p0 : 0.41) ( 1 O 2p0 : 0.43)
( 2 H 1s : -0.17) ( 2 H 2s : -0.51) ( 3 H 1s : -0.17)
( 3 H 2s : -0.51)
Molecular Orbital No. 7:
--------------------------
Occupation: 0.000 Energy: 0.03740 a.u.
( 1 O 1p-1: -0.43) ( 1 O 2p-1: -0.48) ( 2 H 1s : 0.18)
( 2 H 2s : 0.53) ( 3 H 1s : -0.18) ( 3 H 2s : -0.53)
Molecular Orbital No. 8:
--------------------------
Occupation: 0.000 Energy: 0.93316 a.u.
( 1 O 3s : 0.25) ( 1 O 1p0 : -0.22) ( 1 O 2p0 : 0.34)
( 2 H 1s : 0.88) ( 2 H 2s : -0.83) ( 3 H 1s : 0.88)
( 3 H 2s : -0.83)
Molecular Orbital No. 9:
--------------------------
Occupation: 0.000 Energy: 1.00695 a.u.
( 1 O 1p-1: 0.21) ( 1 O 2p-1: -0.42) ( 2 H 1s : -0.91)
( 2 H 2s : 0.95) ( 3 H 1s : 0.91) ( 3 H 2s : -0.95)
Molecular Orbital No. 10:
--------------------------
Occupation: 0.000 Energy: 1.15709 a.u.
( 1 O 1p+1: 0.94) ( 1 O 2p+1: -1.06)
Ground State Dipole Moment
----------------------------
X : 0.000000 a.u. 0.000000 Debye
Y : -0.000000 a.u. -0.000000 Debye
Z : 0.912284 a.u. 2.318794 Debye
Total : 0.912284 a.u. 2.318794 Debye
We then define an orbital space object to store the starting orbitals and define the type of CI expansion:
import multipsi as mtp
space=mtp.OrbSpace(water,scfdrv.mol_orbs)
print(space)
Active space definition:
------------------------
Number of inactive (occupied) orbitals: 5
Number of active orbitals: 0
Number of virtual orbitals: 8
This is a restricted Hartree-Fock wavefunction
By default, OrbSpace defines a restricted Hartree-Fock wavefunction (or restricted open-shell if the molecule is not a singlet). To define the full CI, we simply use the command fci().
space.fci()
print(space)
Active space definition:
------------------------
Number of inactive (occupied) orbitals: 0
Number of active orbitals: 13
Number of virtual orbitals: 0
This is a CASSCF wavefunction: CAS(10,13)
As we can see, internally the code assumes a complete active space, but since all orbitals are included in the active space, it is indeed a full CI.
We can also choose to freeze the core orbitals, as the correlation for core electrons tend to be constant and thus cancel out in relative energies.
space.fci(n_frozen=1)
print(space)
Active space definition:
------------------------
Number of inactive (occupied) orbitals: 1
Number of active orbitals: 12
Number of virtual orbitals: 0
This is a CASSCF wavefunction: CAS(8,12)
This reduces significantly the number of determinants and thus the calculation time. It is sometimes useful to know in advance the number of determinants to see if the calculation is even feasible. For this we can create a CIExpansion object:
print(mtp.CIExpansion(space))
CI expansion:
-------------
Number of determinants: 245025
While expansions as large as hundreds of billions of determinants can in principle be run using a large enough supercomputer, any expansion beyond a billion will be very expensive.
We can now run the configuration interaction:
CIdrv=mtp.CIDriver()
ci_results = CIdrv.compute(water,basis,space)
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 CASSCF wavefunction: CAS(8,12)
CI expansion:
-------------
Number of determinants: 245025
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
Max. iterations : 40
Initial diagonalization : 201
Convergence thresholds:
- Energy change : 1e-08
- Residual square norm : 1e-08
Davidson solver
- Max subspace size : 10
- Min subspace size : 6
- Standard Davidson step
CI Iterations
----------------
Iter. | Average Energy | E. Change | Grad. Norm | Subs. size | Time
----------------------------------------------------------------------
1 -75.798712443 8.5e-14 2.9e-01 1 0:00:00
2 -75.891062197 -9.2e-02 3.7e-02 2 0:00:00
3 -75.899886320 -8.8e-03 6.6e-03 3 0:00:00
4 -75.901675904 -1.8e-03 1.3e-03 4 0:00:00
5 -75.902040583 -3.6e-04 2.6e-04 5 0:00:00
6 -75.902112783 -7.2e-05 6.6e-05 6 0:00:00
7 -75.902134198 -2.1e-05 2.9e-05 7 0:00:00
8 -75.902147466 -1.3e-05 1.5e-05 8 0:00:00
9 -75.902153008 -5.5e-06 5.3e-06 9 0:00:00
10 -75.902154497 -1.5e-06 1.4e-06 10 0:00:00
11 -75.902154933 -4.4e-07 4.2e-07 6 0:00:00
12 -75.902155099 -1.7e-07 2.3e-07 7 0:00:00
13 -75.902155197 -9.8e-08 1.1e-07 8 0:00:00
14 -75.902155229 -3.2e-08 3.7e-08 9 0:00:00
15 -75.902155239 -9.8e-09 1.0e-08 10 0:00:00
16 -75.902155242 -2.6e-09 2.4e-09 6 0:00:00
** Convergence reached in 16 iterations
Final results
-------------
* State 1
- S^2 : 0.00 (multiplicity = 1.0 )
- Energy : -75.90215524166567
- Natural orbitals
1.98787 1.98052 1.72828 1.69442 0.30428 0.26589 0.00035 0.00026 0.01798 0.00528 0.00509 0.00977
If we are interested in several states and not just the ground state, we simply provide the number of states to the compute function:
ci_results = CIdrv.compute(water,basis,space, n_states = 3) #We want 3 states now
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 CASSCF wavefunction: CAS(8,12)
CI expansion:
-------------
Number of determinants: 245025
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
Max. iterations : 40
Initial diagonalization : 203
Convergence thresholds:
- Energy change : 1e-08
- Residual square norm : 1e-08
Davidson solver
- Max subspace size : 30
- Min subspace size : 18
- Standard Davidson step
CI Iterations
----------------
Iter. | Average Energy | E. Change | Grad. Norm | Subs. size | Time
----------------------------------------------------------------------
1 -75.753160108 2.8e-14 3.0e-01 3 0:00:00
2 -75.847512855 -9.4e-02 3.9e-02 6 0:00:00
3 -75.856952388 -9.4e-03 6.6e-03 9 0:00:00
4 -75.858411849 -1.5e-03 1.4e-03 12 0:00:00
5 -75.858787899 -3.8e-04 4.3e-04 15 0:00:00
6 -75.858885366 -9.7e-05 1.2e-04 18 0:00:00
7 -75.858913187 -2.8e-05 5.5e-05 21 0:00:00
8 -75.858926997 -1.4e-05 4.0e-05 24 0:00:00
9 -75.858934626 -7.6e-06 1.5e-05 27 0:00:00
10 -75.858936904 -2.3e-06 3.4e-06 30 0:00:00
11 -75.858937562 -6.6e-07 6.6e-07 18 0:00:00
12 -75.858937787 -2.3e-07 4.7e-07 21 0:00:00
13 -75.858937971 -1.8e-07 8.2e-07 24 0:00:00
14 -75.858938284 -3.1e-07 2.4e-06 27 0:00:00
15 -75.858938912 -6.3e-07 4.6e-06 30 0:00:00
16 -75.858939420 -5.1e-07 1.6e-06 18 0:00:00
17 -75.858939715 -2.9e-07 1.8e-06 20 0:00:00
18 -75.858940203 -4.9e-07 5.0e-06 22 0:00:00
19 -75.858948175 -8.0e-06 2.9e-04 23 0:00:00
20 -75.861855572 -2.9e-03 1.0e-02 24 0:00:00
21 -75.862946904 -1.1e-03 2.7e-03 25 0:00:00
22 -75.863340339 -3.9e-04 1.6e-03 26 0:00:00
23 -75.863577793 -2.4e-04 9.3e-04 27 0:00:00
24 -75.863690332 -1.1e-04 3.1e-04 28 0:00:00
25 -75.863726009 -3.6e-05 1.1e-04 29 0:00:00
26 -75.863742183 -1.6e-05 5.6e-05 30 0:00:00
27 -75.863750555 -8.4e-06 3.3e-05 18 0:00:00
28 -75.863756521 -6.0e-06 3.2e-05 19 0:00:00
29 -75.863761621 -5.1e-06 2.4e-05 20 0:00:00
30 -75.863764000 -2.4e-06 7.6e-06 21 0:00:00
31 -75.863764626 -6.3e-07 2.0e-06 22 0:00:00
32 -75.863764809 -1.8e-07 4.9e-07 23 0:00:00
33 -75.863764852 -4.3e-08 1.1e-07 24 0:00:00
34 -75.863764864 -1.2e-08 3.6e-08 25 0:00:00
35 -75.863764869 -5.0e-09 2.0e-08 26 0:00:00
36 -75.863764873 -3.6e-09 1.8e-08 27 0:00:00
37 -75.863764876 -3.1e-09 1.4e-08 28 0:00:00
38 -75.863764878 -2.0e-09 6.8e-09 29 0:00:00
** Convergence reached in 38 iterations
Final results
-------------
* State 1
- S^2 : 0.00 (multiplicity = 1.0 )
- Energy : -75.90215524166568
- Natural orbitals
1.98787 1.98052 1.72828 1.69442 0.30428 0.26589 0.00035 0.00026 0.01798 0.00528 0.00509 0.00977
* State 2
- S^2 : 0.00 (multiplicity = 1.0 )
- Energy : -75.8490306824275
- Natural orbitals
1.98897 0.99427 1.68092 1.95899 1.00328 0.33533 0.00093 0.00022 0.00674 0.01568 0.00573 0.00894
* State 3
- S^2 : 2.00 (multiplicity = 3.0 )
- Energy : -75.84010870967239
- Natural orbitals
1.98803 0.99430 1.75098 1.75237 0.51335 0.96101 0.00036 0.00110 0.00677 0.01142 0.01196 0.00835
At the end of the calculation, the CI Driver prints a few informations, but the key results can also be accessed directly from the return dictionary.
print("Ground state energy:",ci_results["energies"][0]) #Get the energy of a specific state, here state number 0
print("Array of all state energies:",ci_results["energies"]) #Get a python list of energies
print()
print("Ground state natural orbitals occupation numbers:",ci_results["natural_occupations"][0])
Ground state energy: -75.90215524166568
Array of all state energies: [-75.90215524 -75.84903068 -75.84010871]
Ground state natural orbitals occupation numbers: [1.98787234e+00 1.98051697e+00 1.72827796e+00 1.69442049e+00
3.04282587e-01 2.65893737e-01 3.50022813e-04 2.62155960e-04
1.79768836e-02 5.28475859e-03 5.09315296e-03 9.76894891e-03]
Truncated CI#
To define a truncated CI, one simply needs to define the truncated expansion in the OrbSpace object.
CIS and CISD expansions have a specific command, and higher order truncations can be done by using CI(n) with n the excitation order.
space=mtp.OrbSpace(water,scfdrv.mol_orbs)
space.cisd(n_frozen=1) #Equivalent to space.ci(2,n_frozen=1)
ci_results = CIdrv.compute(water,basis,space)
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 6 8
2 12 8 8
CI expansion:
-------------
Number of determinants: 1425
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
Max. iterations : 40
Initial diagonalization : 201
Convergence thresholds:
- Energy change : 1e-08
- Residual square norm : 1e-08
Davidson solver
- Max subspace size : 10
- Min subspace size : 6
- Standard Davidson step
CI Iterations
----------------
Iter. | Average Energy | E. Change | Grad. Norm | Subs. size | Time
----------------------------------------------------------------------
1 -75.789203149 4.3e-14 2.2e-01 1 0:00:00
2 -75.852063627 -6.3e-02 1.2e-02 2 0:00:00
3 -75.859491109 -7.4e-03 1.9e-03 3 0:00:00
4 -75.860663601 -1.2e-03 1.9e-04 4 0:00:00
5 -75.860816402 -1.5e-04 6.8e-05 5 0:00:00
6 -75.860899549 -8.3e-05 2.6e-05 6 0:00:00
7 -75.860922945 -2.3e-05 4.1e-06 7 0:00:00
8 -75.860926336 -3.4e-06 5.8e-07 8 0:00:00
9 -75.860926930 -5.9e-07 1.7e-07 9 0:00:00
10 -75.860927072 -1.4e-07 3.2e-08 10 0:00:00
11 -75.860927095 -2.3e-08 2.8e-09 6 0:00:00
12 -75.860927096 -1.5e-09 1.8e-10 7 0:00:00
** Convergence reached in 12 iterations
Final results
-------------
* State 1
- S^2 : 0.00 (multiplicity = 1.0 )
- Energy : -75.86092709649164
- Natural orbitals
1.99082 1.98675 1.83937 1.82537 0.17488 0.15640 0.00028 0.00024 0.01214 0.00338 0.00336 0.00702