Calculating the RHF dipole moment#

In this example, we will calculate the (permanent) electronic dipole moment for H2 with an STO-3G basisset at the RHF level of theory.

Setting up the molecular Hamiltonian#

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

import gqcpy

We start by creating the molecule and an associated (restricted) spinor basis:

molecule = gqcpy.Molecule.ReadXYZ("../../gqcp/tests/data/h2_szabo.xyz" , 0)  # create a neutral molecule
N = molecule.numberOfElectrons()

spinor_basis = gqcpy.RSpinOrbitalBasis_d(molecule, "STO-3G")

We can easily find out the overlap matrix of the spinor basis:

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

[[0.99999999 0.65931816]
 [0.65931816 0.99999999]]

We can use the following one-liner to calculate the molecular second-quantized Hamiltonian:

sq_hamiltonian = spinor_basis.quantize(gqcpy.FQMolecularHamiltonian(molecule))  # 'sq' for 'second-quantized'

We can access the core Hamiltonian and the two-electron integrals relatively easily:

print(sq_hamiltonian.core().parameters())

[[-1.12040896 -0.95837989]
 [-0.95837989 -1.12040896]]

print(sq_hamiltonian.twoElectron().parameters())

[[[[0.77460593 0.44410762]
   [0.44410762 0.56967589]]

  [[0.44410762 0.2970285 ]
   [0.2970285  0.44410762]]]


 [[[0.44410762 0.2970285 ]
   [0.2970285  0.44410762]]

  [[0.56967589 0.44410762]
   [0.44410762 0.77460593]]]]

Solving the RHF SCF equations#

In order to solve the RHF SCF equations, we have to set up a solver and its associated environment.

environment = gqcpy.RHFSCFEnvironment_d.WithCoreGuess(N, sq_hamiltonian, S)
solver = gqcpy.RHFSCFSolver_d.DIIS()

The specification of a QCMethod requires something more than just the solution that a solver produces. In order to really confirm that the electronic structure model’s parameters are ‘optimal’, an objective has to be defined.

objective = gqcpy.DiagonalRHFFockMatrixObjective_d(sq_hamiltonian)  # use the default threshold of 1.0e-08

Note that, since we have chosen to use a ‘diagonal Fock matrix objective’, we expect the optimal RHF parameters to represent the canonical spinors.

We then combine objective, solver and environment into the optimize method of the RHF QCMethod, which returns a QCStructure, containing the optimized RHF parameters that satisfy the objective.

rhf_parameters = gqcpy.RHF_d.optimize(objective, solver, environment).groundStateParameters()

C = rhf_parameters.expansion()
print(C.matrix())

[[-0.54893405 -1.21146402]
 [-0.54893405  1.21146402]]

Calculating the RHF dipole moment in the scalar/AO basis#

Since the RHF parameters are variationally determined, we may calculate the electronic dipole moment as the expectation value of the electronic dipole operator. The associated second-quantized operator can be obtained as follows:

dipole_op = spinor_basis.quantize(gqcpy.ElectronicDipoleOperator())  # 'op' for 'operator'

We can access the bare integrals using its method .allParameters(), which returns a list that contains the integrals associated with the Cartesian x-, y- and z-components.

integrals = dipole_op.allParameters()
for component in integrals:
    print(component)

[[0. 0.]
 [0. 0.]]
[[0. 0.]
 [0. 0.]]
[[ 0.         -0.46152275]
 [-0.46152275 -1.40000009]]

In the scalar/AO basis, the RHF 1-DM can be calculated as:

D = rhf_parameters.calculateScalarBasis1DM()
print(D.matrix())

[[0.60265718 0.60265718]
 [0.60265718 0.60265718]]

print(dipole_op.calculateExpectationValue(D))

[0.0, 0.0, -1.4000001021623736]

Calculating the RHF dipole moment in the canonical RHF basis#

We might as well have calculated the second-quantized representation of the dipole operator in the canonical RHF spinor basis. In order to do so, we first transform the spinor basis and only then quantize the dipole operator.

spinor_basis.transform(C)
dipole_op_canonical = spinor_basis.quantize(gqcpy.ElectronicDipoleOperator())  # 'op' for 'operator'

dipole_integrals_canonical = dipole_op_canonical.allParameters()
for component in dipole_integrals_canonical:
    print(component)

[[0. 0.]
 [0. 0.]]
[[0. 0.]
 [0. 0.]]
[[-0.70000005  0.93101945]
 [ 0.93101945 -0.70000005]]

D_canonical = rhf_parameters.calculateOrthonormalBasis1DM()
print(D_canonical.matrix())

[[2. 0.]
 [0. 0.]]

print(dipole_op_canonical.calculateExpectationValue(D_canonical))

[0.0, 0.0, -1.4000001021623738]