Overlaps between R/U/GHF#
# Force the local gqcpy to be imported
import sys
sys.path.insert(0, '../../build/gqcpy/')
import gqcpy
import numpy as np
import numpy.random as rand
np.set_printoptions(precision=3, linewidth=120)
In this notebook, we will explore the multi-configurational behaviour of UHF and GHF, by projecting the UHF and GHF ONVs onto the full ONV basis expressed with respect to the RHF orbitals.
Projecting UHF#
H2#
As an introductory example, let’s start by examining H2. We’ll set up an unrestricted spin-orbital basis, where the alpha- and beta-coefficient matrices are equal. Then, we expect the projection of the UHF ONV onto the ONV basis expressed with respect to the RHF spin-orbitals to yield a linear combination consisting of exactly one 1, and the rest 0.
molecule = gqcpy.Molecule.HChain(2, 1.0)
N = molecule.numberOfElectrons()
r_spinor_basis = gqcpy.RSpinOrbitalBasis_d(molecule, "STO-3G")
K = r_spinor_basis.numberOfSpatialOrbitals()
S = r_spinor_basis.quantize(gqcpy.OverlapOperator())
sq_hamiltonian = r_spinor_basis.quantize(gqcpy.FQMolecularHamiltonian(molecule)) # 'sq' for 'second-quantized'
rhf_environment = gqcpy.RHFSCFEnvironment_d.WithCoreGuess(N, sq_hamiltonian, S)
rhf_solver = gqcpy.RHFSCFSolver_d.Plain()
rhf_objective = gqcpy.DiagonalRHFFockMatrixObjective_d(sq_hamiltonian) # use the default threshold of 1.0e-08
rhf_parameters = gqcpy.RHF_d.optimize(rhf_objective, rhf_solver, rhf_environment).groundStateParameters()
After having done the RHF calculation, we are able to construct restricted and unrestricted spin-orbital bases.
r_spinor_basis.transform(rhf_parameters.expansion())
u_spinor_basis = gqcpy.USpinOrbitalBasis_d.FromRestricted(r_spinor_basis)
We are now in the position to calculate the expansion of the UHF determinant in the full spin-resolved ONV basis expressed with respect to the RHF spin-orbitals.
uhf_determinant = gqcpy.SpinResolvedONV.UHF(2, 1, 1)
Projecting this determinant is done by constructing a LinearExpansion from an ONV projection:
expansion = gqcpy.LinearExpansion_SpinResolved.FromONVProjection(uhf_determinant, r_spinor_basis, u_spinor_basis)
def print_function(coefficient, onv):
print("{:7.4f}".format(coefficient), onv)
expansion.forEach(print_function)
1.0000 10|10
-0.0000 10|01
-0.0000 01|10
0.0000 01|01
H4#
As a more interesting example, we’ll project the H4 UHF triplet state onto the basis of ONVs expressed with respect to the RHF spin-orbitals.
molecule = gqcpy.Molecule.HRingFromDistance(4, 1.0, 0)
N = molecule.numberOfElectrons()
N_P = N//2
N_alpha = N//2
N_beta = N//2
internuclear_repulsion_energy = gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value()
Using Xeno’s code, we find the following coefficient matrices for the RHF solution:
C = [[-0.27745359, -0.8505133, 0.85051937, 2.02075317],
[-0.27745362, -0.85051937, -0.8505133, -2.02075317],
[-0.27745359, 0.8505133, -0.85051937, 2.02075317],
[-0.27745362, 0.85051937, 0.8505133, -2.02075317]]
C = np.array(C)
and the UHF solution:
C_alpha = [[-1.75646828e-01, -1.20606646e-06, 1.20281173e+00, 2.03213486e+00],
[-3.78560533e-01, -1.20281173e+00, -1.20606647e-06, -2.00427438e+00],
[-1.75646828e-01, 1.20606646e-06, -1.20281173e+00, 2.03213486e+00],
[-3.78560533e-01, 1.20281173e+00, 1.20606646e-06, -2.00427438e+00]]
C_alpha = np.array(C_alpha)
C_beta = [[-3.78560533e-01, 1.20281173e+00, 1.21724557e-06, 2.00427438e+00],
[-1.75646828e-01, 1.21724558e-06, -1.20281173e+00, -2.03213486e+00],
[-3.78560533e-01, -1.20281173e+00, -1.21724558e-06, 2.00427438e+00],
[-1.75646828e-01, -1.21724558e-06, 1.20281173e+00, -2.03213486e+00]]
C_beta = np.array(C_beta)
Let’s now set up the corresponding restricted and unrestricted spin-orbital bases.
r_spinor_basis = gqcpy.RSpinOrbitalBasis_d(molecule, "STO-3G")
r_spinor_basis.transform(gqcpy.RTransformation_d(C))
u_spinor_basis = gqcpy.USpinOrbitalBasis_d(molecule, "STO-3G")
u_spinor_basis.transform(gqcpy.UTransformation_d(gqcpy.UTransformationComponent_d(C_alpha), gqcpy.UTransformationComponent_d(C_beta)))
We can now construct the UHF spin-resolved ONV:
uhf_determinant = gqcpy.SpinResolvedONV.UHF(4, 2, 2)
and project it onto the spin-resolved basis of ONVs expressed with respect to the orthonormal restricted spin-orbitals:
expansion = gqcpy.LinearExpansion_SpinResolved.FromONVProjection(uhf_determinant, r_spinor_basis, u_spinor_basis)
expansion.forEach(print_function)
-0.4987 1100|1100
0.4987 1100|1010
0.0000 1100|0110
0.0000 1100|1001
-0.0251 1100|0101
0.0251 1100|0011
-0.4987 1010|1100
0.4987 1010|1010
0.0000 1010|0110
0.0000 1010|1001
-0.0251 1010|0101
0.0251 1010|0011
-0.0000 0110|1100
0.0000 0110|1010
0.0000 0110|0110
0.0000 0110|1001
-0.0000 0110|0101
0.0000 0110|0011
0.0000 1001|1100
-0.0000 1001|1010
-0.0000 1001|0110
-0.0000 1001|1001
0.0000 1001|0101
-0.0000 1001|0011
0.0251 0101|1100
-0.0251 0101|1010
-0.0000 0101|0110
-0.0000 0101|1001
0.0013 0101|0101
-0.0013 0101|0011
0.0251 0011|1100
-0.0251 0011|1010
-0.0000 0011|0110
-0.0000 0011|1001
0.0013 0011|0101
-0.0013 0011|0011
GHF on UHF projection#
Let us start by performing an Unrestricted Hartree-Fock and a Generalized Hartree-Fock calculation on H3. Both should yield different results as H3 posseses a fully spin-symmetry broken solution.
molecule = gqcpy.Molecule.HRingFromDistance(3, 1.889, 0) # create a neutral molecule
N = molecule.numberOfElectrons()
N_alpha = N//2 + 1
N_beta = N//2
spinor_basis = gqcpy.USpinOrbitalBasis_d(molecule, "STO-3G")
S = spinor_basis.quantize(gqcpy.OverlapOperator())
sq_hamiltonian = spinor_basis.quantize(gqcpy.FQMolecularHamiltonian(molecule)) # 'sq' for 'second-quantized'
UHF_environment = gqcpy.UHFSCFEnvironment_d.WithCoreGuess(N_alpha, N_beta, sq_hamiltonian, S)
UHF_solver = gqcpy.UHFSCFSolver_d.Plain(threshold=1.0e-04, maximum_number_of_iterations=1000) # the system is not converging very rapidly
uqc_structure = gqcpy.UHF_d.optimize(UHF_solver, UHF_environment)
print(uqc_structure.groundStateEnergy() + gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value())
uhf_parameters = uqc_structure.groundStateParameters()
-1.3358471555881106
We now have a UHF solution with the wanted parameters saved in a variable. We can now do the same for GHF.
basis = gqcpy.GSpinorBasis_d(molecule, "STO-3G")
gS = basis.quantize(gqcpy.OverlapOperator())
gsq_hamiltonian = basis.quantize(gqcpy.FQMolecularHamiltonian(molecule))
# Generate a random guess in order to find a true GHF solution.
rand.seed(2)
K = spinor_basis.numberOfSpinors()
random_matrix = np.random.rand(K, K)
random_matrix_transpose = random_matrix.T
symmetric_random_matrix = random_matrix + random_matrix_transpose
_, guess = np.linalg.eigh(symmetric_random_matrix)
GHF_environment = gqcpy.GHFSCFEnvironment_d(N, gsq_hamiltonian, gS, gqcpy.GTransformation_d(guess))
GHF_solver = gqcpy.GHFSCFSolver_d.Plain(1.0e-08, 4000)
gqc_structure = gqcpy.GHF_d.optimize(GHF_solver, GHF_environment)
print(gqc_structure.groundStateEnergy() + gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value())
ghf_parameters = gqc_structure.groundStateParameters()
-1.3403026284286064
In order to create a projection, we will need to transform the respective spino(bital) bases to MO basis, as well as transform the UHF basis to a spinor representation.
basis.transform(ghf_parameters.expansion())
spinor_basis.transform(uhf_parameters.expansion())
spinor_basis = gqcpy.GSpinorBasis_d.FromUnrestricted(spinor_basis)
The final ingredient we need is a representation of the GHF determinant.
GHF_determinant = gqcpy.SpinUnresolvedONV.GHF(K, N, ghf_parameters.orbitalEnergies())
We can now create the projection and print it.
expansion = gqcpy.LinearExpansion_SpinUnresolved.FromONVProjection(GHF_determinant, spinor_basis, basis)
expansion.forEach(print_function)
0.0334 111000
-0.2384 110100
-0.6027 101100
-0.0989 011100
0.0976 110010
0.2098 101010
0.0339 011010
0.2638 100110
0.0469 010110
0.0090 001110
0.0780 110001
0.1133 101001
0.0276 011001
0.6000 100101
0.0341 010101
-0.1627 001101
-0.1593 100011
0.0014 010011
0.0583 001011
0.0802 000111