UHF SCF#
# 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=3, linewidth=120)
Molecular setup#
UHF SCF calculations start off in the same way as RHF SCF calculations. We first need a second-quantized Hamiltonian that expresses all the integrals in the AO basis.
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'
print(molecule)
Number of electrons: 3
H (1.09061, 0, 0)
H (-0.545307, 0.9445, 0)
H (-0.545307, -0.9445, 0)
Plain UHF SCF calculations#
Then, we need an UHF SCF solver and an UHF SCF environment.
environment = gqcpy.UHFSCFEnvironment_d.WithCoreGuess(N_alpha, N_beta, sq_hamiltonian, S)
plain_solver = gqcpy.UHFSCFSolver_d.Plain(threshold=1.0e-04, maximum_number_of_iterations=1000) # the system is not converging very rapidly
Using GQCP’s capabilities of ‘injecting’ Python code inside the C++ library, we’ll add an intermediary step that doesn’t modify the environment, but just prints out the current Fock matrices.
This is done by writing a Python function that manipulates the solver’s environment (of type EnvironmentType), and subsequently wrapping that function into a FunctionalStep_EnvironmentType. The instance of type FunctionalStep_EnvironmentType is then subsequently inserted inside the algorithm that should be performed.
In this example, we’ll print out the alpha- and beta- Fock matrices. The Fock matrices are stored as a list of unrestricted Fock-operators. Hence, we have to get alpha/beta component of the last entry and print its parameters.
def print_fock_matrices(environment):
if environment.fock_matrices:
print("Alpha Fock matrix:\n{}\n".format(environment.fock_matrices[-1].alpha.parameters()))
print("Beta Fock matrix:\n{}".format(environment.fock_matrices[-1].beta.parameters()))
else:
print("No Fock matrices in the queue.")
print_step = gqcpy.FunctionalStep_UHFSCFEnvironment_d(print_fock_matrices)
plain_solver.insert(print_step, index=0) # insert at the front
We can use the solver to optimize the UHF parameters.
qc_structure = gqcpy.UHF_d.optimize(plain_solver, environment)
No Fock matrices in the queue.
Alpha Fock matrix:
[[-0.467 -0.42 -0.508]
[-0.42 -0.436 -0.539]
[-0.508 -0.539 -0.348]]
Beta Fock matrix:
[[-0.117 -0.435 -0.46 ]
[-0.435 -0.168 -0.468]
[-0.46 -0.468 -0.311]]
Alpha Fock matrix:
[[-0.488 -0.427 -0.506]
[-0.427 -0.446 -0.535]
[-0.506 -0.535 -0.316]]
Beta Fock matrix:
[[-0.09 -0.415 -0.462]
[-0.415 -0.154 -0.478]
[-0.462 -0.478 -0.346]]
Alpha Fock matrix:
[[-0.5 -0.434 -0.502]
[-0.434 -0.458 -0.531]
[-0.502 -0.531 -0.293]]
Beta Fock matrix:
[[-0.076 -0.403 -0.463]
[-0.403 -0.141 -0.48 ]
[-0.463 -0.48 -0.368]]
Alpha Fock matrix:
[[-0.506 -0.439 -0.499]
[-0.439 -0.467 -0.527]
[-0.499 -0.527 -0.277]]
Beta Fock matrix:
[[-0.068 -0.396 -0.463]
[-0.396 -0.13 -0.481]
[-0.463 -0.481 -0.383]]
Alpha Fock matrix:
[[-0.51 -0.443 -0.497]
[-0.443 -0.474 -0.524]
[-0.497 -0.524 -0.265]]
Beta Fock matrix:
[[-0.064 -0.391 -0.464]
[-0.391 -0.121 -0.481]
[-0.464 -0.481 -0.393]]
Alpha Fock matrix:
[[-0.512 -0.445 -0.497]
[-0.445 -0.48 -0.522]
[-0.497 -0.522 -0.257]]
Beta Fock matrix:
[[-0.062 -0.387 -0.465]
[-0.387 -0.114 -0.48 ]
[-0.465 -0.48 -0.4 ]]
Alpha Fock matrix:
[[-0.513 -0.446 -0.497]
[-0.446 -0.484 -0.52 ]
[-0.497 -0.52 -0.251]]
Beta Fock matrix:
[[-0.061 -0.384 -0.466]
[-0.384 -0.108 -0.479]
[-0.466 -0.479 -0.405]]
Alpha Fock matrix:
[[-0.513 -0.447 -0.497]
[-0.447 -0.488 -0.518]
[-0.497 -0.518 -0.247]]
Beta Fock matrix:
[[-0.061 -0.383 -0.467]
[-0.383 -0.103 -0.479]
[-0.467 -0.479 -0.409]]
Alpha Fock matrix:
[[-0.513 -0.448 -0.497]
[-0.448 -0.49 -0.516]
[-0.497 -0.516 -0.244]]
Beta Fock matrix:
[[-0.062 -0.381 -0.467]
[-0.381 -0.099 -0.478]
[-0.467 -0.478 -0.412]]
Alpha Fock matrix:
[[-0.513 -0.449 -0.498]
[-0.449 -0.493 -0.515]
[-0.498 -0.515 -0.242]]
Beta Fock matrix:
[[-0.062 -0.38 -0.468]
[-0.38 -0.096 -0.477]
[-0.468 -0.477 -0.414]]
Alpha Fock matrix:
[[-0.512 -0.449 -0.498]
[-0.449 -0.495 -0.514]
[-0.498 -0.514 -0.24 ]]
Beta Fock matrix:
[[-0.063 -0.379 -0.468]
[-0.379 -0.093 -0.477]
[-0.468 -0.477 -0.415]]
Alpha Fock matrix:
[[-0.512 -0.449 -0.499]
[-0.449 -0.496 -0.513]
[-0.499 -0.513 -0.239]]
Beta Fock matrix:
[[-0.064 -0.379 -0.469]
[-0.379 -0.09 -0.476]
[-0.469 -0.476 -0.416]]
Alpha Fock matrix:
[[-0.511 -0.449 -0.5 ]
[-0.449 -0.497 -0.512]
[-0.5 -0.512 -0.238]]
Beta Fock matrix:
[[-0.065 -0.378 -0.469]
[-0.378 -0.088 -0.476]
[-0.469 -0.476 -0.417]]
Alpha Fock matrix:
[[-0.511 -0.449 -0.5 ]
[-0.449 -0.499 -0.511]
[-0.5 -0.511 -0.237]]
Beta Fock matrix:
[[-0.066 -0.378 -0.47 ]
[-0.378 -0.086 -0.476]
[-0.47 -0.476 -0.418]]
Alpha Fock matrix:
[[-0.51 -0.449 -0.501]
[-0.449 -0.5 -0.51 ]
[-0.501 -0.51 -0.237]]
Beta Fock matrix:
[[-0.067 -0.378 -0.47 ]
[-0.378 -0.085 -0.475]
[-0.47 -0.475 -0.418]]
Alpha Fock matrix:
[[-0.51 -0.45 -0.501]
[-0.45 -0.5 -0.51 ]
[-0.501 -0.51 -0.236]]
Beta Fock matrix:
[[-0.068 -0.378 -0.47 ]
[-0.378 -0.084 -0.475]
[-0.47 -0.475 -0.418]]
Alpha Fock matrix:
[[-0.509 -0.45 -0.502]
[-0.45 -0.501 -0.509]
[-0.502 -0.509 -0.236]]
Beta Fock matrix:
[[-0.068 -0.378 -0.471]
[-0.378 -0.082 -0.475]
[-0.471 -0.475 -0.419]]
Alpha Fock matrix:
[[-0.509 -0.45 -0.502]
[-0.45 -0.502 -0.509]
[-0.502 -0.509 -0.236]]
Beta Fock matrix:
[[-0.069 -0.377 -0.471]
[-0.377 -0.081 -0.474]
[-0.471 -0.474 -0.419]]
Alpha Fock matrix:
[[-0.509 -0.45 -0.502]
[-0.45 -0.502 -0.508]
[-0.502 -0.508 -0.235]]
Beta Fock matrix:
[[-0.07 -0.377 -0.471]
[-0.377 -0.081 -0.474]
[-0.471 -0.474 -0.419]]
Alpha Fock matrix:
[[-0.508 -0.45 -0.503]
[-0.45 -0.503 -0.508]
[-0.503 -0.508 -0.235]]
Beta Fock matrix:
[[-0.07 -0.377 -0.471]
[-0.377 -0.08 -0.474]
[-0.471 -0.474 -0.419]]
Alpha Fock matrix:
[[-0.508 -0.45 -0.503]
[-0.45 -0.503 -0.508]
[-0.503 -0.508 -0.235]]
Beta Fock matrix:
[[-0.071 -0.377 -0.471]
[-0.377 -0.079 -0.474]
[-0.471 -0.474 -0.419]]
Alpha Fock matrix:
[[-0.508 -0.45 -0.503]
[-0.45 -0.503 -0.507]
[-0.503 -0.507 -0.235]]
Beta Fock matrix:
[[-0.071 -0.377 -0.472]
[-0.377 -0.079 -0.474]
[-0.472 -0.474 -0.419]]
Alpha Fock matrix:
[[-0.508 -0.45 -0.504]
[-0.45 -0.504 -0.507]
[-0.504 -0.507 -0.235]]
Beta Fock matrix:
[[-0.072 -0.377 -0.472]
[-0.377 -0.078 -0.474]
[-0.472 -0.474 -0.42 ]]
Alpha Fock matrix:
[[-0.507 -0.45 -0.504]
[-0.45 -0.504 -0.507]
[-0.504 -0.507 -0.235]]
Beta Fock matrix:
[[-0.072 -0.377 -0.472]
[-0.377 -0.078 -0.473]
[-0.472 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.507 -0.45 -0.504]
[-0.45 -0.504 -0.507]
[-0.504 -0.507 -0.235]]
Beta Fock matrix:
[[-0.072 -0.377 -0.472]
[-0.377 -0.077 -0.473]
[-0.472 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.507 -0.45 -0.504]
[-0.45 -0.504 -0.506]
[-0.504 -0.506 -0.235]]
Beta Fock matrix:
[[-0.073 -0.377 -0.472]
[-0.377 -0.077 -0.473]
[-0.472 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.507 -0.45 -0.504]
[-0.45 -0.504 -0.506]
[-0.504 -0.506 -0.235]]
Beta Fock matrix:
[[-0.073 -0.377 -0.472]
[-0.377 -0.077 -0.473]
[-0.472 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.507 -0.45 -0.504]
[-0.45 -0.505 -0.506]
[-0.504 -0.506 -0.235]]
Beta Fock matrix:
[[-0.073 -0.377 -0.472]
[-0.377 -0.077 -0.473]
[-0.472 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.507 -0.45 -0.504]
[-0.45 -0.505 -0.506]
[-0.504 -0.506 -0.235]]
Beta Fock matrix:
[[-0.073 -0.377 -0.472]
[-0.377 -0.076 -0.473]
[-0.472 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.505 -0.506]
[-0.505 -0.506 -0.235]]
Beta Fock matrix:
[[-0.073 -0.377 -0.472]
[-0.377 -0.076 -0.473]
[-0.472 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.505 -0.506]
[-0.505 -0.506 -0.235]]
Beta Fock matrix:
[[-0.074 -0.377 -0.472]
[-0.377 -0.076 -0.473]
[-0.472 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.505 -0.506]
[-0.505 -0.506 -0.235]]
Beta Fock matrix:
[[-0.074 -0.377 -0.472]
[-0.377 -0.076 -0.473]
[-0.472 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.505 -0.506]
[-0.505 -0.506 -0.235]]
Beta Fock matrix:
[[-0.074 -0.377 -0.472]
[-0.377 -0.076 -0.473]
[-0.472 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.505 -0.506]
[-0.505 -0.506 -0.235]]
Beta Fock matrix:
[[-0.074 -0.377 -0.472]
[-0.377 -0.076 -0.473]
[-0.472 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.505 -0.506]
[-0.505 -0.506 -0.235]]
Beta Fock matrix:
[[-0.074 -0.377 -0.472]
[-0.377 -0.076 -0.473]
[-0.472 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.505 -0.506]
[-0.505 -0.506 -0.235]]
Beta Fock matrix:
[[-0.074 -0.377 -0.472]
[-0.377 -0.075 -0.473]
[-0.472 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.505 -0.506]
[-0.505 -0.506 -0.235]]
Beta Fock matrix:
[[-0.074 -0.377 -0.472]
[-0.377 -0.075 -0.473]
[-0.472 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.505 -0.506]
[-0.505 -0.506 -0.235]]
Beta Fock matrix:
[[-0.074 -0.377 -0.472]
[-0.377 -0.075 -0.473]
[-0.472 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.505 -0.506]
[-0.505 -0.506 -0.235]]
Beta Fock matrix:
[[-0.074 -0.377 -0.472]
[-0.377 -0.075 -0.473]
[-0.472 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.505 -0.505]
[-0.505 -0.505 -0.235]]
Beta Fock matrix:
[[-0.074 -0.377 -0.473]
[-0.377 -0.075 -0.473]
[-0.473 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.505 -0.505]
[-0.505 -0.505 -0.235]]
Beta Fock matrix:
[[-0.074 -0.377 -0.473]
[-0.377 -0.075 -0.473]
[-0.473 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.506 -0.505]
[-0.505 -0.505 -0.235]]
Beta Fock matrix:
[[-0.075 -0.377 -0.473]
[-0.377 -0.075 -0.473]
[-0.473 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.506 -0.505]
[-0.505 -0.505 -0.235]]
Beta Fock matrix:
[[-0.075 -0.377 -0.473]
[-0.377 -0.075 -0.473]
[-0.473 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.506 -0.505]
[-0.505 -0.505 -0.235]]
Beta Fock matrix:
[[-0.075 -0.377 -0.473]
[-0.377 -0.075 -0.473]
[-0.473 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.506 -0.505]
[-0.505 -0.505 -0.235]]
Beta Fock matrix:
[[-0.075 -0.377 -0.473]
[-0.377 -0.075 -0.473]
[-0.473 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.506 -0.505]
[-0.505 -0.505 -0.235]]
Beta Fock matrix:
[[-0.075 -0.377 -0.473]
[-0.377 -0.075 -0.473]
[-0.473 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.506 -0.505]
[-0.505 -0.505 -0.235]]
Beta Fock matrix:
[[-0.075 -0.377 -0.473]
[-0.377 -0.075 -0.473]
[-0.473 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.506 -0.505]
[-0.505 -0.505 -0.235]]
Beta Fock matrix:
[[-0.075 -0.377 -0.473]
[-0.377 -0.075 -0.473]
[-0.473 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.506 -0.505]
[-0.505 -0.505 -0.235]]
Beta Fock matrix:
[[-0.075 -0.377 -0.473]
[-0.377 -0.075 -0.473]
[-0.473 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.506 -0.505]
[-0.505 -0.505 -0.235]]
Beta Fock matrix:
[[-0.075 -0.377 -0.473]
[-0.377 -0.075 -0.473]
[-0.473 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.506 -0.505]
[-0.505 -0.505 -0.235]]
Beta Fock matrix:
[[-0.075 -0.377 -0.473]
[-0.377 -0.075 -0.473]
[-0.473 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.506 -0.505]
[-0.505 -0.505 -0.235]]
Beta Fock matrix:
[[-0.075 -0.377 -0.473]
[-0.377 -0.075 -0.473]
[-0.473 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.506 -0.505]
[-0.505 -0.505 -0.235]]
Beta Fock matrix:
[[-0.075 -0.377 -0.473]
[-0.377 -0.075 -0.473]
[-0.473 -0.473 -0.42 ]]
Alpha Fock matrix:
[[-0.506 -0.45 -0.505]
[-0.45 -0.506 -0.505]
[-0.505 -0.505 -0.235]]
Beta Fock matrix:
[[-0.075 -0.377 -0.473]
[-0.377 -0.075 -0.473]
[-0.473 -0.473 -0.42 ]]
print(qc_structure.groundStateEnergy() + gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value())
-1.3358471555881106
uhf_parameters = qc_structure.groundStateParameters()
print(uhf_parameters.orbitalEnergies().alpha)
[-0.706 -0.111 0.404]
print(uhf_parameters.orbitalEnergies().beta)
[-0.581 0.441 0.601]
print(plain_solver.numberOfIterations())
55
DIIS UHF SCF calculations#
environment = gqcpy.UHFSCFEnvironment_d.WithCoreGuess(N_alpha, N_beta, sq_hamiltonian, S)
diis_solver = gqcpy.UHFSCFSolver_d.DIIS(threshold=1.0e-04, maximum_number_of_iterations=1000, minimum_subspace_dimension=6, maximum_subspace_dimension=6) # the system is not converging very rapidly
qc_structure = gqcpy.UHF_d.optimize(diis_solver, environment)
print(qc_structure.groundStateEnergy() + gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value())
-1.335847157827339
uhf_parameters = qc_structure.groundStateParameters()
print(uhf_parameters.orbitalEnergies().alpha)
[-0.706 -0.111 0.404]
print(uhf_parameters.orbitalEnergies().beta)
[-0.581 0.441 0.601]
print(diis_solver.numberOfIterations())
136
Complex UHF SCF calculations#
The same approach can be used for complex UHF calculations. We just replace the “_d” with “_cd” indicating the use of complex doubles.
complex_spinor_basis = gqcpy.USpinOrbitalBasis_cd(molecule, "STO-3G")
complex_S = complex_spinor_basis.quantize(gqcpy.OverlapOperator())
complex_sq_hamiltonian = complex_spinor_basis.quantize(gqcpy.FQMolecularHamiltonian(molecule)) # 'sq' for 'second-quantized'
complex_environment = gqcpy.UHFSCFEnvironment_cd.WithComplexlyTransformedCoreGuess(N_alpha, N_beta, complex_sq_hamiltonian, complex_S)
complex_plain_solver = gqcpy.UHFSCFSolver_cd.Plain(threshold=1.0e-04, maximum_number_of_iterations=1000) # the system is not converging very rapidly
complex_qc_structure = gqcpy.UHF_cd.optimize(complex_plain_solver, complex_environment)
print(complex_qc_structure.groundStateEnergy().real + gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value())
-1.3358471547992559