Mean-field calculations in magnetic fields

Mean-field calculations in magnetic fields#

In this interactive Notebook, we’ll explore how GQCP can be used to perform electronic structure calculations in the presence of strong magnetic fields.

Since the presence of the magnetic field governs the physical interactions present in the system, we provide two types of magnetic Hamiltonians:

FQMolecularMagneticHamiltonian: The magnetic Hamiltonian that contains the orbital Zeeman and diamagnetic operators.
FQMolecularPauliHamiltonian: The magnetic Hamiltonian that contains the orbital Zeeman, diamagnetic and spin Zeeman operators.

import sys
sys.path.insert(0, '../../build/gqcpy/')

import gqcpy
import numpy as np

In this Notebook, we’ll explore the dissociation of the simple H₂ molecule in the presence of a magnetic field perpendicular on the bond axis.

R_values = np.arange(0.7, 8, 0.1)
H_left = gqcpy.Nucleus(1, 0.0,0.0,0.0)

molecules = []
for R in R_values:
    H_right = gqcpy.Nucleus(1, R,0.0,0.0)  # The molecule is aligned on the x-axis.
    molecules.append(gqcpy.Molecule([H_left, H_right]))

# Apply a magnetic field in the z-axis.
B = gqcpy.HomogeneousMagneticField([0,0,-0.1])

basis_set = "6-31G"

RHF in a magnetic field#

Since the expectation value of the spin operator for RHF is identically zero, the spin Zeeman operator does not contribute to the RHF energy. Therefore, we can use the FQMolecularMagneticHamiltonian.

RHF_energies = []
for molecule in molecules:
    
    N = molecule.numberOfElectrons()

    # Set up the molecular magnetic Hamiltonian in the AO basis.
    spin_orbital_basis = gqcpy.LondonRSpinOrbitalBasis(molecule, basis_set, B)
    hamiltonian = spin_orbital_basis.quantize(gqcpy.FQMolecularMagneticHamiltonian(molecule, B))
    
    # Do an RHF calculation and calculate the total energy.
    S = spin_orbital_basis.quantize(gqcpy.OverlapOperator())
    objective = gqcpy.DiagonalRHFFockMatrixObjective_cd(hamiltonian, 1.0e-03)
    environment = gqcpy.RHFSCFEnvironment_cd.WithCoreGuess(N, hamiltonian, S)
    solver = gqcpy.RHFSCFSolver_cd.Plain(threshold=1.0e-04, maximum_number_of_iterations=1000)
    
    qc_structure = gqcpy.RHF_cd.optimize(objective, solver, environment)
    
    electronic_energy = qc_structure.groundStateEnergy()
    nuclear_repulsion = gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value()

    RHF_energies.append((electronic_energy + nuclear_repulsion).real)  # Omit the imaginary component, since it's zero.

UHF in a magnetic field#

In the UHF case, we can do a singlet calculation (i.e. $S_z = 0$ or $N_\alpha = N_\beta$) or a triplet calculation (i.e. $S_z = 1$ or $N_\alpha = 2$ with $N_\beta = 0$). Since the contribution from the spin Zeeman term is constant, we can do the UHF calculation using a FQMolecularMagneticHamiltonian.

UHF_0_energies = []
UHF_1_energies = []
for molecule in molecules:
    
    # Set up the molecular magnetic Hamiltonian in the AO basis.
    spin_orbital_basis = gqcpy.LondonUSpinOrbitalBasis(molecule, basis_set, B)
    hamiltonian = spin_orbital_basis.quantize(gqcpy.FQMolecularMagneticHamiltonian(molecule, B))
    
    S = spin_orbital_basis.quantize(gqcpy.OverlapOperator())

    
    # Do an UHF calculation for S = 1, calculate the constant spin Zeeman contribution and calculate the total energy.
    N = molecule.numberOfElectrons()
    N_alpha = N
    N_beta = 0
    
    environment = gqcpy.UHFSCFEnvironment_cd.WithCoreGuess(N_alpha, N_beta, hamiltonian, S)
    solver = gqcpy.UHFSCFSolver_cd.Plain(threshold=1.0e-04, maximum_number_of_iterations=1000)
    
    qc_structure = gqcpy.UHF_cd.optimize(solver, environment)
    
    electronic_energy = qc_structure.groundStateEnergy()
    nuclear_repulsion = gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value()
    spin_zeeman_contribution = 0.5 * B.strength()[2] * (N_alpha - N_beta)  # This corresponds to S dot B.

    UHF_1_energies.append((electronic_energy + nuclear_repulsion).real + spin_zeeman_contribution)  # Omit the imaginary component, since it's zero.


    # Do an UHF calculation for S = 0 and calculate the total energy. In this case, the spin Zeeman contribution vanishes.
    N_P = molecule.numberOfElectronPairs()
    N_alpha = N_P
    N_beta = N_P
    
    C_initial = qc_structure.groundStateParameters().expansion()
    environment = gqcpy.UHFSCFEnvironment_cd(N_alpha, N_beta, hamiltonian, S, C_initial)
    solver = gqcpy.UHFSCFSolver_cd.Plain(threshold=1.0e-04, maximum_number_of_iterations=1000)
    
    qc_structure = gqcpy.UHF_cd.optimize(solver, environment)
    
    electronic_energy = qc_structure.groundStateEnergy()
    nuclear_repulsion = gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value()

    UHF_0_energies.append((electronic_energy + nuclear_repulsion).real)  # Omit the imaginary component, since it's zero.

GHF in a magnetic field#

Since GHF is a spin non-collinear theory, its spin Zeeman contribution is not necessarily constant. Therefore, to allow the most flexible description in the presence of an external magnetic field, we use a FQMolecularPauliHamiltonian.

GHF_energies = []
for molecule in molecules:

    N = molecule.numberOfElectrons()

    # Set up the molecular magnetic Pauli Hamiltonian in the AO basis.
    spinor_basis = gqcpy.LondonGSpinorBasis(molecule, basis_set, B)
    hamiltonian = spinor_basis.quantize(gqcpy.FQMolecularPauliHamiltonian(molecule, B))
    
    # Do a GHF calculation and calculate the total energy.
    S = spinor_basis.quantize(gqcpy.OverlapOperator())

    # Initialize the GHF SCF environment with an initial guess that does not consist of spin-orbitals, but consists of true spinors that are not eigenvectors of S_z.
    environment = gqcpy.GHFSCFEnvironment_cd.WithComplexlyTransformedCoreGuess(N, hamiltonian, S)
    solver = gqcpy.GHFSCFSolver_cd.Plain(threshold=1.0e-04, maximum_number_of_iterations=10000)
    
    qc_structure = gqcpy.GHF_cd.optimize(solver, environment)
    
    electronic_energy = qc_structure.groundStateEnergy()
    nuclear_repulsion = gqcpy.NuclearRepulsionOperator(molecule.nuclearFramework()).value()

    GHF_energies.append((electronic_energy + nuclear_repulsion).real)  # Omit the imaginary component, since it's zero.

Visualizations#

import plotly
import plotly.graph_objects as go

figure = go.Figure()

figure.update_layout(
    autosize = False,
    width = 800,
    height = 800
)

# Plot the total energies for RHF, UHF and GHF.
figure.add_trace(
    go.Scatter(
        x = R_values,
        y = RHF_energies,
        name = "RHF",
    )
)

figure.add_trace(
    go.Scatter(
        x = R_values,
        y = UHF_0_energies,
        name = "UHF (S=0)",
    )
)

figure.add_trace(
    go.Scatter(
        x = R_values,
        y = UHF_1_energies,
        name = "UHF (S=1)",
    )
)

figure.add_trace(
    go.Scatter(
        x = R_values,
        y = GHF_energies,
        name = "GHF",
    )
)


# Set the x- and y-axis titles.
figure.update_xaxes(
    title = {
        'text': "Internuclear distance (a.u.)",
        'font': {'size': 18}
    }
)

figure.update_yaxes(
    title = {
        'text': "Energy (a.u.)",
        'font': {'size': 18}
    }
)


# Update the legend's font.
figure.update_layout(
    legend = {
        'font': {'size': 18}
    }
)


figure.show()