PATH INTEGRAL COMPUTATIONAL METHODS
FOR CONDENSED PHASE SYSTEMS: FROM
CLASSICAL TO QUANTUM
A Dissertation
Presented to the Faculty of the Graduate School
of Cornell University
in Partial Fulfillment of the Requirements for the Degree of
Doctor of Philosophy
by
Elliot Christopher Eklund
December 2022
© 2022 Elliot Christopher Eklund
ALL RIGHTS RESERVED
PATH INTEGRAL COMPUTATIONAL METHODS FOR CONDENSED
PHASE SYSTEMS: FROM CLASSICAL TO QUANTUM
Elliot Christopher Eklund, Ph.D.
Cornell University 2022
The path integral formulation of quantum mechanics provides a powerful and
vast set of tools for physicists and chemists to compute thermodynamic and
dynamic quantities of condensed phase systems [1]. Unfortunately, path inte-
gral calculations for large systems have high computational cost and in some
cases converge poorly owing to the infamous sign problem and other numeri-
cal instabilities. Developing new methods to more efficiently perform path in-
tegral calculations in the condensed phase and exploring their numerical stabil-
ity and convergence properties will be useful to many scientists who practice
these methods. In this dissertation we discuss three methods to boost the per-
formance of various path integral calculations, two of which rely on classical
computing while the third utilizes quantum computing.
First, we demonstrate how the ring polymer molecular dynamics (RPMD)
framework for condensed phase systems can be made more efficient by using
mixed time slicing (MTS) to reduce the Trotter dimension of the sub-system
describing the environment [2]. In doing so, we derive a consistent method
for introducing fictitious mass within the RPMD framework for a Hamiltonian
containing multiple sub-systems.
Next, we study mapping variable ring polymer molecular dynamics (MV-
RPMD), a path integral method for non-adiabatic dynamics, and apply MTS
to reduce the Trotter dimension of electronic or nuclear variables. With the
ability to selectively reduce the dimensionality of the MV-RPMD Hamilto-
nian, we study the ability to compute position auto-correlation functions. We
also undertake the first study of the emergence of stable and unstable mani-
folds of a Hamiltonian describing non-adiabatic dynamics using the method
of Lagrangian descriptors[3, 4]. All MV-RPMD calculations are executed using
MAVARIC, a software package developed within the Ananth group. A tutorial
for using MAVARIC is provided toward the end of this dissertation.
Finally, we develop a hybrid quantum-classical algorithm for running path
integral Monte Carlo (PIMC) calculations in which complex time propagator
matrix elements are computed in a position basis by calling a quantum com-
puter. This method can be extended to system-bath models via the use of an in-
fluence functional. We compute the thermal position auto-correlation function
of a simple harmonic oscillator using quantum simulators as proof-of-principle
and compare against exact results. Additionally, we show how our hybrid algo-
rithm can be modified for studying the dynamics of spin systems.
BIOGRAPHICAL SKETCH
Elliot grew up in central Wisconsin in a rural farming town called Manawa, the
self proclaimed Rodeo Capital of the Midwest. As a teenager he spent much of
his free time playing his Fender Stratocaster and getting into trouble with the
neighbor kid, Dylan. Unfortunately, Elliot’s childhood dreams of becoming a
professional basketball player were literally cut short when he grew to a final
height of 5’8”. Accepting this shortcoming, he set his sights on the humble goal
of becoming a famous rock star. Eventually he would come to his senses and at-
tend the University of Wisconsin - Madison with the aspirations of becoming a
doctor of medicine. This pursuit would only last through freshman year though
as Elliot realized he could not stand the constant complaining from the insuf-
ferable pre-med students he was surrounded by. During sophomore year, his
friend David introduced him to Tianqi who was looking for an undergrad to do
labor for her in the Zanni group, and the rest was history. He spent the next
three years learning about ultra fast, two-dimensional infrared spectroscopy
techniques and applying these skills to study protein folding in alpha-B crys-
tallin. It was during his time in the Zanni group that Elliot discovered his love
for physical chemistry and decided to further his education by attending gradu-
ate school to complete a PhD. He attended Cornell University upon graduation
and spent the next six years trying to do useful things with path integrals. Dur-
ing his free time he still enjoys playing guitar as well as talking to his friends,
Jordan and Devika, and playing with his dog Buddy.
iii
For my dog, Buddy:
”Woof, woof”
iv
ACKNOWLEDGEMENTS
A wise man once told me the most important part of a dissertation is the ac-
knowledgements, so I’ll do my best to make this section count.
My academic journey, and general existence, would not have been possible
without my parents, Dennis and Tracy. As a music teacher, my mother un-
derstood the value of a musical education and enrolled my brothers and I in
music lessons from a young age. She has always encouraged my creative and
academic endeavors. My father spent 20 years working as a mail carrier for
the United States Postal Service. He would wake up at 5 AM every weekday,
commute an hour to and from work, and spend all day running from house to
house delivering mail, but he always had enough energy to play basketball or
catch with us when he got home. Thank you both for giving me the freedom
and support to study and pursue the things I enjoy.
Going to high school in rural, small town Wisconsin leaves much to be de-
sired, but there are a few spectacular STEM teachers that deserve thanks. Carey
Celske who taught biology, human anatomy, and any life science the school de-
cided to throw his way. Thank you for providing a rigorous education and not
watering down material - it paid off in the long run. Carmen O’Brien who was
my first of many chemistry teachers. Thank you for offering Advanced Chem-
istry my senior year even though I was literally the only student enrolled. I
will never forget the ad-hoc electrochemistry lab that went horribly wrong and
ended up frying your giant battery. Jackie Hanson who taught AP Calculus.
Thank you for always being enthusiastic and entertaining while teaching math
and allowing me to ask you questions relentlessly throughout the school day.
During undergrad I had the incredible fortunte of winding up in the Zanni
group where I did research for three years. Here, I made lifelong friends and re-
v
alized that I wanted to enter a PhD program after graduation. Thank you Tianqi
Zhang, my grad student mentor, for being patient with me and teaching me tons
of skills for protein expression and purification. Additionally, thank you for not
killing me when I decided to switch projects and abandon biochemistry all to-
gether! Thank you Marty for all your words of wisdom and embodying the
work hard, play hard mentality. Thank you to everyone in the group who was
there while I was an undergrad: Jessi, Ariel, Andy, Josh, Short Kacie, Justin,
Nick, Michal, Jia-Jung, Tom, Arnaldo, Kaarin, Ayanjeet, Huong, Randy, Jack,
Tracey.
I was very lucky to be surrounded by amazing scientists while at Cornell.
In particular, I am very thankful for the members of my graduate committee.
Thank you Professor John Marohn for being the leading authority on all things
spin-system related and engaging conversations about the banjo. Thank you
Professor Greg Ezra for lending me several books through the years (all of which
I hope have been returned), providing insightful discussion regarding some of
the more mathematically complex problems I’ve come across, and conversa-
tions about drumming. I would also like to thank Professors Roger Loring and
Rob DiStasio for insightful discussion during journal club and around the de-
partment.
Thank you to everyone in the Ananth group who I overlapped with during
my time at Cornell: Britta, Srinath, Nathan, Ken, Shreyas, Tim, Matt, Sadrach,
Eric, Jess, Sujata, Siyu, Sutanuka, Rachel, and Pallavi. You have all helped make
this experience as painless as possible! I’ve learned a lot from all of you through
conversations at the white board, discussions during group meetings, and lis-
tening to journal club presentations. I will always cherish group lunches, going
to The Big Red Bardn, rock climbing, and bringing Buddy into the office to bark
vi
at Ken.
Finally, to my graduate advisor Professor Nandini Ananth, thank you for
everything over the past six years. The first time I met Nandini was after a talk
she gave at Wisconsin. I remember being excited by the talk, not because of any
particular thing she said, but from all the equations and math floating around
that I did not yet understand. The idea that I could pretend to be a chemist, but
do math all day was really what led me to Cornell. Thank you for always being
patient and understanding, especially when life outside of lab was rocky. I have
always felt supported and unafraid to come to you for advice on any number of
issues. Thank you for treating us grad students with kindness and respect.
The work in this dissertation was funded by National Science Foundation
EAGER (QAC-QSA), grant number 2038005, and National Science Foundation
CAREER, grant number 1555205.
vii
TABLE OF CONTENTS
Biographical Sketch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v
Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi
1 Introduction 1
2 Ring Polymer Molecular Dynamics 5
2.1 Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Condensed Phase RPMD . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3 Mixed Time Slicing . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4 Handling 1/N . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3 Mapping Variable - Ring Polymer Molecular Dynamics 24
3.1 Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2 Mixed Time Slicing . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.3 Constants of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.4 Model Systems and Numerical Simulation Details . . . . . . . . . 35
3.4.1 Model Systems . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.4.2 Nuclear Position Auto-Correlation Function . . . . . . . . 36
3.4.3 Analysis of Auto-Correlation Functions . . . . . . . . . . . 38
3.5 Lagrangian Descriptors . . . . . . . . . . . . . . . . . . . . . . . . . 41
4 Quantum Computing with Path Integrals: Continuous Systems 47
4.1 Quantum Computing Basics . . . . . . . . . . . . . . . . . . . . . . 49
4.2 Correlation Functions with Path Integrals . . . . . . . . . . . . . . 53
4.3 Hadamard Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.4 Probabalistic Imaginary Time Evolution . . . . . . . . . . . . . . . 58
4.5 Full Circuit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.6 Potential Energy Circuit Design . . . . . . . . . . . . . . . . . . . . 63
4.7 Kinetic Energy Circuit Design . . . . . . . . . . . . . . . . . . . . . 65
4.8 Harmonic Oscillator Data . . . . . . . . . . . . . . . . . . . . . . . 67
4.9 Simulation Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.9.1 Partition Function Trick . . . . . . . . . . . . . . . . . . . . 69
4.9.2 Wavefunction Normalization . . . . . . . . . . . . . . . . . 70
4.9.3 Hashing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
viii
5 Quantum Computing with Spin Systems 72
5.1 Spin Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.2 Hybrid Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.3 Matrix Elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6 MAVARIC 83
6.1 Tutorial Model System and Getting Started . . . . . . . . . . . . . 85
6.1.1 Tutorial Model System . . . . . . . . . . . . . . . . . . . . . 85
6.1.2 Input Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
6.1.3 Compiling and Running MAVARIC . . . . . . . . . . . . . 86
6.2 Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.3 Sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
6.4 Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
6.5 Coding Potentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
6.5.1 Macros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
6.5.2 Valarray . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
6.5.3 State Independent Potential . . . . . . . . . . . . . . . . . . 108
6.5.4 Gradient of State Independent Potential . . . . . . . . . . . 108
6.5.5 State Specific Potential Energy Surfaces . . . . . . . . . . . 109
6.5.6 Gradient of State Specific Potential Energy Surfaces . . . . 110
6.5.7 Potential Energy Coupling . . . . . . . . . . . . . . . . . . 111
6.5.8 Derivative of Potential Energy Coupling . . . . . . . . . . 112
6.5.9 Complicated Model . . . . . . . . . . . . . . . . . . . . . . 113
6.6 Derivative Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 118
A Approximate PITE Expansion 123
ix
LIST OF TABLES
3.1 Parameters for Models I-III provided in atomic units. For all
models, M = β = γ = 1 a.u. Simulations where we investigate
the effect of changing γ use values 1, 10, and 100 a.u. . . . . . . . 36
x
LIST OF FIGURES
2.1 Depiction of ring polymer consisting of 8 beads [22]. . . . . . . . 8
3.1 Depiction of MV-RPMD Hamiltonian [4]. Red beads in the cen-
ter represent electronic variables; grey beads on the outside rep-
resent nuclear beads. Squiggly lines connecting nuclear beads
represent harmonic spring interactions. Straight lines represent
interaction with electronic beads. . . . . . . . . . . . . . . . . . . . 28
3.2 (a) Cartoon depiction of HQN. Since Case 1 has Nn > Ne, only cer-
tain nuclear bead numbers have corresponding electronic bead
numbers. (b) Cartoon depiction of HQE. A single classical nuclear
bead interacts with multiple electronic beads [4]. . . . . . . . . . 32
3.3 Nuclear position auto-correlation function for Model I. The num-
ber of electronic and nuclear beads is indicated by the format
Ne : Nn and corresponds to the following lines: (a) 1:6=red tri-
angles. 3:6=green squares. 6:6=blue circles. Exact=black dashes.
(b) 1:1=blue circles. 4:1=magenta diamonds. Exact=black dashes
[4]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.4 Nuclear position auto-correlation function for Model II. The
number of electronic and nuclear beads is indicated by the for-
mat Ne : Nn and corresponds to the following lines: (a) 1:4=red
triangles. 2:4=green squares. 4:4=blue circles. Exact=black
dashes. (b) 1:1=blue circles. 4:1=magenta diamonds. Ex-
act=black dashes [4]. . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.5 Nuclear position auto-correlation function for Model III. The
number of electronic and nuclear beads is indicated by the for-
mat Ne : Nn and corresponds to the following lines: (a) 1:4=red
triangles. 2:4=green squares. 4:4=blue circles. Exact=black
dashes. (b) 1:1=blue circles. 4:1=magenta diamonds. Ex-
act=black dashes [4]. . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.6 (a) LD calculated over a grid in nuclear phase space shows a sta-
ble manifold (Ws) and an unstable manifold (Wu). The crossing
point lies at P = 0 and Q = −0.13. (b) Cross section at P = −0.2
showing sharp changes in the LD corresponding to the location
of stable and unstable manifolds. (c) We calculate the semiclassi-
cal population estimator using the electronic variables obtained
through the SGD simulation and find that the LDs shown here
correspond to state |1⟩ being mostly occupied and (d) state |2⟩
being mostly unoccupied [4]. . . . . . . . . . . . . . . . . . . . . . 44
xi
3.7 (a) LD calculated over a grid in nuclear phase space shows a sta-
ble manifold (Ws) and an unstable manifold (Wu). The crossing
point lies at P = 0 and Q = 0.2. (b) Cross section at P = 0.2
showing sharp changes in the LD corresponding to the location
of stable and unstable manifolds. (c) We calculate the semiclassi-
cal population estimator using the electronic variables obtained
through SGD simulation and find that the LDs shown here cor-
respond to |2⟩ being mostly occupied and (d) |1⟩ being mostly
unoccupied [4]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.8 (a) Nuclear phase space plot showing HMV where we choose the
electronic variables such that only |1⟩ is occupied. (b) Nuclear
phase space plot showing HMV when only |2⟩ is occupied [4]. . . 46
4.1 Quantum circuit diagram illustrating single qubit passing
through Hadamard gate. Measurement occurs at the end. . . . . 51
4.2 Quantum circuit illustrating two qubits passing through
Hadamard and CNOT gates. Both qubits are measured at the
end. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.3 Hadamard Test circuit diagram. . . . . . . . . . . . . . . . . . . . 56
4.4 Hadamard Test circuit with real and imaginary time evolution
gates explicitly shown. . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.5 Exact PITE circuit diagram. . . . . . . . . . . . . . . . . . . . . . . 59
4.6 Approximate PITE circuit diagram. . . . . . . . . . . . . . . . . . 61
4.7 Full circuit for position basis complex time evolution. Here,
CPIT E is the approximate PITE circuit. . . . . . . . . . . . . . . . . 62
4.8 Approximate PITE circuit with additional ancilla. . . . . . . . . . 62
4.9 Pixel function centered over origin. . . . . . . . . . . . . . . . . . 66
4.10 Harmonic oscillator position auto-correlation function. Here,
N = 3 for the first give points, and N = 4 for the remaining six. . . 68
5.1 XX-chain real time evolution circuit. . . . . . . . . . . . . . . . . . 79
5.2 Spin expectation value for XX chain with two spins. . . . . . . . . 81
5.3 Spin expectation value for XX chain with four spins. . . . . . . . 82
6.1 Average Energy throughout a short Monte Carlo simulation. . . 90
6.2 Monte Carlo simulation start from random initial state. . . . . . . 91
6.3 Average Energy throughout a short Monte Carlo simulation. . . 93
6.4 Monte Carlo simulation starting from saved initial state. . . . . . 94
6.5 Histogram using 300 bins with 105 Trajectories. . . . . . . . . . . 98
6.6 Position Centroid auto-correlation function using 104 trajectories. 103
xii
CHAPTER 1
INTRODUCTION
Understanding how charge and energy transfer occurs in condensed phase sys-
tems is a challenging problem with many important applications [5, 6, 7, 8]. For
example, photosynthesis in plants and photosynthetic bacteria relies on a com-
plicated series of charge transfer and energy transfer reactions [9]. Studying
this process could aid in the future design of bioinspired solar cells and other
technologies. Understanding the mechanisms behind these processes requires a
detailed description of the system at the quantum mechanical level as it evolves
in time. Unfortunately, most time-resolved spectroscopic techniques fail to re-
solve these prcoesses because they occur on a sub-femtosecond time scale. The-
oretical techniques offer an alternate route to gaining insight into such ultrafast
processes. However, they are challenged by the fact that electronic and nuclear
degrees of freedom are often coupled during these reactions, preventing one
from using the Born-Oppenheimer approximation to simplify calculations.
Several methods have been developed to study such non-adiabatic dynam-
ics, each with its own limitations. Multi-Configurational Time Dependent
Hartree (MCTDH) is regarded as the state of the art and attempts to solve the
full time dependent Schrodinger equation head on [10, 11]. As a result, MCTDH
is very accurate, but severely limited to systems of only a few atoms. The need
to simulate larger systems has spurred the development of approximate quan-
tum methods which sacrifice accuracy in exchange for improved scaling cost
with system size. For example, semiclassical techniques work by averaging over
a large number of classical trajectories each weighted by its classical action and
a complex prefactor that accounts for paths close to classical trajectories. Un-
1
fortunately, these methods also struggle to handle large systems and are often
plagued by convergence issues owing the sign problem [12].
One class of methods derived from ring polymer molecular dynamics
(RPMD) [2, 13, 14, 15] has seen early success and could provide a path toward
simulating non-adiabatic dynamics of large systems. Like RPMD, these meth-
ods use classical equations of motion to approximate the Kubo-Transformed
correlation function, but extend to multiple potential energy surfaces [16, 3, 17,
18]. Because they use imaginary time dynamics to approximate real time dy-
namics, they avoid the sign the problem. Additionally, using classical equations
of motion allows them to enjoy linear scaling in computational cost with in-
creasing system size. One particular method known as mapping variable ring
polymer molecular dynamics (MV-RPMD) stands out as it has been shown to
conserve the quantum Boltzmann distribution. MV-RPMD has been success-
fully used to study photo-dissociation dynamics and reaction rates in proton
coupled electron transfer models.
Even with linear scaling offered by RPMD and MV-RPMD, there is always a
desire to improve computational efficiency while retaining accuracy. In this dis-
sertation we examine several new methods designed to achieve this goal. The
first optimization is based on applying mixed time slicing (MTS), in which dif-
ferent degrees of freedom are quantized independently, to RPMD. MTS allows
one to control the level of desired accuracy for different degrees of freedom
which in turn sets how much computational effort is used. We find this to be
particularly useful for system-bath models where we are primarily interested in
treating the system with high accuracy, and the harmonic oscillator bath modes
with lower accuracy. However, doing so requires careful re-scaling of the ficti-
2
tious mass term which we examine.
We also apply MTS to MV-RPMD, in which case it is now the nuclear and
electronic degrees of freedom which are quantized to different extents. This
allows us to study whether accurate treatment of the nuclear or electronic de-
grees of freedom is more important to the overall quality of a simulation. While
studying MV-RPMD, we explore the emergence of stable and unstable mani-
folds using the method of Lagrangian Descriptors. Due to the high dimensional
nature of the MV-RPMD Hamiltonian, we also develop a new non-deterministic
method for calculating Lagrangian Descriptors [4].
In recent years, there has been increasing interest in exploring how quantum
computers might be used to solve problems in chemistry [19, 20]. In this dis-
sertation we introduce a new algorithm which uses hybrid classical-quantum
computing to simulate the dynamics of a condensed phase system. Our algo-
rithm allows us to circumvent the most expensive part of the calculation, matrix
diagonalization, by calling on a quantum computer to directly compute position
matrix elements of a complex time evolution operator. As a result, simulations
of large systems that were previous too expensive for classical methods could
become feasible with our hybrid approach.
This dissertation is organized as follows. In Chapter 2, we review RPMD and
show how to correctly implement MTS for system-bath models. In Chapter 3,
we review MV-RPMD and show MTS can be used to treat nuclear and electronic
degrees of freedom on different footing. Chapter 4 is dedicated to the develop-
ment of a hybrid classical-quantum algorithm for path integral simulations of
condensed phase systems. We examine the success of this method by computing
the position auto-correlation function of harmonic oscillator. Chapter 5 contin-
3
ues Chapter 4, but focuses on spin-systems. Finally, Chapter 6 is a tutorial to
simulating MV-RPMD dynamis with the software package MAVARIC.
4
CHAPTER 2
RING POLYMER MOLECULAR DYNAMICS
In this chapter we introduce ring polymer molecular dynamics (RPMD) and
discuss how it can be used for dynamics in the condensed phase. We demon-
strate how to use mixed time slicing (MTS) to optimize the number of beads
used for system and bath degrees of freedom (DOF), ultimately reducing com-
putational overhead. We highlight that running RPMD with MTS for a system
coupled to a bath is not straightforward and in particular one needs to be careful
when handling the fictitious mass.
2.1 Review
Recall that the partition function for a system at thermal equilibrium is de-
scribed by the canonical ensemble and is given as
( )
Z = Tr e−βĤ (2.1)
where β = 1/kBT , T is temperature, and kB is Boltzmann’s constant. The operator
e−βĤ is referred to as the Boltzmann operator. We use a hat to indicate a quantum
operator, while the lack of a hat indicates a classical operator. The Hamiltonian
operator for a one-dimensional system is defined as
p̂2
Ĥ = + V(x̂) (2.2)
2m
where m is the mass of the system, p̂ and x̂ are momentum and position opera-
tors, respectively, and V is an arbitrary potential which only depends on posi-
tion. Though this review focuses on systems with one physical dimension, the
generalization to two and three dimensions is straightforward. To evaluate the
5
trace in Eq. (2.1) we work in the p∫osition basis
Z = dx ⟨x|e−β(T̂+V̂)|x⟩ (2.3)
2
where the kinetic energy operator is T̂ p̂= 2m and the potential energy operator is
V̂ = V(x̂). Unless explicit bounds are given, all integrals are assumed to be from
negative to positive infinity. Because the kinetic and potential energy operators
do not commute, we must employ the Trotter theorem [21] in order to factor the
exponential: ( )N
e−β(T̂+V̂)
β β
= lim e− N V̂e− N T̂ (2.4)
N→∞
This version of the theorem is referred to as asymmetric Trotter splitting and it
gives the same results as the symmetric version in cases where the trace is taken.
Applying this to Eq. (2.3) gives ∫ ( )
β β N
Z = lim dx ⟨x| e− N V̂e− N T̂ |x⟩ (2.5)
N→∞
β
Let Ω̂ = e− N V̂e−
β
N T̂ so that the integrand of Eq. (2.5) can be expressed as
⟨x|Ω̂ . . . Ω̂|x⟩, i.e. a position matrix element enclosing a string of N copies of Ω̂
N is referred to as the Trotter dimension and will eventually be the number of
“beads” comprising the ring polymer. If we∫insert N − 1 copies of the identity
operator expressed in the position basis, Î = |x⟩⟨x|, inbetween each Ω̂, then Eq.
(2.5) becomes ∫ ∏N
Z = lim dx ⟨x |Ω̂|x
∫ ∫ ∫ →∞ j j+1
⟩ (2.6)
N
j=1
where dx = dx1 . . . dxN , and xN+1 = x1. The fact that xN+1 = x1 is a result of
performing the trace in the position basis.
Our task now is to evaluate ⟨x j|Ω̂|x j+1⟩. Because the potential energy operator
β
is diagonal in position space,we can easily evaluate terms of the form e− 2 V̂ :
⟨x j|Ω̂|x ⟩ = e−βNV(x j)j+1 ⟨x j|e−βN T̂ |x j+1⟩ (2.7)
6
where βN = β/N. To evaluate the remaining matrix element we first insert iden-
tity in the momentum basis ∫to get
dp ⟨x |e−βN T̂j |p⟩⟨p|x j+1⟩ (2.8)
( ) 1 −i
The inner product between position and momentum bases is ⟨p|x⟩ 1 2= e ℏ xp2 .πℏ
Using this relation along with the fact that the kinetic energy operator is diago-
nal in momentum space leads to∫
1 − p2 idp e βN 2m e ℏ (x j−x j+1)p (2.9)
2πℏ
This integral can be evaluated with contour integration or completing the
square [22]. Doing so gives ( ) 1
m 2 − m 2x −x
e 2 ( j j+1)βNℏ2 (2.10)
2πβ ℏ2N
We arrive at our final expression for the partition function by combining Eqs.
(2.10), (2.7), and (2.6): ( ) N ∫
m 2 − ∑β NN j=1 (V(x j)− m 2Z lim dx e 2β2 2 (x j−x j+1) )= Nℏ
N→∞ 2πβ ℏ2
(2.11)
N
What we have succeeded to do thus far is map the partition function of a
quantum system onto an extended classical system consisting of N particles
coupled together through harmonic spring terms with spring constant k = mω2N ,
and ωN = 1 . This spring term appears as the second term in the exponentialβNℏ
of Eq. (2.11). The classical particles are still under the influence of the potential
V(x). A cartoon of the resulting structure, called a ring polymer, is shown in Fig
2.1. Each of the N classical particles is referred to as a “bead”. The delocalized
nature of the ring polymer is what allows it capture quantum nuclear effects de-
spite being only a collection of classical particles. Conventionally if one wanted
7
Figure 2.1: Depiction of ring polymer consisting of 8 beads [22].
to compute equilibrium thermodynamic quantities of a quantum system using
the partition function, the eigenvalues and vectors of the corresponding Hamil-
tonian would need to be found. The most straight forward method for obtain-
ing these quantities is by constructing the Hamiltonian matrix in some basis
set and diagonalizing the resulting matrix, a procedure which has O(n3) scaling
for a square matrix of dimensions n × n. The true ugliness of the scaling lies
in the fact that n itself scales with the size of the Hilbert space of the system.
For example, the full Hilbert space of a molecule with f electrons within the
Born-Oppenheimer approximation scales formally as 2 f . This means matrix di-
agonalization really scales exponentially with increasing system size. Although
in practice modern electronic structure methods offer sub-exponential scaling,
they still quickly become infeasible even for moderately complex systems.
With Eq. (2.11), one can instead using Monte Carlo sampling techniques to
compute the integral over an integrand which now depends only on classical
quantities. This technique is referred to as Path Integral Monte Carlo (PIMC)
[22] and has the advantage that the variance of an estimator being computed,
8
which is often used as a proxy for accuracy, decreases as O( √1 ) [22] up to some
s
constant, where s is the number of random samples (Monte Carlo steps) taken.
Often times sampling schemes can be designed using intuition of the system
so that the constant prefactor is relatively small. Additionally, for classical sys-
tems this prefactor is linear in systems size, n. In general, n = dN p, where d
is the number of physical dimensions, and p is the number of particles. Taken
together, this means PIMC offers a substantial reduction in computational over-
head compared to exact methods.
What would happen if Eq. (2.11) was used to compute dynamical quanti-
ties like correlation functions? This is the underlying question which led to the
development of RPMD. To understand why this is an interesting question we
need to further manipulate Eq. (2.11) and make contact with the expression for
a general co(rrelatio)n f(unctio)n. We start be writingN ∫ ∫ p2
m 2 N/2 ∑ 2β jN −β NN ( +V(x j)− m (x j−x j+1) )Z = lim dx dp e j=1 2m f 2β2 ℏ2N (2.12)
N→∞ 2πβ 2Nℏ 2πm f
where we have inserted N copies of identity in the form of a Gaussian integral
over momentum divided by itself. The new term in the exponential is just the
classical kinetic energy with mass m f , called the fictitious mass. This mass is said
to be fictitious because any choice of m f will yield correct result at equilibrium,
and m f does not need to be the same as m. We define the summation in the
exponential as the RPMD Hamiltonian∑N  p2 m ( ) 2
H  j RP(p, x) =  + V(x j) − 2 x j − x j+1  (2.13)2m 2
j 2β ℏ=1 N
In this definition we make the choice to set the fictitious mass equal to the true
mass of the system, m f = m. Using this definition and truncating the N → ∞
9
limit, we define the RPMD p(artitio)n∫functi∫on
1 N
Z = dx dp e−βN HRP(p,x)RP 2 (2.14)2πℏ
ZRP approximates Z and will give equality when N is taken to infinity. Typi-
cally ZRP, and any thermodynamic quantities derived from it, will converge to
a stable value for relatively small N, at which point further increasing N does
not produce a more accurate result. For example, accurately capturing quantum
effects for a proton at a temperature for 300K requires N = 32 [23].
To make contact the correlation function, recall that the thermal expectation
value of an observable which depends only on position is
( )
⟨Â⟩ 1= Tr e−βĤ Â (2.15)
Z
The corresponding path integral expression, which can be found by working
through Eqs. (2.1 - 2.14) and includ∫ing Â∫, is given as
⟨ ⟩ ≈ 1A dp dx e−βN HRP(p,x)A (x) (2.16)
∑ (2πℏ)NZ
N
RP
where A (x) = 1 NN N j=1 A(x j) [2]. If we repeat this derivation with the product of
two position dependent operato∫rs, we m∫ ight expect to find the result
1 ∞⟨AB⟩ ≈ dp dx e−βN HRP(p,x)N AN(x)BN(x) (2.17)(2πℏ) ZRP −∞
But this is not the case! Rather, Eq. (2.17) corresponds to the zero time limit of
the Kubo-transformed [24] c∫orrelation function:
1 β ( )
C̃ (t) = dλTr e−(β−λ)Ĥ Âe−λĤe+iĤt/ℏB̂e−iĤt/ℏAB (2.18)
βZ 0
One can show that the Fourier Transform of the Kubo-transformed correla-
tion function is related to the Fourier Transform of the standard correlation
function by a temperature dependent factor [2]. This realization led Craig
10
and Manolopoulos to wonder what would happen if this expression was used
to compute entire correlation functions despite equality only formally being
achieved in the t → 0 limit? The idea being that if Eq. (2.17) corresponds to
something meaningful when t → 0, then may it also has meaning when t is
non-zero. They proposed the∫correla∫tion function
⟨ ⟩ 1A(0)B(t) −βN HRP(p0,x0)RP = dp0 dx0 e × AN (x0) BN (xt)N (2.19)(2πℏ) ZRP
where x0 and xt are the position coordinates at t = 0 and some later time, respec-
tively. xt is found by propagating Hamilton’s equations of motion for a given
initial condition. Craig and Manolopoulos found that their correlation function
gave accurate results for short times when compared to exact calculations and
was a big improvement over classical results where quantum mechanical zero
point energy and tunneling are neglected. Since its inception RPMD has been
explored by many researchers and is now a standard tool for running Molecular
Dynamics (MD) simulations. For example, RPMD has been used for comput-
ing reaction rates, dipole absorption spectrum, diffusion coefficients, and many
other experimentally observable dynamic quantities [15].
2.2 Condensed Phase RPMD
In order to describe condensed phase systems RPMD needs to be extended to
allow energy dissipation from a system of interest to an environment. Although
there are many ways to model dissipation, we focus on a popular choice called
the Caldeira–Leggett (CL) model or harmonic bath model [25]. Here the envi-
ronment is modeled as a set of harmonic oscillators linearly coupled to a pri-
mary system through position operators. The Hamiltonian operator for the CL
11
model is
ĤCL = Ĥs + Ĥb + Ĥsb (2.20)
where Ĥs describes the system Hamiltonian,
P̂2
Ĥs = + V(x̂) (2.21)2M
Ĥb describes the collection of free ha(rmonic oscillators,∑F )p̂2
Ĥ i
1 2 2
b = + mbωi qi (2.22)2m 2
i=1 b
and Ĥsb describes the coupling term between system and bath,
∑F ∑F c2
Ĥsb = −x̂ ciq̂i + x̂2 i 2 (2.23)
i=1 i=1 2mbωi
In the Eqs. (2.21 - 2.23), x̂ and P̂ are system position and momentum vari-
ables, respectively, and q̂i and p̂i refer to the ith harmonic oscillator mode posi-
tion and momentum, respectively. F is the number of harmonic bath modes. ci
and ωi are the coupling strength and frequency, respectively, of the ith mode and
their values depend on the spectral density of the environment being modeled.
mb is the mass of each mode. Although each mode has the same mass in our
expression, it is possible to allow each mode to have a unique mass. A popular
model for the environment spectral density is the Debye spectral [26] density
which has the form
η γω
J(ω) = 2 2 (2.24)γ + ω
where ω is a given frequency, γ is the bath cutoff frequency, and η is the cou-
pling strength. Eq. (2.24) is a continuous function of frequency and needs to be
discretized since our model has a finite number, F, of bath modes. There are
several methods for discretizing spectral densities [27], all of which reproduce
12
the exact spectrum in the F → ∞ limit. We chose a discretization [27] that leads
to ( )
π j
ω j = γ tan (2.25)2 F + 1
for the frequencies, and √
mbη γc j = ω j (2.26)F + 1
for the coupling strengths throughout this chapter.
In order to get the corresponding RPMD Hamiltonian for the CL model we
need to repeat much of the math of Section 2.1. Our partition function is now
(
−β(Ĥ′
)
Z = Tr e s+Ĥb) (2.27)
where Ĥ′s = Ĥs + Ĥsb includes all terms involving system variables. Once again,
we evaluate the trace in the position basis, this time over system and bath posi-
tion variables: ∫ ∫
Z = dx dq ⟨x,q|e−β(Ĥ′s+Ĥb)∫ ∫ ∫ |x,q⟩ (2.28)
where q = (q1, . . . , qF), and dq = dq1 . . . dqF . Multiplying and dividing by
N and inserting N − 1 copies of identity in the position basis between each copy
of the Boltzmann operator gi∫ves ∫ ∏N
Z = lim dx d{q} ⟨x j,q j|Ω̂|x j,q→∞ j⟩ (2.29)N
∫ j=1 ∫ ∫ ∫
where Ω̂ = e−βN (Ĥ′s+Ĥb), dx is an N dimensional integral and {q} = q1 . . . qN
is an F × N dimensional integral. We use the notation q j = (q j,1, . . . , q j,F) where
q j,k is the kth bath mode of the jth bead. To evaluate the position basis matrix
element of Ω̂, we employ the Trotter theorem to write the Boltzmann operator
as a product:
⟨x ,q |Ω̂|x ,q ⟩ ≈ ⟨x ,q |e−β Ĥ′N sj j j j j j e−βN Ĥb |x j+1,q j+1⟩ (2.30)
13
Ĥ′s is diagonal in q and Ĥb does not depend on system position, so we can write
Eq. (2.30) as
⟨x |e−β ′N Ĥs(q j)j |x j+1⟩⟨q j|e−βN Ĥb |q j+1⟩ (2.31)
Both matrix elements are now in the same form as Eq. (2.7) and we can insert the
result from Eq. (2.11) to get our final expression. For the bath position matrix
element, we get ( ) F  ∑F ⟨ −β Ĥ m 2b 1 2 2 mq N b b 2 j|e |q j+1⟩ = exp 2 −βN mbωkq j,k − 2 (q2 2 2 2 j,k − q )πβ ℏ β ℏ j+1,k  (2.32)N k=1 N
where k is used to index a particular bath mode. The system position matrix
element is ( ) 1
2
⟨x j|
′
e−βN Ĥs(q j)|x j+1⟩
M
=  2 × (2.32πβNℏ ∑F ∑F 
3)
c M 
exp −β N V(x j) − x j ckq j,k − x2 k 2 j 2m ω2 − 2 (x j − x j+1)2 k=1 k 1 b k 2β= Nℏ
The resulting expressi(on for th)e pa(rtition f)un∫ction∫isF×N N
m 2 M 2
Z = lim b dx d{q}e−βNU(x,{q})2 2 (2.34)N→∞ 2πβNℏ 2πβNℏ
where
∑ U(x, {q}) = U
′
s(x, {q}) +∑Ub(({q}), ) (2.35)N  M F c 
U′(x, {q}) = V(x ) − (x − x 2 2 k s j 2 j j+1) + x jckq j,k − x j 2  (2.36)
j 1 2β ℏ
2 2m ω
= N k=1 b k
and ∑N ∑F [ ]1
U ({q}) = m ω2q2 − mbb (q − q 22 b k j,k 2 2 2 j,k j+1,k
) (2.37)
j 1 k 1 β ℏ= = N
In the above expressions the jth bead of the system ring polymer interacts
with every bath mode through the corresponding jth bead of the bath ring poly-
mer. Finally, we can introduce system and bath momentum in the same way as
14
earlier: ( )
1 N×(F+1)
∫ ∫ ∫ ∫
Z = dP dx d{p} d{q}e−βN HCL(x,P,{q},{p})CL 2 (2.38)2πℏ
where the RPMD expression for the CL Hamiltonian is∑N  P2 ∑Fj p2j,k HCL (x,P, {q}, {p}) = +2m 2m  + U(x, {q}) (2.39)
j=1 s k=1 k
With HCL we can run short time dynamic simulations of a quantum system in-
teracting with a bath of quantum harmonic oscillators for sufficiently large N.
For example, to compute a thermal correlation function we would use Eq. (2.19)
with HCL as our Hamiltonian.
One shortcoming of Eq. (2.39) is that both system and bath sub-systems use
the same number of beads. In general the number beads needed to achieve
convergence is related to how classical or quantum a system is. For example,
modeling an electron requires substantially more beads than for a proton since
the electron is more quantum in nature. This means that if we wish to model an
electron or some other quantum particle coupled to a relatively classical bath,
we are stuck using the same large N needed to converge the electron sub-system
for the bath sub-system as well. In the next section we show how to use Mixed
Time Slicing (MTS) [28] so that the system and bath sub-systems are discribed
by a different number of beads. This increased flexibility leads to a reduction in
the dimensionality of the overall system which in turn reduces computational
overhead.
15
2.3 Mixed Time Slicing
To derive the MTS form of the CL Hamiltonian we begin with the definition of
the partition function and t∫race in∫ the position basis:
Z = dx dq ⟨x,q|e−β(Ĥs′+Ĥb)|x,q⟩ (2.40)
For this derivation we assume that our system has more quantum behavior than
the bath, thus requiring a large Trotter dimension to accurately capture quantum
effects. Our first step is to multiply and divide by Nb in the exponential and split
the Hamiltonian into two pa[rts. (Focusin)g onl(y on the)]matrix element, we haveN
β β b
lim ⟨x,q| exp − Ĥ exp − Ĥ′ |x,q⟩ (2.41)
Nb→∞ N bb N sb
Next, we multiply and divide b(y Ns in )both(expone)ntialscontaining Ĥ′  s
:
lim ⟨x,q| Ns/Nβ β b Nbexp − Ĥb exp − Ĥ′ s  |x,q⟩ (2.42)Nb,Ns→∞ Nb Ns
At this point we can see how MTS allows us to use a different number of beads
for each sub-system. The term involving an exponential of the bath sub-system
Hamiltonian contributes Nb factors and will eventually lead to a ring polymer
with Nb beads. The exponentials containing the system Hamiltonian will result
in a ring polymer with Ns beads.
( ) ( )
To proceed we defined Ω̂b = exp − βN Ĥb , Ω̂s = exp −
β ′ f
b N
Ĥs , and Ω̂ = Ω̂s bΩ̂s
where f = Ns/Nb. This allows us to condense the position basis matrix element:
lim ⟨x,q|Ω̂Nb |x,q⟩ (2.43)
Nb,Ns→∞
We insert Nb − 1 copies of identity in the bath position basis between each Ω̂
operator and group all the int∫egrals over bath position to get
lim d{q} ⟨x|∏  Nb ⟨q j|Ω̂|q ⟩j+1 N ,N →∞  |x⟩ (2.44)b s
j=1
16
Using the fact that Ω̂s is d∫iagonal with respect to q gives∏ Nb 
lim d{q} ⟨x| ⟨q j|Ω̂b|q f →∞ j+1⟩ Ω̂s (q j+1) |x⟩ (2.45)Nb,Ns
j=1
We can evaluate each matrix element ⟨q j|Ω̂b|q j+1⟩ using Eq. (2.32) and pull the
result to the front of the integral:
( ) F×Nb ∫  
N m 2 ∏Nb 
lim b b2 d{q} e
−βUb({q})⟨x| Ω̂ f (q j+1)
N N →∞ 2πβℏ  s,  |x⟩ (2.46)b s
j=1
where ∑Nb ∑F [ 2 ]
UMTS
m N
b ({q}
1 1 b
) b= mbω2 2kq j,k − (q − q )2 (2.47)Nb 2 2 2 2 j,k j+1,kβ ℏj=1 k=1
Although this expression is similar to Eq. (2.37), it is not the same since here
we are dividing the double summation by Nb. To evaluate the remaining matrix
element in the system position basis we explicitly write out the product and
insert copies of the identity operator in the basis of system position:
⟨∫x|Ω̂ fs (q f∏1)Î Ω̂s (q )Î . . . ÎΩ̂
f
2 s (qNb)|x⟩ (2.48)
Nb
= dx ⟨x f( j−1) f+1|Ω̂s (q j) |x j f+1⟩
j=1
We now insert an identity operator in between each pair of adjacent operators
in the remaining prod∫uct to get∏Nb ∏f
dx ⟨x f( j−1) f+l|Ω̂s (q j) |x( j−1) f+l+1⟩ (2.49)
j=1 l=1
We can evaluate the remaining (positio)n matri[x element using Eq. (2.11):1/2 ( )]
⟨x( j−1) f+l|
m N −β
Ω̂s(q j)|x s s( j−1) f+l+1⟩ = [ 2 exp Ve f f (x2 ( j−1) f+l,q j) × (2.50)πβℏ Ns ]
exp −mNs2 (x( j−1) f+l+1 − x
2
2 ( j−1) f+l
)
βℏ
17
where ∑F ∑F c2
Ve f f (x,q) = V(x) − x c 2 kkqk + x 2 (2.51)
k=1 k=1 2mbωk
The entire product now becomes( )
m N Ns/2  ∑Nb ∑f ( )s s exp 2 −β Ve f f (x2 N ( j−1) f+l,q j)πβℏ s j=1 l=1 ∑Nb ∑f
− mNs2 (x
2
2 ( j−1) f+l+1
− x( j−1) f+l) (2.52)
βℏ
k=1 l=1
We are now in a position to arrive at our final expression for the MTS parti-
tion function. To get t(here, w)e co(mbine )Eqs. (2.46) and (2.52):Nb×F Ns ∫ ∫
m 2bNb m 2Z s
Ns
= dx d{q}e−βUMTS (x,{q})MTS 2 2 (2.53)2πβℏ 2πβℏ
where
U MTSMTS (x, {q}) = Ub∑({(q}) + U
MTS
s (x) (+ U
MTS
bs (x), {)q}), (2.54)1 Ns 2
UMTS
2
s (x) = V(x j) −
MNs
2 x j − xN 2 j+1 (2.55)s β ℏ2j=1
and
1 ∑Nb ∑f ∑F ( ( ))c
UMTSbs (x, {q})
k
= x 2( j−1) f+lckq j,k + x( j−1)r+l 2 (2.56)Ns j=1 l=1 k=1 2mbωi
This final set of equation gives us a way to treat our system and bath sub-
systems on different footing. We us a ring polymer with Ns beads for the system
and a ring polymer with Nb beads for each of the F bath modes so that our total
system will have Ns + F × Nb total DOF. Because we decided from the outset of
the derivation to use a larger Trotter dimension for the system than the bath,
we ultimately have the restriction that Ns > Nb and Ns/Nb be an integer. If this
last condition was not true, then the middle sum in UMTSbs would run over a
non-integer value.
18
Our last step is to pull out the corresponding RPMD version of the MTS
CL Hamiltonian by introducing momentum. We can do this by following the
procedure in e(arlier)sections:F×Nb+Ns ∫ ∫ ∫ ∫
1 2 MTS
Z = dx dP d{q} d{p} e−βHCL (x,P,{q},{p})MTS (2.57)2πℏ2
where ∑Ns P2 ∑Nb ∑F p2
HMTSCL (x,P, {q} { }
1 j 1p ) j,k, = + + U
N 2M N 2m MTS
(x, {q}) (2.58)
s j=1 b j=1 k=1 b
Normally at this point we would use HMTSCL to compute a correlation function
using the standard R∫PMD∫proto∫col: ∫
⟨ 1A(x̂)B(x̂)⟩ ≈ dx dP d{q} d{p} e−βHMTS (x0,P0,{q0},{p0})Ā(x0)B̄(xZ t) (2.59)MTS
However, if we were to use this expression we would find our correlation func-
tion would be wrong compared to exact results, even for short times where
RPMD performs well. The reason for this is because instead of grouping 1/N
with β like we had done previously, we were forced to group it with the UMTS .
We had to do so because there was no common factor of 1/N to be pulled out
— in fact, we had a 1/Ns and a 1/Nb. In the next section we explain why the
decision to group 1/N with β or the Hamiltonian matters and the proper way to
go between the two.
2.4 Handling 1/N
To understand why grouping 1/N with β or the Hamiltonian makes a differ-
ence in the outcome of a dynamics simulation we consider a simple example.
Suppose we have a one dimensional Hamiltonian:
p2
H(x, p) = + V(x) (2.60)
2m
19
We wish to compute a position de∫pend∫ent observable in the canonical ensemble:
⟨ 1A(xt)⟩ = dx dp e−βN H(x0,p0)A(xZ t) (2.61)
In this version of the calculation we are grouping 1/N with β, giving rise to βN in
the exponential. Our interpretation of this grouping is that the system is being
sampled at an artificially high temperature since βN = 1/(NTkB). From a pro-
cedural standpoint, the first step in computing Eq. (2.61) is running a Monte
Carlo simulation to sample the position coordinate on the distribution e−βNV(x).
p2
Momentum is directl√y sampled from the Gaussian distribution e−βN 2m with stan-
dard deviation σ = m/βN . The last step then is to integrate Hamilton’s EOM
to propagate the sampled initial conditions:
∂V(x)
ṗ = − (2.62)
∂x
p
ẋ =
m
Now let’s consider a second scenario in which we group 1/N with the Hamilto-
nian instead of β. We define the Hamiltonian
p2 1
HN = + V(x) (2.63)N2m N
The expression for a position dep∫enden∫t observable is now
⟨ 1A(x )⟩ = dx dp e−βHN (x0,p0)t N A(xt) (2.64)Z
Although Eqs. (2.61) and (2.64) are formally the same, they lead to different re-
sults when using Monte Carlo simulation. The first step in computing Eq. (2.64)
is to using Monte Carlo to sample the position coordinate over the distribution
−β V(x) p
2
e N . Next, we sample momentum from the Gaussian distribution e−β N2m . So
far everything in the procedure to compute Eqs. (2.61) and (2.64): the posi-
tion and momentum distributions generated by the sampling procedure will be
20
identical. But now consider the EOM that HN generates:
1 ∂V(x)
ṗ = − (2.65)
N ∂x
p
ẋ =
Nm
Because Eqs. (2.62) and (2.65) are not the same, but the underlying distributions
over which < A(xt) > and < A(xt) >N are computed is the same, then Eqs. (2.61)
and (2.64) cannot yield the same results. Our RPMD correlation functions will
therefore depend on whether or not we group 1/N with β or the Hamiltonian.
Clearly this situation is not desirable and we now demonstrate how it is
remedied. To do so, we repeat the exercise above, but this time instead of work-
ing with position and momentum, we use position and velocity. We will even-
tually recover momentum by performing a Lagrange transform on velocity. Our
expression for the observable A∫is ∫
1 − V(x0)⟨A(x )⟩ = dx dv e β N e−β 1 m 2t N 2 N v0 A(xt) (2.66)Z
where v = ẋ is the velocity of our system and we have decided not to group
1/N with β. Once again, the position distribution is generated by Monte Carlo
V(x)
sampling e−β N , leading to the same distribution as befor√e. Velocity is sampled1 m
from a Gaussian, e−β v22 N , with standard deviation σ = N/(βm). To get xt we
solve Lagrange’s EOM using the Lagrangian
1 ′ 2 1LN = m v − V(x) (2.67)2 N
where m′ = m/N. Our EOM generated by LN is
′ − 1 ∂V(x)m ẍ = (2.68)
N ∂x
Or,
∂V(x)
mẍ = − (2.69)
∂x
21
since each side of the equation as a 1/N. The fact that we can cancel N is crucial
to recovering consistent dynamics independent of how we group 1/N. Consider
what would happen if we had grouped 1/N with β in the Lagrangian version.
The distributions we sample position and momentum over would be the same,
but our Lagrangian would end up being
1
L = mv2 − V(x) (2.70)
2
which clearly generates the same EOM as LN . This means that computing dy-
namic observables using the Lagrangian framework is independent of how we
group 1/N.
We want to find a way to have the same consistency working within the
Hamiltonian formulation. The correct way to achieve this is to start with the La-
grangian of our system then Legendre transform to get our Hamiltonian. Start-
ing with LN , we find the canonical momentum defined as
∂L
p N= (2.71)
∂v
p = m′v
Our Hamiltonian is then the Legendre transform of LN :
HN = pv − LN (2.72)
p2 1
HN = m′
+ V(x)
N
p2 1
HN = N + V(x)m N
Comparing this result to Eq. (2.63) we see that the correct expression for
the Hamiltonian should involving multiplying the kinetic energy term by N
rather than dividing it. Put another away, the mass should be divided by a
22
factor of N rather than multiplied. Using this interpretation, momentum should
p2
b√e sampled on the Gaussian distribution e−β 2m′ , with standard√deviation σ =
m′/β. Using HN along with sampling momentum using σ√= m′/β generates
the same results as H and sampling momentum using σ = m/(βN).
Returning to Eq. (2.58), the correct expression for the Hamiltonian should be
∑Ns P2 ∑Nb ∑F p2
HMTSCL (x,P, {q}, {p})
j j,k
= ′ + ′ + UMTS (x, {q}) (2.73)2M 2m
j=1 j=1 k=1 b
where M′ = M/Ns and √m′b = mb/N. System momentum should be sampled from
a Gaussian with σ √= M′/β and bath momentum should be sampled from a
Gaussian with σ = m′b/β. With these rescaled masses, one can now safely run
RPMD dynamics using MTS on the CL Hamiltonian.
23
CHAPTER 3
MAPPING VARIABLE - RING POLYMER MOLECULAR DYNAMICS
Coupled electron-nuclear dynamics, or non-adiabatic dynamics, plays a cen-
tral role in a wide array of reactions of complex chemical and biological systems
[29, 30, 31]. However, developing theoretical simulation techniques that accu-
rately describe quantum dynamic effects in high-dimensional systems is already
challenging, and non-adiabatic processes additionally require that we treat elec-
tronic and nuclear degrees of freedom within a uniform dynamic framework.
In this chapter we discuss one such method designed for non-adiabatic dy-
namics, mapping variable ring polymer molecular dynamics (MV-RPMD), a
mapping variables, path integral technique for approximate dynamics of multi-
level systems. Although the MV-RPMD Hamiltonian has been successfully
used for the simulation of both adiabatic and non-adiabatic dynamics in several
model multi-level systems [3, 32, 33, 34, 35], the properties of this Hamiltonian
and others like it remain poorly understood.
We begin with a brief review and derivation of MV-RPMD theory. This
derivation differes from the original [3] in that we allow for a free parameter,
γ, when evaluating the SEO wavefunction rather than implicitly setting it equal
to one. Next, we derive two MTS versions of the MV-RPMD Hamiltonian. We
then prove that each version of the MV-RPMD Hamiltonian has a constant of
motion corresponding to the population of the nth state. Finally, we analyze the
performance of these new Hamiltonians by computing position auto-correlation
functions and comparing against exact results. We also study their numerical
stability by applying a modified stochastic Lagrangian Descriptors technique
developed within the Ananth group [4].
24
3.1 Review
MV-RPMD extends RPMD to include multiple potential energy surfaces by
mapping discrete electronic states to continuous Cartesian variables [3]. Con-
sider a Hamiltonian for a general k-level quantum system written as
Ĥ = T (P̂) + V0(Q̂) + Ve(Q̂) (3.1)
where the nuclear kinetic energy is T (P̂) = P̂T P̂/2M, and V0(Q̂) is the state-
independent potential operator. P and Q are momentum and position vectors,
respectively, and M is the mass. The electronic potential is
∑k
Ve(Q̂) = Vnm(Q̂) |ψn⟩ ⟨ψm| (3.2)
n,m=1
where Vnm(Q̂) are elements of the diabatic potential energy matrix and |ψn⟩ rep-
resent the electronic states. The k electronic states are mapped to k indepen-
dent, singly excited oscillator (SEO) states following the Meyer – Miller – Stock
– Thoss mapping protocol [36, 37]:
|ψ ⟩ ⟨ψ | → a†n m nam (3.3)
|ψn⟩ → |01, . . . , 1n, . . . , 0k⟩ ≡ |n⟩ (3.4)
where a†n and am are b[oson c]reation and annihilation operators that obey the
commutation relation a†n, am = δnm. In Eq. (3.3), we use the notation |n⟩ to
represent the nth SEO state that is a product of k 1D oscillator eigenfunctions
with k − 1 oscillators in the ground state and a single oscillator (the nth) in the
first excited state.
Path integral discretization of the canonical partition function is performed
using continuous Cartesian variables for both nuclear and electronic degrees
25
of freedom. This is achieved by inserting N − 1 copies of the identity operator
defined as ∫ ∫
Î = dx dQ|x,Q⟩⟨x,Q|P̂ (3.5)
where x is a k-dimensional continuous electronic state variable, and the projec-
tion operator, defined as ∑k
P̂ = |n⟩⟨n| (3.6)
n
is necessary to constrain |x⟩ to the SEO subspace [38]. Inserting N copies of the
identity in Eq. (3.5), w∫e obtain∫ ∏N 〈 ∣∣ ∣∣ 〉
Z = lim d {Qα} d {xα} xα,Qα ∣∣P̂e−βN ĤP̂∣∣ x→∞ α,Qα (3.7)N
α=1
Following the derivation for RPMD in Chapter 2, the Trotter approximation can
be used to evaluate the n∫uclear only high-temperature matrix elements:∏N
Z = lim d {Q } e−βN HRP({Pα},{Qα})I
→∞ α el
({Qα} , {xα}) (3.8)
N
α=1
where ∫ ∏N 〈 ∣ ∣ 〉
I = d { ∣x } x ∣∣P̂e−βNVe(Q̂α)el α α P̂∣∣∣ xα (3.9)
α=1
Recognizing that each of the N k-dimensional electronic variable integrals can
be written in terms of N trace operations, and replacing the operator traces by
phase space inte∫grals ov∫er the co∫rresponding Wigner functions, we obtain
1
Iel = N d {xα} d {p〈α} d {∆x }(2 απ) ∏N ∣ ∣ 〉
x − ∆xα ∣∣P̂ ∆x Te−βNVe(Q̂α)α P̂∣∣x αα + eipα∆xα (3.10)2 2
α=1
Finally, using the definition of the( SE) O wavefunction,
γ N/4 γ⟨x | ⟩ Tn = (2γ)1/2[x] e−n 2 x x (3.11)
π
26
where [·]n indicates the nth element of the enclosed vector, we analytically eval-
uate the electronic integral.
The form of Eq. (3.11), including γ, comes directly from the expression of the
harmonic oscillator eigenfunctions. When MV-RPMD was originally derived γ
was implicitly set equal to one [1]. However, there is no clear reason why this is
necessarily the correct choice of γ, or why any one choice should be better than
another. For these reasons, we leave γ as a general variable in this derivation,
allowing one to explore how different choices impact the dynamics.
Substituting the previous two expressions back into Eq. (3.8), we obtain our
final expre∫ssion for∫the canon∫ical par∫tition function:
Z ∝ lim d {Q } d {P } d {x } d {p } e−βN HMV({Qα},{Pα},{xα},{pα})α α α α sgn(Θ) (3.12)
N→∞
where sgn(·) returns the mathematical sign of the argument. In Eq. (2.17), the
MV-RPMD Hamiltonian is defined as
1 1
HMV ({Qα} , {Pα} , {xα} , {pα}) = HRP ({Qα} , {Pα}) + G − ln|Θ| (3.13)
βN βN
where HRP is defined in Eq. (2.13), and
∑N 1
G = γxTαxα + p
T
αpα (3.14)γ
α=1
We provide a graphical representation of HMV in Fig. 3.1. We can clearly
see the original RPMD Hamiltonian structure emerge from the harmonic spring
interactions among the grey nuclear beads. The web of straight lines in the
center of the red electronic beads shows that all electronic variables are coupled
to all other electronic variables. This is due to the complicated Γ matrix which
mixes DOF together. This term also couples all nuclear beads to all electronic
27
Figure 3.1: Depiction of MV-RPMD Hamiltonian [4]. Red beads in the cen-
ter represent electronic variables; grey beads on the outside represent nuclear
beads. Squiggly lines connecting nuclear beads represent harmonic spring in-
teractions. Straight lines represent interaction with electronic beads.
beads, however, we have omitted drawing this interaction in Fig. 3.1 for visual
clarity. The electron-nuclear interaction term in Eq. (3.13) is defined as:
Θ = Re(Tr[Γ]) (3.15)
where ∏N
Γ = C (xα,pα)M (Qα) (3.16)
α=1
The electronic matrix, C, i(s defined as ) ( )T
C √ i √ i 1(xα,pα) = γxα + √ pα ⊗ γxα − √ pα − I (3.17)
γ γ 2
where I is the k×k-dimensional identity matrix and the nuclear matrix elements
in the high temperature limit aree−βNVnn(Qα) if n = mMnm (Qα) = lim  (3.18)N→∞ −βNV (Q ) e−βNVnn(Qα)nm α if n , m
28
In deriving the RPMD Hamiltonian in Eq. (2.13), we introduced nuclear mo-
mentum through normalized Gaussians. We note that this procedure cannot be
used to generate a phase space representation in the electronic variables since
the SEO subspace corresponds to a constrained phase space that must be ob-
tained using Wigner transforms with the projection operators [3].
Real-time quantum thermal correlation functions are approximated within
the MV-RPMD framewo∫rk by ∫ ∫ ∫
CAB(t) = d {Qα} d {Pα} d {xα} d {pα} (3.19)
× e−βN HMV ({Qα(0)},{Pα(0)},{xα(0)},{pα(0)})Ā(0)B̄(t) sgn(Θ) (3.20)
where trajectories are propagated to time t under HMV .
3.2 Mixed Time Slicing
As discussed in the previous chapter, path integral based methods converge
as a function of the number of beads, N, and in general the more significant the
quantum effects in a system, the more beads are required to achieve this conver-
gence. In addition, for a given system, it is frequently true that different degrees
of freedom require different numbers of beads to reach convergence. For exam-
ple, describing a tunneling proton requires a larger number of beads than simu-
lating the dynamics of a heavier carbon atom at a given temperature. Since the
computational cost associated with a path integral simulation increases with N,
there is significant motivation to optimize the number of beads in a simulation.
In this section we derive MV-RPMD Hamiltonians for two cases, each using
differing numbers of nuclear beads (Nn) and electronic beads (Ne). The first case,
29
called quantized nuclei (QN), has the condition Nn ⩾ Ne ⩾ 1. The second case,
called quantized electronic (QE), has the condition Nn = 1 and Ne > 1.
Case 1: Starting with the canonical partition function for a general multi-level
system, we use a symmetric Trotter discretization [21] and the identity operator
defined in∫Eq. (3∫.5) to obtain[ ]
β Ne
Z = dQ∫ dx∫⟨Q, x|P̂ e− N Ĥe [ P̂|Q, x⟩ ] (3.21)
β β β Ne
= lim dQ dx⟨Q, x|P̂ e− 2N (T̂+V̂e 0)e− N V̂e ee− 2Ne (T̂+V̂0) P̂|Q, x⟩ (3.22)
Ne→∞
Because the aim of Case 1 is to use more nuclear than electronic beads, we fur-
ther split∫nucle∫ar DOF so tha[t(Nn ⩾ Ne: )
β β β Nn/2Ne
lim dQ dx ⟨Q, x|P̂ e− 2N V̂0e− N T̂ −n n 2N V̂n 0
Ne,Nn→∞ ( e ) ] (3.23)Ne
− β V̂ − β V̂ β β
Nn/2Ne
× e Ne e e 2Nn 0e− N T̂ − V̂n e 2Nn 0 P̂|Q, x⟩ (3.24)
Inserting multiple copies of identity and evaluating matrix elements as de-
scribed earlier, we obt∫ain
ZQN ≈ d {Qα} {
( )
d Pα} d {xα} d {pα} e−βNn HQN sgn ΘQN (3.25)
where βNn = β/N∑n. (The corresponding MV-RPMD Hamiltonian isthenNn PT P M 
H α αQN = + 2 (Qα+1 −Qα)
T (Q 1 −Q ) + V0 (Q )α+ α α  (3.26)2M 2β
α=1 Nn ( )
1 ∑Ne 1 1 ∣
γxT x T ∣ ∣∣+ α α + pαpα − ln ∣ΘQN∣ (3.27)βNn γ βα=1 Nn
and  ∏N ( )e Θ = Re Tr C (x ,p )M Q QN α α MTS (α−1)Nn/Ne+1  (3.28)
α=1
C is the same as in Eq. (3.17), and MMTS has matrix elementse−β Ne
Vnn(Q)
 if n = mMMTSnm(Q) =  (3.29)−βNeV −βN Vnn(Q)nm(Q)e e if n , m
30
where βNe = β/Ne.
The key difference between HMV and HQN can be seen by comparing Θ with
ΘQN. First, we note that Eq. (3.29) uses βNe where Eq. (3.18) uses βN . Further,
in Θ, every bead has both nuclear and electronic state information, whereas in
ΘQN only certain nuclear beads contain electronic information. These beads are
indicated by the index (α−1)Nn/Ne +1 in Eq. (3.28). This is demonstrated by the
illustration in Fig. 3.2.
Case 2: Starting with the multi-level system partition function again, we per-
form a symm∫ etric∫Trotter split[ting w]here we split Ve about T and V0:
β Nn
Z = dQ∫ dx∫⟨Q, x|P̂ e− N Ĥn [ P̂|Q, x⟩ ]
β β N
(3.30)
= lim dQ dx⟨Q, x|P̂ e− 2N V̂n ee− Nn ( ) −
β nT̂+V̂0 e 2N V̂n eP̂ |Q, x⟩
Nn→∞
In this case, we further split electronic DOF Ne times. In the limit where nuclear
DOF are tre∫ated c∫lassically, we can set Nn = 1 to obtain( )
β Ne/2
Z ≈ lim dQ dx⟨Q, x|P̂ e− N V̂e e ×
Ne→∞ ( )
− β β β β
Ne/2
e 2N V̂0e− N T̂ e− 2N V̂n n n 0 e− N V̂e e P̂|Q, x⟩ (3.31)
Evaluating the matrix elements as outlined in the previous section, we arrive at
our final expression fo∫r the partition function:
Z ≈ d {Q } d {P } d {x } d {p } e−βH ( )QEQE α α α α sgn ΘQE (3.32)
where βNn = β/Nn. Since we have fixed Nn = 1, Eq. (3.32) is not, in general, exact
even in the limit Ne → ∞. This is unlike Eq. (3.25) where ZQN can be made exact
in the limit Ne → ∞ and Nn → ∞. For Ca(se 2, the Hamil)tonian is
PT P 1 ∑Ne ∣ ∣
H = + V (Q) + γxT
1 T 1
QE x + p p − ln ∣∣ ∣Θ ∣ (3.33)2M 0 β α α α QEγ α β
α=1
31
and
Θ = ReTr ∏N e QE C (x α,pα)MMTS(Q) (3.34)
α=1
C is defined the same as Eq. (3.17), and MMTS is defined the same as Eq. (3.29).
Examining the form of HQE, we see that a single nuclear bead will inter-
act with Ne electronic beads through ΘQE as pictured in Fig. 3.2. We can see
that each Hamiltonian offers reduced dimensionality compared with the orig-
inal MV-RPMD Hamiltonian. This serves to increase computational efficiency
while also offering a simper interpretation of the system’s behavior. In Section
Figure 3.2: (a) Cartoon depiction of HQN. Since Case 1 has Nn > Ne, only certain
nuclear bead numbers have corresponding electronic bead numbers. (b) Car-
toon depiction of HQE. A single classical nuclear bead interacts with multiple
electronic beads [4].
3.5, we examine the accuracy of HQN and HQE for the simulation of dynamics in
model systems.
3.3 Constants of Motion
Constants of motion provide a valuable tool for assessing the accuracy of dy-
namic simulations. By computing the quantity corresponding to a constant of
32
motion throughout a simulation and verifying it remains constant, one gains
confidence in the correctness of their code. Usually, if there is something wrong
- a bug in the code, or unconverged parameters - the constant will drift. A com-
mon example is checking that individual trajectories conserve energy up to a
specified tolerance for Hamiltonian dynamics.
In the MV-RPMD Hamiltonian the mapping variables report on the popula-
tions of the discrete electronic states. This(allows us to calc)ulate the population
of the nth state of the αth α 1 2
[ ]
bead as Pn = 2 [xα]
2
n + pα n − 1 ; the so-called semi-
classical population estimator [17, 34, 35]. The semiclassical mapping Hamil-
tonian preserves the total population∑[39](and we[ w]o)uld therefore expect the
MV-RPMD Hamiltonian to preserve k [x ]2 2n=1 α n + pα n for each bead, indicat-
ing the total population of a given bead is a constant of motion. As a corollary,
the electronic G-term defined in Eq. (3.14) will also be a constant of motion for
dynamics under the MV-RPMD Hamiltonian. In this section we prove that this
is indeed the case.
Proof : Recall that in classical mechanics a quantity A is a constant of motion
of Hamiltonian H if the followi∑ng h{olds: }
{ ∂A ∂H ∂A ∂HA,H} = − = 0 (3.35)
∂qi ∂pi ∂pi i ∂qi
where {·, ·} den∑otes the Poisson bracket of two variables. In the present case
we have A = k γx2 1 2n=1 α,n + pα,n , which depends only on electronic DOF. Ourγ
Hamiltonian is HMV.
The first term of HMV defined in Eq. (3.13) is independent of xα and pα and
depends only on the nuclear degrees of freedom. The Poisson bracket of A with
this terms w∑ill trivially be zero. The second term can be written as G = 1N (A + F)
where F = ′N T 1j=1 γxαxα + pTαpα with the prime in the sum indicating that theγ
33
summation is over all j , α. The Poisson bracket of {A,G} is then 1N {A, F} +
1
N {A, A} = 0. { }
The final term to evaluate is A,− 1 ln|Θ| = 1 {ln|Θ|, A}. Using linearity of the
βN βN
Poisson bracket, we can split A∑into two terms givingk { } 1 { }{ln|Θ|, A} = γ ln|Θ|, x2α,k + ln|Θ|, p2α,k (3.36)γ
n=1
To evaluate these terms we substitute the definition of the Poisson bracket.
Terms involving partial derivatives with respect to the nuclear variables Qα and
Pα can be ignor{ed since }they give zero contr(ibutio)n: ( [ ])
2 ∂ ln|Θ| 2x ∂Γ
{ln|Θ|, xα, j} = −2x
α, j
α, j =(− ) R( e Tr , (3.37)∂pα, j Θ [ ∂pα],)j
| | 2 ∂ ln|Θ|
2p ∂Γ
lnΘ , pα, j = 2p
α, j
α, j = Re Tr (3.38)
∂xα, j Θ ∂, xα, j
where Γ can be expanded as
Γ = C (x1,p1)M (Q1) . . .C (xα,pα)M (Qα) . . .C (xN ,pN)M (QN) (3.39)
and both C and M have been previously defined in Eq. (3.17) and Eq. (3.18),
respectively. Taking the partial derivative and introducing two matrices, Φα
and Λα, where
Φα=1 = I∏ (3.40)α−1
Φα,1 = C (xi,pi)M (Qi) (3.41)
i=1
and
Λα=N =M (QN) ∏N (3.42)
Λα,N =M (Qα) C (xi,pi)M (Qi)
i=α+1
we obtain ( ) k ( )
∂Γ ∑ ∂C (x
Tr ( ) α
,pα)
= Φα nl (Λ ) (3.43)∂xα, j ∂x
α mn
n,l,m=1 α, j lm
34
Similarly, we have ( ) ∑k ( )∂Γ ∂C (xα,p )Tr = (Φα) αnl (Λα)mn (3.44)∂pα, j ∂pn,l,m=1 α, j lm
Evaluating th( e derivativ)es of t(he matrix C in E)q. (3.17) gives∂C (x ( )α,pα)
( ) = γ xα,lδm, j + xα,mδl, j + i pα,lδx ( ) ( m, j
− pα,mδl, j (3.45)
∂ α, j lm
∂C (x )α,pα) 1
= pα,lδm, j + pα,mδl, j + i xα,mδl, j − xα,lδm, j (3.46)
∂pα, j lm γ
Substituting these equations back into the Poisson bracket in Eq. (3.36) and
evaluating the Kronec∑ker delta functions, we find2 k [ (( ){ln|Θ|, A} = Re (Φα)nl(Λα)mn xα,l pα, jδm, j − xα, j pα,mδl, j
Θ
j,l,m,n(=1 ))] (3.47)
+ xα,m pα, jδl, j − xα, j pα,lδm, j = 0
Thus, {A,HMV} = 0, making A a constant of motion. As a result, the dynamics
conserves total population along a trajectory and G must therefore be a constant
of motion.
3.4 Model Systems and Numerical Simulation Details
3.4.1 Model Systems
Here we investigate the accuracy of the various MV-RPMD Hamiltonians de-
rived earlier through numerical simulations for a series of model systems with
two electronic states coupled to a single nuclear coordinate.
The state-independent potential for all models is
1
V0(Q) = Q2 (3.48)2
35
Diagonal elements of the diabatic potential matrix are
V11(Q) = Q + c, (3.49)
V22(Q) = −Q (3.50)
where c is a constant that determines the asymmetry of the system. The off-
diagonal elements are V12 = V21 = ∆. Parameters used for various model sys-
tems are defined in Table 3.1. These three models allow us to explore a range
Model c ∆
I 0 0.1
II 3 0.1
III 0 0.01
Table 3.1: Parameters for Models I-III provided in atomic units. For all models,
M = β = γ = 1 a.u. Simulations where we investigate the effect of changing γ
use values 1, 10, and 100 a.u.
of different regimes of dynamics. In Model I, dynamics are in the moderate
coupling non-adiabatic regime with zero driving force (symmetric). Model II
uses the same coupling regime, but has a significant driving force. Model III is
a symmetric system with lower coupling making non-adiabatic dynamic effects
more pronounced (βV12 ≪ 1).
3.4.2 Nuclear Position Auto-Correlation Function
We investigate the accuracy of dynamics generated by the various MTS MV-
RPMD Hamiltonian by computing nuclear position auto-correlation functions
and bench-marking against exact quantum results. The thermal real-time posi-
tion auto-correlation function within the MV-RPMD approximation is written
36
as ∫ ∫ ∫ ∫
C (t) = d {Q } d {P } d {x } d {p } e−βN HMVQQ α α α α Q̄(0)Q̄(t) sgn(Θ) (3.51)
∑
where Q̄ = 1 NN α=1 Qα [3].
Initial conditions for all trajectories are sampled using Monte Carlo methods
with Metropolis Hastings acceptance/rejection criteria. Monte Carlo moves are
proposed by randomly selecting to move either nuclear or electronic DOF sep-
arately. We note that, for the case of symmetric, double-well systems, a modest
improvement in convergences is made by introducing “flip” moves that deter-
mine which state of the system is occupied. This helps prevent the system from
becoming trapped in one state.
The time-evolved positions Q(t) are obtained by integrating the MV-RPMD
equations of motion using a 4th-order Adams–Bashforth–Moulton predictor-
corrector integrator with a time step dt = 10−4 a.u. Each trajectory is integrated
for 20 a.u. We use ≈ 93,000 trajectories for each model. Trajectories that fail to
conserve energy to within a specified tolerance are removed from the sampled
ensemble. The time step is chosen to ensure that < 10% of sampled trajectories
fail to conserve energy.
Correlation functions are computed using Eq. (3.51). For the MTS imple-
mentations, we use either HQN or HQE depending on the case being studied.
Error bars for each simulation are within the symbol size plotted. Exact results
for each model are obtained using a discrete-variable representation (DVR) grid
calculation [40].
37
3.4.3 Analysis of Auto-Correlation Functions
We begin our numerical investigations with a study of the effect of changing the
width parameter, γ, in Eq. (3.11) when mapping the discrete electronic states to
SEO states. Since the mapping is formally exact, exact quantum dynamics under
a mapping Hamiltonian will be rigorously independent of γ. However, this is
not necessarily true for approximate dynamics such as those generated within
the MV-RPMD framework.
We find that numerically computing real-time correlation functions gener-
ated by Model I to be identical across three values: γ = 1, 10, 100 a.u. Addition-
ally, we tested and confirmed that this result is true not only for the standard
MV-RPMD Hamiltonian, but also for the Hamiltonians generated by the MTS
derivation.
This finding suggests that MV-RPMD is numerically independent with re-
spect to setting γ. In fact, we argue that the best choice is γ = 1 so that Eq. (3.14)
is symmetric in x and p. This makes numerical instabilities less likely to occur.
Next, we investigate the performance of our MTS Hamiltonians for three
model systems defined in the previous section. For Model I shown in Fig. 3.3,
we find that the fully quantum treatment of both nuclear and electronic DOF
within the MV-RPMD framework produces results closest to exact, as expected.
Unexpectedly, a fully classical treatment where NE = NN = 1 seems to outper-
form all mixed quantization dynamics. Dynamics with electronic variables in
the classical limit (NE = 1 and NN = 6 ) is the least accurate, failing to repro-
duce even the first oscillation. However, the MTS case NE = 3 and NN = 6 is
perhaps more concerning as the correlation function at later times grows be-
38
yond the time zero value, indicating the unphysical nature of the dynamics. We
find that treating the nuclear DOF classically but keeping electronic variables
quantum mechanical leads to some slight improvement in the accuracy of the
correlation function over the classical limit. In general, MV-RPMD performs
Figure 3.3: Nuclear position auto-correlation function for Model I. The number
of electronic and nuclear beads is indicated by the format Ne : Nn and corre-
sponds to the following lines: (a) 1:6=red triangles. 3:6=green squares. 6:6=blue
circles. Exact=black dashes. (b) 1:1=blue circles. 4:1=magenta diamonds. Ex-
act=black dashes [4].
better in simulations of Model II as shown in Fig. 3.4 Again, we find a fully
quantum treatment performs best. However, now the intermediate quantiza-
tion Ne = 2 and Nn = 4 performs surprisingly well. A fully classical treatment
appears to match the exact oscillation frequency, but gradually decays in ampli-
tude. Again, we find that, when the electronic DOF are treated in the classical
limit, the auto-correlation function exhibits unphysical oscillations, while the
classical-limit nuclear simulations are better behaved but still less accurate than
the fully classical simulation. These results suggest that in the asymmetric case,
where the reaction has a significant driving force, the nuclei must be treated
more quantum mechanically than the electronic DOF.
Finally, the nuclear position auto-correlation function for Model III exhibits
39
Figure 3.4: Nuclear position auto-correlation function for Model II. The number
of electronic and nuclear beads is indicated by the format Ne : Nn and corre-
sponds to the following lines: (a) 1:4=red triangles. 2:4=green squares. 4:4=blue
circles. Exact=black dashes. (b) 1:1=blue circles. 4:1=magenta diamonds. Ex-
act=black dashes [4].
a trend very similar to Model I as shown in Fig. 3.5. Like the other models, we
find that a fully quantum treatment produces results closest to exact. For this
model, a fully classical treatment outperforms all intermediate quantizations.
When Ne = 1 and Nn = 4, the amplitude is too negative within the first oscilla-
tion. In the case Ne = 2 and Nn = 4, the function grows past the point at time
zero. Quantizing electronic DOF only, with Ne = 4 and Nn = 1, improves very
slightly upon the classical treatment.
Taken together, these results suggest that for dynamics in all three regimes
it is best to opt for either the fully quantum MV-RPMD Hamiltonian dynamics
or the computationally efficient classical-limit methods. The MTS results, while
very good in some cases (Model II), show unphysical oscillations suggesting
that the dynamics is unstable and its accuracy system-specific. This observation
contrasts with other path integral based methods that have had success treat-
ing nuclei in the classical limit through linearization [18] while quantizing elec-
tronic mapping variables. We note also that the model systems studied here are
40
Figure 3.5: Nuclear position auto-correlation function for Model III. The num-
ber of electronic and nuclear beads is indicated by the format Ne : Nn and corre-
sponds to the following lines: (a) 1:4=red triangles. 2:4=green squares. 4:4=blue
circles. Exact=black dashes. (b) 1:1=blue circles. 4:1=magenta diamonds. Ex-
act=black dashes [4].
likely most challenging for MTS simulations since electron-nuclear dynamics
are strongly coupled; several of these MTS results can be reasonably expected
to be more accurate in system-bath models where dynamics is less strongly cou-
pled. For coupling to a bath we expect oscillations to dampen quickly, in which
case using a MTS Hamiltonian might be a good option since they retain accu-
racy for short times.
3.5 Lagrangian Descriptors
Lagrangian descriptors (LDs) have been shown to be a useful tool to charac-
terize the stability of a Hamiltonian dynamics. Following recent work [41], we
employ a p-norm LD definition ∑n ∫ t0+τ
Mp (w0, t0, τ) = |v pi (w (t; w0) , t)| dt (3.52)
i=1 t0−τ
where p ∈ (0, 1],w is a phase space vector, vi(w, t) = dwi/dt, and τ is the integra-
tion time. It has been shown that evolution forward in time will reveal stable
41
manifolds, while backward integration reports on unstable manifolds [41, 42,
43].
For low-dimensional systems with relatively simple Hamiltonians, LDs of
a constant energy manifold are typically computed by fixing a subset of coor-
dinates and inverting the Hamiltonian to determine values for the remaining
coordinates given a fixed energy. However, even for the model systems studied
here in the classical limit (Ne = Nn = 1) we are unable to invert the MV-RPMD
Hamiltonian, HMV. Rather, we employ a stochastic gradient descent (SGD) al-
gorithm [44] modified to produce a constant energy manifold when we fix the
nuclear phase space variables. Standard SGD approaches find stationary points
on multidimensional surfaces; we modify the algorithm to find values for the
electronic phase space variables that yield the correct total energy at a particular
nuclear phase space grid point.
Here we compute LDs for Model III with β = 3 and explore the constant
energy manifold E = 1.2 a.u. We scan over a grid of nuclear position and
momentum and use our modified SGD algorithm to find x and p such that
HMV(Q, P, x,p) = 1.2 a.u. The a∣lgorithm wo∣ rks by first supplying an initial phase
space point v0. From here, ∣∣E − H (v)0∣MV ∣ < δ is calculated. If the inequality
holds, then the algorithm terminates because a viable phase space point has
been found. Otherwise, v is updated by
( ( ))  
vi+1 ← vi j j + η E − HMV vi sign ∂HMV i ∂v  (3.53)j
Here, vi+1j is the ith update of the jth element of v, and η is a positive factor that
scales the overall magnitude of each update. At each update, j is randomly
chosen among the elements which are not held constant. In our case, Q and P
are held constant for a given phase space point, and v j is randomly chosen from
42
∣∣ ( )∣∣
x and p. The update is repeated until ∣∣E − H vi ∣MV ∣ is sufficiently small. In our
simulations, we use η = 0.1 and tolerance δ = 10−7.
Below, we utilize this method to compute p-norm LDs in which τ = 4 a.u.
and p = 1/2 and study saddle points of HMV with NE = NN = 1. We compute
LDs over the P-Q phase space grid and we do not find it necessary to define an
interaction region for our calculations. This study represents the first investiga-
tion of non-adiabatic Hamiltonian dynamics by identifying stable and unstable
manifolds in the phase space of the classical-limit 6D MV-RPMD Hamiltonian.
The nuclear phase space LDs for Model III are shown in Figs. 3.6 and 3.7
below. Each exhibits a characteristic crossover for a reactive system in which
we identify the arms of the ”X” forming a negative slope as corresponding to a
stable manifold, and the arms forming a line with positive slope as an unstable
manifold. In Fig. 3.6 the point at which the stable and unstable manifolds cross
is at (P = 0,Q = −0.13), but this point is not fixed. In Fig. 3.7 we find the crossing
point occurs at (P = 0,Q = 0.2). The changing location of the crossing point can
be attributed to the relative populations of the electronic states calculated from
the electronic variable values obtained by the modified SGD solver for each
nuclear grid point.
In the case of Fig. 3.6, we show that state |1⟩ is mostly occupied, while |2⟩
is mostly unoccupied. Fig. 3.7 is the result when these relative populations are
swapped with |2⟩ being the mostly occupied state and |1⟩ mostly unoccupied.
Cross section plots demonstrate the points of inflection for the LD as a function
of nuclear position with fixed momentum. These results allow us to identify the
top and bottom enclosed sectors of the LD plot as regions of phase space with
reactive trajectories, while the left and right sectors correspond to nonreactive
43
Figure 3.6: (a) LD calculated over a grid in nuclear phase space shows a stable
manifold (Ws) and an unstable manifold (Wu). The crossing point lies at P = 0
and Q = −0.13. (b) Cross section at P = −0.2 showing sharp changes in the LD
corresponding to the location of stable and unstable manifolds. (c) We calculate
the semiclassical population estimator using the electronic variables obtained
through the SGD simulation and find that the LDs shown here correspond to
state |1⟩ being mostly occupied and (d) state |2⟩ being mostly unoccupied [4].
trajectories.
To further understand the location of the crossing point - important because
this point lies on the normally hyperbolic invariant manifold (NHIMs) - we plot
the MV-RPMD Hamiltonian on a nuclear phase space grid with electronic vari-
ables fixed to values corresponding to one state being fully occupied and the
other state being unoccupied. The Fig. 3.8 show the barrier region in the nu-
clear Q direction shifting in concert with the crossing point we observe in the
44
Figure 3.7: (a) LD calculated over a grid in nuclear phase space shows a stable
manifold (Ws) and an unstable manifold (Wu). The crossing point lies at P = 0
and Q = 0.2. (b) Cross section at P = 0.2 showing sharp changes in the LD
corresponding to the location of stable and unstable manifolds. (c) We calculate
the semiclassical population estimator using the electronic variables obtained
through SGD simulation and find that the LDs shown here correspond to |2⟩
being mostly occupied and (d) |1⟩ being mostly unoccupied [4].
the LDs as the electronic state populations are modulated.
45
Figure 3.8: (a) Nuclear phase space plot showing HMV where we choose the
electronic variables such that only |1⟩ is occupied. (b) Nuclear phase space plot
showing HMV when only |2⟩ is occupied [4].
46
CHAPTER 4
QUANTUM COMPUTING WITH PATH INTEGRALS: CONTINUOUS
SYSTEMS
Academic and industrial research efforts aimed at advancing the field of
quantum computing have seen an explosion of growth in recent years [45].
Much of this growth has been spurred on by recent advances in hardware as
well as Google’s claim to have recently achieved quantum supremacy [46]. As
research into quantum hardware continues to advance the challenge becomes
identifying problems that quantum computers will be able to solve and con-
structing algorithms to do so. Several problems in chemistry are expected ben-
efit substantial from the aid of quantum computing [19, 20].
Quantum computing algorithms are often discussed as two classes: far-term
algorithms which assume access to fault tolerant machines, such as Quantum
Fourier Transform (QFT)[47] and phase estimation[48], and near-term algo-
rithms which anticipate noisy qubits and error accumulation, such as Varia-
tional Quantum Eigensolver (VQE) [49]. In this thesis we limit ourselves to
discussion of near-term algorithms since fault tolerant machines do not yet ex-
ist. Near-term algorithms run on Noisy Intermediate Scale Quantum (NISQ)
machines [50]. This name comes from the fact that modern quantum computers
contain a small number - on the order of 10’s - of qubits. The qubits are noisy
in the sense that they will eventually decohere as they couple to their environ-
ment, causing erroneous results. Although error correcting schemes to protect
against incorrect results exists, they require on the order of 1000’s of qubits to
implement [51].
Instead of using quantum computers to execute algorithms in their entirety,
47
NISQ era algorithms often use hybrid quantum-classical schemes to limit the
accumulation of error [50]. In general, hybrid algorithms stem from a classi-
cal algorithm, but are modified to query a quantum computer only in certain
instances. Because a large part of the algorithm is computed using classical re-
sources, deep quantum circuits are avoided which mitigates error accumulation
and noisy results. A popular example of a hybrid algorithm is VQE, which is
used for ground state estimation in a chemistry context [49]. In this algorithm
an easily prepared variationally parameterized wavefunction is initialized on
a quantum computer and individual terms of the molecular Hamiltonian are
measured. The ground state energy is then estimated by reconstructing the
measured terms of the Hamiltonian on a classical computer. Based on this en-
ergy a classical optimization algorithm then proposes a new parameterization
for the wavefunction which should move the trial state incrementally closer to
the exact ground state eignefunction of the Hamiltonian. This new wavefunc-
tion is then prepared on the quantum computer and the loop repeats until some
convergence criteria is met.
In this chapter we discuss a new hybrid algorithm which uses quantum com-
puting to increase the performance of a classical Path Integral method to com-
pute correlation functions. We begin with a brief overview of quantum comput-
ing basics so the reader has enough background to understand quantum circuit
diagrams that appear throughout the chapter. We then outline the classical Path
Integral algorithm and identify where its performance can be improved with
quantum computing, specifically the task of evaluating position basis matrix
elements of a complex-time evolution operator. The next three sections are ded-
icated to building a generic quantum circuit which can be used to compute the
desired matrix elements. We then discuss a specific implementation for the case
48
of a system described by harmonic oscillator and how to curb numerical difficul-
ties that arise from propagation under the kinetic energy operator. Finally, we
examine the accuracy of our method by computing the position auto-correlation
function for a harmonic oscillator and compare against exact results. The final
section discusses helpful simulation notes for anyone who might wish to repro-
duce this experiment in the future.
4.1 Quantum Computing Basics
The basic unit of information in quantum computing is the qubit. As implied
by the name, a qubit is simply any two-level (binary) quantum system [52]. The
most generic state a qubit can occupy is
|ψ⟩ = c0|0⟩ + c1|1⟩ (4.1)
Here, both c0 and c1 can be complex numbers. There are many physical systems
that are used to realize qubits and currently no clear best system has been de-
cided upon. Popular systems include trapped ions[53], superconducting quan-
tum circuits[54], and quantum dots[55]. It should also be mentioned quantum
computing formalisms based on general n-level systems known as qudits exist.
We join several qubits together by taking the tensor product of each individ-
ual qubit. For a system with n qubits, we have
[ ] [ ] [ ]
|ψ⟩ = c0,qn−1 |0⟩ + c1,qn−1 |1⟩ ⊗ c0,qn−2 |0⟩ + c1,qn−2 |1⟩ . . . ⊗ c0,q0 |0⟩ + c1,q0 |1⟩ (4.2)
where q j indexes the jth qubit. By expanding each term in the product, we can
write |ψ⟩ in a more compact summation notation:
∑1 ∑1
|ψ⟩ = . . . cqn−1,...,q0 |qn−1, . . . , q0⟩ (4.3)
qn−1=0 q0=0
49
Here, cqn−1,...,q0 is the coefficient corresponding to the state |qn−1, . . . , q0⟩, and we
drop explicit reference to the tensor product.
The string qn−1, . . . , q0 encodes a binary number. Often times it will be con-
venient to switch between binary and decimal representation of our quantum
system. The equivalent representation in decimal is
2∑n−1
|ψ⟩ = ck|k⟩ (4.4)
k=0
This is easily seen since the binary string qn−1, . . . , q0 can represent all integers
from 0 to 2n − 1, inclusive. As an example, if we have n = 2 qubits, the two
representations are:
|ψ⟩ = c00|00⟩ + c01|01⟩ + c10|10⟩ + c11|11⟩ (4.5)
|ψ⟩ = c0|0⟩ + c1|1⟩ + c2|2⟩ + c3|3⟩ (4.6)
The set of basis states in either representation are referred to as the computa-
tional basis. Eq. (4.4) captures the power of information storage using qubits.
For the price of n qubits, we are able to store 2n pieces of information. This ex-
ponential scaling of information representation does not exist in classical com-
puting. In the classical case n bits store only n pieces of information.
In any model of computing we need a way to store information and a
method to perform operations and manipulate the information. Qubits pro-
vide our model for information storage; quantum circuit diagrams will be our
model for information manipulation.
In Fig. 4.1 is an example of a simple quantum circuit diagram. Time is un-
derstood to flow from left to right as we go across the circuit. Kets on the left
side represent initial states fed into the algorithm. Blocks with writing represent
50
|0⟩ H
Figure 4.1: Quantum circuit diagram illustrating single qubit passing through
Hadamard gate. Measurement occurs at the end.
quantum gates which act on the state feed into them. The meter at the end of the
circuit, or right side, indicates measurement. The double wire exiting the mea-
surement gate represents flow of classical information; single wires represent
flow of quantum information. A double wire will always follow after a mea-
surement because doing so collapses the quantum wavefunction to a classical
state.
In this specific circuit we being with a single qubit in the state |0⟩. Next a
Hadamard gate, H, which represents theopera 
tor,

1 
H = √ 2 1 1 
 (4.7)1 −1
is applied. Directly after applying H our qubit is in the state |ψ⟩ = √1 [|0⟩ + |1⟩].
2
Finally, the qubit is measured. In this example the qubit will be left in either |0⟩
or |1⟩ with equal probability.
Let’s consider a more complicated example in Fig. 4.2 involving two qubits.
In this circuit, we first put the top qubit into a superposition state by applying
the Hadamard gate:
| H 100⟩ −→ √ [|0⟩ + |1⟩] ⊗ |0⟩ (4.8)
2
The next operation in the circuit is a controlled-X, or controlled-NOT, gate.
The filled in circle on the top wire is a controlling gate that is conditional on the
state of the top qubit. In this case it will apply X to the second qubit only if the
51
|0⟩ H
|0⟩ X
Figure 4.2: Quantum circuit illustrating two qubits passing through Hadamard
and CNOT gates. Both qubits are measured at the end.
first qubit is |1⟩. X is referred to as a NOT gate because it flips the state from |0⟩
to |1⟩, or vice versa. This can be seen from the definition of X given as0 1 

X =  (4.9)1 0
After applying the controlled-NOT gate, our state is:
1 | ⟩ | ⟩ ⊗ | ⟩ c√ [ 0 + 1 ] 0 −−
N−−O−→T 1√ [|00⟩ + |11⟩] (4.10)
2 2
Finally, we measure the system and collapse the wavefunction. In this case we
will find the system to be in either state |00⟩ or |11⟩ with equal probability.
In quantum computing answers are obtained by collecting measurements
and analyzing the resulting statistics. For example, suppose we asked: How
much does the state |01⟩ contribute to |ψ⟩ after we run the circuit in Fig. 4.2 ? To
answer this question we would need to run the circuit many times and record
each measured final state. Only after collecting a large amount of data points
could we safely conclude that |01⟩ does not contribute to the final state.
In the rest of this chapter we explain the function of gates and further prop-
erties of quantum circuits as they arise. For more extensive review of the basics
of quantum computing, see Ref [52].
52
4.2 Correlation Functions with Path Integrals
In this section we outline the classical algorithm for computing correlation func-
tions using Path Integrals and discuss how to rephrase it as a hybrid algorithm.
Suppose we are interested in the thermal correlation function of two operators,
 and B̂, that evolve according to the Hamiltonian, Ĥ, in the canonical ensemble.
That is
1 ( )
CAB(t) = Tr e−βĤ ÂeiĤt/ℏB̂e−iĤt/ℏ (4.11)Z
where Z is the system partition function and β = 1/kbT . An alternative form
known as the symmetrized correlation(function is defi)ned by [56]:1 ∗
GAB(t) = Tr ÂeiĤτ /ℏB̂e−iĤτ/ℏ (4.12)Z
with complex time τ = t − iℏβ/2. The two are related through a Fourier Trans-
form:
C̃ (ω) = e1/2βℏωAB G̃AB(ω) (4.13)
where we define the Fourier Transform a∫s ∞
G̃AB(ω) ≡ (2π)−1/2 GAB(t)e−iωtdt (4.14)
−∞
and similarly for C̃(t)AB. Typically, it is easier to compute GAB(t) and take the
Fourier Transform to arrive at CAB(t). In this work we assume  and B̂ are op-
erators that are functions only of position. Expressing Eq. (4.12) in a position
basis gives ∫ ∫
1 ∗
GAB(t) = dx dx′A(x)⟨x|eiĤτ /ℏ|x′⟩B(x′)⟨x′|e−iĤτ/ℏ|x⟩Z
Next, we multiply and∫divid∫e by N in eac(h expo)nential:1 ∗ N
GAB(t) = dx dx′ A(x) ⟨x| eiĤ∆τ /ℏ ( |x
′⟩ ) (4.15)Z N
× B(x′) ⟨x′| e−iĤ∆τ/ℏ |x⟩ (4.16)
53
where ∆τ = τ/N. Finally, we insert N − 1 copies of identity in the position basis
between both expande∫d expo∫nentials, yiel
1 ∏
ding
N−1
G (t) = dx+ dx− A(x−
∗
AB N) ⟨x− iĤ∆τ /ℏ −Z j+1|e |x j ⟩ (4.17)
j=0 ∏N−1
× B(x+ ) ⟨x+ |e−iĤ∆τ/ℏN j+1 |x+j ⟩
j=0
Here, x+ = (x+0 , . . . x
+
n ) is a vector defining the forward path and x− = (x−0 , . . . x
−
n )
defines the backward path.∏Let Θ(x
+, x−) be
N−1 ∏N−1
∗
Θ(x+, x−) = ⟨x− |eiĤ∆τ /ℏ|x−⟩ ⟨x+ −iĤ∆τ/ℏ +j+1 j j+1|e |x j ⟩ (4.18)
j=0 j=0
We can re-write Eq. (4.17) as ∫ ∫
1
G + − − +AB(t) = dx dx A(xN)B(xN) sgnΘ |Θ| (4.19)Z
Our goal is to come up with a sampling function so that we can use Monte
Carlo methods to compute the integrals in Eq. (4.19). To do this we define, F, a
normalization co∫nstant∫to |Θ| ∣∣∣a∏sN−1 ∏N−1∗
F = dx+ dx− ∣∣∣∣ ⟨x− |eiĤ∆τ /ℏj 1 |x−j ⟩ ⟨x+ −iĤ∆τ/ℏ +j 1|e |x j ⟩∣
∣
+ + ∣∣∣∣∣ (4.20)j=0 j=0
We can now write GAB(t) in a way that lends itself to Monte Carlo sampling by
multiplying and diving throu∫gh wit∫h F: ( )
1 |Θ|
GAB(t) = F dx+ dx− A(x−N)B(x
+ ) sgnΘ (4.21)
Z
F 〈 〉 N F
= A(x−N)B(x
+
N) sgnΘ W (4.22)Z
where W = |θ|/F.
In every step of the algorithm we need to evaluate matrix elements of the
generic form ⟨x |e−iĤ∆τ/ℏj+1 |x j⟩. This is normally done by diagonalizing the Hamil-
tonian matrix represented in some basis set and then inserting identity in the
54
eighenvalue basis. The current state of the art is to work in a Discrete Variable
Representation (DVR) for which first quantized Hamiltonian matrix elements
have known expressions, saving one from having to explicitly compute any in-
tegrals [57, 40]. However, it is still necessary to perform matrix diagonalization
once the Hamiltonian matrix is constructed.
Suppose we have a d dimensional system in which each dimension requires
M grid points. Our Hamiltonian matrix will then be of size Md × Md. Even
for small grids where M is on the order of 10 this matrix becomes intractable
to diagonalize for a system with anything greater than a few dimensions. Grid
based methods are often limited to low-dimensional systems because of this
exponential scaling.
Quantum computers are able to store information exponentially with in-
creasing number of qubits. This means that if we are able to map the problem of
evaluating matrix elements onto a quantum computer, we can side-step step the
classical exponential scaling issue and open the door to computing correlation
functions of systems beyond a few dimensions. We demonstrate how this can
be accomplished throughout the next few sections.
4.3 Hadamard Test
The simplest quantum algorithm to compute matrix elements of an operator  is
the Hadamard Test [52]. Fig. 4.3 illustrates the circuit executing the Hadamard
Test for the matrix element ⟨ϕ|V̂ |ψ⟩. The circuit contains n qubits initialized to the
state |0⟩⊗, where |0⟩⊗ is shorthand for |0⟩ ⊗ |0⟩ . . .⊗ |0⟩, and one ancialla qubit. We
now work through the action of the Hadamard Test circuit step-by-step. First,
55
|0⟩ H H
|0⟩⊗n U(ψ) U(ϕ) V̂
Figure 4.3: Hadamard Test circuit diagram.
we encounter a Hadamard gate on the ancilla qubit:
| ⟩ ⊗ | ⟩⊗n −→H 10 0 √ [|0⟩ + |1⟩] ⊗ |0⟩⊗n (4.23)
2
Next, we encounter a controlled gate which performs the operation U(ψ) only if
the ancialla state is |1⟩. U(ψ) is meant to represent a generic unitary gate which
takes an arbitrary input state and outputs the state |ψ⟩. The action of the con-
trolled gate is then:
1 | ⟩ | ⟩ ⊗ | ⟩⊗n −cU−−(→ψ) 1 [√ [ 0 + 1 ] 0 √ |0⟩ ⊗ |0⟩⊗n + |1⟩ ⊗ |ψ⟩] (4.24)
2 2
The next gate is also a controlled gate, except the open circle indicates that |0⟩ is
now the controlling state instead of |1⟩. Hence, when the ancilla state is |0⟩, the
gate U(ϕ) is applied to the bottom wire and the state |ϕ⟩ is prepared:
1 [
√ |0⟩ ⊗ |0⟩⊗n + |1⟩ ⊗ |ψ⟩
] −x−c−U−(→ϕ) 1 [√ |0⟩ ⊗ |ϕ⟩ + |1⟩ ⊗ | ]ψ⟩ (4.25)
2 2
Here, xcU(ϕ) is used indicate that the controlling gate is an open circle. The next
controlling gate simply applies the unitary operator V̂ to the appropriate state:
1 [ ] cV 1 [ ]
√ |0⟩ ⊗ |ϕ⟩ + |1⟩ ⊗ |ψ⟩ −→ √ |0⟩ ⊗ |ϕ⟩ + |1⟩ ⊗ V̂ |ψ⟩ (4.26)
2 2
The final gate before we measure the system is another Hadamard gate applied
to the ancilla:
1 [ ] [ ]
√ |0⟩ ⊗ |ϕ⟩ + |1⟩ ⊗ V̂ |ψ⟩ −→
H 1 |0⟩ ⊗ (|ϕ⟩ + V̂ |ψ⟩) + |1⟩(|ϕ⟩ − V̂ |ψ⟩) (4.27)
2 2
56
We now ask what is the probability of measuring the ancialla to be |0⟩? This is
given by the squared modulus of the coefficient,
1 [ ] [ ] [ ]
p = ⟨ϕ| ⟨ | † | ⟩ | ⟩ 1+ ψ V̂ ϕ + V̂ ψ = 2 + ⟨ϕ|V̂ |ψ⟩ + ⟨ϕ|V̂ |ψ⟩∗0 (4.28)4 4
where we used the fact that V̂ is unitary. Similarly, the probability to measure
the ancialla in |1⟩ is
1 [ ] [ ] [ ]
p1 = ⟨ϕ|−⟨ψ|V̂† |ϕ⟩ − V̂ |ψ⟩
1
= 2 − (⟨ϕ|V̂ |ψ⟩ + ⟨ϕ|V̂ |ψ⟩∗) (4.29)
4 4
Taking the difference of the probabilities we are able to extract the real part of
the matrix element:
p0 − p1 = Re⟨ϕ|V̂ |ψ⟩ (4.30)
The probabilities are determined experimentally by running the circuit many
times and counting each time |0⟩ and |1⟩ are observed. The experimentally prob-
abilities are then
state 0 count state 1 count
p0 = , p = (4.31)total count 1 total count
To determine the imaginary part of the matrix element we replace the first
Hadamard gate with Rx(θ = π/2), where Rx(θ) is a rotation about the x-axis of
the Bloch sphere:   cos θ −isin θ 2 2

Rx(θ) =  (4.32)
−isin θ cos θ 2 2
Working through the circuit we find that p0 − p1 = Im⟨ϕ|V̂ |ψ⟩. Again, we can de-
termined the imaginary part of the matrix element by running the correspond-
ing circuit several times and accumulating statistics to experimentally deter-
mine p0 and p1. ⟨ϕ|V̂ |ψ⟩ is of course given by ⟨ϕ|V̂ |ψ⟩ = Re⟨ϕ|V̂ |ψ⟩ + i Im⟨ϕ|V̂ |ψ⟩.
We are interested in the matrix element ⟨x j+1|e−iĤ∆τℏ|x j⟩. The corresponding
Hadamard Test circuit is shown in Fig. 4.4 Here, A is either the Hadamard or
57
|0⟩ A H
|0⟩⊗n U(x j) U(x j+1) e−∆βĤ e−iĤ∆tℏ
Figure 4.4: Hadamard Test circuit with real and imaginary time evolution gates
explicitly shown.
x-rotation gate. Unfortunately, life is not so simple and the circuit in Fig 4.4
cannot be constructed. This is because quantum computers are only capable
of performing unitary evolution, but our circuit includes e−∆βĤ which is a non-
unitary operator. In the next section we show how to work around this problem.
4.4 Probabalistic Imaginary Time Evolution
Imaginary time evolution is accomplished on a quantum computer is by em-
bedding the non-unitary operator in a larger operator which is unitary. Do-
ing so requires adding additional ancialla qubits to the qubits representing our
wavefunction, resulting in an extended system [58, 59]. Probabalistic Imaginary
Time Evolution (PITE) is a method for accomplishing this which uses the mini-
mal number of additional anciall qubits, i.e. one, to create the extended system
[60]. In this section we introduce the PITE circuit and explain how it performs
imaginary time evolution.
√
Let M = m e−∆βĤ0 where we impose the conditions 0 < m0 < 1 and m0 , 1/ 2
so that we may safely perform a Taylor series expansion later. We define the
Hermitian operator √
≡ M + 1 − M
2
Θ arccos √ (4.33)
2 √
which acts on the n qubits systems. If we define κ = sgn(m0−1/ 2), then cosΘ =
58
√ √ √ √
(M + 1 − M2)/ 2 and sinΘ = κ(M − 1 − M2)/ 2.
|0⟩ H W W†
| ⊗nψ⟩ eiκΘ e−iκΘ
Figure 4.5: Exact PITE circuit diagram.
The circuit which performs Exact PITE is given in Fig. 4.5. Here, W is given
by
1 1 −i
W = √   (4.34)2 1 i
Directly before measurement, the overall state prepared by the circuit is
√
|Ψ⟩ = |0⟩ ⊗ M|ψ⟩ + |1⟩ ⊗ 1 − M2|ψ⟩ (4.35)
The state |ψ⟩ undergoes imaginary time evolution when the ancialla qubit is
measured to be in the state |0⟩. In other words, the imaginary time evolved state
which we are interested in obtaining will only be prepared with some probabil-
ity, hence the name Probabalistic Imaginary Time Evolution.
We now walk through the Exact PITE circuit step-by-step. The first two gates
act only on the ancilla qubit:
| ⟩ ⊗ | ⟩ −→H 10 ψ √ [|0⟩ + |1⟩
W 1 [ ]
] ⊗ |ψ⟩ −→ (1 − i)|0⟩ ⊗ |ψ⟩ + (1 + i)|1⟩ ⊗ |ψ⟩ (4.36)
2 2
Next, we apply two controlling gates:
−xc−e
iκΘ 1 [ ]−→ [(1 − i)|0⟩ ⊗ e
iκΘ|ψ⟩ + (1 + i)|1⟩ ⊗ |ψ⟩ ] (4.37)2−iκΘ
−ce−−→ 1 (1 − i)|0⟩ ⊗ eiκΘ|ψ⟩ + (1 + i)|1⟩ ⊗ e−iκΘ|ψ⟩ (4.38)
2
59
The last gate is W† which is applied to the ancilla qubit:
W†−→ 1
[ ]
√ (1 − i)(|0⟩ + i|1⟩)eiκΘ ⊗ |ψ⟩ + (1 + i)(|0⟩ − i|1⟩)e−iκθ ⊗ |ψ⟩ (4.39)
2 2
Since the imaginary time evolved state which we are interested in is associated
with the ancialla state |0⟩, we will focus on that term. Rearranging the above
equation gives,
[ ]
|Ψ⟩ 1= √ |0⟩ ⊗ (1 − i)eiκΘ|ψ⟩ + (1 + i)e−iκΘ|ψ⟩ + . . . (4.40)
2 2
1 [( ) ( )]
= √ |0⟩ ⊗ e+iκΘ + e−iκΘ − i e+iκΘ − e−iκΘ |ψ⟩ (4.41)
2 2
1
= |0⟩ ⊗ [cosκΘ + sinκΘ] |ψ⟩ (4.42)
2
1 ( √ √ )
= |0⟩ ⊗ M + 1 − M2 + M − 1 + M2 |ψ⟩ = |0⟩ ⊗ M|ψ⟩ (4.43)
2
This demonstrates that the Exact PITE algorithm does in fact prepare the imagi-
nary time evolved state M|ψ⟩whenever the ancilla qubit is measure to be in state
|0⟩.
Unfortunately, it is extremely difficult to create a gate which performs the
operation eiκΘ. In practice, we perform a Taylor series in κΘ out to first order in
∆β instead. The details of this expansion are given in Appendix A. This gives us
κΘ = θ0 − Ĥs1∆β + O(∆β2) (4.44)
where [ √ √ ]
θ0 = κarccos (m + 1 − m20 0)/ 2 (4.45)
and √
s1 = m0/ 1 − m20 (4.46)
60
|0⟩ H W Rz(−2θ0) W†
|ψ⟩ ⊗n U (∆β) U (∆β)†RT E RT E
Figure 4.6: Approximate PITE circuit diagram.
We now have eiκΘ ≈ eiΘ0e−iĤs1∆β. The corresponding circuit that executes Ap-
proximate PITE is then given in Fig 4.6. Here, URT E is a Real Time Evolution
(RTE) gate given by e−iĤγ, where γ is a real number. In our case, γ = ∆β. The Rz
gate performs a rotation about the z-axis of the Bloch sphere and is given by
R (θ) = z e−i
θ
 2 0 

 (4.47)0 e+i θ2
The final state of the entire system produced by Approximate PITE right before
measurement is
√ √ 2 m 
 ( )
|Ψ⟩ = |0⟩ ⊗ m (1 − H∆β)|ψ⟩ + |1⟩ ⊗  1 − m2 0 0 0 + H∆β |ψ⟩ + O ∆β2 (4.48)1 − m20
We now consider the probability of measuring the ancialla qubit to be in state
|0⟩. Taking the inner product of the corresponding coefficient and keeping only
terms up to O(∆β) gives ( )
p ∝ m20(1 − 2∆β⟨ψ|Ĥ|ψ⟩) + O ∆β2 (4.49)
up to a normalization constant. What is important to note here is that p increases
as ∆β increases. In other words, working with a smaller ∆β is beneficial in terms
of Approximate PITE. This is good news for us because Path Integral techniques
also depend on ∆β being small so that the Hamiltonian can be Trotterized. Be-
cause of this, Approximate PITE does not impose any further limitations to the
classical Path Integral algorithm. In the next secion we combine Approximate
PITE with the Hadamard Test and develope our final circuit.
61
4.5 Full Circuit
To construct the full circuit we add an extra ancialla qubit to the Hadamard
Test circuit in Fig. 4.3 and replace the e−iĤ∆τ/ℏ gate with the Approximate PITE
circuit. Fig. 4.7 illustrates the full circuit. Here, the gate CPIT E denotes the circuit
for Approximate PITE adapted for two ancialla qubits, which is shown in Fig.
4.8.
|0⟩ A H
|0⟩ CPIT E(∆β)
|0⟩⊗n U(x j) U(x j+1) URT E(∆t)
Figure 4.7: Full circuit for position basis complex time evolution. Here, CPIT E is
the approximate PITE circuit.
One of the key advantages to our method for computing complex-time ma-
trix elements is that initial state preparation is simple. This is because we di-
rectly map position space grid points onto the computational basis, as will be
explained in the next section. Because of this, U(x j) and U(x j+1) can always be
implemented with a set of NOT gates. This is important because initial state
preparation is a major challenge in quantum computing. In general, there is no
|0⟩
|0⟩ H W Rz(−2θ0) W†
| ⊗nψ⟩ U †RT E(∆β) URT E(∆β)
Figure 4.8: Approximate PITE circuit with additional ancilla.
62
guarantee that an efficient algorithm exists to prepare an arbitrary wavefunc-
tion.
A second advantage is that our circuit inherently uses a small time step be-
cause it is designed for use in a Path Integral setting. Small time steps are
good for quantum circuit design because they prevent the need for deep circuits
which risk heavy error accumulation. This is because the way longer time steps
are achieved is by successive application of small time step evolution operators;
each additional applications extends the depth of the circuit.
A final advantage is that URT E can be treated as a black box. This is because
the construction of the final circuit did not depend on the specific implementa-
tion details of URT E. This means different implementations of URT E can easily be
swapped in.
On that note, we still need to decide what the implementation details of URT E
will be in our specific use case. The simplest scheme for first quantized systems
is the split operator approach [61]:
URT E(γ) ≈ e−iT̂γe−iV̂γ (4.50)
This splitting is accurate to first order in γ and involves evolution first under the
potential energy operator then the kinetic energy operator. In the next couple
sections we discus the circuits needed to implement each operator.
4.6 Potential Energy Circuit Design
As an initial test system, we work with an harmonic oscillator potential, V(x) =
1mω2x22 . Our challenge is to design a circuit which allows us to implement
63
e−iV(x)γ. To accomplish this, we first need to decide on a mapping between our
discrete position grid and the computational basis.
Suppose we have a one-dimensional grid that spans −d < x < d and n qubits
at our disposal [62]. If we discretize the grid into evenly spaced pieces and map
each grid point onto the 2n available computational basis states, then the spacing
between each grid point is ∆x = 2d/2n. If we index the lth grid point as xl, then
xl = −d + l ∗ ∆x. Now we substitute the binary expansion for l,
∑n−1
l = k jj 2 (4.51)
− ∑
j=0
to give xl = d + ∆x n−1j=0 k 2 jj . Finally, we pull −d inside the sum to give
∑n−1
xl = ∆x (k jj 2 + α) (4.52)
j=1
where α = −d+∆x/2( x) n . We can now insert this expression for xl into the potential∆
energy evolution operator:
∏n−1
−i 1 mω2 x2e 2 γ = e−iδ(k 2
j
j +α)(k pp 2 +α) (4.53)
j,p=0
where δ = 1mω22 γ. Next, we expand the exponential and collect terms by order
1
of k. This gives e−i 2 mω2 x2γ = T1T2T3, where
2
T1 = e−iδα∏ (4.54)n−1
T = e−i2δα k j2
j
2 (4.55)
∏j=0n−1
j p
T3 = e−iδα k j2 ·kp2 (4.56)
j,p=0
T1 can be implemented with a single phase gate and T2 can be implemented
with n phase gates. The most expensive gate is T3 which requires n2 controlled
64
phase gates. This means the overall resource requirements to implement simple
harmonic oscillator scales as n2. The corresponding quantum circuit diagram
can be found in Fig. S1 of the Supporting Materials of Ref. [63]. The red layer
implements T1, the blue layer implements T2, and the green layer implements
T3.
4.7 Kinetic Energy Circuit Design
Evolution under the kinetic energy operator can be achieved using the Quantum
Fourier Transform (QFT) algorithm [61]:
e−iTγ = QFT†e−iγp
2/2mQFT (4.57)
QFT is a canonical quantum computing algorithm and its implementation can
be found in Ref [47]. This method works by first transforming the computa-
tional basis from a position space representation to momentum space. In the
new basis e−iγp2/2m is a simple diagonal operator which means it can be imple-
mented with a set of phase gates. After applying −iγp2e /2m, we transform back
into the position basis.
Unfortunately, taking the Fourier Transform of eigenstates in the position ba-
sis leads to numerical issues. This results from the fact that position basis states
are Dirac delta functions and their Fourier Transforms are therefore infinitely
delocalized in momentum space. Of course we cannot represent an infinite do-
main using a finite grid and are forced to truncate the Fourier Transformed func-
tion in momentum space, leading to numerical difficulties.
One way to circumvent this problem is to work in a basis of ”smeared” delta
65
functions, called pixel functions, where each function is centered about a unique
grid point [64]. The jth pixel function center about the grid point x j is given by
the following definition( : ) √ n ( )
iπ(x − x ) 1 2∑−1j π(2 j − 1)(x − x )P j(x) = exp − jn−1 cos (4.58)2d 2 d 2d
j=1
Ref. [64] details how the pixel function is derived. Fig. 4.9 is an example of a
pixel function centered over the origin. Because each pixel function is sharply
peaked about one grid point, we can approximate our position basis set by a
pixel basis set, x j → P j, and map the pixel basis onto our computational basis,
P j → | j⟩, without having to make changes to our previous work.
Figure 4.9: Pixel function centered over origin.
Pixel functions have the property that they map to momentum eigenstates
in momentum space, hereafter called k-space [64]. In particular, the kth compu-
tational basis state in k-space encodes the fun(ction )
1 iπkx
ϕk(x) = (2d)− 2 exp (4.59)d
where k is an integer. Each state is mapped to the computational basis in the in
k-space, ϕk → |k⟩. After applying the unitary kinetic energy operator, this state
is transformed according to
|k⟩ → e−iℏ2π2k2γ/dm|k⟩ (4.60)
66
where m is the mass of the particle. This means we only need to apply a phase
gate to each qubit register in k-space to implemented unitary kinetic energy
evolution.
In the beginning of this chapter we said we would only concern ourselves
with near-term algorithms, of which QFT is not. Then we designed a circuit to
perform evolution under the kinetic energy operator which uses QFT. This is
because there is currently no good way to avoid using QFT when dealing with
first quantized methods for dynamics. Thankfully, our circuit does not require
many repeated applications of QFT and it may be possible to execute our circuit
on a real quantum computer in the near future. We acknowledge, however,
that reliable, noise reduced results will most likely not be obtainable until fault
tolerant quantum computing is available.
We now have all the ingredient to simulate a simple harmonic oscillator. In
the next section we evaluate initial results of the method.
4.8 Harmonic Oscillator Data
As an initial test case we compute the position auto-correlation function, Cxx(t),
for harmonic oscillator using a quantum simulator and the procedure detailed
in this chapter. The results are shown in Fig. 4.10. We work in atomic units
with β = m = ω = 1 a.u., and use a grid with d = 10 discretized over n = 7
qubits. Each point was generated using on the order of 106 Monte Carlo steps
and error bars were collected by repeating the numerical experiment 10 times.
For 0 ≤ t < 2 a.u., a Trotter dimension of three was used; for 2 < t ≤ 5 a.u., a
Trotter dimension of four was used.
67
Figure 4.10: Harmonic oscillator position auto-correlation function. Here, N = 3
for the first give points, and N = 4 for the remaining six.
Overall, the simulations is in good agreement with the exact results out to 4
a.u. The last two points are far from the exact results, but moving qualitatively
in the right direction. We suspect that better results could be reached for the last
two points by using a larger N and increasing the number of Monte Carlo point.
To our knowledge this marks the first time a Path Integral calculation has
been successfully performed with the aid of quantum computing. This is impor-
tant because it offers a new and distinct avenue to explore hybrid algorithms for
chemical dynamics and is considerably different than current techniques based
on the variational formulation of the time-dependent Schrodinger equation. We
hope to perform similar calculations using a system bath model and expect our
method to excel for condensed phase systems.
68
4.9 Simulation Notes
In this section we discuss a few tips we found helpful for running simulations.
4.9.1 Partition Function Trick
We can avoid computing the partition function in Eq. (4.21) with a handy trick.
First, we write the position auto-corre(lation function∗ ) as
Tr (x̂eiĤτ/ℏ x̂e−iĤτ /ℏCxx(t) = ) (4.61)
Tr ÎeiĤτ/ℏ Îe−iĤτ∗/ℏ
Notice that if we add the terms in the exponentials of the denominator to-
gether we recover Z. This is because the real time components cancel and leave
behind −β. This means we can think of the partition function as correlating two
identity operators. We can expand the numerator and denominator in a position
bases to get∫ ∫
d∫x+ d∫x− x−∏∏N−1 − iĤ ∏ +∆τ∗/ℏ − + N−1 + −iĤ∆τ/ℏ|x jN j=0 ⟨x j 1|e |x+ j ⟩xN j=0 ⟨x j 1|e ⟩+
dx+ dx− N−1⟨x− |eiĤ∆τ∗/ℏ| −⟩ ∏x N−1 (4.62)⟨x+ |e−iĤ∆τ/ℏj=0 j 1 j j=0 j 1 |x+j ⟩+ +
Following the same logic in Section 4.2 we introduce a sampling function by
multiplying and dividing b∫y F ∫
F d∫x+ d∫x−x−N x+N sgnΘ(|Θ|/F) (4.63)
F dx+ dx−sgnΘ(|Θ|/F)
We define our sampling function W =〈|Θ|/F and〉cancel the F’s in front to give
x− +N〈xNsgn〉ΘC (t) Wxx = (4.64)
〈 〉 sgnΘ W
We can compute sgnΘ W during the same Monte Carlo simulation used to com-
pute the numerator by simply adding each sgn together and dividing by the
69
total number of Monte Carlo simulation steps taken. This trick allows us to
avoid explicitly computing the partition function, saving the need to perform
an additional independent Monte Carlo simulation.
4.9.2 Wavefunction Normalization
Computing the correct real and imaginary components using the Hadamard
Test in conjunction with PITE is a bit tricky owing to the fact that the operator
M does not leave |ψ⟩ normalized. However, the overall wavefunction for the
extended wavefunction will remain normalized throughout the circuit. Here
we outline the correct normalization.
Let Ns be the number of times the success state is measured and N0 and N1
be the number of times the first ancilla qubit is measured in state |0⟩ and |2⟩,
respectively. The ratios needed to compute p0 and p1 are given by
N0 Np 10 = , p1 = (4.65)Ns Ns
The specific matrix element that will be returned by the quantum computer is
⟨x |m e−iĤ∆τj 0 |xk⟩. Dividing by m0 then gives the correct matrix element to be used
in the calculation.
4.9.3 Hashing
Once a matrix element is computed using the quantum computer, it should be
stored in a table using a standard hashing algorithm. This way if the same
matrix element is needed again later in the calculation, we can look it up in
70
the table rather than make another call to the quantum computer. This limits
redundant calls to the quantum computer and helps improve efficiency.
71
CHAPTER 5
QUANTUM COMPUTING WITH SPIN SYSTEMS
In the previous chapter we presented a hybrid algorithm for computing corre-
lation functions derived from a classical Path Integral algorithm. The hybrid
algorithm was designed with first quantized systems in mind and we tested it
by computing the position auto-correlation function of a harmonic oscillator.
In this chapter we show that we can easily adapt the algorithm for correlation
functions of spin systems.
Spin systems provide an important test bed for quantum algorithms because
their Hamiltonians can typically be easily constructed on a quantum computer.
In fact, some models of quantum computing naturally use spin 1/2 systems as
qubits, for example quantum dot qubits [55]. Studying spin systems is also of
interest in the quantum machine learning community because several optimiza-
tion problems can be mapped directly to spin Hamiltonians [65]. In this context,
finding low energy states of the Hamiltonian is the same as finding optimized
solutions to a given problem.
Several hybrid algorithms have been developed to compute dynamics of
spin systems, but none are derived from a Path Integral framework. These
methods are mostly derived from the Time Dependent Variational Schrodinger
(TDVS) equation and use an initial ansatz to parameterize the initial states that
the quantum circuit can prepare [66]. Such ansatz are used because quantum
computers cannot efficiently prepare arbitrary intial states. TDVS algorithms
thus provide the best approximation of the evolved state at a specified time for
a given ansatz. This adds additional error to the calculation because the systems
is now restricted to a sub-space of the full Hilbert space which is spanned by the
72
ansatz. For example, we might limit ourselves to initial states which can be pre-
pared by using only nearest neighbor two-qubit rotations along the x direction.
One key advantage to working with Path Integrals is that we avoid the issue
of initial state preparation altogether. This is because we are only concerned
with computing matrix elements between states which map directly onto the
computational basis. As a results we have no need to work with an ansatz and
do not limit ourselves to a sub-space of Hilbert space.
In the next section we discuss spin systems and why they are important
models for quantum computing. We then derive a hybrid algorithm to compute
correlation functions by following much of the same steps we saw in Chapter
4, Section 2. Next, we discuss how to build a quantum circuit to measure ma-
trix elements of the time-evolution operator. Finally, we present some prelimi-
nary results and conclude by discussing how the calculation might be improved
upon in the future.
5.1 Spin Systems
Consider a system of n spins each of which can either be in a spin-up or spin-
down state. We assume the spins exists on a one-dimensional grid and are
equally spaced from their nearest neighbors. If we only allow for local coupling
and an external magnetic field, the system is described by the transverse-field
Ising model [67]: ∑n ( )
Ĥ = −J σzσz + gσxj j 1 j (5.1)+
j=1
73
where σz and σx and the z and x Pauli matrices, respectively. Both J and g are
real and positive. The first term corresponds to interaction among the spins;
paired spins of the same orientation lower the overall energy while spins of
opposite orientation raise the energy. We use periodic boundary conditions so
that σn+1 = σ1. The second term accounts for the interaction of each spin with
the external magnetic field.
It is worth taking a moment to discuss how quantum computers are used to
compute energy expectation values by examining this simple model. Suppose
we first tried implementing a classical algorithm using diagonalization. The
Hilbert space of our system is 2n because each spin in the chain has two states.
This means the exact representation matrix representation of H is size 2n × 2n.
Classical diagonalization will not always be feasible since the size of matrix
grows exponentially in n.
We can avoid this exponential scaling by using a hybrid algorithm in which
each term is individually measured on a quantum computer, then the total ex-
pected energy is added on a classical computer. We are able to do this because
each Pauli matrix or product of Pauli matrices is a unitary matrix and can thus
be implemented on a quantum circuit. This fact is a crucial reason why spin
systems are the ideal test beds for quantum computers. In general, it is not true
that an arbitrary Hermitian matrix can be decomposed into a sum of unitary
matrices. Even in instances where this is the case, finding an efficient decompo-
sition of the unitary to elements of a universal gate set is not easy. In fact, this
task usually involves search for a decomposition of the unitary into a product
of Pauli matrices. The beauty of spin systems is that they are already written as
products of Pauli matrices!
74
If we measure each term of the Hamiltonian individually on a quantum com-
puter, there will be 2n different circuits to construct: n from the first term, and
n from the second. Each circuit corresponding to the first term will require two
gates and each circuit corresponding to the second term will require one gate.
Hence, the total resource scaling is 2n + n. This means we can compute the
energy using a hybrid algorithm that scales linearly with increasing n for our
nearest neighbor coupling model. Clearly, this is a huge advantage over the
exponentially scaling classical algorithm.
It is important to note that this method does not necessarily give the ground
state energy. If we use the Hadamard Test from Chapter 4, Section 3 to measure
expectation values, then this method simply takes a state given to it and returns
its energy expectation value. The ground state energy will only be returned if
the ground state was the state given to the method. However, we usually do not
know what the ground state is. This is where a classical optimization routine is
used to search the Hilbert space for the ground state [49].
A more complicated spin systems is the XX chain model described by the
following Hamiltonian [68]: ∑n ( )
Ĥ = −J σxσx y yj j+1 + σ jσ j (5.2)+1
i=1
where σy is the y Pauli matrix. This model allows for local spin coupling in
both x and y orientation. We will work with the XX chain model later on when
running simulations to compute correlation functions.
In dynamics we are of course interested in the time evolution operator cor-
responding to a H given as U = e−iHt/ℏ. In the next section we discuss how to
compute U on a quantum computer and a hybrid scheme to evaluate correlation
functions.
75
5.2 Hybrid Algorithm
The steps we take to derive a hybrid algorithm for spin systems are very similar
to those we saw in Chapter 4, Section 2, except we will be working in a finite dis-
crete basis rather than a continuous one. The symmetrized correlation function
between two observables  and B̂ in the canonical ensemble is given as
1 ( )
C −iĤτ +iĤτ
∗
AB(t) = Tr Âe B̂e (5.3)Z
where τ = t−iβ/2 and Z is the canonical ensemble partition function. We work in
atomic units so that ℏ = 1. We compute the trace in the basis of spin eigenstates.
For a single spin we have |σ⟩ where |σ = ±1⟩. For n spins we have |σ1⟩⊗ . . .⊗ |σn⟩
where each |σ j = ±1⟩. We compress the notation by defining |σ j⟩ = |σ1⟩⊗. . .⊗|σn⟩,
where j = f (σ1) . . . f (σn) and 



0 if σ = −1
f (σ) =  (5.4)1 if σ = 1
In other words, j is a binary representation of the state of the entire system.
For example, the state |1⟩ ⊗ |1⟩ ⊗ |−1⟩ = |σ110⟩. Using the equivalent decimal
representation for j, identity is given as
∑M
I = |σ j⟩⟨σ j| (5.5)
j=0
where M = 2n − 1. Our correlation function is now given as:
1 ∑M ∗
C (t) = ⟨σ |Âe−iĤτB̂e+iĤτAB j |σ j⟩ (5.6)Z
j=0
Next, we insert identity:
1 ∑M
C −iĤτ +iĤτ
∗
AB(t) = A(σ j)B(σk)⟨σ |e |σZ j k⟩⟨σk|e |σ j⟩ (5.7)
j,k=0
76
where we assume both observables are Hermitian and only functions of σ. We
now multiply and divide by N in both exponentials:
1 ∑M ∗ ∗
CAB(t) = A(σ j)B(σk)⟨σ j|e−iĤ∆τ . . . e−iĤ∆τ|σ +iĤ∆τk⟩⟨σk|e . . . e+iĤ∆τ |σ j⟩ (5.8)Z
j,k=0
where we explicitly write out each for the N exponential factors. Inserting N − 1
factors of identity between each of the resulting exponential factors yeilds
1 ∏N ∑M ∗
CAB(t) = A(σ+j )B(σ
− )⟨σ+ −iĤ∆τ + − +iĤ∆τ −
Z N lN j
|e |σ j ⟩⟨σl |e |σ ⟩ (5.9)k+1 k k+1 lk
k=1 jk ,lk=0
where we use σ+ for forward paths and σ− for backward paths. The index jk
represents the jth state and kth time step. The same is true for lk.
Just like in the continuous case, we use Monte Carlo sampling methods to
evaluate the sums in Eq. (5.9). To uncover our sampling function first define
F(σ⃗) = ⟨σ+ |e−iĤ∆τ|σ+ ⟩⟨σ− |e+iĤ∆τ∗j j l |σ−k+1 k k+1 l ⟩ (5.10)k
and ∏N ∑M
N = |F(σ⃗)| (5.11)
k=1 jk ,lk=0
Here, σ⃗ is a vector of all states for all time steps of both forward and backward
paths. We can write Eq. (5.9) as
1 ∏N ∑MN + |F|CAB(t) = A(σ j )B(σ− )sgn(F) (5.12)Z N lN N
k=1 jk ,lk=0
Our sampling function is then W = |F|/N . Finally,
N 〈 〉
CAB(t) = A(σ+j )B(σ
−
l )sgn(F) (5.13)Z N N W
Just like in Chapter 4, Section 9, we can avoid computing the normalization
constant and partition function if we write Z as a correlation function of identity
77
operators. We do not go into the details here as the derivation is very similar to
that in Chapter 4, Section 9.
The format of our hybrid scheme will again be to call upon the quantum
computer to evaluate matrix elements of the form ⟨σ |e−iĤ∆τk |σ j⟩ and use a clas-
sical computer to execute the Monte Carlo subroutine. Because each state |σ j⟩
maps directly onto the computational basis we can prepare each states using
only NOT gates. In the next section we discuss how to obtain these matrix ele-
ments using a quantum computer.
5.3 Matrix Elements
We are after matrix elements of the non-unitary operator e−iĤ∆τ. As discussed in
Chapter 4, Section 3, this presents an issue since quantum computers are only
capable of simulating unitary evolution. Methods exist for performing non-
unitary of discrete systems on quantum computers, such as the Quantum Imag-
inary Time Evolution (QITE) algorithm [58]. Unfortunately, we did not have
time to explore QITE and how it might be used in our Path Integral scheme.
Instead we limited ourselves to purely real time correlation functions by setting
β = 0. Under this condition the time evolution operator is of the form e−iĤ∆t.
Substituting the XX chainHamiltonian from Eq. (5.2) gives ∑n
exp  i tJ x x y y + ∆ σ σ + σ σ j j+1 j j 1 (5.14)+
i=1
We would like to write the exponential sum as a product of exponentials, how-
ever σx yj and σ j do not commute for the same index. We instead approximate the
78
product using a Trotter splitting to first order in ∆t:
∏n ( ) ( )
exp (−iH∆t) ≈ exp +i∆tJσx x y yjσ j+1 exp +i∆tJσ jσ j 1 (5.15)+
j=1
Each term in the product is a unitary operator and acts only on two spins at
a time. From a quantum circuit perspective, each unitary will be decomposed
to a two qubit gate and we will have 2n such gates. Fig. 5.1 provides a circuit
x x y y
for executing Eq. (5.15). Here U x e+i∆tJσ jσ j 1 and Uy e+i∆tJσ σ+ j j+1j, j+1 = j, j 1 = .+
|σ⟩1
U x y1,2 U1,2
|σ2⟩
U x y2,3 U2,3
|σ3⟩
U x Uy3,4 3,4
|σ4⟩
U x y4,5 U4,5
|σ5⟩
U x5,6 U
y
5,6
|σ6⟩
Figure 5.1: XX-chain real time evolution circuit.
The circuit consists of two layers: the blue layer performs real time evolution
under coupled Pauli x operators while the red layer performs real time evolu-
tion under coupled Pauli y operators. We assume that the state passed into the
algorithm has been prepared elsewhere.
In the next section we use the our hybrid algorithm to compute correlation
functions for the XX chain models.
79
5.4 Results
In this section we present preliminary results using our method to compute
correlation functions of the XX chain Hamiltonian. We first consider a two spin
systems with J = 1 and ∆t = 0.1. The observable we compute is the expectation
value of the orientation of the first spin along the z direction, ⟨σz⟩. Although this
quantity is not a correlation function, it is useful for evaluating the dynamics of
spin systems. The systems is intially prepared in a domain wall state given by
|↑↓⟩.
Fig. 5.2 compares exact results using matrix diagonalization shown as yel-
low crosses with our hybrid Path Integral Monte Carlo (PIMC) method shown
as blue connected dots. It is clear from these results that the hybrid PIMC
method does not perform well. The only point in agreement with the exact
results is at t = 0 and the error bars are very large. We can make a substantial
improvement by using the fact that the XX chain Hamiltonian commutes with
the total spin along the z direction:
∑n
M = σzi (5.16)
i=1
This means M is a conserved quantity. Instead of sampling the entire Hilbert
space, we can designed our Monte Carlo sampling algorithm to only consider
states that conserve M. We call this method PIMC cons. and represent it by black
connected dots in Fig. 5.2. Now we see our results are in excellent agreement
with the exact results and our errors are too small to see.
Next, we consider a four spin system and again compute the expectation
value of the spin along the z direction of the first spin. Results are shown in
Fig. 5.3. The system is initial prepared in the domain wall configuration |↑↑↓↓⟩.
80
Figure 5.2: Spin expectation value for XX chain with two spins.
Exact results are represented by yellow crosses and results obtained with the
hybrid PIMC cons. method are represented with black connected dots. We see
good agreement with the exact method until around 0.3 a.u. Although we see
that the exact results remain within the error bars of PIMC cons. throughout the
calculation, this is partly due to how large the error bars are.
The reason for the large disagreement with exact results in the four spin
systems is most likely due to the sign problem. Without coupling to a bath or
contribution from the Boltzmann operator, there is nothing damping the com-
plex phase oscillations and achieving convergence becomes increasingly diffi-
cult with time. This is why we are only able to run our calculation for a short
time.
One direction in which this method may see success is coupling to a bath
represented by a collection of spins. Here the bath would help rid us of the sign
problem and allow for longer simulations. Additionally, if each bath ”mode”
81
Figure 5.3: Spin expectation value for XX chain with four spins.
is a spin then the additional cost of the bath is one qubit per bath spin. We
can thus readily expand our system to incorporate a bath without substantially
increasing the resource requirements of the quantum computer.
82
CHAPTER 6
MAVARIC
This chapter is dedicated to MApping VARiable Integration Code, or
MAVARIC, a code developed originally for MV-RPMD simulation but has
been extended to other methods that are derived from the RPMD framework.
The link to download MAVARIC is: https://github.com/ElliotEklund/
MAVARIC.git. At the end of the chapter we outline an improved algorithm for
computing gradients of Γ from the MV-RPMD Hamiltonian.
In this dissertation we only focus on using MAVARIC for MV-RPMD simula-
tions, but tutorials on other simulations can be found in the user manual avail-
able for download via the link above. MAVARIC is designed so that a molecular
dynamics calculation can be split into three manageable phases. These phases
are Monte Carlo, Sampling, and Dynamics. All three phases can be run sequen-
tially, or broken up into pieces depending on the user’s needs. In brief, we
outline what each phase is meant to accomplish in the overall computation.
The first phase, Monte Carlo, is used to establish convergence parameters at
equilibrium. This includes performing bead convergence and energy estimator
convergence. Monte Carlo performs a standard Monte Carlo simulation using
the Metropolis-Hastings algorithm for acceptance/rejection criteria.
The Sampling phase is used after convergence of equilibrium parameters has
been established. During Sampling the user samples decorrelated trajectories
that will eventually be used to perform dynamics calculations. These samples
are generated via Monte Carlo simulation. However, unlike the Monte Carlo
phase, Sampling does not compute the energy estimator. As a result, Sampling
83
is much faster at running a Monte Carlo simulation because it is able to skip the
energy estimator calculation. The sampled trajectories can be histogramed to
check the position centroid distribution.
The final phase is Dynamics. Once a sufficient number of trajectories have
been sampled, we are ready to run a dynamics calculation. During this phase,
we can monitor how well our trajectories are conserving energy. This is use-
ful for deciding an appropriate time step to be used by the integrator during
the simulation. We are also able to calculate position centroid auto-correlation
functions during the Dynamics phase. As mentioned above, the equations of
motion (EOM) are propogated with a fourth-order Adams-Bashforth-Moulton
(ABM) integrator, and the initial steps are found with a fourth-order Runge-
Kutta integrator.
In the following sections, we describe each of these phases in detail. This
is done as a tutorial that works through all the steps necessary to calculate the
position centroid auto-correlation function of a model system. At the end we
provide an algorithm for computing gradients of Γ which has linear scaling in
the number of nuclear beads.
84
6.1 Tutorial Model System and Getting Started
6.1.1 Tutorial Model System
The model system we use in the tutorial is Model 1 of [3]. This is a two state
system with state specific potential energy surfaces
V11(Q) = aQ + c (6.1)
V22(Q) = −aQ (6.2)
and off-diagonal coupling elements V12 = V21 = ∆. The state independent po-
tential is
1
V0(R) = kQ2 (6.3)2
The values of the parameters are c = 0, k = a = 1, and ∆ = 0.1. Other parameters
that are needed to fully define the Hamiltonian are M = β = 1 and N = 4. All
units are in a.u. except for N, which is unitless. Throughout the rest of this
chapter we refer to this model as Model 1.
Upon downloading MAVARIC the potential energy surfaces and coupling
will be configured to Model 1. The input files will also be configured for these
parameters, however, you will be asked to change these parameters throughout
the tutorial. At the end of the tutorial, in Section 5, we will explain how to
change the potential energy surfaces to your system of interest.
85
6.1.2 Input Files
MAVARIC contains five input files through which calculations are requested
and certain model parameters can be specified. These files are found in the di-
rectory InputFiles. It is crucial that the number of lines in each file never change.
Adding or deleting a line will cause the program to crash. In addition, the for-
mat of each line should never change. The format should always be a descrip-
tive phrase, followed by a colon (:), followed by a positive number. There are
only primitive error handling measures in place to detect deviations from this
format - proceed with caution!
Many of the parameters used to request a calculation take binary input.
These parameters use 1 to request the specified calculation, and 0 otherwise.
Here is a list of all binary parameters used by MAVARIC. The parameters are
listed under the file they can be found in.
MonteCarlo Sampling Dynamics
Run MC Run Sampling Run Dynamics
Save PSV Save Sampled Trajectories Read Trajectories
Save MC Data Histogram Positions Check Energy Conservation
Read PSV Read PSV
Read MC Data
6.1.3 Compiling and Running MAVARIC
MAVARIC is intended for Linux and Mac OS X operating systems. MAVARIC is writ-
ten in C++ and will only work with compilers compatible with C++ 11 or higher.
86
MAVARIC provides a makefile for compiling the code which can be found in the di-
rectory MAVARIC. The default compiler in the makefile is set to Intel’s distribution of
the C++ compiler. To change the compiler, modify the second line of the makefile. Just
make sure it is C++ 11 compatible! To compile the code, enter make into the terminal.
6.2 Monte Carlo
The Monte Carlo phase is used to determine equilibrium convergence parameters. You
can request calculations for the Monte Carlo phase using the input file MonteCarlo,
located in the InputFiles directory. After downloading MAVARIC, the file MonteCarlo
should contain exactly seven lines in the following format:
MonteCarlo
Run MC:0
MC Moves:1e5
Estimator Rate:1e2
Save PSV:0
Save MC Data:0
Read PSV:0
Read MC Data:0
In the directory InputFiles, you should also find a file called SystemParameters and
a file called ElecParameters. SystemParameters contains parameters for specifying the
system model. ElecParameters contains parameters that specify the electronic states of
our model. After download, SystemParameters should contain exactly four lines in the
following format:
SystemParameters
Mass:1
87
Beads:4
Temperature:1
MC Step Size:1
ElecParameters should contain exactly two lines in the following format:
ElecParameters
States:2
MC step size:1
The parameters Mass, Beads, and Temperature correspond to M, N, and T in
β = 1/T of the MV-RPMD Hamiltonian. Both M and T are in a.u. States is the number
of electronic states in our model. Both files contain a parameter called MC Step Size
which will be discussed below. Its units are also given in atomic units.
To begin the tutorial, change Run MC:0 to Run MC:1 in the MonteCarlo file. Run the
code after you have made this change. If done successfully, you will see the following
message:
Calculating...
Running Monte Carlo calculations ...
MonteCarlo Results:
System Acceptance Ratio: 25.8661
Electronic Acceptance Ratio: 22.9134
Average Energy: 0.619438
Successfully wrote energy estimator file to Results.
Monte Carlo Simulation Run Time: 0.163655
88
Finished Calculations
Because Monte Carlo simulations are a stochastic process, your System Accep-
tance Ratio, Electronic Acceptance Ratio, and Average Energy will dif-
fer slightly from the example above. Additionally, Monte Carlo Simulation Run
Time will be dependent on the machine executing MAVARIC.
The lines System Acceptance Ratio and Electronic Acceptance Ratio
tell us the percentage of moves that were successfully accepted for the system and elec-
tronic degrees of freedom, respectively. In this case, roughly 26% of the time, our system
moves were accepted, and roughly 23% of the time, our electronic moves were accepted.
Typically, we want both acceptance ratios to be around 50%. If we lower the parameter
MC Step Size in the SystemParameters file, we should see the System Acceptance Ra-
tio increase. Similarly, if we lower the parameter MC Step Size in the ElecParameters
file, we should see the Electronic Acceptance Ratio increase.
Navigate to the SystemParameters file and change MC Step Size to 0.55. Simi-
larly, in ElecParameters, change MC Step Size to 0.48. After making these changes,
run MAVARIC. You should now find acceptance ratios close to those reported below:
System Acceptance Ratio: 49.7047
Electronic Acceptance Ratio: 49.4914
Once we have found good Monte Carlo step sizes, we can work on converging en-
ergy. In the directory Results is a file called energy estimator. This file contains the energy
estimator computed throughout the most recently requested Monte Carlo simulation.
Every time the Monte Carlo phase is executed, energy estimator is overwritten with the
latest results. In 6.1 we plot energy estimator.
89
Figure 6.1: Average Energy throughout a short Monte Carlo simulation.
As one would expect, the average energy fluctuates rapidly at the beginning of the
simulation, and gradually flattens toward the end as our system approaches equilib-
rium. Model 1, and most other models we will work with, take far more than 105
Monte Carlo steps for energy to converge. To change the number of Monte Carlo steps
in a given simulation, modify MC Moves in the file MonteCarlo. For now, increase the
value of MC Moves to 107. The next time we run MAVARIC, the Monte Carlo simula-
tion will take 107 steps.
Typically, we do not need to write the average energy to file after every Monte Carlo
step. For example, writing the average energy after every 103 steps is probably high
enough resolution to understand how the average energy is converging. We can control
how often the average energy is written to file with the parameter Estimator Rate.
Change the value of Estimator Rate to 103. Now when we run the calculation using
MC steps:1e7 and Estimator Rate:1e3, we will write out the average energy
after every 103 steps, giving us a total of 104 points being written to file. After making
these changes, re-run MAVARIC and plot energy estimator. It should look similar to Fig
6.2 below.
90
Figure 6.2: Monte Carlo simulation start from random initial state.
Notice that Fig 6.2 is much flatter than Fig 6.2 toward the end of the simulation.
Now we can be more confident that energy is converging.
After we run a Monte Carlo simulation, our system should be in a state that is much
closer to equilibrium in comparison to its initial state. Unfortunately, every time we
re-run a Monte Carlo simulation, our initial state is randomly generated, and we lose
information about the equilibrium state gained from the previous simulation. Instead
of starting from scratch each time a Monte Carlo simulation is run, MAVARIC allows
the user to save the final state of the system following a Monte Carlo simulation, and
use it as a starting point for later simulations.
To save the information from your Monte Carlo calculation, open the file MonteCarlo
and change Save PSV to 1 and Save MV Data to 1 as well. Make these changes, then
re-run the simulation. You should find the following print statement:
Calculating...
91
Running Monte Carlo calculations ...
MonteCarlo Results:
System Acceptance Ratio: 49.9421
Electronic Acceptance Ratio: 49.9112
Average Energy: 0.575902
Successfully wrote energy estimator file to Results.
Successfully saved PSV to Results.
Successfully saved MC data to Results.
Monte Carlo Simulation Run Time: 14.4884
Finished Calculations
This indicates that the phase space variables (PSV) and other Monte Carlo data were
successfully saved after the simulation. The directory Results will now contain two new
files called PSV and mc data. These files contain the saved data. Fig 6.3 below is a plot
of the average energy for the simulation we just ran. We will say more about it in a
moment.
To read in data that has been saved, change Read PSV to 1 and Read MC Data
to 1 in the MonteCarlo file. Make these changes, then re-run the MAVARIC. If done
successfully, you should see the following message:
Calculating...
MonteCarlo Results:
92
Figure 6.3: Average Energy throughout a short Monte Carlo simulation.
Successfully read PSV file from Results.
Successfully read MC datat from Results.
System Acceptance Ratio: 49.9243
Electronic Acceptance Ratio: 49.9717
Average Energy: 0.550436
Successfully wrote energy estimator file to Results.
Successfully saved PSV to Results.
Successfully saved MC data to Results.
Monte Carlo Simulation Run Time: 14.7865
This indicates that the files were successfully read. Let’s look at the energy estimator
plot for the simulation that we just ran in Fig. 6.4 below.
The fluctuations in Fig 6.4 are much smaller than Fig 6.3 and we no longer see mas-
sive fluctuations at the start of the simulation. This is because after running the simula-
93
Figure 6.4: Monte Carlo simulation starting from saved initial state.
tion that generated Fig 6.3, we saved the final phase space variables and energy. Then,
when we ran the simulation that produced Fig 6.4, we first read in those saved phase
space variables and energy, then started the simulation. This means our starting point
was much closer to equilibrium for the simulation that produced Fig 6.4 in comparison
to Fig 6.3.
This demonstrates how running a Monte Carlo simulation can be broken up into
chunks. If you find that your Monte Carlo simulation was not long enough to converge
energy, you can simply start another simulation at the point where you left off, running
the calculation for longer this time. Be careful when reading in files. If no files were
saved previously, then trying to read in PSV and mc data will cause an error.
Once you are satisfied with how energy is converging, you can increase the number
of beads in the SystemParameters file and repeat the process outlined in this section. This
is typically done for a range of bead numbers to establish bead convergence.
94
6.3 Sampling
The Sampling phase is used to generate a distribution of trajectories that will eventually
be used in the Dynamics Phase. The trajectories are sampled using a Monte Carlo sim-
ulation, however, the energy estimator is not calculated during Sampling because it is
assumed that the average energy has already been converged. The step size of system
and electronic moves are set by the parameter MC Step Size in the input files Sys-
temParameters and ElecParameters, respectively. Requesting calculations for Sampling
is done via the input file Sampling. After installation, the Sampling file should contain
exactly seven lines in the following format:
Sampling
Run Sampling:0
Number of Trajectories:1e4
Decorrelation Length:1e2
Save Sampled Trajectories:0
Histogram Positions:0
Number of Bins:300
Read PSV:0
To use Sampling, simply change Run Sampling to 1. In the current example, 104
trajectories will be sampled. Between each new trajectory that is sampled, we perform
an additional 102 Monte Carlo steps to ensure the trajectories are sufficiently decorre-
lated. To change the number of sampled trajectories or decorrelation length, modify
the parameters Number of Trajectories and Decorrelation Length, respec-
tively.
Let’s try running a Sampling calculation. In the MonteCarlo file, configure the pa-
rameters in the following way:
95
MonteCarlo
Run MC:1
MC Moves:1e6
Estimator Rate:1e3
Save PSV:0
Save MC Data:0
Read PSV:0
Read MC Data:0
In the Sampling file, only change Run Sampling to 1. After making these changes,
run the code. If everything worked correctly, you should see the following output:
Calculating...
Running Monte Carlo calculations ...
MonteCarlo Results:
System Acceptance Ratio: 50.0492
Electronic Acceptance Ratio: 49.9594
Average Energy: 0.6922
Successfully wrote energy estimator file to Results.
Monte Carlo Simulation Run Time: 1.4755
Running Sampling ...
Sampling Run Time: 0.509911
96
Finished Calculations
MAVARIC began the program by running a Monte Carlo simulation, producing
the corresponding Monte Carlo output messages when it finished. Then, it entered a
Sampling calculation, also producing an output message upon completion.
Let’s run another Sampling calculation, but this time, we will generate a histogram
of the sampled position centroids. To do this, open the Sampling file and change the
Histogram Positions parameters to 1. Additionally, change Number of Tra-
jectories to 1e5. We can control the number of bins our histogram uses by mod-
ifying the Number of Bins parameter. For now, we will leave it at 300 bins. After
making these changes, re-run the simulation. If done correctly, the Sampling output
message should now include a histogram statistics report that looks like the following:
--- Histogram Statistics ---
Mode: -0.342636
Average: -0.00551759
Standard Deviation: 1.24748
Skew: 0.0284287
To view the histogram, plot the file pos histogram in the Results directory. Fig. 6.5
below is an example of what the histogram should look like. Of course, increasing the
number of trajectories will produce a smoother histogram.
The Sampling calculation we have done so far have been have been after first run-
ning a Monte Carlo calculation. The final state of our system after the Monte Carlo
calculation is what Sampling uses as its initial state. If equilibrium was successfully es-
97
Figure 6.5: Histogram using 300 bins with 105 Trajectories.
tablished during Monte Carlo, than all our sampled trajectories produced during Sam-
pling should follow the equilibrium Boltzmann distribution.
Instead of starting with a Monte Carlo simulation, we can read in phase space vari-
ables that have been saved to the PSV file in the Results directory. This way, we do
not need to run Sampling immediately after Monte Carlo. If we have saved our phase
space variables after a Monte Carlo simulation, we can read them back in at a later time
and start MAVARIC at Sampling. To read in saved phase space variables, change Read
PSV to 1.
In addition, we can save trajectories that we have sampled. To do this, change Save
Sampled Trajectories to 1. Now, after we run a sampling calculation, all the tra-
jectories will be saved to the directory Results/Trajectories.
As a demonstration, run a Monte Carlo simulation and select to save both PSV and
mc data. Do this without running Sampling; i.e., set Run Sampling to 0. If done
successfully, you should see the following message:
98
Calculating...
Running Monte Carlo calculations ...
MonteCarlo Results:
System Acceptance Ratio: 49.8436
Electronic Acceptance Ratio: 50.1265
Average Energy: 0.622538
Successfully wrote energy estimator file to
Results.
Successfully saved PSV to Results.
Successfully saved MC data to Results.
Monte Carlo Simulation Run Time: 1.46195
Finished Calculations
Now, run Sampling without doing a Monte Carlo simulation first; i.e. set Run MC
to 0 and Run Sampling to 1. In addition, set Save Sampled Trajectories and
Read PSV both to 1. If done successfully, you will see the following message:
Calculating...
Running Sampling ...
Successfully read PSV from Results.
--- Histogram Statistics ---
99
Mode: 0.0924216
Average: 0.000271159
Standard Deviation: 1.24538
Skew: -0.014498
Successfully saved sampled trajectories to Results/Trajectories.
Sampling Run Time: 11.1057
Finished Calculations
The directory Results/Trajectories should now contain four files titled P, Q, xelec, and
pelec. These are the sampled trajectories that you have just saved. Saving sampled
trajectories allows the user to read them in later when working with the Dynamics
phase. That way, work can be done in chunks and one does not need to go directly
from Sampling to Dynamics.
6.4 Dynamics
The Dynamics phase calculates the position centroid auto-correlation function within
the MV-RPMD framework. EOM are propogated using a 4th order ABM predictor-
corrector integration scheme. Requesting calculations for Dynamics is done via the
input file Dynamics. After download, the Dynamics file should contain exactly six lines
in the following format:
Dynamics
Run Dynamics:0
Run Time:5
Time Step:0.01
100
Read Trajectories:0
Check Energy Conservation:0
Conservation Tolerance:0.1
To run Dynamics, change Run Dynamics to 1. The parameter Run Time specifies
the total time in atomic units for Dynamics to run. Time Step specifies the size of the
time step, dt , used by the integrator.
As a first example, configure the MonteCarlo, Sampling, and Dynamics input files to
match the following.
Monte Carlo
Run MC:1
MC Moves:1e7
Estimator Rate:1e4
Save PSV:0
Save MC Data:0
Read PSV:0
Read MC Data:0
Sampling
Run Sampling:1
Number of Trajectories:1e4
Decorrelation Length:1e2
Save Sampled Trajectories:0
Histogram Positions:0
Number of Bins:300
Read PSV:0
101
Dynamics
Run Dynamics:1
Run Time:5
Time Step:0.005
Read Trajectories:0
Check Energy Conservation:0
Conservation Tolerance:0.1
Run MAVARIC with this set of configurations. You should see the following output
after MAVARIC finishes:
Calculating...
Running Monte Carlo calculations ...
MonteCarlo Results:
System Acceptance Ratio: 49.9658
Electronic Acceptance Ratio: 50.0111
Average Energy: 0.562923
Successfully wrote energy estimator file to Results.
Monte Carlo Simulation Run Time: 16.0427
Running Sampling ...
Sampling Run Time: 0.570696
102
Running Dynamics ...
Successfully wrote pos auto corr to Results.
Dynamics Run Time: 65.9003
Finished Calculations
After the calculation has finished, the Results directory will contain a file called
pos auto corr. This file contains the position centroid auto-correlation function. In Fig.
6.6 is a plot of the auto-correlation function for the simulation we just ran.
Figure 6.6: Position Centroid auto-correlation function using 104 trajectories.
To run a calculation using more or less trajectories, change the parameter Number
of Trajectories in Sampling.
MAVARIC provides a way to monitor how well trajectories are conserving energy
throughout the course of a calculation. If we change Test Energy Conservation
to 1 in the Dynamics file, MAVARIC will now check how well energy is being con-
served. The parameter Conservation Tolerance allows us to specify the tolerance
103
to which we want each trajectory to conserve energy. For example, when Conserva-
tion Tolerance is set to 0.1, this means a trajectory will be considered ”broken”
(not conserving energy) if its initial energy and current energy differ by more than 10%.
Once a trajectory is broken, MAVARIC will stop propagating it and move on to the next
trajectory. At the end of the Dynamics calculation, the percentage of broken trajectories
is reported.
As a demonstration, run MAVARIC using the same input file configurations as pre-
viously, but with Test Energy Conservation set to 1. If done correctly, the Dy-
namics output message should be similar to the following:
Running Dynamics ...
Percentage of Trajectories Broken: 30
Dynamics Run Time: 53.209
This indicates that 30% of our trajectories failed to conserve energy to within 10%
of their initial energy. A proper Dynamics simulation should aim for a much lower
percentage of broken trajectories. Lowering Time Step will decrease the number of
trajectories that fail to conserve energy.
Finally, if sampled trajectories were saved during the Sampling phase, they can be
read in during the Dynamics phase. To do this, change Read Trajectories to 1.
This will read in trajectories stored in the directors Results/Trajectories. If no trajectories
have been saved, then trying to read them in will cause an error. Being able to read
in trajectories means we do not need to run MonteCarlo, Sampling, and Dynamics in
sequence. If sampled trajectories have been saved, then we can skip MonteCarlo and
Sampling, and jump straight to Dynamics.
104
6.5 Coding Potentials
So far, this tutorial has only worked with one model that was pre-configured upon
download. To change that model, you will need to change the potential energy surfaces
and coupling. Making these changes is done in the file Potentials.h which is found
in the directory Potentials. Potentials.h is a special file type called a C++ header file.
This is the only C++ file a user who wishes to use MAVARIC as a black box device will
have to worry about modifying.
Potentials.h makes use of six functions that need to be edited to fully specify a
model. Throughout this section, we will discuss those functions in detail using Model
I as an example. At the end of the section, we will consider a more complicated model
and the changes that would need to be made to implement it.
Before discussing the different functions, we need to become somewhat familiar
with two important C++ features that Potentials.h makes use of. These are macros
and valarrays.
6.5.1 Macros
Typically, our potential energy surfaces and couplings contain parameters that need to
be specified by the user. For example, if our system independent potential energy is the
following,
1
V (Q) = kQ20 (6.4)2
then we need to be able to specify the parameter k in Potentials.h. We specify pa-
rameters using a C++ macro. The best way to understand a macro is to see an example.
Toward the top of Potentials.h are the following three lines:
105
#define DELTA 0.1
#define A 1.0
#define D 0.0
Each of these defines a macro. Any new macro that you define must be placed
before the line struct Potentials{ . The general format for defining a macro is:
#define NAME value.
It is common practice for the name of the macro to be in all caps. Once the macro
is defined, you can use it anywhere in Potentials.h by typing its name. You should
think of using a macro the same way you would a regular variable.
The only parameters that do not need to be defined as macros are the number of
beads, number of electronic states, the product of the number of beads with number of
electronic states, and mass. This is because Potentials.h always has access to these
variables. They are represented by the following variables:
int num beads
int num states
int elec size
double mass
These variables are defined in the lines immediately following struct Poten-
tials{ . They can be used in any of the functions inside Potentials.h.
106
6.5.2 Valarray
The standard vector and array classes in C++ lack support for many vector features that
are convenient when working with mathematical operations. There is a lesser known
C++ vector class called valarray which does support these features. Potentials.h
makes use of valarray when working with a vector of bead positions. We will now
provide a brief discussion of how to work with a valarray, however, more information
can be found at http://en.cppreference.com/w/cpp/numeric/valarray.
Without using valarray, mathematical operations on vectors in C++ would have to
be performed with for loops. For example, element-wise multiplication between two
vectors v and w would involve looping over all the elements in v and w and multiplying
them one-by-one. If the two vectors are valarrays, we can make use of the element-wise
multiplication function that valarray supports. The syntax for element-wise multiplica-
tion is v * w. As another example, let’s consider summing all the elements in a vector.
The valarray function for this operation is called sum(). To sum the elements of a valar-
ray v, we simply type v.sum().
In general, if there is a mathematical operation you want to perform over all the
elements of a vector, valarray should have a function for doing so. As a final example,
let’s compute the dot product between v and w. To do this, we need to combine two
operations, element-wise multiplication and summation. The syntax we would use is
(v * w).sum().
Vallary also supports array slicing - the ability to access chunks of an array at one
time. For example, if we have an valarray x with ten elements, we can access the first
five elements using the syntax x[std::slice(0,5,1)]. This syntax means starting
at position 0, grab the first 5 elements in intervals of 1. In general, the syntax is
x[std::slice(a,b,n)]. This means starting at position a, return the first b ele-
ments in intervals of n.
107
Now that we have covered the two main data types used by Potentials.h, we
will turn to discussing the six functions that Potentials.h uses.
6.5.3 State Independent Potential
The state independent potential energy of a model is defined in the function V0. After
download, V0 should look like the following:
inline double V0(const std::valarray<double> Q){
return 0.5*(Q*Q).sum();
}
∑
Mathematically, V0 corresponds to the term Nα=1 V0(Qα) in the MV-RPMD Hamiltonian.
V0 takes a valarray of b∑ead positions as its argument, Q, and returns a single value, the
energy corr∑esponding Nα=1 V0(Qα). In Model I, V0 is a harmonic oscillator, which means
V0 returns N 1 2α=1 2 Qα.
6.5.4 Gradient of State Independent Potential
The gradient of the state independent potential energy with respect to bead position is
defined in the function dV0 dQ. After download, dV0 dQ should look like the follow-
ing:
inline std::valarray<double> dV0 dQ
(const std::valarray<double> Q){
return Q;
}
108
Mathematically, dV0 dQ corresponds to ∇QV0. We need this term to calculate the
forces used in the integration scheme. dV0 dQ takes a valarray of bead positions as
its argument, Q, and returns a valarray corresponding to ∇QV0. The elements of the
returned valarray are (∂V0(Q1) , ∂V0(Q2) , . . . , ∂V0(QN )∂Q ∂Q ∂Q ).1 2 N
In Model I, V0 is a harmonic oscillator. Mathematically, its derivative with respect
to a particular bead’s position is ∂V0(Qα)∂Q = Qα. The valarray returned by dV0 dQ will beα
(Q1,Q2, . . . ,QN).
6.5.5 State Specific Potential Energy Surfaces
The state specific potential energy surfaces are defined in the the function Velec. After
download, Velec should look like the following:
inline void Velec(const std::valarray<double> &Q,
std::valarray<double> &Vout){
//state 1
Vout[std::slice(0,num beads,1)] = A*Q + D;
//state 2
Vout[std::slice(num beads,num beads,1)] = -A*Q;
}
Mathematically, Velec corresponds to the diagonal Vnn(Qα) elements of the M Ma-
trix in Eq. (3.18). Velec takes two valarrays as arguments. The first argument, Q, is a
valarray of bead positions: (Q1,Q2, . . . ,QN). The second argument, Vout, corresponds
109
to the evaluation of Vnn over Q. Vout is a valarray of length elec size. In Model I, the
state specific potential energy surfaces are
V11(Qα) = aQα + d (6.5)
V22(Qα) = −aQα (6.6)
Velec needs to evaluate both V11 and V22 over every element of (Q1,Q2, . . . ,QN). Let’s
think about how many computations we will need to do during this process. The num-
ber of beads in our system is num beads, and the number of states is num states. This
means we need to perform num beads * num states number of computations. V11
will be computed num beads times for each Qα, and V22 will be computed num beads
times for each Qα as well.
The first num beads calculations will be to evaluate V11(Qα) over all num beads
positions and store the result in Vout. This means the first num beads elements of
Vout are now (V11(Q1),V11(Q2), . . . ,V11(QN)).
Next, we evaluate state V22(Qα) over all num beads positions and store
the result in Vout. Similarly, the second num beads elements of Vout are
now (V22(Q1),V22(Q2), . . . ,V22(QN)). In total, Vout will contain the elements
(V11(Q1),V11(Q2), . . . ,V11(QN),V22(Q1),V22(Q2), . . . ,V22(QN)).
Velec uses valarray’s slicing features to store elements in Vout. To access the first
num beads elements of Vout, we write Vout[std::slice(0,num beads,1)]. To
access the next num beads elements. we write Vout[std::slice(num beads,num beads,1)].
6.5.6 Gradient of State Specific Potential Energy Surfaces
The gradient of the state specific potential energy surfaces with respect to bead position
is defined in the function dVelec. After download, dVelec should look like the fol-
lowing:
110
inline void dVelec(const std::valarray<double> &Q,
std::valarray<double> &Vout){
//state 1
Vout[std::slice(0,num beads,1)] = A;
//state 1
Vout[std::slice(num beads,num beads,1)] = -A;
}
Mathematically, dVelec corresponds to (∇QV11,∇QV22). The layout of the function
is identical to Velec, except the elements of Vout are now ∂Vnn(Qα)∂Q . For Model I, weα
have
∂V11(Qα)
= a (6.7)
∂Qα
∂V22(Qα)
= −a (6.8)
∂Qα
The elements of Vout are (V ′ ′ ′ ′ ′ ′11(Q1),V11(Q2), . . . ,V11(QN),V22(Q1),V22(Q2), . . . ,V22(QN)).
6.5.7 Potential Energy Coupling
The potential energy coupling is defined in the function V couple. After download,
V couple should look like the following:
inline double V couple(const double &Q, int state1, int state2) {
return DELTA;
}
111
Mathematically, V couple corresponds to the off-diagonal potential energy cou-
pling, Vn,m(Qα), where n , m. This function takes three arguments: Q, state1, and
state2. The first argument, Q, is the position of a single bead, Qα. Note that this is
not a valarray. The next two arguments indicate which two states we are calculating the
coupling element between. For example, state1 = 0 and state2 = 1 would corre-
spond to the element V1,2(Qα). Notice that the C++ syntax for the states and the indices
of the mathematical description differ by one. This is because C++ starts its indexing at
zero! Keep this in mind when working with Potentials.h.
In Model I, we have a two state system with constant coupling. As a consequence,
we do not need to worry about which states we are computing the coupling between
since they all have the same coupling. Below, we will consider a case where coupling
is not constant and does depend on which states we are considering. In Model I, the
coupling term DELTA is defined as a macro.
6.5.8 Derivative of Potential Energy Coupling
The derivative of the potential energy coupling with respect to bead position is defined
in the function dV couple. Following download, dV couple should look like the code
below:
inline double dV couple(const double &Q, int state1, int state2) {
return 0;
}
Mathematically, dV couple corresponds to ∂Vn,m(Qα)∂Q , where n , m. This functionα
takes three arguments, which are the same as arguments as V couple: Q, state1, and
112
state2. The layout of dV couple is identical to V couple, except now we are dealing
with derivatives. In Model I, which uses constant coupling, dV couple returns 0.
It is important to recognize that dV couple is returning a number and not a valar-
ray. For reasons that are not important to using MAVARIC, it is easier to work with the
derivative of the potential energy coupling with respect to a single bead rather than the
entire gradient.
6.5.9 Complicated Model
Let’s walk through one more example using a significantly more complicated model.
We will use a three state model with non-constant coupling. Our state independent
potential is
V0(Qα) = c 4 21Qα + c2Qα (6.9)
The three state dependent terms are
(
V (Q ) = D 1 − e−B1(Qα−R )1) 211 α 1( (6.10)
V (Q ) = D 1 − e−B2(Qα−R2))222 α 2( ) (6.11)
V33(Qα) = D3 1 − e−B3(Qα−R3) 2 (6.12)
Finally, the coupling terms are
2
V12(Qα) = V21(Qα) =A e−g12(Qα−Q12)12 (6.13)
2
V13(Qα) = V (Q ) =A e−g13(Qα−Q13)31 α 13 (6.14)
V (Q ) = V (Q ) =A e−g23(Qα−Q
2
23)
23 α 32 α 23 (6.15)
Before touching any code, we need to compute the necessary derivatives. The deriva-
tive of the state independent potential is
∂
V (Q ) = 4c Q3 + 2c Q (6.16)
∂Q 0 α 1 α 2 αα
113
The derivatives of the state dependent terms are
∂
V (Q ) = −2B D eB1(R1−Qα) (eB1(R1−Q)11 α 1 1 − 1) (6.17)
∂Qα
∂
V (Q ) = −2B D eB2(R2−Qα)22 α 2 2 (eB2(R2−Q) − 1) (6.18)
∂Qα
∂
V (Q ) = −2B D eB3(R3−Qα) (eB3(R3−Q)33 α 3 3 − 1) (6.19)
∂Qα
Finally, the derivatives of the coupling terms are
∂ ∂ 2
V12(Qα) = V21(Qα) = − 2g12A12(Qα − Q12)e−g12(Qα−Q12) (6.20)
∂Qα ∂Qα
∂ ∂ 2
V (Q ) = V −g13(Qα−Q13)
Q 13 α Q 31
(Qα) = − 2g13A13(Qα − Q13)e (6.21)
∂ α ∂ α
∂ ∂
V (Q ) = V (Q ) = − 22g A (Q − Q )e−g23(Qα−Q23)23 α 32 α 23 23 α 23 (6.22)
∂Qα ∂Qα
A good place to get started coding up our model is to define all our macros - there
are a lot in this model.
#define C1 1.0
#define C2 2.0
#define D1 0.5
#define D2 0.3
#define D3 0.2
#define B1 10.0
#define B2 20.0
#define B3 30.0
#define R1 1.0
#define R2 1.0
#define R3 1.0
114
#define A12 17.0
#define A13 20.0
#define A23 23.0
#define G12 1.0
#define G13 1.0
#define G23 1.0
#define Q12 3.0
#define Q13 3.0
#define Q23 3.0
Remember, these need to be defined before the line struct Potentials{ .
Please note that the values we’ve assigned to the macros above are completely arbi-
trary and only for the sake of demonstration.
To code the state independent potential energy surface, we could do the following:
inline double V0(const std::valarray<double> Q){
return B1*(Q*Q*Q*Q).sum() + B2*(Q*Q).sum();
}
For the gradient of the state independent potential energy surface, we would write,
inline std::valarray<double> dV0 dQ
(const std::valarray<double> Q){
return 4.0*C1*(Q*Q*Q) + 2.0*C2*(Q);
115
}
Now we will concentrate on the state dependent potential energy surfaces. The
code should look like the following:
inline void Velec(const std::valarray<double> &Q,
std::valarray<double> &Vout){
//state 1
Vout[std::slice(0,num beads,1)] =
D1*pow(1.0 - exp(-B1*(Q - R1)) , 2);
//state 2
Vout[std::slice(num beads,num beads,1)] =
D2*pow(1.0 - exp(-B2*(Q - R2) , 2);
//state 3
Vout[std::slice(2*num beads,num beads,1)] =
D3*pow(1.0 - exp(-B3*(Q - R3) , 2);
}
The gradient of the state dependent potential energy surfaces will be very similar to
the code we just wrote:
inline void dVelec(const std::valarray<double> &Q,
std::valarray<double> &Vout){
//state 1
116
Vout[std::slice(0,num beads,1)] =
-2.0*B1*D1*(exp(B1*(R1-Q)))*(exp(B1*(R1 - Q)) - 1.0);
//state 2
Vout[std::slice(num beads,num beads,1)] =
-2.0*B2*D2*(exp(B2*(R2-Q)))*(exp(B2*(R2 - Q)) - 1.0);
//state 3
Vout[std::slice(2*num beads,num beads,1)] =
-2.0*B3*D3*(exp(B3*(R3-Q)))*(exp(B3*(R3 - Q)) - 1.0);
}
Finally, we code up the coupling terms. Now that we have non-constant coupling,
we need to use if statements to decide which two states we are evaluating the coupling
between. Here is an example of one way we can achieve this.
inline double V couple(const double &Q, int state1, int state2) {
if((state1==0 state2==1) || (state1==1 state2==0)){
return A12*exp(-G12*pow(Q - Q12 , 2));
}
else if((state1==0 state2==2) || (state1==2 state2==0)){
return A13*exp(-G13*pow(Q - Q13 , 2));
}
else {
return A23*exp(-G23*pow(Q - Q23 , 2));
}
}
117
For the derivative of the coupling terms, we follow the same format as in V couple,
but substitute in the derivatives.
inline double dV couple(const double &Q, int state1, int state2) {
if((state1==0 state2==1) || (state1==1 state2==0)){
return -2.0*G12*A12*(Q - Q12)*exp(-G12*pow(Q - Q12 , 2));
}
else if((state1==0 state2==2) || (state1==2 state2==0)){
return -2.0*G13*A13*(Q - Q13)*exp(-G13*pow(Q - Q13 , 2));
}
else {
return -2.0*G23*A23*(Q - Q23)*exp(-G23*pow(Q - Q23 , 2));
}
}
Remember that after you make any changes to Potentials.h, you need to re-
compile MAVARIC. Please note that the code demonstrations above were written for
the sake of demonstration and are not well optimized. A more efficient implementa-
tion would be to define a single macro any time a series of constants are multiplied in a
function.
6.6 Derivative Algorithm
In this section we introduce an algorithm which speeds up the calculation of deriva-
tives of the nuclear - electronic coupling term, Γ, in the MV-RPMD Hamiltonian from
O(N2) to O(N), where N is the number of beads. This algorithm works for both deriva-
tives with respect to nuclear beads and electronic beads. Throughout this section we
118
assume systems are one-dimensional with respect to nuclear spatial coordinates. We
also assume the number of nuclear and electronic beads is the same.
To begin, we consider the standard protocol for computing derivatives of Γ with
respect to nuclear beads and show this method scales as O(N2). Recall the definition of
Γ: ∏N
Γ = C (xα,pα)M (Qα) (6.23)
α=1
where C andM are defined in Eqs. (3.17) and (3.18), respectively. Both C andM are k×k
matrices, where k is the number of electronic states. We are interested in the gradient of
Γ with respect to nuclear beads: ( )
∇ ∂Γ ∂ΓQΓ = , . . . , (6.24)
∂Q1 ∂QN
Our first step in computing ∇QΓ i to(compute the gr)adient of M:
∇ M ∂M1 ∂MNQ = , . . . , (6.25)
∂Q1 ∂QN
where we use Mα = M(Qα) as well as C(xα,pα) for short. Each term in Eq. (6.24) can
now be computed in the following way:
∂Γ ( )
= C1 ∇QM 1 . . .CPMN (6.26)∂Q1
... ( (6.27)∂Γ
= C1M1 . . .CN ∇QM
)
N (6.28)∂QN
It is easy to verify that each term requires 2N − 1 matrix multiplications to compute.
Since there are N terms in the gradient that we need to compute, we have in total N(2N−
1) matrix multiplications. Finally, the cost of each matrix multiplication is k3, giving us
a cost of Cquad = N(2N − 1)k3. For a given k, we say this method scales as O(N2).
Now we consider a method which uses a recursive algorithm to compute strings of
matrix products that are eventually used in the final calculation of ∇QΓ. Once again, we
119
start by computing ∇QM as in Eq. (6.25). Next, we define two vectors recursively called
left string array, l, and right string array, r. Left string array is defined by:
l1 = C1 (6.29)
l j = l j−1M j−1C j (6.30)
Right string array is defined by:
rP = I (6.31)
r j = C jM jr j+1 (6.32)
The jth element of ∇QΓ can now be computed using
(∇ ) ( )QΓ j = l j ∇QM j r j (6.33)
Each term in l, except l1, requires two matrix multiplications. Hence, l requires 2(N − 1)
matrix multiplications to construct. The same is true for r, except rN is the term which
does not require matrix multiplications. Combined, l and r use 4(N − 1) matrix mul-
tiplication. Finally, each term in ∇QΓ requires two matrix multiplications to compute,
giving us a total of 4(N − 1) + 2N matrix multiplications. The cost the algorithm is
Clin = (4(N − 1) + 2N)k3. For a given k, we say this method scales as O(N).
We have shown how to reduce the problem of computing derivatives of Γ with re-
spect to nuclear beads from quadratic to linear scaling. The key to doing so involved
recycling matrix products by defining l and r recursively. For example, the product
CN MN is only computed once and then stored in rN−1. In the quadratic scaling algo-
rithm, this product is re-computed for each new ∂Γ∂Q , except for j = N. Although it isj
custom to drop prefactors when giving complexity scaling, it is interesting to note that
both methods of a k3 dependence in their prefactors.
This method can also be used to computed derivatives with respect to electronic
120
variables. Consider ∇xΓ given as
  ∂Γ x . . .
∂Γ

11 x1k
 . . 

∇xΓ = 
.. .. 
 (6.34)
∂Γ
x . . .
∂Γ
N1 xNk
and ∇pΓ given as  
 ∂Γ . . . ∂Γ p11 p1k
∇ 
.. .. pΓ =  . .∂Γ . . . ∂Γ 
 (6.35)
pN1 pNk
Each gradient has k × N terms and each term equires 2N − 1 matrix multiplications to
compute, hence there are 2kN(2N − 1) matrix multiplications required to construct both
gradients. Each matrix multiplication costs k3 giving us a final cost of 2kN(2N − 1)k3 =
k4(4N2 − 2N). Again, the standard algorithm has O(N2) scaling, but this time with a
prefactor of k4.
For the optimized algorithm we once again define a left string array and right string
array:
l′1 = I (6.36)
l′j = l
′
j−1C j−1M j−1 (6.37)
and
r′P =MN (6.38)
r′j =M ′jC j+1r j+1 (6.39)
Notice that l′ and r′ have different definitions than l and r. We can compute electronic
gradients using
(∇xΓ) ′jk = l j (∇xC) ′jk r j (6.40)
and ( ) ( )
∇ ′ ′pΓ = l ∇ C r (6.41)jk j p jk j
121
( )
where (∇ C) ∂C jkx jk = ∂x and ∇
∂C jk ′
pC = ∂p . Because we resue l and r′ when computingjk jk jk
the gradients of Γ with respect to electronic position and momentum, we only pay the
cost to construct l′ and r′ once. This cost requires 4(2N − 1) matrix multiplications.
Each term in either gradient requires two matrix multiplications, and there are kN such
terms. In total we have 2kN + 4(2N − 1) matrix multiplications. Finally, our overall cost
2k4N + 4k3(2N − 1). Once again, we have succeeded in optimizing our algorithm from
O(N2) to O(N). We still retain the O(k4) prefactor.
122
APPENDIX A
APPROXIMATE PITE EXPANSION
Recall that √
≡ M + I − M
2
Θ arccos √ (A.1)
2
where
M = m e−∆βĤ0 (A.2)
√
with the conditions 0 < m0 < 1 and m0 , 1/ 2. Our goal is derive an approximation for
κΘ to first order in ∆β. Start by expanding M to first order in ∆β gives M ≈ m0(I − ∆βĤ).
Substituting this into Eq. (A.1) yields

arccos 
√
  2 2 m0(I − ∆βĤ) + I − m (I − ∆βĤ)

0 √  (A.3)2
Expanding the quadratic term under the radi√cal and throwing out terms with ∆β
2 gives


2
m0(I − ∆βĤ) + I − m0(I − 2∆βĤ)arccos √ 
2  (A.4)
√
Next we expand the radical using A + bX ≈ A1/2 + 1 A−12 bX. In our case, A = (1 − m20)I,
b = 2m20∆β, and X = Ĥ. This leads to √ − − 2 ∆ βm2m (I ∆βĤ) + 1 m I + 0 Ĥ 0 0 1−m2
arccos  √ 0  (A.5)2
where we used A1/2 = (1 − m20)I and A−1 = (1 − m−20 )I. Grouping terms in powers of ∆β
yields
 √m 1 − m2   0 + 1  m2  
arccos  √ 0 I + √  0 2 − m 0 ∆βĤ2 2 1 − m  (A.6)0
Finally, we expand arccos using the approximation arccos(A + bX) ≈ arccos(A) − bX(I −
A2)−1/2. This yields
κΘ ≈ θ0 − κ √ b ∆βĤ (A.7)
1 − a2
123
where  √ 2
θ0 = κ arccos m0 + 1 − m

0 √  (A.8)2
1
√ b =  m2 0 − 2 m0 (A.9)2 1 − m0
and √
m0 + 1 − m2
√ 0a = (A.10)
2
To simplify the final term, notice that
1 ( √ )
1 − a2 = 1 −( 1 + 2√m 20 1 −)m0 (A.11)21
= (1 − 2m0 √1 − m20 ) (A.12)21
= ( 1 − 2m 1 − m2 2 2( 0 √ 0 +))m0 − m0 (A.13)2 1 2
= √ m0 − 1 − m20 (A.14)2
This leaves us with  
√ √
 √  m0 − 12 
b 1 −( √
m0
= m0 (A.15)
1 − a2 )2
m0 − 1 − m20
√
where we canceled factors of 1/ 2. We need to be careful when taking the square root.
√ √
Note that 1 − a2 will always be positive since a takes on values [1/ 2, 1], where m0 =
√
1/ 2 gives the singular value of a = 1. However, if we simply take the squ√are root of the
denominator on the right hand side of Eq. (A.15) we are left with m0 − 1 − m20 which
can take both positive and negative values. To correct for this, we need to multiply by
κ = sgn(m0 − 1/2). This gives us  ( √ 1 )  √ m  0κ − 1− − 2 −√2 
m0 (A.16)m0 1 m0  1 m0 
1
=κ ( √ ) m √− 1 − m2  0 0 

m − 1 − m2 1 − m2 
m0 (A.17)
0 0 0
= κs1 (A.18)
124
√
where s = m / 1 − m21 0 0. Now we have
κΘ ≈ θ − κ20 s1∆βĤ (A.19)
κΘ ≈ θ0 − ∆βĤ, (A.20)
where we use the fact κ2 = 1. This is the expression we set out to derive.
125
BIBLIOGRAPHY
[1] Richard P. Feynman. “Quantum Mechanics and Path Integrals”. In: (1965).
[2] Ian R. Craig and David E. Manolopoulos. “Quantum statistics and classi-
cal mechanics: Real time correlation functions from ring polymer molec-
ular dynamics”. In: The Journal of Chemical Physics 121.8 (2004), pp. 3368–
3373. DOI: 10.1063/1.1777575.
[3] Nandini Ananth. “Mapping variable ring polymer molecular dynamics:
A path-integral based method for nonadiabatic processes”. In: The Journal
of Chemical Physics 139.12 (2013), p. 124102. DOI: 10.1063/1.4821590.
[4] Elliot C. Eklund and Ananth. Nandini. “Investigating the Stability and
Accuracy of a Classical Mapping Variable Hamiltonian for Nonadiabatic
Quantum Dynamics”. In: Regular and Chaotic Dynamics 26 (2021), pp. 131–
146. DOI: 10.1134/S1560354721020039.
[5] M D Newton and N Sutin. “Electron Transfer Reactions in Condensed
Phases”. In: Annual Review of Physical Chemistry 35.1 (1984), pp. 437–480.
DOI: 10.1146/annurev.pc.35.100184.002253.
[6] Harry B. Gray and Jay R. Winkler. “ELECTRON TRANSFER IN PRO-
TEINS”. In: Annual Review of Biochemistry 65.1 (1996). PMID: 8811189,
pp. 537–561. DOI: 10.1146/annurev.bi.65.070196.002541.
[7] Robert I. Cukier and Daniel G. Nocera. “PROTON-COUPLED ELEC-
TRON TRANSFER”. In: Annual Review of Physical Chemistry 49.1 (1998).
PMID: 9933908, pp. 337–369. DOI: 10.1146/annurev.physchem.49.
1.337.
[8] Michael Grätzel. “Photoelectrochemical cells”. In: Nature 414.6861 (2001),
pp. 338–344. DOI: 10.1038/35104607.
126
[9] Md. Wahadoszamen et al. “The role of charge-transfer states in energy
transfer and dissipation within natural and artificial bacteriochlorophyll
proteins”. In: Nature Communications 5.1 (2014), p. 5287. DOI: 10.1038/
ncomms6287.
[10] “The multi-configurational time-dependent Hartree approach”. In: Chem-
ical Physics Letters 165.1 (1990), pp. 73–78. ISSN: 0009-2614. DOI: https:
//doi.org/10.1016/0009-2614(90)87014-I.
[11] “The multiconfiguration time-dependent Hartree (MCTDH) method: a
highly efficient algorithm for propagating wavepackets”. In: Physics Re-
ports 324.1 (2000), pp. 1–105. ISSN: 0370-1573. DOI: https://doi.org/
10.1016/S0370-1573(99)00047-2.
[12] Michael Thoss and Haobin Wang. “SEMICLASSICAL DESCRIPTION OF
MOLECULAR DYNAMICS BASED ON INITIAL-VALUE REPRESENTA-
TION METHODS”. In: Annual Review of Physical Chemistry 55.1 (2004).
PMID: 15117255, pp. 299–332. DOI: 10.1146/annurev.physchem.
55.091602.094429.
[13] Ian R Craig and David E Manolopoulos. “Chemical reaction rates from
ring polymer molecular dynamics”. In: The Journal of chemical physics 122.8
(2005), p. 084106.
[14] Ian R Craig and David E Manolopoulos. “A refined ring polymer molecu-
lar dynamics theory of chemical reaction rates”. In: The Journal of chemical
physics 123.3 (2005), p. 034102.
[15] Scott Habershon et al. “Ring-Polymer Molecular Dynamics: Quantum Ef-
fects in Chemical Dynamics from Classical Trajectories in an Extended
Phase Space”. In: Annual Review of Physical Chemistry 64.1 (2013). PMID:
127
23298242, pp. 387–413. DOI: 10.1146/annurev-physchem-040412-
110122.
[16] Nandini Ananth. “Path Integrals for Nonadiabatic Dynamics: Multistate
Ring Polymer Molecular Dynamics”. In: Annual Review of Physical Chem-
istry 73 (2022), pp. 299–322.
[17] Jeremy O. Richardson and Michael Thoss. “Communication: Nonadi-
abatic ring-polymer molecular dynamics”. In: The Journal of Chemical
Physics 139.3 (2013), p. 031102. DOI: 10.1063/1.4816124.
[18] Pengfei Huo and David F. Coker. “Communication: Partial linearized
density matrix dynamics for dissipative, non-adiabatic quantum evolu-
tion”. In: The Journal of Chemical Physics 135.20 (2011), p. 201101. DOI: 10.
1063/1.3664763.
[19] Yudong Cao et al. “Quantum Chemistry in the Age of Quantum Com-
puting”. In: Chemical Reviews 119.19 (Oct. 2019), pp. 10856–10915. DOI:
10.1021/acs.chemrev.8b00803.
[20] Bela Bauer et al. “Quantum Algorithms for Quantum Chemistry and
Quantum Materials Science”. In: Chemical Reviews 120.22 (Nov. 2020),
pp. 12685–12717. DOI: 10.1021/acs.chemrev.9b00829.
[21] H. F. Trotter. “On the Product of Semi-Groups of Operators”. In: Proceed-
ings of the American Mathematical Society 10.4 (1959), pp. 545–551. ISSN:
00029939, 10886826. (Visited on 06/26/2022).
[22] Mark Tuckerman. “Statistical Mechanics: Theory And Molecular Simula-
tion”. In: (Jan. 2001).
128
[23] Nicholas Boekelheide, Romelia Salomón-Ferrer, and Thomas F. Miller.
“Dynamics and dissipation in enzyme catalysis”. In: Proceedings of the Na-
tional Academy of Sciences 108.39 (2011), pp. 16159–16163. DOI: 10.1073/
pnas.1106397108.
[24] Ryogo Kubo. “Statistical-Mechanical Theory of Irreversible Processes. I.
General Theory and Simple Applications to Magnetic and Conduction
Problems”. In: Journal of the Physical Society of Japan 12.6 (1957), pp. 570–
586. DOI: 10.1143/JPSJ.12.570.
[25] A. O. Caldeira and A. J. Leggett. “Influence of Dissipation on Quantum
Tunneling in Macroscopic Systems”. In: Phys. Rev. Lett. 46 (4 Jan. 1981),
pp. 211–214.
[26] Michael Thoss, Haobin Wang, and William H. Miller. “Self-consistent hy-
brid approach for complex systems: Application to the spin-boson model
with Debye spectral density”. In: The Journal of Chemical Physics 115.7
(2001), pp. 2991–3005. DOI: 10.1063/1.1385562.
[27] Ian R. Craig, Michael Thoss, and Haobin Wang. “Proton transfer reac-
tions in model condensed-phase environments: Accurate quantum dy-
namics using the multilayer multiconfiguration time-dependent Hartree
approach”. In: The Journal of Chemical Physics 127.14 (2007), p. 144503. DOI:
10.1063/1.2772265.
[28] Ryan P. Steele et al. “Mixed time slicing in path integral simulations”. In:
The Journal of Chemical Physics 134.7 (2011), p. 074112. DOI: 10.1063/1.
3518714.
[29] Wolfgang Domcke and David R. Yarkony. “Role of Conical Intersections
in Molecular Spectroscopy and Photoinduced Chemical Dynamics”. In:
129
Annual Review of Physical Chemistry 63.1 (2012). PMID: 22475338, pp. 325–
352. DOI: 10.1146/annurev-physchem-032210-103522.
[30] Jianshu Cao et al. “Quantum biology revisited”. In: Science Advances 6.14
(2020), eaaz4888. DOI: 10.1126/sciadv.aaz4888.
[31] Laurie J. Butler. “CHEMICAL REACTION DYNAMICS BEYOND THE
BORN-OPPENHEIMER APPROXIMATION”. In: Annual Review of Phys-
ical Chemistry 49.1 (1998), pp. 125–171. DOI: 10 . 1146 / annurev .
physchem.49.1.125.
[32] Jessica R. Duke and Nandini Ananth. “Simulating Excited State Dynam-
ics in Systems with Multiple Avoided Crossings Using Mapping Variable
Ring Polymer Molecular Dynamics”. In: The Journal of Physical Chemistry
Letters 6.21 (2015). PMID: 26722962, pp. 4219–4223. DOI: 10.1021/acs.
jpclett.5b01957.
[33] Jessica Ryan Duke and Nandini Ananth. “Mean field ring polymer molec-
ular dynamics for electronically nonadiabatic reaction rates”. In: Faraday
Discuss. 195 (0 2016), pp. 253–268. DOI: 10.1039/C6FD00123H.
[34] Timothy J. H. Hele and Nandini Ananth. “Deriving the exact nonadiabatic
quantum propagator in the mapping variable representation”. In: Faraday
Discuss. 195 (0 2016), pp. 269–289. DOI: 10.1039/C6FD00106H.
[35] Sadrach Pierre et al. “A mapping variable ring polymer molecular dy-
namics study of condensed phase proton-coupled electron transfer”. In:
The Journal of Chemical Physics 147.23 (2017), p. 234103. DOI: 10.1063/1.
4986517.
[36] Hans-Dieter Meyera) and William H. Miller. “A classical analog for elec-
tronic degrees of freedom in nonadiabatic collision processes”. In: The
130
Journal of Chemical Physics 70.7 (1979), pp. 3214–3223. DOI: 10.1063/1.
437910.
[37] Gerhard Stock and Michael Thoss. “Semiclassical Description of Nonadi-
abatic Quantum Dynamics”. In: Phys. Rev. Lett. 78 (4 Jan. 1997), pp. 578–
581. DOI: 10.1103/PhysRevLett.78.578.
[38] Nandini Ananth and Thomas F. Miller. “Exact quantum statistics for elec-
tronically nonadiabatic systems using continuous path variables”. In: The
Journal of Chemical Physics 133.23 (2010), p. 234103. DOI: 10.1063/1.
3511700.
[39] Michael Thoss and Gerhard Stock. “Mapping approach to the semiclassi-
cal description of nonadiabatic quantum dynamics”. In: Phys. Rev. A 59 (1
Jan. 1999), pp. 64–79. DOI: 10.1103/PhysRevA.59.64.
[40] Daniel T. Colbert and William H. Miller. “A novel discrete variable rep-
resentation for quantum mechanical reactive scattering via the S-matrix
Kohn method”. In: The Journal of Chemical Physics 96.3 (1992), pp. 1982–
1991. DOI: 10.1063/1.462100.
[41] C. Lopesino et al. “A Theoretical Framework for Lagrangian Descriptors”.
In: International Journal of Bifurcation and Chaos 27.01 (2017), p. 1730001.
DOI: 10.1142/S0218127417300014.
[42] Atanasiu Stefan Demian and Stephen Wiggins. “Detection of Periodic Or-
bits in Hamiltonian Systems Using Lagrangian Descriptors”. In: Interna-
tional Journal of Bifurcation and Chaos 27.14 (2017), p. 1750225. DOI: 10.
1142/S021812741750225X.
[43] Shibabrat Naik, Vı́ctor J. Garcı́a-Garrido, and Stephen Wiggins. “Finding
NHIM: Identifying high dimensional phase space structures in reaction
131
dynamics using Lagrangian descriptors”. In: Communications in Nonlinear
Science and Numerical Simulation 79 (2019), p. 104907. ISSN: 1007-5704. DOI:
https://doi.org/10.1016/j.cnsns.2019.104907.
[44] Sebastian Nowozin Suvrit Sra and Stephen J. Wright. “Optimization for
Machine Learning”. In: (2012).
[45] Evan R. MacQuarrie et al. In: Nature Reviews Physics 2.11 (2020). DOI: 10.
1038/s42254-020-00247-5.
[46] Frank Arute et al. “Quantum supremacy using a programmable super-
conducting processor”. In: Nature 574.7779 (2019), pp. 505–510. DOI: 10.
1038/s41586-019-1666-5.
[47] Y. S. Weinstein et al. “Implementation of the Quantum Fourier Trans-
form”. In: Phys. Rev. Lett. 86 (9 Feb. 2001), pp. 1889–1891. DOI: 10.1103/
PhysRevLett.86.1889.
[48] Miroslav Dob šı ček et al. “Arbitrary accuracy iterative quantum phase
estimation algorithm using a single ancillary qubit: A two-qubit bench-
mark”. In: Phys. Rev. A 76 (3 Sept. 2007), p. 030306. DOI: 10 . 1103 /
PhysRevA.76.030306.
[49] Alberto Peruzzo et al. “A variational eigenvalue solver on a photonic
quantum processor”. In: Nature Communications 5.1 (2014), p. 4213. DOI:
10.1038/ncomms5213.
[50] John Preskill. “Quantum Computing in the NISQ era and beyond”. In:
Quantum 2 (Aug. 2018), p. 79. ISSN: 2521-327X. DOI: 10.22331/q-2018-
08-06-79.
132
[51] Christian Kraglund Andersen et al. “Repeated quantum error detection
in a surface code”. In: Nature Physics 16.8 (2020), pp. 875–880. DOI: 10.
1038/s41567-020-0920-y.
[52] Michael A. Nielsen and Isaac L. Chuang. Quantum Computation and Quan-
tum Information: 10th Anniversary Edition. Cambridge University Press,
2010. DOI: 10.1017/CBO9780511976667.
[53] Colin D. Bruzewicz et al. “Trapped-ion quantum computing: Progress and
challenges”. In: Applied Physics Reviews 6.2 (2019), p. 021314. DOI: 10.
1063/1.5088164.
[54] M. H. Devoret and R. J. Schoelkopf. “Superconducting Circuits for Quan-
tum Information: An Outlook”. In: Science 339.6124 (2013), pp. 1169–1174.
DOI: 10.1126/science.1231930.
[55] Ke Wang et al. “Spin manipulation in semiconductor quantum dots
qubit”. In: Chinese Physics B 27.9 (Sept. 2018), p. 090308. DOI: 10.1088/
1674-1056/27/9/090308.
[56] S. Bonella et al. “Path integral based calculations of symmetrized time
correlation functions. II”. In: The Journal of Chemical Physics 133.16 (2010),
p. 164105. DOI: 10.1063/1.3493449.
[57] Nancy Makri. “Quantum Dissipative Dynamics: A Numerically Exact
Methodology”. In: The Journal of Physical Chemistry A 102.24 (June 1998),
pp. 4414–4427. DOI: 10.1021/jp980359y.
[58] Mario Motta et al. “Determining eigenstates and thermal states on a
quantum computer using quantum imaginary time evolution”. In: Nature
Physics 16.2 (2020), pp. 205–210. DOI: 10.1038/s41567-019-0704-4.
133
[59] Shi-Ning Sun et al. “Quantum Computation of Finite-Temperature Static
and Dynamical Properties of Spin Systems Using Quantum Imaginary
Time Evolution”. In: PRX Quantum 2 (1 Feb. 2021), p. 010317. DOI: 10.
1103/PRXQuantum.2.010317.
[60] Taichi Kosugi, Yusuke Nishiya, and Yu-ichiro Matsushita. Probabilistic
imaginary-time evolution by using forward and backward real-time evolution
with a single ancilla: first-quantized eigensolver of quantum chemistry for
ground states. 2021. DOI: 10.48550/ARXIV.2111.12471.
[61] Ivan Kassal et al. “Polynomial-time quantum algorithm for the simulation
of chemical dynamics”. In: Proceedings of the National Academy of Sciences
105.48 (2008), pp. 18681–18686. DOI: 10.1073/pnas.0808245105.
[62] Giuliano Benenti and Giuliano Strini. “Quantum simulation of the single-
particle Schrödinger equation”. In: American Journal of Physics 76.7 (2008),
pp. 657–662. DOI: 10.1119/1.2894532.
[63] Pauline J. Ollitrault, Guglielmo Mazzola, and Ivano Tavernelli. “Nona-
diabatic Molecular Quantum Dynamics with Quantum Computers”. In:
Physical Review Letters 125 (26 2020), p. 260511. ISSN: 10797114. DOI: 10.
1103/PhysRevLett.125.260511.
[64] Hans Hon Sang Chan et al. Grid-based methods for chemistry simulations on
a quantum computer. 2022. DOI: 10.48550/ARXIV.2202.05864.
[65] Osvaldo Simeone. An Introduction to Quantum Machine Learning for Engi-
neers. 2022. DOI: 10.48550/ARXIV.2205.09510. URL: https://
arxiv.org/abs/2205.09510.
134
[66] Alexander Miessen, Pauline J. Ollitrault, and Ivano Tavernelli. “Quan-
tum algorithms for quantum dynamics: a performance study on the spin-
boson model”. In: (2021), pp. 1–12.
[67] Glen Bigan Mbeng, Angelo Russomanno, and Giuseppe E. Santoro. The
quantum Ising chain for beginners. 2020. DOI: 10.48550/ARXIV.2009.
09208.
[68] Adam Smith et al. “Simulating quantum many-body dynamics on a cur-
rent digital quantum computer”. In: npj Quantum Information 5 (1 Dec.
2019). ISSN: 20566387. DOI: 10.1038/s41534-019-0217-0.
135