OO-DOCI calculations

DOCI is, in essence, just another type of CI, so doing DOCI calculations is very analogous to doing FCI calculations. In the OO-DOCI case, we enhance the regular DOCI method with the use of an orbital optimizer.

# 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_cristina.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())
nuc_rep = gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value()
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 OO-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)

We now have to declare our orbital optimizer.

optimizer = gqcpy.DOCINewtonOrbitalOptimizer(onv_basis, solver, environment)

The orbital optimizer can now do the required optimization given the spinor basis and the second quantized Hamiltonian.

optimizer.optimize(spinor_basis, sq_hamiltonian)

print(optimizer.eigenvalue() + nuc_rep)

-1.137263337693952

optimized_orbitals = optimizer.eigenvector()
energies = optimizer.eigenvalues()

print(optimized_orbitals)

[-0.99360078  0.11294904]