Custom algorithms showcase - CUHF#
# 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=6, linewidth=120)
So, recently, there has been a lot of research in our group concerned with CUHF: ‘constrained’ UHF. This is trying to modify UHF SCF with a constraint such that it should convergence to ROHF.
Since the nature of the proposed algorithms is inherently related to modifying an existing UHF SCF algorithm (plain, DIIS, …) but with a small extra step or two, I think this is an excellent showcase for GQCP’s ‘injection’ capabilities`.
Molecular setup#
Let us start by defining the molecular system under consideration. We’ll take an H3-ring, 1 a.u. apart and in the STO-3G basisset.
molecule = gqcpy.Molecule.HRingFromDistance(3, 1.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 (0.57735, 0, 0)
H (-0.288675, 0.5, 0)
H (-0.288675, -0.5, 0)
Number of alpha electrons: 2
Number of beta electrons: 1
In order to start any calculations, we need to calculate the molecular integrals. In GQCP, this is done by creating a restricted spin-orbital basis:
spinor_basis = gqcpy.USpinOrbitalBasis_d(molecule, "STO-3G") # in AO basis; atomic spin-orbitals
S = spinor_basis.overlap() # in AO basis
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:
sq_hamiltonian = spinor_basis.quantize(gqcpy.FQMolecularHamiltonian(molecule)) # in AO basis
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 electronic energy: {:.8f}\n".format(uhf_energy))
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 electronic energy: -3.63114631
UHF C_alpha:
[[-0.331961 0.002746 -1.815423]
[-0.371811 -1.569163 0.89747 ]
[-0.371698 1.56648 0.902192]]
UHF C_beta:
[[ 0.471022 -1.784404 0.002791]
[ 0.300995 0.923314 -1.56935 ]
[ 0.301463 0.928342 1.566292]]
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 # in AO basis
D_beta = environment.density_matrices[-1].beta # 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.matrix() @ np.linalg.inv(X.T)
natural_occupation_numbers, V = np.linalg.eigh(P_MO)
print(natural_occupation_numbers)
# 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_beta_constrained = F_bb - Lambda_AO
Fock_matrices = gqcpy.ScalarUSQOneElectronOperator_d(gqcpy.ScalarUSQOneElectronOperatorComponent_d(F_alpha_constrained), gqcpy.ScalarUSQOneElectronOperatorComponent_d(F_beta_constrained))
environment.replace_current_fock_matrix(Fock_matrices)
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-04, 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-04.
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()
[4.439377e-17 5.000000e-01 1.000000e+00]
[1.715894e-07 5.000000e-01 9.999998e-01]
[8.026167e-05 5.000000e-01 9.999197e-01]
[9.968217e-05 5.000000e-01 9.999003e-01]
[9.398111e-05 5.000000e-01 9.999060e-01]
[8.313397e-05 5.000000e-01 9.999169e-01]
[7.068445e-05 5.000000e-01 9.999293e-01]
[5.783872e-05 5.000000e-01 9.999422e-01]
[4.553746e-05 5.000000e-01 9.999545e-01]
[3.455296e-05 5.000000e-01 9.999654e-01]
[2.536143e-05 5.000000e-01 9.999746e-01]
[1.809853e-05 5.000000e-01 9.999819e-01]
[1.262744e-05 5.000000e-01 9.999874e-01]
[8.659465e-06 5.000000e-01 9.999913e-01]
[5.863353e-06 5.000000e-01 9.999941e-01]
[3.934205e-06 5.000000e-01 9.999961e-01]
[2.623099e-06 5.000000e-01 9.999974e-01]
[1.741339e-06 5.000000e-01 9.999983e-01]
[1.152584e-06 5.000000e-01 9.999988e-01]
[7.613817e-07 5.000000e-01 9.999992e-01]
[5.022954e-07 5.000000e-01 9.999995e-01]
[3.310818e-07 5.000000e-01 9.999997e-01]
[2.181019e-07 5.000000e-01 9.999998e-01]
[1.436207e-07 5.000000e-01 9.999999e-01]
[9.455066e-08 5.000000e-01 9.999999e-01]
[6.223572e-08 5.000000e-01 9.999999e-01]
[4.096068e-08 5.000000e-01 1.000000e+00]
[2.695648e-08 5.000000e-01 1.000000e+00]
[1.773938e-08 5.000000e-01 1.000000e+00]
[1.167347e-08 5.000000e-01 1.000000e+00]
[7.681621e-09 5.000000e-01 1.000000e+00]
[5.054749e-09 5.000000e-01 1.000000e+00]
[3.326156e-09 5.000000e-01 1.000000e+00]
[2.188684e-09 5.000000e-01 1.000000e+00]
[1.440196e-09 5.000000e-01 1.000000e+00]
[9.476747e-10 5.000000e-01 1.000000e+00]
[6.235856e-10 5.000000e-01 1.000000e+00]
[4.103292e-10 5.000000e-01 1.000000e+00]
[2.700028e-10 5.000000e-01 1.000000e+00]
[1.77666e-10 5.00000e-01 1.00000e+00]
[1.169069e-10 5.000000e-01 1.000000e+00]
[7.692639e-11 5.000000e-01 1.000000e+00]
[5.061874e-11 5.000000e-01 1.000000e+00]
[3.330782e-11 5.000000e-01 1.000000e+00]
[2.191711e-11 5.000000e-01 1.000000e+00]
[1.442164e-11 5.000000e-01 1.000000e+00]
[9.489689e-12 5.000000e-01 1.000000e+00]
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): -3.63052195
CUHF C_alpha (thesis):
[[-3.470008e-01 9.096012e-01 -1.567860e+00]
[-3.815161e-01 -1.805661e+00 7.964018e-05]
[-3.469808e-01 9.097352e-01 1.567786e+00]]
CUHF C_beta (thesis):
[[ 4.064815e-01 -1.567828e+00 8.846784e-01]
[ 2.610739e-01 8.720879e-06 -1.826967e+00]
[ 4.064861e-01 1.567818e+00 8.846941e-01]]
The CUHF (Bultinck) algorithm#
A different CUHF algorithm was formulated by P. Bultinck. We’ll use the same injection approach.
Again, we’ve adapted the algorithm in the ghf.
def constrain_function_bultinck(environment):
"""Replace the alpha- and beta- UHF coefficient matrices by their constrained counterparts, as explained in the CUHF algorithm devised by Patrick Bultinck."""
# 'Extract' the required ingredients from the environment.
N_beta = environment.N.beta
C_alpha = environment.coefficient_matrices[-1].alpha
C_a = C_alpha.matrix()
C_beta = environment.coefficient_matrices[-1].beta
C_b = C_beta.matrix()
# Apply the 'constraint'.
C_b[:, :N_beta] = C_a[:, :N_beta] # make sure the occupied beta spin-orbitals are the same as the lowest alpha spin-orbitals
new_coefficients = gqcpy.UTransformation_d(gqcpy.UTransformationComponent_d(C_a), gqcpy.UTransformationComponent_d(C_b))
# Overwrite the current coefficient matrices with their constraint versions.
# We haven't modified C_alpha, so we only have to replace the current beta coefficient matrix.
environment.replace_current_coefficient_matrix(new_coefficients)
constrain_step_bultinck = gqcpy.FunctionalStep_UHFSCFEnvironment_d(constrain_function_bultinck, "Replace the alpha- and beta- UHF coefficient matrices by their constrained counterparts, as explained in the CUHF algorithm devised by Patrick Bultinck.")
We now modify the plain UHF SCF algorithm with this constrain step, which we insert right before the next iteration starts, so at the end of the current algorithm.
CUHF_algorithm_bultinck = gqcpy.UHFSCFSolver_d.Plain(threshold = 1.0e-04, maximum_number_of_iterations = 2)
CUHF_algorithm_bultinck.insert(constrain_step_bultinck, 4) # at the end
print(CUHF_algorithm_bultinck.description())
An iterative algorithm (with a maximum of 2 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. 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.
5. Replace the alpha- and beta- UHF coefficient matrices by their constrained counterparts, as explained in the CUHF algorithm devised by Patrick Bultinck.
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.
We then proceed as usual, by creating an environment and optimizing the CUHF model parameters.
environment = gqcpy.UHFSCFEnvironment_d.WithCoreGuess(N_alpha, N_beta, sq_hamiltonian, S)
cuhf_qc_structure_bultinck = gqcpy.UHF_d.optimize(CUHF_algorithm_bultinck, environment)
cuhf_energy_bultinck = cuhf_qc_structure_bultinck.groundStateEnergy()
cuhf_parameters_bultinck = cuhf_qc_structure_bultinck.groundStateParameters()
print("CUHF electronic energy (P. Bultinck): {:.8f}\n".format(cuhf_energy_bultinck))
print("CUHF C_alpha (P. Bultinck):\n{}\n".format(cuhf_parameters_bultinck.expansion().alpha.matrix()))
print("CUHF C_beta (P. Bultinck):\n{}\n".format(cuhf_parameters_bultinck.expansion().beta.matrix()))
CUHF electronic energy (P. Bultinck): -3.62680102
CUHF C_alpha (P. Bultinck):
[[-0.358528 0.354135 -1.775391]
[-0.358528 -1.714601 0.581006]
[-0.358528 1.360467 1.194385]]
CUHF C_beta (P. Bultinck):
[[-0.358528 -1.758272 0.365065]
[-0.358528 0.573311 -1.72797 ]
[-0.358528 1.221002 1.339733]]
Bonus: changing the thesis algorithm#
Something seemed peculiar to me about the transformations of the density matrices and operators in the ‘thesis’-inspired CUHF modification. I proceeded by ‘fixing’ these transformation formulas in the following modification.
def constrain_function_thesis_fixed(environment):
"""Replace the alpha- and beta- UHF Fock matrices by their constrained counterparts, as explained in the CUHF 'thesis' algorithm."""
# Read the relevant intermediates from the environment.
F = environment.fock_matrices[-1]
D = environment.density_matrices[-1]
N = environment.N # Number of alpha and beta electrons.
# Form the closed-shell Fock matrix and the UHF modification in the AO basis.
F_cs = (F.alpha + F.beta) / 2.0
Delta_UHF = (F.alpha - F.beta) / 2.0
# Form the 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. For the orthonormal basis, we use the Löwdin basis of the alpha basis functions.
X = spinor_basis.lowdinOrthonormalization().alpha # X transforms from the AO basis to the Löwdin basis.
P_AO = (D.alpha + D.beta) / 2.0
P_MO = P_AO.transformed(X)
natural_occupation_numbers, V = np.linalg.eigh(P_MO.matrix()) # Use NumPy for the diagonalization.
# Construct the Langrange multipliers to add them to the 'constrained' Fock matrices. They
# should be constructed in the basis of the natural occupations.
V = gqcpy.UTransformationComponent_d(V) # V transforms from the Löwdin basis to the natural occupations.
Delta_UHF_NO = Delta_UHF.transformed(X).transformed(V) # In the natural occupation (NO) basis.
Delta_UHF_NO = Delta_UHF_NO.parameters()
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_NO = gqcpy.ScalarUSQOneElectronOperatorComponent_d(Lambda_NO)
Lambda_AO = Lambda_NO.transformed(V.inverse()).transformed(X.inverse()) # In the AO basis.
# Overwrite the most recent UHF Fock matrices with the CUHF modifications.
F_alpha_constrained = F_aa + Lambda_AO
F_beta_constrained = F_bb - Lambda_AO
new_fock_matrices = gqcpy.ScalarUSQOneElectronOperator_d(F_alpha_constrained, F_beta_constrained)
environment.replace_current_fock_matrix(new_fock_matrices)
constrain_step_thesis_fixed = gqcpy.FunctionalStep_UHFSCFEnvironment_d(constrain_function_thesis_fixed, description="Replace the alpha- and beta- UHF Fock matrices by their constrained counterparts, as explained in the CUHF 'thesis' algorithm, but with 'fixed' transformation formulas.")
We proceed similarly to before, inserting the modifying step in the plain UHF SCF algorithm.
CUHF_algorithm_thesis_fixed = gqcpy.UHFSCFSolver_d.Plain(threshold = 1.0e-04, maximum_number_of_iterations = 1000)
CUHF_algorithm_thesis_fixed.insert(constrain_step_thesis_fixed, 2)
print(CUHF_algorithm_thesis_fixed.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, but with 'fixed' transformation formulas.
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-04.
environment = gqcpy.UHFSCFEnvironment_d.WithCoreGuess(N_alpha, N_beta, sq_hamiltonian, S)
cuhf_qc_structure_thesis_fixed = gqcpy.UHF_d.optimize(CUHF_algorithm_thesis_fixed, environment)
cuhf_energy_thesis_fixed = cuhf_qc_structure_thesis_fixed.groundStateEnergy()
cuhf_parameters_thesis_fixed = cuhf_qc_structure_thesis_fixed.groundStateParameters()
print(environment.density_matrices[-1].alpha.matrix())
[[ 0.946746 -1.508967 0.947898]
[-1.508967 3.405967 -1.511373]
[ 0.947898 -1.511373 0.949052]]
print("CUHF electronic energy (thesis): {:.8f}\n".format(cuhf_energy_thesis_fixed))
print("CUHF C_alpha (thesis):\n{}\n".format(cuhf_parameters_thesis_fixed.expansion().alpha.matrix()))
print("CUHF C_beta (thesis):\n{}\n".format(cuhf_parameters_thesis_fixed.expansion().beta.matrix()))
CUHF electronic energy (thesis): -3.63052195
CUHF C_alpha (thesis):
[[-3.470438e-01 9.090350e-01 -1.568179e+00]
[-3.815161e-01 -1.805661e+00 7.420704e-04]
[-3.469378e-01 9.103012e-01 1.567467e+00]]
CUHF C_beta (thesis):
[[ 4.064937e-01 -1.568129e+00 8.841394e-01]
[ 2.610740e-01 6.388414e-04 -1.826967e+00]
[ 4.064738e-01 1.567517e+00 8.852330e-01]]