Configuration interaction
Contents
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)
* Warning * Environment variable OMP_NUM_THREADS not set.
* Warning * Setting OMP_NUM_THREADS to 8.
* Info * Reading basis set from file: /opt/anaconda3/envs/vlxenv/lib/python3.9/site-packages/veloxchem/basis/6-31G
Molecular Basis (Atomic Basis)
================================
Basis: 6-31G
Atom Contracted GTOs Primitive GTOs
O (3S,2P) (10S,4P)
H (2S) (4S)
Contracted Basis Functions : 13
Primitive Basis Functions : 30
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 Scheme : Cauchy Schwarz + Density
ERI Screening Mode : Dynamic
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 * SAD initial guess 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.02 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.0
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)
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: 122760
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: 122760
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
Max. iterations : 40
Initial diagonalization : 200
Max subspace size : 10
Convergence thresholds:
- Energy change : 1e-08
- Residual square norm: 1e-08
Standard Davidson step
CI Iterations
-------------
Iter. | Average Energy | E. Change | Grad. Norm | Time
----------------------------------------------------------
1 -75.815781515 0.0e+00 2.6e-01 0:00:00
2 -75.894292645 -7.9e-02 2.9e-02 0:00:00
3 -75.901150187 -6.9e-03 3.5e-03 0:00:00
4 -75.901963550 -8.1e-04 6.4e-04 0:00:00
5 -75.902129581 -1.7e-04 8.7e-05 0:00:00
6 -75.902150481 -2.1e-05 1.4e-05 0:00:00
7 -75.902154386 -3.9e-06 2.6e-06 0:00:00
8 -75.902155098 -7.1e-07 4.9e-07 0:00:00
9 -75.902155213 -1.1e-07 9.5e-08 0:00:00
10 -75.902155237 -2.5e-08 1.5e-08 0:00:00
11 -75.902155242 -4.2e-09 3.1e-09 0:00:00
** Convergence reached in 11 iterations
Final results
-------------
* State 1
- Energy: -75.90215524152757
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.98787 1.98052 1.72827 1.69442 0.30428 0.26591 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: 122760
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
Max. iterations : 40
Initial diagonalization : 200
Max subspace size : 30
Convergence thresholds:
- Energy change : 1e-08
- Residual square norm: 1e-08
Standard Davidson step
CI Iterations
-------------
Iter. | Average Energy | E. Change | Grad. Norm | Time
----------------------------------------------------------
1 -75.769165226 0.0e+00 2.6e-01 0:00:00
2 -75.850094458 -8.3e-02 3.2e-02 0:00:00
3 -75.857867208 -8.5e-03 3.7e-03 0:00:00
4 -75.858732659 -9.6e-04 7.6e-04 0:00:00
5 -75.858903324 -1.9e-04 1.4e-04 0:00:00
6 -75.858929884 -3.2e-05 2.6e-05 0:00:00
7 -75.858935632 -7.0e-06 6.3e-06 0:00:00
8 -75.858937230 -2.2e-06 2.3e-06 0:00:00
9 -75.858937671 -7.4e-07 8.9e-07 0:00:00
10 -75.858937816 -2.7e-07 2.4e-07 0:00:00
11 -75.858937853 -8.3e-08 8.5e-08 0:00:00
12 -75.858937864 -2.9e-08 3.0e-08 0:00:00
13 -75.858937868 -8.9e-09 7.3e-09 0:00:00
** Convergence reached in 13 iterations
Final results
-------------
* State 1
- Energy: -75.90215524245849
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.98787 1.98052 1.72828 1.69443 0.30427 0.26590 0.00035 0.00026 0.01798 0.00528 0.00509 0.00977
* State 2
- Energy: -75.84903068223876
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.98897 0.99427 1.68096 1.95900 1.00328 0.33528 0.00093 0.00022 0.00674 0.01568 0.00573 0.00894
* State 3
- Energy: -75.82562767830413
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.98928 0.99437 1.94571 1.62697 0.40629 1.00029 0.00027 0.00075 0.00666 0.00601 0.01415 0.00925
At the end of the calculation, the CI Driver prints a few informations, but the key results can also be accessed directly in python.
print("Return dictionary:",ci_results)
print()
print("Ground state energy:",CIdrv.get_energy(0)) #Get the energy of a specific state, here state number 0
print("Array of all state energies:",CIdrv.get_energies()) #Get a python list of energies
print()
print("Ground state natural orbitals occupation numbers:",CIdrv.get_natural_occupations(0))
Return dictionary: {'orbital_space': <multipsi.OrbSpace.OrbSpace object at 0x10d44d940>, 'converged': True, 'n_iterations': 13, 'energies': array([-75.90215524, -75.84903068, -75.82562768]), 'ci_vectors': [<multipsi.CIVector.CIVector object at 0x17b3b5f70>, <multipsi.CIVector.CIVector object at 0x10d3fccd0>, <multipsi.CIVector.CIVector object at 0x10d44d970>]}
Ground state energy: -75.90215524245849
Array of all state energies: [-75.90215524245849, -75.84903068223876, -75.82562767830413]
Ground state natural orbitals occupation numbers: [1.98787204e+00 1.98051728e+00 1.72827832e+00 1.69442749e+00
3.04273735e-01 2.65895102e-01 1.79768537e-02 9.76908009e-03
5.28474000e-03 5.09315765e-03 3.50038350e-04 2.62162451e-04]
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: 729
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
Max. iterations : 40
Initial diagonalization : 200
Max subspace size : 10
Convergence thresholds:
- Energy change : 1e-08
- Residual square norm: 1e-08
Standard Davidson step
CI Iterations
-------------
Iter. | Average Energy | E. Change | Grad. Norm | Time
----------------------------------------------------------
1 -75.811173597 0.0e+00 1.6e-01 0:00:00
2 -75.856674958 -4.6e-02 6.9e-03 0:00:00
3 -75.860829076 -4.2e-03 2.4e-04 0:00:00
4 -75.860924620 -9.6e-05 5.4e-06 0:00:00
5 -75.860927029 -2.4e-06 1.7e-07 0:00:00
6 -75.860927094 -6.5e-08 5.3e-09 0:00:00
7 -75.860927097 -2.3e-09 1.2e-10 0:00:00
** Convergence reached in 7 iterations
Final results
-------------
* State 1
- Energy: -75.86092709655223
- S^2 : -0.00 (multiplicity = 1.0 )
- 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