Gauge-origin independence for London kinetic integrals#
# 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=8, linewidth=120)
nuclei = [gqcpy.Nucleus(1, 0.0, 0.0, 0.0), gqcpy.Nucleus(1, 0.0, 0.0, 1.0)]
molecule = gqcpy.Molecule(nuclei)
Kinetic energy integrals#
Let’s prepare two London-spinor bases, one for each gauge origin.
B1 = gqcpy.HomogeneousMagneticField([0.0, 0.0, 0.1], [0.0, 0.0, 0.0]) # Gauge origin at the origin.
B2 = gqcpy.HomogeneousMagneticField([0.0, 0.0, 0.1], [25.0, -3.2, 0.6]) # Gauge origin at a random point.
spinor_basis1 = gqcpy.LondonGSpinorBasis(molecule, "STO-3G", B1)
spinor_basis2 = gqcpy.LondonGSpinorBasis(molecule, "STO-3G", B2)
From the formulas, it is immediately clear that the overlap integrals are gauge origin independent. Let’s confirm this.
S1 = spinor_basis1.quantize(gqcpy.OverlapOperator()).parameters()
S2 = spinor_basis2.quantize(gqcpy.OverlapOperator()).parameters()
print(S1)
[[0.99999999+0.j 0.79658829+0.j 0. +0.j 0. +0.j]
[0.79658829+0.j 0.99999999+0.j 0. +0.j 0. +0.j]
[0. +0.j 0. +0.j 0.99999999+0.j 0.79658829+0.j]
[0. +0.j 0. +0.j 0.79658829+0.j 0.99999999+0.j]]
np.allclose(S1, S2)
True
Furthermore, the nuclear attraction and two-electron repulsion integrals are also gauge origin independent.
V1 = spinor_basis1.quantize(gqcpy.NuclearAttractionOperator(molecule.nuclearFramework())).parameters()
V2 = spinor_basis2.quantize(gqcpy.NuclearAttractionOperator(molecule.nuclearFramework())).parameters()
print(V1)
[[-2.03852055+0.j -1.60241665+0.j 0. +0.j 0. +0.j]
[-1.60241665+0.j -2.03852055+0.j 0. +0.j 0. +0.j]
[ 0. +0.j 0. +0.j -2.03852055+0.j -1.60241665+0.j]
[ 0. +0.j 0. +0.j -1.60241665+0.j -2.03852055+0.j]]
np.allclose(V1, V2)
True
g1 = spinor_basis1.quantize(gqcpy.CoulombRepulsionOperator()).parameters()
g2 = spinor_basis2.quantize(gqcpy.CoulombRepulsionOperator()).parameters()
#print(g1)
np.allclose(g1, g2)
True
Now, let’s prepare the kinetic energy integrals in both spinor basis. (Note that we are calculating only the term of the magnetic kinetic energy operator that is zero-th order in the magnetic field.)
T1 = spinor_basis1.quantize(gqcpy.KineticOperator()).parameters()
T2 = spinor_basis2.quantize(gqcpy.KineticOperator()).parameters()
print(T1)
[[0.76003188+0.j 0.38325367+0.j 0. +0.j 0. +0.j]
[0.38325367+0.j 0.76003188+0.j 0. +0.j 0. +0.j]
[0. +0.j 0. +0.j 0.76003188+0.j 0.38325367+0.j]
[0. +0.j 0. +0.j 0.38325367+0.j 0.76003188+0.j]]
print(T2)
[[1.55408187+0.j 1.0157846 +0.j 0. +0.j 0. +0.j]
[1.0157846 +0.j 1.55408187+0.j 0. +0.j 0. +0.j]
[0. +0.j 0. +0.j 1.55408187+0.j 1.0157846 +0.j]
[0. +0.j 0. +0.j 1.0157846 +0.j 1.55408187+0.j]]
Apparently, the matrix representation of the operator $-\dfrac{1}{2} \nabla^2$ is gauge-origin dependent! Let’s take a further look at the orbital Zeeman (first order in the magnetic field) and diamagnetic (second order in the magnetic field) operators.
P1 = spinor_basis1.quantize(gqcpy.OrbitalZeemanOperator(B1, B1.gaugeOrigin())).parameters()
P2 = spinor_basis2.quantize(gqcpy.OrbitalZeemanOperator(B2, B2.gaugeOrigin())).parameters()
print(P1)
[[0.+0.j 0.+0.j 0.+0.j 0.+0.j]
[0.+0.j 0.+0.j 0.+0.j 0.+0.j]
[0.+0.j 0.+0.j 0.+0.j 0.+0.j]
[0.+0.j 0.+0.j 0.+0.j 0.+0.j]]
print(P2)
[[-1.58809999+0.j -1.26506187+0.j 0. +0.j 0. +0.j]
[-1.26506187+0.j -1.58809999+0.j 0. +0.j 0. +0.j]
[ 0. +0.j 0. +0.j -1.58809999+0.j -1.26506187+0.j]
[ 0. +0.j 0. +0.j -1.26506187+0.j -1.58809999+0.j]]
D1 = spinor_basis1.quantize(gqcpy.DiamagneticOperator(B1, B1.gaugeOrigin())).parameters()
D2 = spinor_basis2.quantize(gqcpy.DiamagneticOperator(B2, B2.gaugeOrigin())).parameters()
print(D1)
[[0.00162381+0.j 0.00140007+0.j 0. +0.j 0. +0.j]
[0.00140007+0.j 0.00162381+0.j 0. +0.j 0. +0.j]
[0. +0.j 0. +0.j 0.00162381+0.j 0.00140007+0.j]
[0. +0.j 0. +0.j 0.00140007+0.j 0.00162381+0.j]]
print(D2)
[[0.7956738 +0.j 0.63393101+0.j 0. +0.j 0. +0.j]
[0.63393101+0.j 0.7956738 +0.j 0. +0.j 0. +0.j]
[0. +0.j 0. +0.j 0.7956738 +0.j 0.63393101+0.j]
[0. +0.j 0. +0.j 0.63393101+0.j 0.7956738 +0.j]]
We can see that these operators are also gauge origin dependent!
However, by including the orbital Zeeman and diamagnetic terms, the magnetic kinetic energy becomes gauge invariant, as shown below.
K1 = T1 + P1 + D1
K2 = T2 + P2 + D2
print(K1)
[[0.76165569+0.j 0.38465374+0.j 0. +0.j 0. +0.j]
[0.38465374+0.j 0.76165569+0.j 0. +0.j 0. +0.j]
[0. +0.j 0. +0.j 0.76165569+0.j 0.38465374+0.j]
[0. +0.j 0. +0.j 0.38465374+0.j 0.76165569+0.j]]
np.allclose(K1, K2)
True
We can also include the Spin-Zeeman operator to the kinetic energy to obtain the kinetic energy when minimal coupling is applied at the Pauli-Hamiltonian level.
SZ1 = spinor_basis1.quantize(gqcpy.SpinZeemanOperator(B1)).parameters()
SZ2 = spinor_basis2.quantize(gqcpy.SpinZeemanOperator(B2)).parameters()
print(SZ1, SZ2, sep='\n\n')
[[ 0.05 +0.j 0.03982941+0.j 0. +0.j 0. +0.j]
[ 0.03982941+0.j 0.05 +0.j 0. +0.j 0. +0.j]
[ 0. +0.j 0. +0.j -0.05 +0.j -0.03982941+0.j]
[ 0. +0.j 0. +0.j -0.03982941+0.j -0.05 +0.j]]
[[ 0.05 +0.j 0.03982941+0.j 0. +0.j 0. +0.j]
[ 0.03982941+0.j 0.05 +0.j 0. +0.j 0. +0.j]
[ 0. +0.j 0. +0.j -0.05 +0.j -0.03982941+0.j]
[ 0. +0.j 0. +0.j -0.03982941+0.j -0.05 +0.j]]