Presentation for the NKS-meeting 2021 by Audun Skau Hansen ✉️ at the Department of Chemistry, The Hylleraas Centre for Quantum Molecular Sciences, University of Oslo (2021)
import numpy as np # Numpy for numerics
import sympy as sp # We need sympy for defining our functions
import matplotlib.pyplot as plt # Matplotlib for visualization
import braketlab as bk # <- pip install braketlab (if you don't have it)
sp.init_printing() # <- pretty display of mathematical expressions
%matplotlib notebook
plt.figure.max_open_warning = False # <- surpress warnings from Matplotlib
Bra-ket notation
Students in physical chemistry at UiO are early on introduced to Dirac notation, named after Paul Dirac who introduced the formalism in a surprisingly short paper back in 1939. This notation allows for a very compact and convenient mathematical treatment of quantum theory.
It is built around the concept of a Hilbert space, in which a certain linear vector space is combined with a specific inner product. This space is spanned by so-called kets:
which is the same thing as a vector
A feature of the kets that many students may find unfamiliar is that they can be constructed from functions.
Where you would normally have Euclidean vectors such as:
a = bk.ket([1.0, 2], name = "a") # <- just like np.array([1.0, 2])
b = bk.ket([0, -2], name = "b")
bk.show(a,b) # <- visualize
the kets allows you to define something like
x = sp.symbols("x") # <- symbolic 'x' from sympy
a = bk.ket( x*sp.exp(-x**2), name = "a" )
b = bk.ket( sp.exp(-2*x**2), name = "b" )
bk.show(a,b)
Of course, the functions can be of arbitrary dimensionality, such as two:
x,y = bk.get_default_variables(1,2)
a = bk.ket( x*sp.exp(-.2*(x**2 + y**2) ), name = "a") # A 2D ket
bk.show(a)
...three...
x,y,z = sp.symbols("x y z")
a = bk.ket( x*sp.exp(-.01*(x**2 + y**2 + z**2)), name = "a") # A 3D ket
bk.show(a) #visualize with https://github.com/K3D-tools/K3D-jupyter
/opt/anaconda3/lib/python3.7/site-packages/traittypes/traittypes.py:101: UserWarning: Given trait value dtype "float64" does not match required type "float32". A coerced copy has been created.
np.dtype(self.dtype).name))
Output()
...or five (which is hard to visualize):
x,y,z,p,q = bk.get_default_variables(1,5)
a = bk.ket( sp.exp(-(2*x**2 + y**2 + z**2 + p**2 + q**2)), name = "5D")
a
\(\vert 5D \rangle\)
You may furthermore evaluate these functions in points:
a(.1,0,2,.1,.1)
\(\displaystyle 0.01759747241562341\)
...or over ranges and grids:
a(np.linspace(-1,1,100), 0,0,0,0)
array([0.13533528, 0.14660577, 0.15855578, 0.17120015, 0.18455134,
0.19861924, 0.21341081, 0.22892992, 0.24517699, 0.26214881,
0.27983825, 0.2982341 , 0.31732079, 0.33707829, 0.3574819 ,
0.37850215, 0.40010472, 0.42225034, 0.44489483, 0.4679891 ,
0.49147919, 0.51530642, 0.53940751, 0.56371481, 0.58815653,
0.61265702, 0.63713716, 0.66151466, 0.68570454, 0.7096196 ,
0.73317087, 0.75626817, 0.77882065, 0.8007374 , 0.82192803,
0.84230328, 0.86177563, 0.88025996, 0.89767412, 0.91393958,
0.92898197, 0.94273168, 0.9551244 , 0.9661016 , 0.97561098,
0.98360693, 0.99005084, 0.99491147, 0.99816514, 0.99979596,
0.99979596, 0.99816514, 0.99491147, 0.99005084, 0.98360693,
0.97561098, 0.9661016 , 0.9551244 , 0.94273168, 0.92898197,
0.91393958, 0.89767412, 0.88025996, 0.86177563, 0.84230328,
0.82192803, 0.8007374 , 0.77882065, 0.75626817, 0.73317087,
0.7096196 , 0.68570454, 0.66151466, 0.63713716, 0.61265702,
0.58815653, 0.56371481, 0.53940751, 0.51530642, 0.49147919,
0.4679891 , 0.44489483, 0.42225034, 0.40010472, 0.37850215,
0.3574819 , 0.33707829, 0.31732079, 0.2982341 , 0.27983825,
0.26214881, 0.24517699, 0.22892992, 0.21341081, 0.19861924,
0.18455134, 0.17120015, 0.15855578, 0.14660577, 0.13533528])
(which in this case corresponds to computing the integral \(\int_{\mathbb{R}^n} \delta(\mathbf{r}-\mathbf{r}') a(\mathbf{r}') d\mathbf{r}'\))
Bra-ket algebra
Just as you can do all sorts of algebra on Euclidean vectors:
a = bk.ket( [ 2,4], name = "a")
b = bk.ket( [-2,3], name = "b")
bk.show(a, b, a-b, a+b)
the same can be done to the more abstract kets:
a = bk.ket( y*sp.exp(-2*(x**2 + y**2))/(1+.1*x**2) )
b = bk.ket( .1*(y**2+x**2)*sp.sin(x**2)*sp.exp(-.5*y**2 -.5*x**2))
bk.show(.1*a+2*b)
Students are likely familiar with the inner product in the context of Euclidean vectors, the dot-product:
a = bk.ket([1,0])
b = bk.ket([0,1])
a.bra@b # <- dot product (these are orthogonal so should be zero)
\(\displaystyle 0.0\)
...but may not have seen this way of expressing it before. For any ket:
a = bk.ket(y*sp.exp(-2*(x**2 + y**2))/(1+.1*x**2), name = "a")
a
\(\vert a \rangle\)
we may define a bra:
a.bra
\(\langle a \vert\)
Whenever a bra runs into a ket on its right, they form a so called bra-ket - the inner product:
a.bra@a
\(\displaystyle 0.09571343617285585\)
The inner product is thus a mapping of two kets onto a complex number:
Hilbert space is spanned by all square normalizeable functions in \(\mathbb{R}^3\),
provided the inner product
Orthonormality - normality and orthogonality
Quantum mechanics deals with distributions evolving in time according to a wave equation. As would be the case with most waves, the distributions will then at some point oscillate between positive and negative values, thus calling for a different probability interpretation than the one we normally use.
Rather than having the distribution \(P\) normalized to one over \(\mathbb{R}^3\):
we require our wavefunction, represented by \(\vert \Psi \rangle\), to be square-normalized to one, meaning that
Using BraketLab, we satistfy this condition by having
in the following way:
psi = bk.ket( 2*sp.cos(2*x) * sp.exp(-.2*x**2), name = "$\\psi$" ) # <- some ket in Hilbert space
psi_normalized = (psi.bra@psi)**-.5*psi #normalization
psi_normalized.__name__ = "$\\psi_n$"
bk.show(psi, psi_normalized)
print("<psi | psi > = ", psi.bra@psi)
print("<psi_n | psi_n> = ", psi_normalized.bra@psi_normalized)
<psi | psi > = 5.605245682605472
<psi_n | psi_n> = 1.0
Being orthogonal to another function does not neccessarily mean that they do not have regions where they are both non-zero:
psi_a = bk.ket( 5*sp.cos(2*x) * sp.exp(-.1*x**2), name = "$\\psi_a$" )
psi_b = bk.ket( 5*sp.sin(1*x) * sp.exp(-.1*x**2), name = "$\\psi_b$" )
bk.show(psi_a, psi_b)
print("<psi_a | psi_b > = ", psi_a.bra@psi_b, "(orthogonal)")
<psi_a | psi_b > = 0.0 (orthogonal)
yet, we call integrals like the one above overlap integrals. (Students may be shocked to learn that quantum chemistry is more often than not carried out in a non-orthogonal basis where these are non-zero.)
Behind the scenes
Inner products may be numerically estimated using a range of methods depending on the basis and the dimensionality.
For 1D integrals, BraketLab in general utilize quadrature methods. In higher dimensions, the integrals are computed using a combination of zero-variance Monte-Carlo methods [1], including (1) importance sampling, interpolated control variates and coordinate transforms. Repeated calculations of identical quantities is avoided by using the lru_cache decorator.
In the long run, BraketLab will have basis-specific integrators, analytical and numerical, with fallback to multipurpose numerical methods where required.
[1] Caffarel, M. (2019). Evaluating two-electron-repulsion integrals over arbitrary orbitals using zero variance Monte Carlo: Application to full configuration interaction calculations with Slater-type orbitals. The Journal of Chemical Physics, 151(6), 064101.
Outer products and operators
Kets can be combined in many ways. For instance, the (cartesian) product of two kets is a new ket in a space often called Liouville space:
x,y,z = bk.get_default_variables(1,3)
psi_a = bk.ket( sp.exp(-.1*x**2), name = "a")
psi_b = bk.ket( x*sp.exp(-.2*x**2), name = "b")
bk.show(psi_a, psi_b)
ab = psi_a@psi_b
ab
\(\vert ab \rangle\)
bk.show(ab)
ab.ket_sympy_expression
\(\displaystyle 1.0 x_{2; 0} e^{- 0.1 x_{1; 0}^{2}} e^{- 0.2 x_{2; 0}^{2}}\)
In addition to operations such as addition, subtraction, division and multiplication, kets may be acted upon by several other operators. The students encounter such an operator if they happen to multiply a ket with a bra to its right, the so-called mapping operator:
a = bk.ket( sp.exp(-4*x**2), name = "a")
b = bk.ket( x*sp.exp(-6*x**2), name = "b")
a@b.bra
The trace of the mapping operator is called a projection operator:
bk.trace(a@b.bra)
Operators in general act on kets, resulting in new kets:
Many other operators are of interest in quantum mechanics, such as translation, differentiation, multiplication by a variable.
Certain operators called Hermitian operators represent observable quantities in quantum mechanics. Operators are thus a cruical part of the algebra og quantum theory. Several operators are available in BraketLab, let's look at some of them:
psi = bk.ket( sp.exp(-4*(x+3)**2))
T = bk.get_translation_operator(np.array([2.1]))
# Apply the translation operator one, two and three times
Tpsi = T*psi
TTpsi = T*Tpsi
TTTpsi = T*TTpsi
bk.show(psi, Tpsi, TTpsi, TTTpsi)
a = bk.ket( sp.exp(-x**2), name = "a(x)")
D = bk.get_differential_operator(order = [1])
Da = D*a
Da.__name__ = "$\\frac{d}{dx} a(x)$"
bk.show(a, Da)
It generalizes to many dimensions:
a = bk.ket( x*sp.exp(-(x**2 + y**2)**.5))
D = bk.get_differential_operator(order = [1,1])
D2a = D*a
bk.show(a)
bk.show(D2a)
The Hamiltonian operator of the electronic Hydrogen-equation is the Hermitian operator representing the energy of the Hydrogen system. It consists of the potential energy term:
V = bk.onebody_coulomb_operator()
V
$ -\frac{1}{\vert \mathbf{r} - (0.000000, 0.000000, 0.000000) \vert }$
and the kinetic term
T = bk.kinetic_operator()
T
$ -\frac{1}{2} \nabla^2 $
Together, they describe the total energy of the electron in the Hydrogen atom, which in textbooks is introduced as the time-independent Schrödinger equation. Operators, like matrices, have eigenfunctions and eigenvalues. For Hermitian operators, the eigenfunctions and eigenvalues constitute the possible outcomes of measurements.
The eigenfunctions of the Hydrogen-atom are the so-called spherical harmonics multiplied with a radial factor:
n,l,m = 1,0,0
psi_0 = bk.basisbank.get_hydrogen_function(n,l,m)
psi_0.get_ket_sympy_expression()
\(\displaystyle \frac{1.0 e^{- \sqrt{x_{0; 0}^{2} + x_{0; 1}^{2} + x_{0; 2}^{2}}}}{\sqrt{\pi}}\)
We can of course verify that they are eigenfunctions, in which case they should be proportional to eachother:
x_,y_,z_ = np.random.uniform(-1,1,(3,30))
Hpsi_0 = V*psi_0 + T*psi_0
Hpsi_0(x_,y_,z_)/psi_0(x_,y_,z_)
array([-0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
-0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
-0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5])
...and the factor of proportionality is the eigenvalue. The eigenvalue is in this case the ground state energy, confirmed by computing $\langle H \rangle = \langle \psi_0 \vert \hat{H} \vert \psi_0 \rangle $
print("E_0 = %.2f Hartrees" % (psi_0.bra@(V*psi_0) + psi_0.bra@(T*psi_0)))
E_0 = -0.50 Hartrees
Other functions do not satisfy this proportionality for Hydrogen, albeit the energy may not be too far off:
x,y,z = sp.symbols("x y z")
psi_trial = bk.ket( sp.exp(-.5*(x**2 +y**2 + z**2)) )
psi_trial = (psi_trial.bra@psi_trial)**-.5*psi_trial #normalize
print("E_0 = %.2f Hartrees" % (psi_trial.bra@(V*psi_0) + psi_trial.bra@(T*psi_0)))
Hpsi_trial = V*psi_trial + T*psi_trial
Hpsi_trial(x_,y_,z_)/psi_trial(x_,y_,z_)
E_0 = -0.47 Hartrees
array([-0.89210527, -0.1883901 , -0.00925463, -0.01432623, -0.04769066,
-0.18410504, -0.15937876, -0.01014062, -0.01446683, -1.57452854,
-0.01552328, -0.15282104, -0.00246156, -0.00742671, -0.07203871,
-0.01400314, -0.00217867, -0.02712906, -0.01107158, -0.01054015,
-0.02023079, -0.01571116, -0.01014115, -0.01996933, -0.09101338,
-0.16793262, -0.1185765 , -0.95278777, -0.26525749, -0.31099361])
Behind the scenes
BraketLab is jointly powered by sympy and numpy. Analytical derivation is easily performed with sympy, whereby integrals are computed numerically. In this sense, BraketLab is simply a wrapping of these two modules into a framework very similar to what you would find in a book on quantum theory.
Stationary states
Eigenfunctions of system-specific Hamiltonians are generally called stationary states, and while their mathematical formulations may be intimidating to students, we mostly look them up in tables (or generate them on the fly) when we need them. For Hydrogen, these are:
Nt = 200
t = np.linspace(-1,1,Nt)*25
n,l,m = 3,2,0 #quantum numbers
psi_nlm = bk.basisbank.get_hydrogen_function(n,l,m)
px = psi_nlm(t[:, None],t[None, :], 0)
p = (np.conjugate(px)*px).real
plt.figure(figsize = (3,3))
plt.imshow(p)
plt.title("n,l,m=%i,%i,%i" % (n,l,m))
plt.show()
Superposition
With a complete basis, we may resolve the identity of any function in this basis by projection.
For instance, the quantum harmonic oscillator has eigenfunctions in the shape of sinusoids enveloped by gaussians:
bk.basisbank.get_harmonic_oscillator_function(0)
\(\vert 0 \rangle\)
bk.show( bk.basisbank.get_harmonic_oscillator_function(0, position = [0]),
bk.basisbank.get_harmonic_oscillator_function(1),
bk.basisbank.get_harmonic_oscillator_function(2),
bk.basisbank.get_harmonic_oscillator_function(3),
)
If we conjure up another function, not present in this basis and slightly offset from the basis centre:
x = sp.symbols("x")
psi_sto = bk.ket( sp.exp( - sp.sqrt((x-1)**2) ) , name = "sto")
psi_sto = (psi_sto.bra@psi_sto)**-.5 * psi_sto # normalize
bk.show(psi_sto)
...the Harmonic oscillator eigenbasis is complete, meaning that with sufficiently many functions, any other L2-function can be represented as a linear combination / superposition of states by projection:
p = bk.basisbank.get_harmonic_oscillator_function(0)
for i in range(1,15):
p += bk.basisbank.get_harmonic_oscillator_function(i)
P = bk.trace(p@p.bra) # <- form a projection operator
psi_ho = P*psi_sto #<- project psi_sto onto the HO-basis
bk.show(psi_sto, psi_ho) #<- compare the original to the projection
(Notice, however, how poorly the gaussians manage to capture the pointy peak (cusp) of the exponential function. This is a well-known issue which has stuck with quantum chemistry ever since the 1950s when Boys introduced Gaussian type orbitals in place of Slater type orbitals, which in turn had been introduced in the 1930s in order to avoid complications with the radial part of the hydrogen-solution.)
Time evolution
As time progresses, the function will evolve as a wave in the potential of the harmonic oscillator. With the known eigenbasis and eigenvalues of the Harmonic Oscillator, we may now watch what happens as time starts to flow:
psi_ho.run(dt = 0.02)
What happens here is that the full solution evolves according to:
as we expect for a time-independent potential.
Measurement
Measuring stuff in quantum mechanics consists of projecting the wavefunction into the basis of the eigenfunctions of the Hermitian operator being observed:
where the square of the coefficients \(\vert c_i \vert^2\) are interpreted as the probability of finding the system in the corresponding state.
BraketLab (do not yet fully) support measurement of the kets, for instance position:
ms = psi_ho.measure(repetitions = 1000)
ns, bins = np.histogram(ms)
bins = .5*(bins[1:] + bins[:-1])
plt.figure()
plt.plot(bins, ns/np.sum(ns), ".")
plt.show()
Behind the scenes
A given Hermitian operator may be measured on a ket by providing it's eigenbasis to the measurement-method of the ket. Note that this is not fully implemented yet. Currently, when calling the measurement method, a Metropolis-Hastings random walk is performed in the probability distribution of choice (discrete in the case of a finite eigenbasis). By default, a measurement corresponds to measuring the position.
Atoms and molecules
From a quantum perspective, many-electron atoms may be seen as pertubed version of the Hydrogen atom. More electrons entering the system furthermore shifts us away from the Hydrogen-like solutions. Beyond atoms, we may consider molecules as perturbed versions of atoms. This motivates us to explore many-body quantum physics through the lens of independent particles, combined with the variational principle, stating that any trial wavefunction yields an upper estimate of the ground state energy of the system in question:
Just keep in mind that the further away from the the hydrogen-case we move, the less reliable estimate is provided from the hydrogenlike, single-particle picture.
Helium is a famous case, which was close to perfectly solved for by Egil Hylleraas at the golden age of quantum physics almost a hundred years ago. The problem illustrates the many facets of the role of electronic correlation in quantum chemistry (see for instance Tew, D. P., Klopper, W., & Helgaker, T. (2007). Electron correlation: The many‐body problem at the heart of chemistry. Journal of computational chemistry, 28(8), 1307-1320.)
The Hamiltonian of Helium is (in atomic units, i.e. \(\hbar = 4 \pi \epsilon_0 = m_e = q_e = 1\)):
The ground state energy has been experimentally determined to be about -2.9 Hartrees.
In the spirit of our current discussion, let's check how close we may get with a product of the Hydrogen-functions:
If we neglect the two-body term, we can compute the remaining contributions with BraketLab in the following way:
psi_0 = bk.basisbank.get_hydrogen_function(1,0,0)
E = 0
# both contributions to the kinetic energy
E_kin = psi_0.bra@(bk.kinetic_operator()*psi_0) * 2
# both contributions to the potential energy from the core
Z = 2
E_nuc = psi_0.bra@(bk.onebody_coulomb_operator()*psi_0) * Z * 2
print("E_0 < %.2f Hartrees" % (E_kin + E_nuc), " ( with inter-electronic repulsion neglected)")
E_0 < -2.99 Hartrees ( with inter-electronic repulsion neglected)
Not surprisingly, our energy is below the experimental result, because we neglected the repulsion between the electrons. It is actually twice what we would find for a ionized Helium atom with only one electron. Let's now add the inter-electronic contribution (but be a bit patient, this is a 6-dimensional integral):
psi_he = psi_0@psi_0 #the product wavefunction
E_vv = psi_he.bra@(bk.twobody_coulomb_operator()*psi_he)
print("E_0 < %.2f Hartrees" % (E_kin + E_nuc + E_vv), " ( with inter-electronic repulsion included)")
E_0 < -2.36 Hartrees ( with inter-electronic repulsion included)
This is not so bad, but Hylleraas did much better with a wavefunction on the form
# Hylleraas-type function
x1,x2,x3 = bk.get_default_variables(0,3) #variables for electron 0
x4,x5,x6 = bk.get_default_variables(1,3) #variables for electron 1
s = (x1**2 + x2**2 + x3**2)**.5 + (x4**2 + x5**2 + x6**2)**.5
t = (x1**2 + x2**2 + x3**2)**.5 - (x4**2 + x5**2 + x6**2)**.5
u = ((x1-x4)**2 + (x2-x5)**2 + (x3-x6)**2)**.5
k = 3.5
beta = 0.35
gamma = 0.15
delta = -0.13
epsilon = 0.01
xi = -0.07
Psi_hylleraas = bk.ket( sp.exp(-.5*k*s)*(1 + beta*u) )
# Add more terms if you like: (1 + beta*u + gamma*t**2 + delta*s + epsilon*s**2 + xi * u**2) )
# normalize
Psi_hylleraas = (Psi_hylleraas.bra@Psi_hylleraas)**-.5 * Psi_hylleraas
# print wavefunction
Psi_hylleraas.ket_sympy_expression
\(\displaystyle 1.17379724508579 \left(0.35 \left(\left(x_{0; 0} - x_{1; 0}\right)^{2} + \left(x_{0; 1} - x_{1; 1}\right)^{2} + \left(x_{0; 2} - x_{1; 2}\right)^{2}\right)^{0.5} + 1\right) e^{- 1.75 \left(x_{0; 0}^{2} + x_{0; 1}^{2} + x_{0; 2}^{2}\right)^{0.5} - 1.75 \left(x_{1; 0}^{2} + x_{1; 1}^{2} + x_{1; 2}^{2}\right)^{0.5}}\)
T0 = bk.get_kinetic_operator(p=0) # kinetic operator for electron 0
T1 = bk.get_kinetic_operator(p=1) # ...and 1
V0 = bk.get_onebody_coulomb_operator(Z=2, p = 0) # electron-nuclear interaction for electron 0
V0.ops[0][0].position = np.zeros(3)
V1 = bk.get_onebody_coulomb_operator(Z=2, p = 1) # electron-nuclear interaction for electron 1
V1.ops[0][0].position = np.zeros(3)
V_ee = bk.operator_expression(bk.twobody_coulomb_operator(0,1)) # electron-electron interaction, 0-1
H = T0 + T1 + V0 + V1 + V_ee # The full Helium-Hamiltonian
# evaluate variance in local energy and energy
xx = np.random.uniform(-1,1,(6,10000))
print("Std in energy:", np.std((H*Psi_hylleraas)(*xx)/Psi_hylleraas(*xx)))
print("< H > = ", Psi_hylleraas.bra@(H*Psi_hylleraas))
Std in energy: 0.24383023133704104
< H > = -2.910636492069882
BraketLab can deal with the electron repulsion integrals implicitly as above, or explicitly as follows:
bk.eri_mci(psi_0, psi_0, psi_0, psi_0)
\(\displaystyle 0.623744497031888\)
In turn, this facilitates a very modular approach to implementing traditional quantum chemistry methods such as Hartree-Fock or DFT (if you have the time to wait for the integrals). It furthermore allows for thinking outside of the box of independent particles, gaussian type orbitals and slater determinants, which is why I have called BraketLab out-of-the-box outside-of-the-box thinking.
It's out-of-the-box since you may simply do
# pip install braketlab
...to get it on your laptop. The outside-of-the-box-thinking-part is up to the students, but there are not many confines here (*) so they are free to explore. I'd love to get your feedback and proposals for features, just keep in mind that this is a work in progress. Thank you for your attention!
(*) Yet, many bugs and missing features.
(By the way, also check out BubbleBox for learning molecular dynamics and thermodynamics)