Checking the Hartree-Fock Stability conditions

Checking the Hartree-Fock Stability conditions#

Force the local gqcpy to be imported#

import sys
sys.path.insert(0, '../../build/gqcpy/')

import gqcpy
import numpy as np

The RHF stability conditions#

We will run an RHF calculation on $H_4$.

def RHF_calculation(molecule, basis_set='sto-3g'):
    N = molecule.numberOfElectrons()
    basis = gqcpy.RSpinOrbitalBasis_d(molecule, basis_set)
    S = basis.quantize(gqcpy.OverlapOperator())
    
    hamiltonian = gqcpy.FQMolecularHamiltonian(molecule)
    rsq_hamiltonian = basis.quantize(hamiltonian)
    objective = gqcpy.DiagonalRHFFockMatrixObjective_d(rsq_hamiltonian, 1.0e-5)  
    
    environment = gqcpy.RHFSCFEnvironment_d.WithCoreGuess(N, rsq_hamiltonian, S) 
    solver = gqcpy.RHFSCFSolver_d.Plain(threshold=1.0e-06, maximum_number_of_iterations=1000)

    qc_structure = gqcpy.RHF_d.optimize(objective, solver, environment)
    nuc_rep = gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value()
    print(qc_structure.groundStateEnergy() + nuc_rep)
    
    return qc_structure

H4 = gqcpy.Molecule.HRingFromDistance(4, 1.889, 0)

RHF_model = RHF_calculation(H4).parameters()

-1.7609380236331802

These parameters contain everything there is to know about the RHF wave function model. We can now ask this model to calculate the stability matrices. For that we need a Hamiltonian in the RHF MO basis.

restricted_basis = gqcpy.RSpinOrbitalBasis_d(H4, 'sto-3g')
fq_hamiltonian = gqcpy.FQMolecularHamiltonian(H4)
restricted_hamiltonian = restricted_basis.quantize(fq_hamiltonian)
rsq_hamiltonian_mo = restricted_hamiltonian.transformed(RHF_model.expansion())
restricted_stability_matrices = RHF_model.calculateStabilityMatrices(rsq_hamiltonian_mo)

Now we can check the stabilities. All the stability checks are done by the stability matrices themselves.

restricted_stability_matrices.isInternallyStable()

True

restricted_stability_matrices.isExternallyStable()

False

Since the wave function model is not externally stable, we can verify which external instability it contains.

restricted_stability_matrices.isTripletStable()

False

restricted_stability_matrices.isComplexConjugateStable()

False

We can also print the stability description. Note that this runs the calculation of diagonalizing the stability matrix.

restricted_stability_matrices.printStabilityDescription()

The real valued RHF wavefunction is internally stable.
The real valued RHF wavefunction contains a real->complex instability.
The real valued RHF wavefunction contains a restricted->unrestricted instability.

Since there is a real->complex instability present in the wavefunction, we can find a lower lying, complex stable state, by re-running the RHF calculation using complex valued parameters.

def complex_RHF_calculation(molecule, basis_set='sto-3g'):
    N = molecule.numberOfElectrons()
    basis = gqcpy.RSpinOrbitalBasis_cd(molecule, basis_set)
    S = basis.quantize(gqcpy.OverlapOperator())
    
    hamiltonian = gqcpy.FQMolecularHamiltonian(molecule)
    rsq_hamiltonian = basis.quantize(hamiltonian)
    objective = gqcpy.DiagonalRHFFockMatrixObjective_cd(rsq_hamiltonian, 1.0e-5)  
    
    environment = gqcpy.RHFSCFEnvironment_cd.WithComplexlyTransformedCoreGuess(N, rsq_hamiltonian, S) 
    solver = gqcpy.RHFSCFSolver_cd.Plain(threshold=1.0e-06, maximum_number_of_iterations=1000)

    qc_structure = gqcpy.RHF_cd.optimize(objective, solver, environment)
    nuc_rep = gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value()
    print(qc_structure.groundStateEnergy().real + nuc_rep)
    
    return qc_structure

complex_RHF_model = complex_RHF_calculation(H4).parameters()

-1.7720378398040086

complex_restricted_basis = gqcpy.RSpinOrbitalBasis_cd(H4, 'sto-3g')
fq_hamiltonian = gqcpy.FQMolecularHamiltonian(H4)
complex_restricted_hamiltonian = complex_restricted_basis.quantize(fq_hamiltonian)
complex_rsq_hamiltonian_mo = complex_restricted_hamiltonian.transformed(complex_RHF_model.expansion())
complex_restricted_stability_matrices = complex_RHF_model.calculateStabilityMatrices(complex_rsq_hamiltonian_mo)

complex_restricted_stability_matrices.printStabilityDescription()

The complex valued RHF wavefunction is internally stable.
The complex valued RHF wavefunction contains a restricted->unrestricted instability.

The UHF stability conditions#

We will run an UHF calculation on $H_3$.

def UHF_calculation(molecule, basis_set='STO-3G'):
    N_a = molecule.numberOfElectronPairs()
    N_b = molecule.numberOfElectrons() - N_a
    basis = gqcpy.USpinOrbitalBasis_d(molecule, basis_set)
    S = basis.quantize(gqcpy.OverlapOperator())

    hamiltonian = gqcpy.FQMolecularHamiltonian(molecule)
    sq_hamiltonian = basis.quantize(hamiltonian) 

    environment = gqcpy.UHFSCFEnvironment_d.WithCoreGuess(N_a, N_b, sq_hamiltonian, S) 
    solver = gqcpy.UHFSCFSolver_d.Plain(1.0e-06, 3000)

    qc_structure = gqcpy.UHF_d.optimize(solver, environment)
    
    return qc_structure

H3 = gqcpy.Molecule.HRingFromDistance(3, 1.889, 0)

UHF_model = UHF_calculation(H3).parameters()

These parameters contain everything there is to know about the UHF wave function model. We can now ask this model to calculate the stability matrices. For that we need a Hamiltonian in the UHF MO basis.

unrestricted_basis = gqcpy.USpinOrbitalBasis_d(H3, 'STO-3G')
usq_hamiltonian_mo = unrestricted_basis.quantize(fq_hamiltonian)
unrestricted_stability_matrices = UHF_model.calculateStabilityMatrices(usq_hamiltonian_mo)

Now we can check the stabilities.

unrestricted_stability_matrices.isInternallyStable()

False

unrestricted_stability_matrices.isExternallyStable()

False

unrestricted_stability_matrices.isSpinUnconservedStable()

False

unrestricted_stability_matrices.isComplexConjugateStable()

True

unrestricted_stability_matrices.printStabilityDescription()

The real valued UHF wavefunction contains an internal instability.
The real valued UHF wavefunction is stable within the real/complex UHF space.
The real valued UHF wavefunction contains an unrestricted->generalized instability.

The GHF stability conditions#

We will run an GHF calculation on $H_3$.

def real_GHF_calculation(molecule, basis_set='STO-3G'):
    N = H3.numberOfElectrons()
    basis = gqcpy.GSpinorBasis_d(H3, basis_set)
    S = basis.quantize(gqcpy.OverlapOperator())

    hamiltonian = gqcpy.FQMolecularHamiltonian(molecule)
    gsq_hamiltonian = basis.quantize(hamiltonian)

    environment = gqcpy.GHFSCFEnvironment_d.WithCoreGuess(N, gsq_hamiltonian, S) 
    solver = gqcpy.GHFSCFSolver_d.Plain(1.0e-08, 4000)
    qc_structure = gqcpy.GHF_d.optimize(solver, environment)
    nuc_rep = gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value()
    print(qc_structure.groundStateEnergy() + nuc_rep)
    
    return qc_structure

GHF_model = real_GHF_calculation(H3).parameters()

-1.3358471594598837

These parameters contain everything there is to know about the GHF wave function model. We can now ask this model to calculate the stability matrices. For that we need a Hamiltonian in the GHF MO basis.

generalized_basis = gqcpy.GSpinorBasis_d(H3, 'STO-3G')
gsq_hamiltonian_mo = generalized_basis.quantize(fq_hamiltonian)
generalized_stability_matrices = GHF_model.calculateStabilityMatrices(gsq_hamiltonian_mo)

generalized_stability_matrices.isInternallyStable()

False

generalized_stability_matrices.isExternallyStable()

False

generalized_stability_matrices.printStabilityDescription()

The real valued GHF wavefunction contains an internal instability.
The real valued GHF wavefunction contains a real->complex external instability.

Since we noticed that there’s a real->complex external instability, let’s try to find a complex GHF solution that is lower in energy.

def complex_GHF_calculation(molecule, basis_set='STO-3G'):
    N = H3.numberOfElectrons()
    basis = gqcpy.GSpinorBasis_cd(H3, basis_set)
    S = basis.quantize(gqcpy.OverlapOperator())

    hamiltonian = gqcpy.FQMolecularHamiltonian(molecule)
    gsq_hamiltonian = basis.quantize(hamiltonian)

    environment = gqcpy.GHFSCFEnvironment_cd.WithComplexlyTransformedCoreGuess(N, gsq_hamiltonian, S) 

    solver = gqcpy.GHFSCFSolver_cd.Plain(threshold=1.0e-08, maximum_number_of_iterations=4000)
    qc_structure = gqcpy.GHF_cd.optimize(solver, environment)
    nuc_rep = gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value()
    print(qc_structure.groundStateEnergy().real + nuc_rep)

    return qc_structure

complex_GHF_model = complex_GHF_calculation(H3).parameters()

-1.340302628428606

generalized_basis_c = gqcpy.GSpinorBasis_cd(H3, "STO-3G")
gsq_hamiltonian_mo = generalized_basis_c.quantize(fq_hamiltonian).transformed(complex_GHF_model.expansion())
generalized_stability_matrices = complex_GHF_model.calculateStabilityMatrices(gsq_hamiltonian_mo)

generalized_stability_matrices.printStabilityDescription()

The complex valued GHF wavefunction is internally stable.