FanCI Quickstart

FanCI Quickstart#

FanCI is a Python library for computing various post-Hartree-Fock methods (of the Flexible Ansätze for N-electron CI, or “FanCI”, type) using the PyCI library. The main objective of FanCI is to provide a simple and flexible framework to implement and test new CI methods. FanCI is not intended to be a high performance library, but rather a tool to test new ideas and methods.

Supported FanCI methods#

FanCI wavefunctions can be based on the open-shell fullci_wfn class, the closed-shell doci_wfn class, or the generalized genci_wfn class of PyCI.

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

Available methods include:

Closed-shell#

Open-shell#

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.

# Import modules
import pyci
from pyci.fanci import AP1roG, pCCDS
# optional
import numpy as np
from pyci.test import datafile

# System information
filename = datafile("lih_sto6g.fcidump")
ham = pyci.hamiltonian(filename)
e_dict = {}

Computing reference (Hatree-Fock) energy#

The following code shows how to compute the reference (Hartree-Fock) energy (minus the constant (nuclear-nuclear repulsion) term) using the PyCI package.

# Get Hartree-Fock energy
hf_wfn = pyci.doci_wfn(ham.nbasis, 2, 2)
hf_wfn.add_hartreefock_det()
hf_op = pyci.sparse_op(ham, hf_wfn)
e_dict["HF"] = hf_op.solve(n=1)[0][0] - ham.ecore

AP1roG#

The following code shows how to create an initial-guess array for the nonlinear FanCI parameters of the wavefunction. In the case of AP1roG, it is usually ok to leave the wavefunction parameters at zero. The last parameter, the energy, is set to the Hartree-Fock energy.

# Initialize AP1roG instance
ap1rog = AP1roG(ham, 2)

# Make initial guess
ap1_params = np.zeros(ap1rog.nparam, dtype=pyci.c_double)
ap1_params[-1] = e_dict["HF"]

One can use the built-in optimize method to optimize the wavefunction parameters using a SciPy optimizer. If one would like to use a different optimizer, the objective function and its Jacobian are available for use as obj.compute_objective(x) and obj.compute_jacobian(x).

Optimize the AP1roG wavefunction:

# Optimize wavefunction
ap1rog_results = ap1rog.optimize(ap1_params, use_jac=True)
e_dict["AP1roG"] = ap1rog_results.x[-1]

pCCD+S#

The following code shows how the optimized parameters from a simpler wavefunction (AP1roG) can be used as the initial guess to a more complicated and more accurate wavefunction (pCCD+S).

Use the AP1roG wave function to generate an initial guess for the pCCD+S wave function; optimize it:

# Initialize pCCD+S instance
pccds = pCCDS(ham, 2, 2)

# Make initial guess from AP1roG params
pccds_params = np.zeros(pccds.nparam, dtype=pyci.c_double)
pccds_params[:pccds.wfn.nocc_up * pccds.wfn.nvir_up] = ap1rog_results.x[:-1]
pccds_params[-1] = ap1rog_results.x[-1]

# Optimize wavefunction
pccds_results = pccds.optimize(pccds_params, use_jac=False)
e_dict["pCCD+S"] = pccds_results.x[-1]

Compare results#

Compare the energies and verify that \(E_\text{HF}\) > \(E_\text{AP1roG}\) > \(E_\text{pCCD+S}\):

# Print energies from various methods
print(f"{'METHOD':>8s}, {'ENERGY':>16s}")
for name, energy in e_dict.items():
    print(f"{name:>8s}, {energy:>16.9e}")

  METHOD,           ENERGY
      HF, -8.947289175e+00
  AP1roG, -8.963531034e+00
  pCCD+S, -8.963613544e+00

Computing reduced density matrices#

PyCI can be used to approximately evaluate the 1- and 2- electron reduced density matrices of the FanCI wavefunction object. For example, here we use a larger, although still much smaller than FullCI, projection space (seniority-(0,2) + excitation-(0,1,2,3) CI) over which to compute the CI coefficients \(c_i = \langle\psi_i\mid\Psi_\text{pCCD+S}\rangle\), which are then used to approximately (but accurately) evaluate the pCCD+S RDMs in a reasonable manner.

# Make larger projection space over which to evaluate RDMs
proj_wfn = pyci.fullci_wfn(ham.nbasis, 2, 2)
# Add seniority- 0 and 2 determinants
pyci.add_seniorities(proj_wfn, 0, 2)
# Add singly, doubly, triply excited determinants
pyci.add_excitations(proj_wfn, 0, 1, 2, 3)

# Evaluate coefficients with optimized pCCD+S
coeffs = pccds.compute_overlap(pccds_results.x[:-1], proj_wfn.to_occ_array())

# Compute RDMs using larger projection space and pCCD+S coefficients
d1, d2 = pyci.compute_rdms(proj_wfn, coeffs)
rdm1, rdm2 = pyci.spinize_rdms(d1, d2)

Validate the RDMs; ensure that \(\text{tr}\left(\gamma\right) \approx n_\text{elec}\) and that \(\frac{1}{2}\sum_{pq}{\Gamma_{pqpq}} \approx n_\text{pair}\):

# Validate the RDMs
print(f"Number of electrons = {np.trace(rdm1):.1f}")
print(f"Number of pairs     = {np.einsum('pqpq', rdm2) / 2.0:.1f}")

Number of electrons = 4.1
Number of pairs     = 6.1