PyCI Quickstart

PyCI Quickstart#

PyCI is a flexible quantum chemistry Configuration Interaction library for Python 3. The main objective of PyCI is to provide a simple and flexible framework to implement and test new CI methods. PyCI is not intended to be a high performance library, but rather a tool to test new ideas and methods.

Supported CI methods#

Through its main classes, PyCI supports the pair-occupied spatial orbital doci_wfn, orthogonal spin-orbital fullci_wfn, and generalized orbital genci_wfn wavefunction symmetries. The package also supports the use of the restricted and generalized Hamiltonians (unrestricted Hamiltonians can be constructed as a special case of generalized Hamiltonian). These features, combined with routines to control the type of configurations added, according to their excitation level or seniority, allows the user to easily program any CI method.

In the following tutorial we will look at the following wavefunction models using Be as our toy model:

  • Full CI

  • CISD

  • GKCI

  • HCI

  • Seniority-truncated CI

Loading data and defining the Hamiltonian#

First, we need to load the integrals data from a quantum chemistry package. PyCI supports the use of files in the FCIDUMP file format. The following code shows how to load the integrals from a FCIDUMP file and define the hamiltonian object.

import pyci
import numpy as np
from pyci.test import datafile  # optional, data file location for tests

# load system information
filename = datafile("be_ccpvdz.fcidump")
occs = (2,2)                        # (alpha, beta) occupation numbers
ham = pyci.hamiltonian(filename)    # initialize second-quantized operator instance
e_dict = {}

Computing reference (Hatree-Fock) energy#

The following code shows how to compute the reference (Hartree-Fock) energy using the PyCI package.

# contruct empty fci wave function class instance from # of basis functions and occupation
wfn0 = pyci.fullci_wfn(ham.nbasis, *occs)
wfn0.add_hartreefock_det()        # add Hartree-Fock determinant to wave function

# initialize sparse matrix operator (hamiltonian into wave function)
op = pyci.sparse_op(ham, wfn0)

# solve for the lowest eigenvalue and eigenvector
e_vals, e_vecs0 = op.solve(n=1, tol=1.0e-9)
e_dict["HF"] = (len(wfn0), e_vals[0])

print("Ground state wave function: HF")
print("Number of determinants: {}".format(len(wfn0)))
print("Energy: {}".format(e_vals[0]))

Ground state wave function: HF
Number of determinants: 1
Energy: -14.572337630953388

Full CI#

Once the integrals are loaded, in the hamiltonian object, PyCI offers a simple and flexible framework to implement any CI method following the next steps:

  1. Define the wavefunction model (the configurations to be included in the wavefunction).

  2. Find the representation of the Hamiltonian in the selected configuration space.

  3. Obtain the desired eigenvalues and eigenvectors of the Hamiltonian.

The Full CI method gives the exact solution of the electronic Schrödinger equation within a given basis set. It is given by the linear combination of all possible configurations in a given basis set. The following code shows how to define a Full CI wavefunction and compute the energy for the lowest singlet state of the beryllium atom.

wfn1 = pyci.fullci_wfn(ham.nbasis, *occs) 
wfn1.add_all_dets()                          # add all determinants to the wave function

# Solve the CI matrix problem
op = pyci.sparse_op(ham, wfn1)    

e_vals, e_vecs1 = op.solve(n=1, tol=1.0e-9)
e_dict["Full-CI"] = (len(wfn1), e_vals[0])

print("Ground state wave function: Full CI")
print(f"Number of determinants: {len(wfn1)}")
print(f"Full CI energy: {e_vals[0]:.8f} a.u.")

Ground state wave function: Full CI
Number of determinants: 8281
Full CI energy: -14.61740951 a.u.

CISD#

The CISD method is a truncated version of the Full CI method that includes all single and double excitations. These calculations can be achieved by truncating the Full CI space to the subset of single and double excitations. The following code shows how to define a CISD wavefunction and compute the energy for the lowest singlet state of the beryllium atom.

# Create a CISD wave function
excitations = (0, 1, 2)     # excitations to include (0 = reference, 1 = single, 2 = double)
wfn2 = pyci.fullci_wfn(ham.nbasis, *occs)
pyci.add_excitations(wfn2, *excitations, ref=None)

op = pyci.sparse_op(ham, wfn2)
e_vals, e_vecs2 = op.solve(n=1, tol=1.0e-9)
e_dict["CISD"] = (len(wfn2), e_vals[0])

print("Ground state wave function: CISD")
print(f"Number of determinants: {len(wfn2)}")
print(f"CISD energy: {e_vals[0]:.8f} a.u.")

Ground state wave function: CISD
Number of determinants: 757
CISD energy: -14.61735579 a.u.

Seniority-truncated CI#

Another way to truncate the Full CI space is by truncating the space according to the seniority quantum number. The seniority quantum number is defined as the number of unpaired electrons in a given configuration. The following code shows how to define a seniority-truncated (0,2) CI wavefunction and compute the energy for the lowest singlet state of the beryllium atom.

# Defining a seniority truncated CI wave function
seniorities = [0, 2]
wfn3 = pyci.fullci_wfn(ham.nbasis, *occs)
pyci.add_seniorities(wfn3, *seniorities)

op = pyci.sparse_op(ham, wfn3)
e_vals, e_vecs3 = op.solve(n=1, tol=1.0e-9)
e_dict["SenTrunc-CI"] = (len(wfn3), e_vals[0])

print("Ground state wave function: Seniority Truncated CI (0,2)")
print(f"Number of determinants: {len(wfn3)}")
print(f"Seniority Truncated CI energy: {e_vals[0]:.8f} a.u.")

Ground state wave function: Seniority Truncated CI (0,2)
Number of determinants: 2275
Seniority Truncated CI energy: -14.61722574 a.u.

Griebel-Knapeck CI (GKCI)#

The Griebel-Knapeck CI method is a CI method which tries to include the determinants such that the CI truncation error is equal to the basis set truncation error. For this, the selection of the determinants to include in the expansion is based on a cost function. The following code shows how to define a GKCI wavefunction and compute the energy for the lowest singlet state of the beryllium atom.

wfn4 = pyci.fullci_wfn(ham.nbasis, *occs)
pyci.add_gkci(wfn4, t=-0.5, p=1.0, mode="cntsp", dim=3, energies=None, width=None)

op = pyci.sparse_op(ham, wfn4)
e_vals, e_vecs4 = op.solve(n=1, tol=1.0e-9)
e_dict["GKCI"] = (len(wfn4), e_vals[0])

print("Ground state wave function: Generalized Seniority Truncated CI (t=-0.5, p=1.0)")
print(f"Number of determinants: {len(wfn4)}")
print(f"Generalized Seniority Truncated CI energy: {e_vals[0]:.8f} a.u.")

Ground state wave function: Generalized Seniority Truncated CI (t=-0.5, p=1.0)
Number of determinants: 169
Generalized Seniority Truncated CI energy: -14.61684259 a.u.

Heat Bath CI#

The Heat Bath CI starts with the Hartree-Fock determinant then use a sampling method to generate determinants connected to the HF determinant with a Hamiltonian matrix element larger than a given threshold. The algorithm is repeated until convergence, each time using the lowest eigenvalue in the selected subspace as the new reference determinant. The following code shows how to define a Heat Bath CI wavefunction and compute the energy for the lowest singlet state of the beryllium atom.

wfn5 = pyci.fullci_wfn(ham.nbasis, *occs)

# Add Hartree-Fock determinant
wfn5.add_hartreefock_det()
dets_added = 1

# Create CI matrix operator and initial Hartree-Fock solution
op = pyci.sparse_op(ham, wfn5)
e_vals, e_vecs5 = op.solve(n=1, tol=1.0e-9)

# Run HCI iterations
niter = 0
eps = 5.0e-4
while dets_added:
    # Add connected determinants to wave function via HCI
    dets_added = pyci.add_hci(ham, wfn5, e_vecs5[0], eps=eps)
    # Update CI matrix operator
    op.update(ham, wfn5)
    # Solve CI matrix problem
    e_vals, e_vecs5 = op.solve(n=1, tol=1.0e-9)
    niter += 1
e_dict["HCI"] = (len(wfn5), e_vals[0])

print("Ground state wave function: HCI")
print(f"Number of determinants: {len(wfn5)}")
print(f"HCI energy: {e_vals[0]:.8f} a.u.")
print(f"Number of iterations used: {niter}")

Ground state wave function: HCI
Number of determinants: 282
HCI energy: -14.61740392 a.u.
Number of iterations used: 4

The following code shows the performance of all the methods discussed in this tutorial in recovering the correlation energy of the beryllium atom.

max_corr = e_dict["Full-CI"][1] - e_dict["HF"][1]
e_ref = e_dict["HF"][1]

print(f"{'Model':<15} {'# Dets':>10} {'Energy [a.u.]':>10} {'E_corr [%]':>10}")
for key in e_dict:
    print(f"{key:<15} {e_dict[key][0]:>10} {e_dict[key][1]:>10.5f} {100*(e_dict[key][1]-e_ref)/max_corr:>10.2f}%")

Model               # Dets Energy [a.u.] E_corr [%]
HF                       1  -14.57234      -0.00%
Full-CI               8281  -14.61741     100.00%
CISD                   757  -14.61736      99.88%
SenTrunc-CI           2275  -14.61723      99.59%
GKCI                   169  -14.61684      98.74%
HCI                    282  -14.61740      99.99%

Overlap#

PyCI also provides a simple framework to compute the overlap between two wavefunction objects. The following code shows how to compute the overlap between the different wavefunction models discussed in this tutorial.

ovl0 = pyci.compute_overlap(wfn0, wfn1, e_vecs0, e_vecs1)
ovl1 = pyci.compute_overlap(wfn1, wfn1, e_vecs1, e_vecs1)
ovl2 = pyci.compute_overlap(wfn2, wfn1, e_vecs2, e_vecs1)
ovl3 = pyci.compute_overlap(wfn3, wfn1, e_vecs3, e_vecs1)
ovl4 = pyci.compute_overlap(wfn4, wfn1, e_vecs4, e_vecs1)
ovl5 = pyci.compute_overlap(wfn5, wfn1, e_vecs5, e_vecs1)

print(f"<HF|FCI>   =  {ovl0:.5f}")
print(f"<FCI|FCI>  = {ovl1:.5f}")
print(f"<CISD|FCI> = {ovl2:.5f}")
print(f"<SCI|FCI>  = {ovl3:.5f}")
print(f"<GKCI|FCI> = {ovl4:.5f}")
print(f"<HCI|FCI>  = {ovl5:.5f}")

<HF|FCI>   =  0.95238
<FCI|FCI>  = 1.00000
<CISD|FCI> = 1.00000
<SCI|FCI>  = -0.99998
<GKCI|FCI> = 0.99995
<HCI|FCI>  = 1.00000

Computing reduced density matrices#

PyCI also provides a way to compute the reduced density matrices for a given wavefunction object. The following code shows how to compute the one- and two-particle reduced density matrices for the Full CI wavefunction of the beryllium atom.

d1, d2 = pyci.compute_rdms(wfn5, e_vecs5[0])
rdm1, rdm2 = pyci.spinize_rdms(d1, d2)

print("RDM1:")
print(f"Shape = {rdm1.shape}")
print(f"Number of electrons = {np.trace(rdm1):.1f}")


print("RDM2:")
print(f"Shape = {rdm2.shape}")
print(f"Number of electron pairs = {np.einsum('pqpq', rdm2) / 2.0:.1f}")

RDM1:
Shape = (28, 28)
Number of electrons = 4.0
RDM2:
Shape = (28, 28, 28, 28)
Number of electron pairs = 6.0