%\documentclass[letterpaper,12pt]{article}
\documentclass[letterpaper]{tufte-handout}
\usepackage{times}
\usepackage[]{hyperref}
%\usepackage[section]{placeins}
\usepackage{graphicx}
\usepackage{amsmath,amsfonts,amssymb}
\usepackage{mathrsfs} % formal script font; use \mathscr{...}
\usepackage[titletoc]{appendix} % adds ``Appendix'' in TOC
\usepackage{datetime} % gives \currenttime
\usdate% keeps standard \LaTeX date format for \today
\usepackage{gensymb}
\usepackage{bbm}% for \mathbbm % blackboard bold
\usepackage{mathtools} % gives \substack
\title{Computing Geodesics in the Landscape Ensemble}
\date{Modified: \today\ at \currenttime}
\author{Vale Cofer-Shabica}
\newcommand{\laeq}[1]{\label{eqn:#1}}
\renewcommand{\refeq}[1]{eq.~\ref{eqn:#1}}
\renewcommand{\vec}[1]{\mathbf{#1}}
\newcommand{\mat}[1]{\;\textrm{\textbf{#1}}\,}
%\newcommand{\mat}[1]{\mbox{\normalsize $\,\mathtt{#1}\,$}}
%\newcommand{\mat}[1]{\,\mathtt{#1}\,}
\newcommand{\trans}[1]{{#1}^{\mathsf{T}}}
\newcommand{\norm}[1]{\left\lVert#1\right\rVert}
\newcommand{\svec}[1]{\vec{#1}}
\newcommand{\cvec}[1]{\vec{#1}}
\renewcommand{\refeq}[1]{eq.~\ref{eqn:#1}}
\newcommand{\notesym}{\textsuperscript{*}}
\begin{document}
\maketitle
\newthought{We previously discussed} how a path-integral formulation of the diffusion equation suggests that geodesics on an energy landscape will be the most important components of the dynamics on that surface. Today we will derive the key equations of the algorithm the group uses to construct geodesics.
\section{Introduction}
Last time we discussed in detail, the origin of the key equation from Wang and Stratt's 2007 paper on geodesics~\cite[-1\baselineskip]{wang:2007:geodesics}:
\begin{equation}\laeq{pathd}
G(\vec{R}_{0} \to \vec{R} , t) = \int_{\vec{R}_{0},0}^{\vec{R},t} \mathscr{D} [\vec{R}(\tau)]\exp\left[ -\frac{1}{4D} \int_{0}^{t} {\left(\frac{d\vec{R}}{d\tau} \right)}^{2} d\tau\right]
\end{equation}
This is a path integral expression for free diffusion and we inferred that the dominant paths when diffusion is slow ($D$ is small) are those that minimize the integral in the exponential---those that minimize the action.
We have not yet specified the nature of the boundary conditions in this space nor do we yet have any machinery to construct such action-minimizing paths.
To situate our thinking within the problem area, a broad-brushed overview of geodesic analysis as typically practiced by the group follows
\begin{enumerate}
\item Construct endpoint pairs \emph{(some combination of molecular dynamics or Monte Carlo methods)}
\item Select a landscape energy \emph{(perhaps an ensemble average, perhaps something else)}
\item \textbf{Construct geodesics at that energy \emph{(subject of this talk)}}
\item Analyze their properties \emph{(many ways, often starting with their length)}
\item Infer something interesting about the system under study
\end{enumerate}
Today's discussion sits at the middle of the process and brings us from the \emph{`why?'} to the \emph{`how?'} of one of the groups lenses for understanding the dynamics of chemical systems. We will treat the problem in a Cartesian coordinate system of $N$ particles with equal masses. References will be provided along the way for deviations from this basic formulation.
\FloatBarrier\
\subsection{Boundary Conditions: The Landscape Ensemble}
\begin{figure}
\includegraphics[width=.95\textwidth]{wang2007}
\caption{\label{fig:pele} The potential energy landscape ensemble takes the landscape energy, $E_L$, as its thermodynamic variable of control. This fixed maximum potential energy plays a similar role to the microcanonical energy or the canonical temperature. In the region depicted, the shaded mountains are forbidden, while the other areas are accessible. \emph{Figure from Wang and Stratt.}}
\end{figure}
We constructed \refeq{pathd} with the goal of understanding the \emph{inherent dynamics}\sidenote[][-2\baselineskip]{The inherent dynamics of a system can be thought of as the system's classical dynamics stripped of noise (\emph{e.g.} high frequency vibrations or rotations) that hamper identification of the essential steps or motions of the system. They describe the dominant thoroughfares of configuration space---the major motions during a chemical or physical transformation.} or most efficient paths of systems under study. There are many ways that we could formulate a notion of most efficient path, but an expedient one comes from combining \refeq{pathd} with the potential energy landscape ensemble\cite{wang:2007:pele}.
The landscape ensemble is a statistical mechanical ensemble like the canonical or microcanonical ensembles. However, instead of a fixed temperature or energy, its thermodynamic variable of control is a fixed maximum potential energy, called the landscape energy, $E_L$. Within this ensemble, all configurations such that the potential is less than the landscape energy,
\[
V(\cvec{R}) \le E_L
\]
are admitted with equal probability\sidenote[][]{Equiprobability is a result of the same entropy maximization arguments that lead to a similar result in the microcanonical ensemble.}. In the thermodynamic limit, this ensemble gives the same configuration space probability densities that all other such ensembles give.
The landscape ensemble partitions configuration space into allowed and dis-allowed regions and provides boundary conditions for \refeq{pathd}---we must minimize the action subject to the landscape energy:
\begin{equation} \laeq{peleS}
S = \int_{0}^{t}d\tau {\left(\frac{d\vec{R}}{d\tau} \right)}^{2} \quad \textrm{s.t.} \quad V(\cvec{R}(\tau)) \le E_L
\end{equation}
This condition takes the form of a non-holonomic constraint\sidenote[][-1\baselineskip]{A holonomic constraint is one that can be expressed as an \emph{equality} relationship between a constant and a function of coordinates and time, \emph{e.g.} $0 = f(q_1, q_2, \ldots, t)$.}. Solutions to problems with non-holonomic are often difficult, but the form of solution is specified by the Kuhn-Tucker theorem\cite{kuhntucker:1951}. The theorem holds that the desired solution will be the union of paths that unconditionally minimize the functional in \refeq{peleS} and those where strict equality
\FloatBarrier\
\section{Unconditional minimization}
In this section, we follow Goldstein~\cite[-2\baselineskip]{goldstein:1980} to find the form of the (shortest) paths, which unconditionally minimize the action. Again, we will work in Cartesian coordinates for an otherwise unconstrained system. See Jacobson and Stratt~\cite{jacobson:2014} for work dealing with free rotations and Frechette and Stratt~\cite{frechette:2016} for a case where a constraint is applied to the free steps.
Using calculus to find the extrema of a function, for instance $h(x)$, is a familiar operation. By equating the derivative of $h(x)$ with $0$, we locate values of $x$ for which $h(x)$ is unchanging or stationary; these locations often correspond to local extrema of the function. The problem treated by the calculus of variations is analogous. We seek to extremize an object that takes a \emph{function} (in our case, a path) as its argument. That is, we find some path such that the object is stationary with respect to variations in that function.
From \refeq{peleS}, we seek to minimize an integral over kinetic energy, which in general is a function of coordinates and their derivatives. Working in one dimension\footnote{This procedure extends to systems with an arbitrary number of degrees of freedom}, we seek to find the path $R(\tau)$ that extremizes:
\begin{equation}\laeq{J}
J = \int_{0}^{t}d\tau f \left(R, \dot{R} \right)
\end{equation}
That is, we seek $R(\tau)$ such that $J$ is unchanged for infinitesimal variations in $R$. We can represent the paths neighboring the correct (extremal) path, $R(\tau)$ by adding an arbitrary function, $\eta(\tau)$ in proportion to a small variable parameter, $\alpha$:
\begin{equation}\laeq{eta}
R(\tau,\alpha) = R(\tau) + \alpha\eta(\tau)
\end{equation}
As we are interested in $R$ such that it connects $R(0)$ to $R(t)$, we require that:
\begin{equation}\laeq{eta:bounds}
\eta(0)=\eta(t)=0
\end{equation}
$J$ in \refeq{J} is now a function of $\alpha$ as well:
\begin{equation}\laeq{Ja}
J(\alpha) = \int_{0}^{t}d\tau f \left(R(\tau,\alpha), \dot{R}(\tau, \alpha) \right)
\end{equation}
We can encode the requirement for stationary $J$ in the following familiar form\sidenote[][]{We require that the derivative be evaluated at $\alpha=0$ because this corresponds to the unperturbed path we seek.}:
\begin{equation}
0 = {\left. \frac{d J}{d \alpha} \right|}_{\alpha=0}
\end{equation}
We can bring the $\alpha$ derivative inside the integral in \refeq{Ja} and applying the chain rule, we have:
\begin{equation*}
\frac{d J}{d \alpha} = \int_{0}^{t}d\tau \left( \frac{\partial f}{\partial R}\frac{\partial R}{\partial \alpha} + \frac{\partial f}{\partial \dot{R}}\frac{\partial \dot{R}}{\partial \alpha} \right)
\end{equation*}
From \refeq{eta}, we have that $\frac{\partial R}{\partial \alpha} = \eta$ and $\frac{\partial \dot{R}}{\partial \alpha} = \dot{\eta}$, the arbitrary function of $\tau$ and its time derivative. This leaves:
\begin{equation}\laeq{J:1}
\frac{d J}{d \alpha} = \int_{0}^{t}d\tau \left( \frac{\partial f}{\partial R}\eta + \frac{\partial f}{\partial \dot{R}} \dot{\eta} \right)
\end{equation}
Consider the second term in the integral:
\begin{equation*}
\int_{0}^{t}d\tau\frac{\partial f}{\partial \dot{R}}\dot{\eta}
\end{equation*}
This expression is susceptible by integration by parts
\sidenote[][-3\baselineskip]{
Recall the general form,
\[
\int_a^b u dv = {\left. uv \right|}_a^b - \int_a^b v du
\]
And identify:
\begin{center}
\vspace{-\baselineskip}
\begin{tabular}{c c c}
$\begin{aligned}
u &= \frac{\partial f}{\partial \dot{R}} \\
dv &= \dot{\eta} d\tau
\end{aligned}$
& \hspace{3em} &
$\begin{aligned}
du &= \frac{d}{d\tau} \left( \frac{\partial f}{\partial \dot{R}} \right) d\tau \\
v &= \eta
\end{aligned}$
\end{tabular}
\end{center}
}
Re-arranging leaves us with:
\begin{equation*}
\int_{0}^{t}d\tau\frac{\partial f}{\partial \dot{R}}\dot{\eta} = {\left. \frac{\partial f}{\partial \dot{R}} \eta \right|}_0^t - \int_{0}^{t}d\tau \frac{d}{d\tau} \left( \frac{\partial f}{\partial \dot{R}} \right) \eta
\end{equation*}
Because of the bounds on $\eta$, $\eta(0)=\eta(t)=0$, the surface terms vanish! Inserting back into \refeq{J:1} and factoring out $\eta$ gives:
\begin{equation}\laeq{J:2}
{\left. \frac{d J}{d \alpha} \right|}_{\alpha=0} =
\int_{0}^{t}d\tau \left[ \frac{\partial f}{\partial R} - \frac{d}{d\tau} \left( \frac{\partial f}{\partial \dot{R}} \right) \right] \eta = 0
\end{equation}
where evaluating the derivative at $\alpha=0$ gives $R$ and $\dot{R}$ as functions of $\tau$ alone. Since $\eta$ is an arbitrary function, the integral is stationary \emph{only} if
\begin{equation}\laeq{G2-11}
\frac{\partial f}{\partial R} - \frac{d}{d\tau}\frac{\partial f}{\partial\dot{R}} = 0
\end{equation}
which gives the conditions on $R$ to minimize $J$\sidenote[][-3\baselineskip]{Notice that inserting the Lagrangian for $f$ gives $S=J$, the action, and yields Lagrange's equations of motion.}.
\subsection{Minimizing The Force-Free Action}
From \refeq{peleS}, we have $f={\left(\frac{d\vec{R}}{d\tau} \right)}^{2}$ in \refeq{G2-11}. Proceeding, we find that the path which minimizes \refeq{peleS} is specified by the differential equation:
\begin{equation}
\ddot{\vec{R}} = 0
\end{equation}
which, when integrated and solved for its boundary values, yields the equation for a straight line (and the geodesic!):
\begin{equation}\laeq{straightLine}
\vec{R}(\tau) = \vec{R}(0) \left( 1-\frac{\tau}{t} \right) + \vec{R}(t)\left( \frac{\tau}{t} \right)
\end{equation}
this amounts to linearly interpolating between $\vec{R}(0)$ and $\vec{R}(t)$ in time $t$ for each of the $3N$ coordinates. To actually compute such paths we discretize\sidenote[][]{Consider $\vec{R}(\tau + \Delta) - \vec{R}(\tau)$.}\refeq{straightLine} as follows:
\begin{equation}\laeq{step:free}
\boxed{
\cvec{R}^{(t+1)} = \cvec{R}^{(t)} + \delta R \frac{\cvec{R} - \cvec{R}^{(t)}}{\left| \cvec{R} - \cvec{R}^{(t)} \right|}
}
\end{equation}
Where $\delta R$ is a small\sidenote[][-3\baselineskip]{\emph{Small} here is relative to features of the landscape. In practice, $\delta R$ should be as large as it can be such that the lengths of computed geodesics are still converged to a maximum.} step-size.
\section{Strict Equality: Paths along boundaries}
So far we have successfully dealt with the first part of the solution suggested by the Kuhn-Tucker theorem, namely the unconditional minimization of the action integral. However, unless we are extraordinarily lucky, the straight-line motion implied by \refeq{step:free} will surely send the system into a region violating the landscape criterion, $V(\cvec{r}) > E_L$. Once this happens, we need a way to proceed.
See Ma and Stratt~\cite[-3\baselineskip]{ma:2015} for an entirely different approach for tracing the boundaries in a hard sphere system. See Cofer-Shabica and Stratt~\cite[-1\baselineskip]{dvcs:2017:jcp} for a method to construct escape steps that preserve the center of mass for systems with hetero-atoms
\subsection{A Newton-Raphson root search}
The Kuhn-Tucker theorem suggests that we must closely trace the boundaries of the obstacles specified by $V(\cvec{R}(\tau)) \le E_L$. We aim, therefore to take the most efficient path out of the obstacle once \refeq{step:free} causes us to stumble into it. That is, we seek the nearest solution to $V(\cvec{R}(\tau)) = E_L$. A classic strategy for solving this problem is a Newton-Raphson root search.
The key insight behind Newton-Raphson is that the gradient points in the direction of most rapid change in the function. We can exploit this by following the gradient to the solution we seek. Here, we expand the potential to first order and then demand that this approximation yield the landscape energy:
\begin{align}\laeq{expansion}
V(\cvec{R}) \approx V(\cvec{R}_0) + \left. \cvec{\nabla}V\right|_{\cvec{R}_0} \cdot (\cvec{R} - \cvec{R}_0) = E_L
\end{align}
where $\cvec{R}_0$ was the first offending step in \refeq{step:free}. There are no bounds on the number of possible solutions to this equation, but we make the \emph{guess} that heading downhill via the steepest path will be optimal\sidenote[][]{In a system with hetero-atoms we could accomplish the same goal without perturbing the center of mass by following the inverse-mass-weighted gradient. You can convince yourself of this by considering the direction the a dynamical system moves under the action of a force from \emph{rest}.}. This suggests a solution of the form:
\begin{equation}\laeq{NR:vanilla}
(\cvec{R} - \cvec{R}_0) = C \left. \cvec{\nabla} V\right|_{\cvec{R}_0} \equiv \Delta\cvec{R}
\end{equation}
Rearranging for the step-length, $C$, immediately yields:
\[
C = \frac{E_L - V(\cvec{R}_0)}{{\left| \left. \cvec{\nabla} V\right|_{\cvec{R}_0} \right|}^2}
\]
Which we insert back into \refeq{NR:vanilla} to find:
\begin{equation}\laeq{dR:noCM}
(\cvec{R} - \cvec{R}_0) = - \frac{V \left(\cvec{R}_0\right) - E_L}{ {\left| {\left. \cvec{\nabla}V \right|}_{\cvec{R}_0} \right|}^2} \; {\left. \cvec{\nabla}V \right|}_{\cvec{R}_0}
\end{equation}
Applying \refeq{dR:noCM} iteratively leads to the relation:
\begin{equation}\laeq{step:escape}
\boxed{
\cvec{R}_{n+1} = \cvec{R}_{n} - \frac{V \left(\cvec{R}_n\right) - E_L}{ {\left| {\left. \cvec{\nabla}V \right|}_{\cvec{R}_n} \right|}^2} \; {\left. \cvec{\nabla}V \right|}_{\cvec{R}_n}
}
\end{equation}
where $\cvec{R}_0$ is the first $\cvec{R}^{(t)}$ with $V\left( \cvec{R}^{(t)} \right) > E_L$ from \refeq{step:free}. This procedure is applied repeatedly until $V(\vec{R}_n) < E_L$. And then free steps (\refeq{step:free}) are resumed with $\vec{R}^{t+1}=\vec{R}_n$. Carried out in conjunction with \refeq{step:free}, the system will repeatedly stray into the forbidden region and then escape via the expression we just derived. This will cause the path to wind around and hug the bounds of obstacles exactly as required by the Kuhn-Tucker theorem.
\section{Optimizing Paths}
Together, \refeq{step:free} for free steps and \refeq{step:escape} for escape steps, provide the machinery to compute paths satisfying the Kuhn-Tucker theorem on arbitrary potentials. However, while the theorem imposes a necessary condition, it is not sufficient on its own---it provides no guarantee that the path will be \emph{shortest}\sidenote[][]{It is worth noting that the free-step portion of our algorithm gives the correct answer in the limit of low-obstacle density.}.
The algorithm described above will generate a path, $\vec{R}_0 \to \vec{R}$. A first step in refinement is to check the reverse path and see if it is shorter, $\vec{R} \to \vec{R}_0$; if so this new path is retained. Then a point along the path is randomly selected, $\vec{R}_i$ and displaced using small single-atom moves to a new point $\vec{R}_i'$. The same algorithm is used to compute geodesics, $\vec{R}_0 \to \vec{R}_i'$ and $\vec{R}_i' \to \vec{R}$. If the length of this new, combined-path is shorter than the original, we keep the new path otherwise we discard it. This process is repeated until the path-lengths are converged\sidenote[][]{Typical length-reductions are of the order of 5\%.}. One of the methods of this manner of shortening, is that it uses the algorithm we already have and so the optimized paths retain any properties of paths produced by the first algorithm, for instance obeying the Kuhn-Tucker theorem.
\section{Conclusion}
We have derived the key equations of the algorithm for computing geodesics on a potential energy surface specified in Cartesian coordinates. Free steps are given by \refeq{step:free} and when the system encounters an obstacle, \refeq{step:escape} will lead it out. An optimization process further refines the paths until they are converged in length.
\bibliographystyle{unsrt}
\bibliography{\jobname}
\end{document}