Molecular setup

Molecular setup#

Let us start by defining the molecular system under consideration. We’ll take an H3-ring, 1 angstrom apart and in the STO-3G basisset.

molecule = gqcpy.Molecule.HRingFromDistance(3, 1.889, 0)
N = molecule.numberOfElectrons()
N_P = molecule.numberOfElectronPairs()  # rounds down
N_alpha = N_P + 1
N_beta = N_P

print(molecule)
print("Number of alpha electrons: {}".format(N_alpha))
print("Number of beta electrons: {}".format(N_beta))

Number of electrons: 3 
H  (1.09061, 0, 0)
H  (-0.545307, 0.9445, 0)
H  (-0.545307, -0.9445, 0)

Number of alpha electrons: 2
Number of beta electrons: 1

In order to start any unrestricted calculations, we need to calculate the molecular integrals. In GQCP, this is done by creating an unrestricted spin-orbital basis:

spinor_basis = gqcpy.USpinOrbitalBasis_d(molecule, "STO-3G")  # in AO basis; atomic spin-orbitals
K = spinor_basis.numberOfSpinors()

S = spinor_basis.quantize(gqcpy.OverlapOperator())

Then, we will quantize the molecular Hamiltonian (in which the interactions terms are the electron’s kinetic energy, the nucleus-electron attraction and the interelectronic repulson) in that spin-orbital basis:

fq_hamiltonian = gqcpy.FQMolecularHamiltonian(molecule)
sq_hamiltonian = spinor_basis.quantize(fq_hamiltonian)

The standard ‘plain’ UHF SCF algorithm#

Let’s establish a base line for the algorithms by examining the standard ‘plain’ UHF SCF algorithm. In GQCP, we have implemented algorithms and environments:

an environment acts as a ‘calculation space’, storing intermediate calculations;
an algorithm and its composing steps continuously modify the environment.

For UHF SCF, we use an UHFSCFEnvironment:

environment = gqcpy.UHFSCFEnvironment_d.WithCoreGuess(N_alpha, N_beta, sq_hamiltonian, S)

where we have initialized the alpha- and beta-coefficient matrices with the (generalized) eigenvectors of the core Hamiltonian (in AO basis).

Let’s now take a look at the plain UHF SCF solver:

plain_algorithm = gqcpy.UHFSCFSolver_d.Plain(threshold=1.0e-04, maximum_number_of_iterations=1000)

print(plain_algorithm.description())

An iterative algorithm (with a maximum of 1000 iterations) consisting of the following steps:
An algorithmic step consisting of 4 algorithmic steps:
	1. Calculate the current UHF alpha and beta density matrices (in the AO basis) and place them in the environment.
	2. Calculate the current UHF Fock matrices (expressed in the scalar/AO basis) and place them in the environment.
	3. Solve the generalized eigenvalue problem for the most recent scalar/AO Fock matrices. Add the associated coefficient matrices and orbital energies to the environment.
	4. Calculate the current electronic UHF energy and place it in the environment.

With the following convergence criterion:
A convergence criterion that checks if the norm of the difference of two iterates (the UHF spin resolved density matrix in AO basis) is converged, with a tolerance of 1.00e-04.

With this algorithm, let’s check the UHF SCF procedure.

uhf_qc_structure = gqcpy.UHF_d.optimize(plain_algorithm, environment)

uhf_energy = uhf_qc_structure.groundStateEnergy()
uhf_parameters = uhf_qc_structure.groundStateParameters()

print("UHF total energy: {:.8f}\n".format(uhf_energy + gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value()))
print("UHF C_alpha:\n{}\n".format(uhf_parameters.expansion().alpha.matrix()))
print("UHF C_beta:\n{}\n".format(uhf_parameters.expansion().beta.matrix()))

UHF total energy: -1.33584716

UHF C_alpha:
[[-4.599757e-01 -9.968298e-01 -5.353406e-01]
 [-4.599315e-01  9.966246e-01 -5.357604e-01]
 [-3.016910e-01  2.559108e-04  1.183562e+00]]

UHF C_beta:
[[ 2.806758e-01 -6.474947e-01  9.968995e-01]
 [ 2.807662e-01 -6.479857e-01 -9.965549e-01]
 [ 6.386637e-01  1.041127e+00 -2.275071e-04]]

The CUHF (Scuseria/thesis) algorithm#

One of the CUHF algorithms described by Scuseria’s group can be found in the PhD thesis by T. Tsuchimochi. Since this algorithm uses a modification of the UHF alpha- and beta-Fock matrices, this is an ideal opportunity to use GQCP’s injection capabilities.

I’ve taken the reference implementation from Xeno’s ghf repository, and adapted it to GQCP.

We start off by defining a Python function that modifies the alpha- and beta- Fock matrices in an UHFSCFEnvironment. We then wrap that Python function into a special class FunctionalStep_UHFSCFEnvironment, which allows for the bridging between Python functions and GQCP’s Algorithm framework.

def constrain_function_thesis(environment):
    """Replace the alpha- and beta- UHF Fock matrices by their constrained counterparts, as explained in the CUHF 'thesis' algorithm."""

    # Prepare some variables to be used in the algorithm.
    X = spinor_basis.lowdinOrthonormalization().alpha.matrix()  # the transformation matrix TO an orthonormal basis
    
    F_alpha = environment.fock_matrices[-1].alpha.parameters()  # in AO basis
    F_beta = environment.fock_matrices[-1].beta.parameters()  # in AO basis
    
    D_alpha = environment.density_matrices[-1].alpha.matrix()  # in AO basis
    D_beta = environment.density_matrices[-1].beta.matrix()  # in AO basis
    
    N_alpha = environment.N.alpha
    N_beta = environment.N.beta
    
    # Form the closed-shell Fock matrix and the UHF modification.
    F_cs = (F_alpha + F_beta) / 2.0  # the closed-shell Fock matrix, in AO basis
    Delta_UHF = (F_alpha - F_beta) / 2.0  # in AO basis

    # Form modified Fock matrices.
    F_aa = F_cs + Delta_UHF
    F_bb = F_cs - Delta_UHF
    
    # Find the natural occupation numbers and vectors by diagonalizing the charge-density matrix in an orthonormal basis.
    P_AO = (D_alpha + D_beta) / 2.0
    P_MO = np.linalg.inv(X) @ P_AO @ np.linalg.inv(X.T)
    
    natural_occupation_numbers, V = np.linalg.eigh(P_MO)
    
    # Construct the Langrange multipliers to add them to the 'constrained' Fock matrices.
    Delta_UHF_NO = np.linalg.inv(V) @ np.linalg.inv(X) @ Delta_UHF @ np.linalg.inv(X.T) @ np.linalg.inv(V.T)  # in the natural occupation (NO) basis
    Lambda_NO = np.zeros(np.shape(Delta_UHF_NO))
    Lambda_NO[:N_beta, N_alpha:] = -Delta_UHF_NO[:N_beta, N_alpha:]
    Lambda_NO[N_alpha:, :N_beta] = -Delta_UHF_NO[N_alpha:, :N_beta]
    Lambda_AO = X @ V @ Lambda_NO @ V.T @ X.T
    

    # Overwrite the most recent UHF Fock matrices with the cUHF ones
    
    F_alpha_constrained = F_aa + Lambda_AO
    F_alpha_constrained = gqcpy.ScalarUSQOneElectronOperatorComponent_d(F_alpha_constrained)
    F_beta_constrained = F_bb - Lambda_AO
    F_beta_constrained = gqcpy.ScalarUSQOneElectronOperatorComponent_d(F_beta_constrained)

    fock_matrix = gqcpy.ScalarUSQOneElectronOperator_d(F_alpha_constrained, F_beta_constrained)
    environment.replace_current_fock_matrix(fock_matrix)

constrain_step_thesis = gqcpy.FunctionalStep_UHFSCFEnvironment_d(constrain_function_thesis, description="Replace the alpha- and beta- UHF Fock matrices by their constrained counterparts, as explained in the CUHF 'thesis' algorithm.")

We’ll now modify the plain UHF SCF algorithm with this modification constrain_step_thesis. From the algorithmic description printed above, we’ll have to insert the Fock matrix modification at index 2.

CUHF_algorithm_thesis = gqcpy.UHFSCFSolver_d.Plain(threshold = 1.0e-12, maximum_number_of_iterations = 1000)

CUHF_algorithm_thesis.insert(constrain_step_thesis, 2)
print(CUHF_algorithm_thesis.description())

An iterative algorithm (with a maximum of 1000 iterations) consisting of the following steps:
An algorithmic step consisting of 5 algorithmic steps:
	1. Calculate the current UHF alpha and beta density matrices (in the AO basis) and place them in the environment.
	2. Calculate the current UHF Fock matrices (expressed in the scalar/AO basis) and place them in the environment.
	3. Replace the alpha- and beta- UHF Fock matrices by their constrained counterparts, as explained in the CUHF 'thesis' algorithm.
	4. Solve the generalized eigenvalue problem for the most recent scalar/AO Fock matrices. Add the associated coefficient matrices and orbital energies to the environment.
	5. Calculate the current electronic UHF energy and place it in the environment.

With the following convergence criterion:
A convergence criterion that checks if the norm of the difference of two iterates (the UHF spin resolved density matrix in AO basis) is converged, with a tolerance of 1.00e-12.

Since we’ve now created the CUHF_algorithm_thesis, which is still an UHF SCF solver, we can set up an environment and optimize the CUHF model parameters.

environment = gqcpy.UHFSCFEnvironment_d.WithCoreGuess(N_alpha, N_beta, sq_hamiltonian, S)

cuhf_qc_structure_thesis = gqcpy.UHF_d.optimize(CUHF_algorithm_thesis, environment)
cuhf_energy_thesis = cuhf_qc_structure_thesis.groundStateEnergy()
cuhf_parameters_thesis = cuhf_qc_structure_thesis.groundStateParameters()

print("CUHF electronic energy (thesis): {:.8f}\n".format(cuhf_energy_thesis))
print("CUHF C_alpha (thesis):\n{}\n".format(cuhf_parameters_thesis.expansion().alpha.matrix()))
print("CUHF C_beta (thesis):\n{}\n".format(cuhf_parameters_thesis.expansion().beta.matrix()))

CUHF electronic energy (thesis): -2.91619904

CUHF C_alpha (thesis):
[[-4.754485e-01 -1.125071e+00 -7.419621e-13]
 [-3.746105e-01  5.983634e-01 -9.967273e-01]
 [-3.746105e-01  5.983634e-01  9.967273e-01]]

CUHF C_beta (thesis):
[[ 2.552886e-01  1.853517e-13 -1.194431e+00]
 [ 4.804962e-01  9.967273e-01  5.171993e-01]
 [ 4.804962e-01 -9.967273e-01  5.171993e-01]]

sq_hamiltonian.transform(cuhf_parameters_thesis.expansion())

Dense CIS calculations#

In order to do a dense CIS calculation, GQCP uses its own framework for mathematical optimization. We’ll have to create:

a dense eigenvalue problem solver
an associated eigenvalue problem environment

onv_basis = gqcpy.SpinResolvedSelectedONVBasis.CIS(K//2, N_alpha, N_beta, True)

solver = gqcpy.EigenproblemSolver.Dense_d()
environment = gqcpy.CIEnvironment.Dense(sq_hamiltonian, onv_basis)

The encompassing QCMethod is CI (for any type of configuration interaction), so we’ll use that to find the optimized and associated wave function model: a LinearExpansion. Internally, QCMethods return a QCStructure, which wraps energies and optimized parameters for the ground state and excited states.

qc_structure = gqcpy.CI(onv_basis, onv_basis.dimension()).optimize(solver, environment)

for i in range(onv_basis.dimension()):
    energy = qc_structure.energy(i) #+ gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value()
    if np.abs(energy - cuhf_energy_thesis) > 1e-8:
        print("excited state (electronic) energy: {}".format(energy))
    else:
        print("ground state (electronic) energy: {}".format(energy))

excited state (electronic) energy: -2.923132380774984
ground state (electronic) energy: -2.91619903663055
ground state (electronic) energy: -2.9161990366305486
excited state (electronic) energy: -2.896197641014147
excited state (electronic) energy: -2.5611054054604567
excited state (electronic) energy: -2.436718810736164
excited state (electronic) energy: -2.3232093422666957
excited state (electronic) energy: -2.176823953117537
excited state (electronic) energy: -2.1768239531175366
excited state (electronic) energy: -2.1392307239340025

Molecular setup

Contents

Molecular setup#

The standard ‘plain’ UHF SCF algorithm#

The CUHF (Scuseria/thesis) algorithm#

Dense CIS calculations#