AP1roG calculations#
# Force the local gqcpy to be imported
import sys
sys.path.insert(0, '../../build/gqcpy/')
import gqcpy
import numpy as np
np.set_printoptions(precision=3, linewidth=120)
Setting up the molecular Hamiltonian in the canonical RHF spin-orbital basis#
The AP1roG geminal pairing scheme is set up in the canonical RHF (spin-orbital) basis, so we’ll have to do an RHF calculation first. There are plenty of other examples that focus on Hartree-Fock calculations in more detail.
molecule = gqcpy.Molecule.ReadXYZ("../../gqcp/tests/data/ch4_crawdad.xyz", 0) # '0': Create a neutral molecule.
N = molecule.numberOfElectrons()
N_P = molecule.numberOfElectronPairs()
spin_orbital_basis = gqcpy.RSpinOrbitalBasis_d(molecule, "STO-3G") # The basis of the atomic spin-orbitals.
K = spin_orbital_basis.numberOfSpatialOrbitals()
S = spin_orbital_basis.quantize(gqcpy.OverlapOperator())
Since the initial spin-orbital basis is the AO basis, the Hamiltonian is expressed in the non-orthogonal AOs. This is exactly what we need to start an RHF calculation.
hamiltonian = gqcpy.FQMolecularHamiltonian(molecule)
sq_hamiltonian = spin_orbital_basis.quantize(hamiltonian)
environment = gqcpy.RHFSCFEnvironment_d.WithCoreGuess(N, sq_hamiltonian, S)
solver = gqcpy.RHFSCFSolver_d.DIIS()
objective = gqcpy.DiagonalRHFFockMatrixObjective_d(sq_hamiltonian)
Using this objective makes sure that the optimized expansion coefficients yield a diagonal Fock matrix, so the solution represents the canonical RHF spinor basis.
rhf_parameters = gqcpy.RHF_d.optimize(objective, solver, environment).groundStateParameters()
C = rhf_parameters.expansion()
Since we have the canonical RHF spinor expansion coefficients now, we can transform the underlying spin-orbital basis and then re-quantize the molecular Hamiltonian, in order to let both instances be in-sync with their basis transformations.
However, gqcpy offers a different approach, transforming the spin-orbital basis and the Hamiltonian with their member APIs
spin_orbital_basis.transform(C)
sq_hamiltonian.transform(C)
Right now, the spin-orbital basis and Hamiltonian are expressed in the canonical RHF spin-orbitals.
AP1roG calculations#
AP1roG has been formulated as a QCModel/QCMethod. Therefore, we should create a solver (which is able to solve the AP1roG PSEs), an associated environment and an objective. For AP1roG, the solver is a non-linear equation solver.
environment = gqcpy.PSEnvironment.AP1roG(sq_hamiltonian, N_P)
solver = gqcpy.NonLinearEquationSolver.Newton()
qc_structure = gqcpy.AP1roG(sq_hamiltonian, N_P).optimize(solver, environment)
print(qc_structure.groundStateEnergy())
-53.25026301675117
print(qc_structure.groundStateParameters().geminalCoefficients().asMatrix())
[[ 1.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 -7.790e-04 -7.788e-04 -7.793e-04 -1.948e-03]
[ 0.000e+00 1.000e+00 0.000e+00 0.000e+00 0.000e+00 -1.416e-02 -1.413e-02 -1.421e-02 -2.450e-02]
[ 0.000e+00 0.000e+00 1.000e+00 0.000e+00 0.000e+00 -1.036e-02 -5.941e-02 -2.103e-02 -1.861e-02]
[ 0.000e+00 0.000e+00 0.000e+00 1.000e+00 0.000e+00 -4.347e-02 -1.461e-02 -3.162e-02 -1.871e-02]
[ 0.000e+00 0.000e+00 0.000e+00 0.000e+00 1.000e+00 -3.626e-02 -1.663e-02 -3.670e-02 -1.872e-02]]
vAP1roG calculations#
We have also extended the AP1roG method to be variationally optimized, resulting in the vAP1roG method. In short, what this method does is analogous to AP1roG, but after determining the optimal geminal coefficients, a set of optimal Lagrange multipliers is also searched for.
Since these Lagrange multipliers are determined through solving a linear equation, we will have to supply a linear equations solver to the vAP1roG method.
non_linear_environment = gqcpy.PSEnvironment.AP1roG(sq_hamiltonian, N_P)
non_linear_solver = gqcpy.NonLinearEquationSolver.Newton()
linear_solver = gqcpy.LinearEquationSolver.ColPivHouseholderQR()
qc_structure = gqcpy.vAP1roG(sq_hamiltonian, N_P).optimize(non_linear_solver, non_linear_environment, linear_solver)
print(qc_structure.groundStateEnergy())
-53.25026301675117
print(qc_structure.groundStateParameters().geminalCoefficients().asMatrix())
[[ 1.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 -7.790e-04 -7.788e-04 -7.793e-04 -1.948e-03]
[ 0.000e+00 1.000e+00 0.000e+00 0.000e+00 0.000e+00 -1.416e-02 -1.413e-02 -1.421e-02 -2.450e-02]
[ 0.000e+00 0.000e+00 1.000e+00 0.000e+00 0.000e+00 -1.036e-02 -5.941e-02 -2.103e-02 -1.861e-02]
[ 0.000e+00 0.000e+00 0.000e+00 1.000e+00 0.000e+00 -4.347e-02 -1.461e-02 -3.162e-02 -1.871e-02]
[ 0.000e+00 0.000e+00 0.000e+00 0.000e+00 1.000e+00 -3.626e-02 -1.663e-02 -3.670e-02 -1.872e-02]]
print(qc_structure.groundStateParameters().lagrangeMultipliers())
[[-0.001 -0.001 -0.001 -0.002]
[-0.014 -0.014 -0.014 -0.025]
[-0.01 -0.059 -0.021 -0.019]
[-0.043 -0.015 -0.032 -0.019]
[-0.036 -0.017 -0.037 -0.019]]
Orbital optimization for vAP1roG#
We’ve also implemented an second-order orbital optimizer that uses a Newton step in every iteration. (The API isn’t quite up-to-par yet.)
The orbital optimize requires an initial guess.
G_initial = qc_structure.groundStateParameters().geminalCoefficients()
optimizer = gqcpy.AP1roGLagrangianNewtonOrbitalOptimizer(G_initial, oo_convergence_threshold=1.0e-04)
optimizer.optimize(spin_orbital_basis, sq_hamiltonian)
We can see that the electronic energy has lowered due to the orbital optimization.
print(optimizer.electronicEnergy())
-53.28557071355452
The converged geminal coefficients and multipliers can be found, too.
print(optimizer.geminalCoefficients().asMatrix())
[[ 1. 0. 0. 0. 0. -0.001 -0.001 -0.001 -0.001]
[ 0. 1. 0. 0. 0. -0.007 -0.007 -0.007 -0.091]
[ 0. 0. 1. 0. 0. -0.007 -0.091 -0.007 -0.007]
[ 0. 0. 0. 1. 0. -0.091 -0.007 -0.007 -0.007]
[ 0. 0. 0. 0. 1. -0.007 -0.007 -0.091 -0.007]]
print(optimizer.multipliers())
[[-0.001 -0.001 -0.001 -0.001]
[-0.007 -0.007 -0.007 -0.091]
[-0.007 -0.091 -0.007 -0.007]
[-0.091 -0.007 -0.007 -0.007]
[-0.007 -0.007 -0.091 -0.007]]