DOCI calculations

DOCI is, in essence, just another type of CI, so doing DOCI calculations is very analogous to doing FCI calculations.

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

import gqcpy
import numpy as np

Preparing the Hamiltonian in the canonical RHF basis

As always, we’ll first set up the Hamiltonian in an orthonormal spinor basis. For this example, we’ll use the canonical RHF basis

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

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

S = spinor_basis.quantize(gqcpy.OverlapOperator())
sq_hamiltonian = spinor_basis.quantize(gqcpy.FQMolecularHamiltonian(molecule))

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

C = rhf_parameters.expansion()
spinor_basis.transform(C)
sq_hamiltonian.transform(C)

Dense DOCI calculations

For FCI, we would create the full spin-resolved ONV basis. For DOCI, we need the full seniority-zero ONV basis.

onv_basis = gqcpy.SeniorityZeroONVBasis(K, N_P)  # number of spatial orbitals, number of electron pairs

The subsequent steps are very analogous. We create a solver and associated environment, and use the CI QCMethod to find the optimal parameters and energies.

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

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

electronic_energy = qc_structure.groundStateEnergy()
print(electronic_energy)

-1.8515616052384543

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

[-0.9936273   0.11271556]

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

-1.1372759430763961

‘Davidson’ DOCI calculations

Here, we’ll use a Davidson solver to find the ground state parameters. Since this is an iterative procedure, we’ll have to supply an initial guess.

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

We’ll then feed this initial guess to the environment, and let the CI method use the solver to find the optimized parameters (expansion coefficients).

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

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

electronic_energy_davidson = qc_structure_davidson.groundStateEnergy()
print(electronic_energy_davidson)

-1.8515616052384543

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

[-0.9936273   0.11271556]