Molecular FCI calculations

Molecular FCI calculations#

In this example, we will calculate the ground state energy for H2 with an STO-3G basisset at the FCI level of theory.

Setting up the molecular Hamiltonian#

# Force the local gqcpy to be imported
import sys
sys.path.insert(0, '../../build/gqcpy/')

import gqcpy
import numpy as np

We start by creating the molecule and an associated spinor basis.

molecule = gqcpy.Molecule.ReadXYZ("../../gqcp/tests/data/h2_szabo.xyz" , 0)  # create a neutral molecule
N = molecule.numberOfElectrons()
N_P = molecule.numberOfElectronPairs()

spin_orbital_basis = gqcpy.RSpinOrbitalBasis_d(molecule, "STO-3G")
K = spin_orbital_basis.numberOfSpatialOrbitals()

Unfortunately, that spinor basis is not orthonormal, which is a requirement for our CI calculations.

print(spin_orbital_basis.quantize(gqcpy.OverlapOperator()).parameters())

[[0.99999999 0.65931816]
 [0.65931816 0.99999999]]

We will proceed by finding the canonical RHF orbitals, which correspond to an orthonormal spinor basis.

fq_hamiltonian = gqcpy.FQMolecularHamiltonian(molecule)
hamiltonian = spin_orbital_basis.quantize(fq_hamiltonian)
S = spin_orbital_basis.quantize(gqcpy.OverlapOperator())

environment = gqcpy.RHFSCFEnvironment_d.WithCoreGuess(N, hamiltonian, S)
solver = gqcpy.RHFSCFSolver_d.Plain()
objective = gqcpy.DiagonalRHFFockMatrixObjective_d(hamiltonian)  # use the default threshold of 1.0e-08
rhf_parameters = gqcpy.RHF_d.optimize(objective, solver, environment).groundStateParameters()

Using this orthonormal spinor basis, we can quantize the molecular Hamiltonian.

spin_orbital_basis.transform(rhf_parameters.expansion())
hamiltonian.transform(rhf_parameters.expansion())

Specifying the ‘F’ in FCI#

FCI, in essence, uses a linear expansion in the full (spin-resolved) Fock space. In GQCP (gqcpy), we have access to a class that represents an ONV basis for that Fock space:

onv_basis = gqcpy.SpinResolvedONVBasis(K, N_P, N_P)  # number of spatial orbitals, number of alpha-electrons, number of beta-electrons

Dense FCI calculations#

In order to do a dense FCI calculations, GQCP uses its own framework for mathematical optimization. We’ll have to create:

a dense eigenvalue problem solver
an associated eigenvalue problem environment

solver = gqcpy.EigenproblemSolver.Dense_d()
environment = gqcpy.CIEnvironment.Dense(hamiltonian, onv_basis)

The encompassing QCMethod is CI (for any type of configuration interaction), so we’ll use that to find the optimized and associated wave function model: a LinearExpansion. Internally, QCMethods return a QCStructure, which wraps energies and optimized parameters for the ground state and excited states.

qc_structure = gqcpy.CI(onv_basis).optimize(solver, environment)

electronic_energy = qc_structure.groundStateEnergy()
print(electronic_energy)

-1.8515616052384547

print(qc_structure.groundStateParameters().coefficients())

[-9.93627296e-01 -1.96580806e-16 -1.28531973e-16  1.12715560e-01]

In order to find the total energy, we have to account for the nuclear repulsion. In GQCP, we use different kinds of first-quantized operators, and NuclearRepulsionOperator is one of them.

energy = electronic_energy + gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value()
print(energy)

-1.1372759430763966

Iterative FCI calculations (with a Davidson algorithm)#

In order to tackle larger problems for which we cannot store the Hamiltonian matrix (in the ONV basis) in memory, we must use an iterative diagonalization procedure. In GQCP, we have implemented Davidson’s algorithm. In order to use it to do FCI calculations, we must again create a solver and associated environment instance.

Since the Davidson algorithm is iterative, we will have to supply an initial guess vector. Since we’re working in the canonical RHF basis, a linear expansion that has a 1 at the position of the RHF determinant is a good initial guess.

x0 = gqcpy.LinearExpansion_SpinResolved.HartreeFock(onv_basis).coefficients()

solver_davidson = gqcpy.EigenproblemSolver.Davidson()
environment_davidson = gqcpy.CIEnvironment.Iterative(hamiltonian, onv_basis, x0)

qc_structure_davidson = gqcpy.CI(onv_basis).optimize(solver_davidson, environment_davidson)

We can now verify that both types of solvers find the same, optimized linear wave function model and associated energy.

electronic_energy_davidson = qc_structure_davidson.groundStateEnergy()
print(electronic_energy_davidson)

-1.8515616052384547

print(qc_structure_davidson.groundStateParameters().coefficients())

[-9.93627296e-01 -1.24364039e-16 -1.42130330e-16  1.12715560e-01]

‘Unrestricted’ FCI calculations#

In the previous examples, we have expressed the Hamiltonian in a restricted spin-orbital basis. Furthermore, CI calculations should be independent of the underlying one-electron (spinor/spin-orbital) basis. Can we reproduce this result, for an unrestricted spin-orbital basis?

We’ll start by creating an unrestricted spin-orbital (AO) basis, and Löwdin-orthonormalizing it.

u_spin_orbital_basis = gqcpy.USpinOrbitalBasis_d(molecule, "STO-3G")
u_spin_orbital_basis.lowdinOrthonormalize()

Now that we have an orthonormal unrestricted spin-orbital basis, we may quantize the molecular Hamiltonian in it.

u_hamiltonian = u_spin_orbital_basis.quantize(fq_hamiltonian)

We can now reuse the full spin-resolved ONV basis, and proceed to initialize an environment and solver.

solver = gqcpy.EigenproblemSolver.Dense_d()
environment = gqcpy.CIEnvironment.Dense(u_hamiltonian, onv_basis)

qc_structure = gqcpy.CI(onv_basis).optimize(solver, environment)

electronic_energy = qc_structure.groundStateEnergy()
print(electronic_energy)

-1.8515616052384558

This is exactly the result we obtain from the restricted case!