Package Guide

Define Grid for Calculations

Start by defining grid used in calculations

using HKQM

ceg = ElementGridSymmetricBox(5u"Å", 4, 24)

ElementGridSymmetricBox 4^3 elements and 24^3 Gauss points per element

This creates a cubic box with with side length of 5 Å that are divided to 4 elements and 24 Gauss-Lagrange points for each element. Resulting in total of (4*24)^3=884736 points.

The grid is also an Array that can be used as one

typeof(ceg) <: AbstractArray

true

size(ceg)

(24, 24, 24, 4, 4, 4)

eltype(ceg)

SVector{3, Float64} (alias for StaticArraysCore.SArray{Tuple{3}, Float64, 1, 3})

The values are x-, y- and z-coordinates of the grid point in bohr.

ceg[1,1,1,3,3,3]

3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 0.005684272564815695
 0.005684272564815695
 0.005684272564815695

Operator algebra

Generate basic operators for the grid

r = position_operator(ceg)
p = momentum_operator(ceg)

Operator 4^3 elements, 24^3 Gauss points per element

Operators have units defined with Unitful.jl package

unit(r)

a₀

unit(p)

ħ a₀^-1

These operator are vector operator and have length defined

length(r)

Individual components can be accessed with indexing

x = r[1]
y = r[2]
z = r[3]

Operators support basic algebra operations

r + r
2 * r
r + [1u"bohr", 2u"Å", 1u"pm"]
r + 1u"bohr"
r / 2
x * y

Units are checked for the operations and operations that do not make sense are prohibited

r + p

Vector operations are supported

r² = r ⋅ r
l  = r × p

Common functions can be used also

sin(1u"bohr^-1" * x)
exp(-1u"bohr^-2" * r²)

Functions require that the input is unitless

exp(-r²)

Quantum States

Quantum states can be created from operators

ψ = QuantumState( exp(-2u"bohr^-2" * r²) )

Quantum state

States can be normalized

normalize!(ψ)

Quantum state

Complex conjugate can be taken with

ϕ = conj(ψ)
conj!(ψ)

Quantum states have linear algebra defined

2ψ - ψ

Quantum state

Inner product can be calculated with bracket function

bracket(ψ, 2ψ)

2.0

Operators can be applied to quantum state by multiplication

x * ψ

Quantum state

Vector operators return arrays of quantum state

r * ψ

3-element Vector{QuantumState{Array{Float64, 6}, Float64}}:
 Quantum state
 Quantum state
 Quantum state

Quantum states have units

unit(ψ)

unit( x * ψ )

a₀

Other Unitful functions like dimension and uconvert are defined also.

Expectational values of operators can be calculated with bracket funtion

bracket(ψ, x, ψ)

-2.2188173684881792e-18 a₀

Slater Determinant

Slater determinant is orthonormal set of quantum states

st = SlaterDeterminant([ψ, (1u"bohr^-1"*x)*ψ])

SlaterDetermiant 2 orbitals

Slater determinant is an array of orbitals represented by quantum states

length(st)

st[1]

Quantum state

Hamilton Operator

Hamilton operator is a special operator that is needed for Helmholtz Greens function. To create it you need to create potential energy operator first

V = -30u"eV" * exp(-0.1u"bohr^-2" * r²)

H = HamiltonOperator(V)

Operator 4^3 elements, 24^3 Gauss points per element

Default mass is one electron mass and it can be customized with m keyword

H_2me = HamiltonOperator(V; m=2u"me_au")

Solving Eigen-States of a Hamiltonian

You need to generate initial state for Hamiltonian that gives negative energy!

ψ = QuantumState( exp(-0.2u"bohr^-2" * r²) )
normalize!(ψ)

bracket(ψ, H, ψ)

-0.4887869883070721 Eₕ

After that Helmholtz Greens function can be used to generate better estimate for the lowest eigen-state

ϕ, E = solve_eigen_states(H, ψ; max_iter=10, rtol=1E-6)

(QuantumState[Quantum state], Quantity{Float64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(Eₕ,), 𝐋^2 𝐌 𝐓^-2, nothing}}[-0.4978561191412073 Eₕ])

You can add more states to the solution by giving more initial states

ψ111 = particle_in_box(ceg, 1,1,1)
ψ112 = particle_in_box(ceg, 1,1,2)

ϕ2, E2 = solve_eigen_states(H, ψ111, ψ112)

(QuantumState[Quantum state, Quantum state], Quantity{Float64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(Eₕ,), 𝐋^2 𝐌 𝐓^-2, nothing}}[-0.4978559490725577 Eₕ, -0.17245088784778714 Eₕ])

Once estimate is self consistent a true solution has been found.

Solving Hartree-Fock equation

Hartree-Fock equation can be solver with scf command.

V = -100u"eV" * exp(-0.1u"bohr^-2" * r²)
H = HamiltonOperator(V)

ψ₁ = QuantumState( exp(-0.2u"bohr^-2" * r²) )
ψ₂ = 1u"bohr^-1"*r[1]*QuantumState( exp(-0.2u"bohr^-2" * r²) )
sd = SlaterDeterminant(ψ₁, ψ₂)

SlaterDetermiant 2 orbitals

To check that all eigen-values are negative calculate Fock matrix and look for diagonal values.

fock_matrix(sd, H)

After that you can solve Hartree-Fock equations

sd1 = scf(sd, H; tol=1E-6, max_iter=10)

tol is maximum chance in orbital overlap until convergence is archived. max_iter is maximum iterations calculated.

Hartree-Fock energy is calculated by calling hf_energy

hf_energy(sd1, H)

Orbital energies can be found from diagonal of Fock matrix

fock_matrix(sd1, H)

Check also that off-diagonal elements are insignificant to make sure the system real solution has been found.

Approximate Nuclear Potential

There is an approximation to nuclear potential defined in J. Chem. Phys. 121, 11587 (2004). It allows approximate electronic structure calculations.

Here is an example for Hydrogen molecule.

Define nuclear positions

r₁ = [0.37, 0., 0.] .* 1u"Å"
r₂ = [-0.37, 0., 0.] .* 1u"Å"

3-element Vector{Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}}:
 -0.37 Å
   0.0 Å
   0.0 Å

After that create nuclear potential

V₁ = nuclear_potential_harrison_approximation(ceg, r₁, "H")
V₂ = nuclear_potential_harrison_approximation(ceg, r₂, "H")

V = V₁ + V₂

Operator 4^3 elements, 24^3 Gauss points per element

and Hamiltonian

H = HamiltonOperator(V)

Operator 4^3 elements, 24^3 Gauss points per element

Create initial orbital

ϕ = particle_in_box(ceg, 1,1,1)
ψ = SlaterDeterminant( ϕ )

SlaterDetermiant 1 orbitals

Solve SCF equations

ψ1 = scf(ψ, H)

Calculations in Magnetic Field

To calculate magnetic field effects on the above hydrogen molecule, start by defining magnetic field

B = [0., 0., 100.0].*u"T"
A = vector_potential(ceg, B...)

Operator 4^3 elements, 24^3 Gauss points per element

After that create a magnetic field Hamiltonian

Hm = HamiltonOperatorMagneticField(V,A)

Operator 4^3 elements, 24^3 Gauss points per element

You can then just solve the system like above

ψm = scf(ψ, Hm)

It is a good idea to start form the no magnetic field case. Because wavefunction is complex in magnetic field and thus more expensive to calculate. By starting from no field case you minimize the iterations needed to solve SCF equations.

ψm = scf(ψ1, Hm)

Magnetic current can be extracted with

# current for specific orbital
j1 = magnetic_current(ψm[1], Hm)

# total magnetic current
j = magnetic_current(ψm, Hm)

# para magnetic current
jp = para_magnetic_current(ψm)

In the future there will be a way to visualize or export the magnetic current!