CCSD calculations#
# 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)
CCSD calculations on H2O // STO-3G#
In this example, we will try to reproduce the crawdad CCSD results.
molecule = gqcpy.Molecule.ReadXYZ("../../gqcp/tests/data/h2o_crawdad.xyz")
N = molecule.numberOfElectrons()
N_P = molecule.numberOfElectronPairs()
print(molecule)
Number of electrons: 10
O (0, -0.143226, 0)
H (1.63804, 1.13655, -0)
H (-1.63804, 1.13655, -0)
The first task is to generate the RHF reference, so we’ll do an RHF SCF calculation.
spin_orbital_basis = gqcpy.RSpinOrbitalBasis_d(molecule, "STO-3G")
K = spin_orbital_basis.numberOfSpatialOrbitals()
S = spin_orbital_basis.quantize(gqcpy.OverlapOperator())
fq_hamiltonian = gqcpy.FQMolecularHamiltonian(molecule)
hamiltonian = spin_orbital_basis.quantize(fq_hamiltonian)
The current Hamiltionian is expressed in the non-orthogonal AO basis. This is exactly what we need to start an RHF calculation.
environment = gqcpy.RHFSCFEnvironment_d.WithCoreGuess(N, hamiltonian, S)
solver = gqcpy.RHFSCFSolver_d.DIIS()
objective = gqcpy.DiagonalRHFFockMatrixObjective_d(hamiltonian) # use the default threshold of 1.0e-08
Using this objective makes sure that the optimized expansion coefficients yield a diagonal Fock matrix, so they will represent the canonical RHF spin-orbital basis.
rhf_qc_structure = gqcpy.RHF_d.optimize(objective, solver, environment)
rhf_energy = rhf_qc_structure.groundStateEnergy() + gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value()
print(rhf_energy)
-74.9420799280923
This is very much in line with crawdad’s results: -74.942079928192.
rhf_parameters = rhf_qc_structure.groundStateParameters()
C = rhf_parameters.expansion()
Now, let’s transform the spin-orbital basis and the second-quantized Hamiltonian to the RHF MOs.
spin_orbital_basis.transform(C) # Now represents a restricted spin-orbital that contains the canonical RHF MOs.
print(spin_orbital_basis.expansion().matrix())
[[-9.944e-01 -2.392e-01 2.405e-16 -9.368e-02 4.245e-31 1.116e-01 1.180e-16]
[-2.410e-02 8.857e-01 -1.159e-15 4.796e-01 -2.150e-30 -6.696e-01 -7.428e-16]
[ 2.047e-18 -2.675e-16 -6.073e-01 -9.061e-16 1.729e-17 -8.882e-16 9.192e-01]
[-3.162e-03 8.590e-02 1.498e-15 -7.474e-01 2.613e-30 -7.385e-01 -4.996e-16]
[ 1.524e-33 1.173e-31 -6.203e-17 3.753e-30 1.000e+00 -3.628e-32 -2.313e-16]
[ 4.594e-03 1.440e-01 -4.530e-01 -3.295e-01 -1.272e-16 7.098e-01 -7.325e-01]
[ 4.594e-03 1.440e-01 4.530e-01 -3.295e-01 1.272e-16 7.098e-01 7.325e-01]]
hamiltonian.transform(C)
print(hamiltonian.core().parameters())
[[-3.256e+01 -5.683e-01 5.022e-16 -2.095e-01 9.840e-31 2.563e-01 3.331e-16]
[-5.683e-01 -7.564e+00 1.141e-15 -5.223e-01 2.425e-30 1.220e+00 1.332e-15]
[ 5.251e-16 1.300e-15 -6.018e+00 1.825e-15 -9.824e-17 6.295e-16 1.767e+00]
[-2.095e-01 -5.223e-01 1.814e-15 -6.610e+00 -3.945e-30 -1.306e+00 2.442e-15]
[ 1.002e-30 2.451e-30 -9.824e-17 -3.977e-30 -7.347e+00 3.776e-30 1.368e-15]
[ 2.563e-01 1.220e+00 4.025e-16 -1.306e+00 3.744e-30 -5.291e+00 2.220e-16]
[ 3.619e-16 1.562e-15 1.767e+00 2.409e-15 1.368e-15 3.246e-16 -5.513e+00]]
Before continuing with the CCSD specifics, we must first prepare the molecular Hamiltonian in the correct spinor basis. Since we have implemented CCSD using general spinors, we should use a GSpinorBasis.
spinor_basis = gqcpy.GSpinorBasis_d.FromRestricted(spin_orbital_basis) # Represents a general spinor basis, based off the restricted canonical RHF MOs.
# We can inspect the two non-zero blocks (top-left and bottom-right) of the coefficient matrix.
print(spinor_basis.expansion().matrix())
[[-9.944e-01 -2.392e-01 2.405e-16 -9.368e-02 4.245e-31 1.116e-01 1.180e-16 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00]
[-2.410e-02 8.857e-01 -1.159e-15 4.796e-01 -2.150e-30 -6.696e-01 -7.428e-16 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00]
[ 2.047e-18 -2.675e-16 -6.073e-01 -9.061e-16 1.729e-17 -8.882e-16 9.192e-01 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00]
[-3.162e-03 8.590e-02 1.498e-15 -7.474e-01 2.613e-30 -7.385e-01 -4.996e-16 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00]
[ 1.524e-33 1.173e-31 -6.203e-17 3.753e-30 1.000e+00 -3.628e-32 -2.313e-16 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00]
[ 4.594e-03 1.440e-01 -4.530e-01 -3.295e-01 -1.272e-16 7.098e-01 -7.325e-01 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00]
[ 4.594e-03 1.440e-01 4.530e-01 -3.295e-01 1.272e-16 7.098e-01 7.325e-01 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00]
[ 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 -9.944e-01 -2.392e-01 2.405e-16
-9.368e-02 4.245e-31 1.116e-01 1.180e-16]
[ 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 -2.410e-02 8.857e-01 -1.159e-15
4.796e-01 -2.150e-30 -6.696e-01 -7.428e-16]
[ 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.047e-18 -2.675e-16 -6.073e-01
-9.061e-16 1.729e-17 -8.882e-16 9.192e-01]
[ 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 -3.162e-03 8.590e-02 1.498e-15
-7.474e-01 2.613e-30 -7.385e-01 -4.996e-16]
[ 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 1.524e-33 1.173e-31 -6.203e-17
3.753e-30 1.000e+00 -3.628e-32 -2.313e-16]
[ 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 4.594e-03 1.440e-01 -4.530e-01
-3.295e-01 -1.272e-16 7.098e-01 -7.325e-01]
[ 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 4.594e-03 1.440e-01 4.530e-01
-3.295e-01 1.272e-16 7.098e-01 7.325e-01]]
We should note that the ordering of the MOs in this spinor basis is not by ascending energy, but they are ordered alphas first and then all betas.
For H2O, we have 10 electrons, so this means we shouldn’t occupy the first 10 spinors (since that would correspond to occupying all alpha-spin-orbitals and then 3 beta-spin-orbitals), but we should occupy the first 5 alpha spin-orbitals and the first 5 beta-spin-orbitals.
This issue is circumvented by constructing an OrbitalSpace, which takes care of setting up the occupied and virtual orbital spaces. We don’t have to construct this OrbitalSpace by hand: GQCP provides an API through the creation of the ‘GHF’ reference ONV.
M = spinor_basis.numberOfSpinors()
ghf_onv = gqcpy.SpinUnresolvedONV.GHF(M, N, rhf_parameters.spinOrbitalEnergiesBlocked())
orbital_space = ghf_onv.orbitalSpace()
print(orbital_space.description())
An orbital space with 14 orbitals.
The occupied orbital indices are: 0 1 2 3 4 7 8 9 10 11
The active orbital indices are:
The virtual orbital indices are: 5 6 12 13
CCSD is also implemented as a QCMethod. This means that we’ll have to set up an environment and a particular solver.
hamiltonian = spinor_basis.quantize(fq_hamiltonian)
solver = gqcpy.CCSDSolver_d.Plain()
environment = gqcpy.CCSDEnvironment_d.PerturbativeCCSD(hamiltonian, orbital_space)
It’s interesting to see that the initial, perturbative, T1 amplitudes are zero, which is a direct consequence of Brillouin’s theorem: the occupied-virtual block of the Fock(ian) matrix is zero!
print(environment.t1_amplitudes[-1].asMatrix())
[[ 2.953e-15 -1.205e-18 -0.000e+00 -0.000e+00]
[ 5.534e-14 8.803e-17 -0.000e+00 -0.000e+00]
[ 1.582e-16 2.386e-13 -0.000e+00 -0.000e+00]
[-2.485e-13 -5.697e-16 -0.000e+00 -0.000e+00]
[ 2.030e-32 -1.052e-16 -0.000e+00 -0.000e+00]
[-0.000e+00 -0.000e+00 2.995e-15 5.041e-18]
[-0.000e+00 -0.000e+00 5.578e-14 3.847e-16]
[-0.000e+00 -0.000e+00 9.113e-17 2.385e-13]
[-0.000e+00 -0.000e+00 -2.496e-13 -1.666e-16]
[-0.000e+00 -0.000e+00 -4.396e-33 -1.052e-16]]
The initial T2-amplitudes are actually the MP2 T2-amplitudes, so the initial CCSD correlation energy is the MP2 correlation energy. Crawdad reports a value of -0.049149636120.
print(environment.correlation_energies[-1])
-0.04914963614943018
Let’s now proceed by actually optimizing the CCSD wave function model. This may take a while if gqcp is built in Debug mode!
ccsd_qc_structure = gqcpy.CCSD_d.optimize(solver, environment)
Crawdad lists a converged CCSD energy as: -0.070680088376.
ccsd_correlation_energy = ccsd_qc_structure.groundStateEnergy()
print(ccsd_correlation_energy)
-0.0706800879370939