Following internal instabilities in GHF

Contents

Following internal instabilities in GHF#

Import gqcpy and numpy#

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

import gqcpy
import numpy as np
from scipy import linalg as la

Define a function to perform a GHF calculation#

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

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

    if guess == None:
        environment = gqcpy.GHFSCFEnvironment_d.WithCoreGuess(N, gsq_hamiltonian, S) 
    else:
        environment = gqcpy.GHFSCFEnvironment_d(N, gsq_hamiltonian, S, guess)

    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

Create a molecule and run the initial GHF calculation#

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

GHF_model = GHF_calculation(H3).parameters()

-1.3358471594598837

Run a stability check#

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

generalized_stability_matrices.printStabilityDescription()

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

Follow the internal instability#

def calculateNewGuess(GHF_stability_matrices, GHF_QC_Model):

    # The real GHF stability matrix consists of the sum of the A and B submatrices.
    Stability_matrix = GHF_stability_matrices.subMatrixA() + GHF_stability_matrices.subMatrixB()

    # diagonalize the stability matrix
    eigenvalues, eigenvectors = np.linalg.eigh(Stability_matrix)

    # Extract the lowest eigenvector
    lowest_eigenvector = eigenvectors[:, 0]

    # Create the rotation matrix that rotates the coefficients to the lowest eigenvector
    occupied_orbitals =  GHF_QC_Model.numberOfElectrons()
    virtual_orbitals = int(len(GHF_QC_Model.expansion().matrix()[0]) - occupied_orbitals)
    
    K = np.zeros(((occupied_orbitals + virtual_orbitals), (occupied_orbitals + virtual_orbitals)))
 
    lowest_eigenvector = lowest_eigenvector.reshape((occupied_orbitals, virtual_orbitals))

    K[occupied_orbitals:, :occupied_orbitals] = -1 * lowest_eigenvector.T.conjugate()
    K[:occupied_orbitals, occupied_orbitals:] = lowest_eigenvector

    rotated_coefficients = GHF_QC_Model.expansion().matrix() @ la.expm(-K)

    return gqcpy.GTransformation_d(rotated_coefficients)

Set up a new guess by rotating the old one to the lowest eigenvector of the Hessian#

new_guess = calculateNewGuess(generalized_stability_matrices, GHF_model)

Run a calculation with the new guess and check stability#

GHF_model_rotated = GHF_calculation(H3, guess=new_guess).parameters()

-1.3403026284286064

generalized_basis_rotated = gqcpy.GSpinorBasis_d(H3, 'STO-3G')
fq_hamiltonian_rotated = gqcpy.FQMolecularHamiltonian(H3)
gsq_hamiltonian_mo_rotated = generalized_basis.quantize(fq_hamiltonian_rotated).transformed(GHF_model_rotated.expansion())
generalized_stability_matrices_rotated = GHF_model_rotated.calculateStabilityMatrices(gsq_hamiltonian_mo_rotated)

generalized_stability_matrices_rotated.printStabilityDescription()

The real valued GHF wavefunction is internally stable.
The real valued GHF wavefunction is externally stable.