The ipsocentric magnetic inducibility

Contents

The ipsocentric magnetic inducibility#

# 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)

Molecular setup#

molecule = gqcpy.Molecule.HChain(2, 1.0)
N = molecule.numberOfElectrons()
print(molecule)

Number of electrons: 2 
H  (0, 0, 0)
H  (0, 0, 1)

Create a grid#

origin = np.array([-2.0, -2.0, -2.0])
steps = [5, 5, 5]
step_sizes = [0.8, 0.8, 0.8]

grid = gqcpy.CubicGrid(origin, steps, step_sizes)

Set up an RHF calculation#

spin_orbital_basis = gqcpy.RSpinOrbitalBasis_d(molecule, "STO-3G")
print("Number of orbitals: ", spin_orbital_basis.numberOfSpatialOrbitals())
hamiltonian = spin_orbital_basis.quantize(gqcpy.FQMolecularHamiltonian(molecule))

Number of orbitals:  2

S = spin_orbital_basis.quantize(gqcpy.OverlapOperator())
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.

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

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

[[-0.527546 -1.567823]
 [-0.527546  1.567823]]

E_RHF_total = qc_structure.groundStateEnergy() + gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value()
print("Total RHF energy: ", E_RHF_total)

Total RHF energy:  -1.0659994621433042

C_complex = gqcpy.RTransformation_cd(C.matrix().astype(complex))
complex_spin_orbital_basis = gqcpy.RSpinOrbitalBasis_cd(molecule, "STO-3G")
complex_spin_orbital_basis.transform(C_complex)

spin_orbital_basis.transform(C)
hamiltonian.transform(C)

orbital_space = rhf_parameters.orbitalSpace()

k_kappa = rhf_parameters.calculateOrbitalHessianForImaginaryResponse(hamiltonian, orbital_space)

print("k_kappa (Response force constant matrix A in Ax=-b)")
print(k_kappa)

k_kappa (Response force constant matrix A in Ax=-b)
[[0.-4.34392j]]

L = complex_spin_orbital_basis.quantize(gqcpy.AngularMomentumOperator())
print("Integrals over the angular momentum operator (x,y,z)")
for L_m in L.allParameters():
    print(L_m)

Integrals over the angular momentum operator (x,y,z)
[[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]]

F_kappa_B = rhf_parameters.calculateMagneticFieldResponseForce(L)

print("F_kappa_B (Response force vector b in Ax=-b) for magnetic field perturbation")
print(F_kappa_B)

F_kappa_B (Response force vector b in Ax=-b) for magnetic field perturbation
[[-0.-0.j -0.-0.j -0.-0.j]]

x_B = np.linalg.solve(k_kappa, -F_kappa_B)
print("x_B (Linear response x in Ax=-b) for magnetic field perturbation")
print(x_B)

x_B (Linear response x in Ax=-b) for magnetic field perturbation
[[0.+0.j 0.+0.j 0.+0.j]]

p = complex_spin_orbital_basis.quantize(gqcpy.LinearMomentumOperator())
print("Integrals over the linear momentum operator (x,y,z)")
for p_m in p.allParameters():
    print(p_m)

Integrals over the linear momentum operator (x,y,z)
[[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.54894j]
 [0.-0.54894j 0.+0.j     ]]

F_kappa_G = rhf_parameters.calculateGaugeOriginTranslationResponseForce(p)
print("F_kappa_G (Response force vector b in Ax=-b) for gauge origin translation perturbation")
print(F_kappa_G)

F_kappa_G (Response force vector b in Ax=-b) for gauge origin translation perturbation
[[ 0.+1.097881j -0.-0.j       -0.-1.097881j  0.+0.j        0.+0.j       -0.-0.j      ]]

x_G = np.linalg.solve(k_kappa, -F_kappa_G)
print("x_G (Linear response x in Ax=-b) for gauge origin translation perturbation:")
print(x_G)

x_G (Linear response x in Ax=-b) for gauge origin translation perturbation:
[[ 0.25274-0.j  0.     +0.j -0.25274+0.j  0.     -0.j  0.     -0.j  0.     +0.j]]

j_op = complex_spin_orbital_basis.quantize(gqcpy.CurrentDensityOperator())

J_field = gqcpy.QCModel_RHF_cd.calculateIpsocentricMagneticInducibility(grid, orbital_space, x_B, x_G, j_op)
J_field = np.array(J_field.values())

print(J_field)

[[[ 0.000000e+00+0.j -1.938959e-06+0.j  0.000000e+00+0.j]
  [ 0.000000e+00+0.j -1.938959e-06+0.j  0.000000e+00+0.j]
  [ 0.000000e+00+0.j  2.801371e-05+0.j  0.000000e+00+0.j]]

 [[ 0.000000e+00+0.j -1.696106e-05+0.j  0.000000e+00+0.j]
  [ 0.000000e+00+0.j -1.696106e-05+0.j  0.000000e+00+0.j]
  [ 0.000000e+00+0.j  8.853624e-05+0.j  0.000000e+00+0.j]]

 [[ 0.000000e+00+0.j -4.354813e-05+0.j  0.000000e+00+0.j]
  [ 0.000000e+00+0.j -4.354813e-05+0.j  0.000000e+00+0.j]
  [ 0.000000e+00+0.j  2.156915e-04+0.j  0.000000e+00+0.j]]

 ...

 [[ 0.000000e+00+0.j -9.110816e-04+0.j  0.000000e+00+0.j]
  [ 0.000000e+00+0.j -9.110816e-04+0.j  0.000000e+00+0.j]
  [ 0.000000e+00+0.j -2.298519e-03+0.j  0.000000e+00+0.j]]

 [[ 0.000000e+00+0.j -1.874316e-04+0.j  0.000000e+00+0.j]
  [ 0.000000e+00+0.j -1.874316e-04+0.j  0.000000e+00+0.j]
  [ 0.000000e+00+0.j -5.545698e-03+0.j  0.000000e+00+0.j]]

 [[ 0.000000e+00+0.j  9.056225e-04+0.j  0.000000e+00+0.j]
  [ 0.000000e+00+0.j  9.056225e-04+0.j  0.000000e+00+0.j]
  [ 0.000000e+00+0.j -3.270967e-03+0.j  0.000000e+00+0.j]]]