Why is the objective function used in Nudged Elastic Band method reasonable? In Nudged Elastic Band (NEB) method, which is used to find reaction pathways when both initial and final states are know, an objective function is first constructed and then minimized to find reaction pathways.
NEB creates a string of replicas (or 'images') of the system between initial and final states, and connects them with springs to represent a path connecting initial and final states. Initially, the images may be interpolated between reactant and product linearly. Then an optimization algorithm is applied to relax the images down towards the Minimum Energy Path (MEP). 
The objective function is constructed as a combination of potential energy of each image and "spring energy" between each image as equation below shows.
$$
S(\vec R_0,...,\vec R_N)=\sum_{i=1}^{N-1}E(\vec R_i)+\sum_{i=1}^N\frac{k}{2}(\vec R_i-\vec R_{i-1})^2
$$
Where $\vec R_i$ is the coordinates of i-th image.
NEB solved "corner-cutting" and "sliding-down" problems in this model by considering vertical and parallel components of potential force and string force. I will not introduce here, you can look at this article for more details.
My question is why the objective function above is reasonable. Does minimization of this objective function sufficiently leads to the MEP? If yes, how to prove it? If no, why is this "chain-of-states" model widely used?
Thank you very much for answering.
 A: I have no background in this matter, but I think some basic intuition is in order.
Sliding down:
Suppose you have chain of rubber balls connected by elastic springs. Hold up the chain and let it dangle. Notice that the energy of the system is exactly as above, with $E(R_n)=mgZ_n$ where $Z_n$ is the height of ball $n$. 
What will it look like? Numbering the ball at the bottom as $n=1$ and the ones above it as $n=2,3...$, you have that the force on ball $n$ from the ones hanging below it is $mg(n-1)$. Thus the stretch distance of the spring below it is $d(n)=mg(n-1)/k$. Defining the height of ball 1 as 0, the height of ball $n$ will be
$$h(n)=\sum_{J=1}^nd(J)=\frac{g m \left(n^2-n\right)}{2 k}$$ which looks like this:
ListLinePlot[
 Table[{0, (g m (-n + n^2))/(2 k) /. {g -> 1, m -> 1, k -> 1}}, {n, 
   10}], Axes -> False, PlotMarkers -> Automatic, AspectRatio -> 4]


As you can see, it gets stretched out at the top because the links at the top are supporting more weight.
Now suppose you drape the rope over a turtle (or some other hill-shaped object). What will it look like? From the previous example, it should be intuitively obvious that at the top of the turtle, the balls will be stretched farther apart, whereas at the bottom edges the springs will be more relaxed and the balls will be closer together. In other words, the resolution is worse at the top of the turtle's shell than it is at the bottom edges.
This is a problem. When you try to find a minimum energy surface, you want to get good resolution at the saddle point region (aka, point of no return), but this method does the exact opposite. That is what is meant when the article says "sliding down": the beads sag towards the unimportant parts, and get stretched thin over the important parts. 
Notice that the force that causes this problem is parallel to the springs, as mentioned in the article.
Cutting corners:
Now suppose you are in a valley between two hills, and that the valley curves in one direction. Drop a rope into the valley, and pull it taut. 
What does it look like?
If you don't pull too hard, it will straighten out and do its best to conform to the shape of the valley, sort of like this:

But, if you pull it too hard, it will start to cut across the corner, like this:

Notice that the force that causes this problem is perpendicular to the springs, as mentioned in the article.
A: Ok, so, this is a pretty late answer, but I was thinking about this and I believe I figured it out - though I'm not even sure whether the original NEB was implemented with this exact justification in mind, seeing how cavalier it seems to be with its notation. But nevertheless, it's a good justification. Strap in, because it's a bit of a ride - that involves quantum mechanics and Feynman path integrals, no less!
So, let's first begin with the notion of 'action'. Action is a quantity defined over a path. If we have a motion in time from $t=0$ to $t=T$ following a path $\mathbf{x}(t)$ then the action is the integral in time of the Lagrangian, $\mathcal{L}$:
$$
S = \int_0^T \mathcal{L}(t)dt= \int_0^T\frac{1}{2}m\dot{\mathbf{x}}^2-V(\mathbf{x})dt
$$
This factors in an interesting property of quantum mechanics. If we ask ourselves "how likely it is for a particle in that can be found at $\mathbf{x}_0$ at $t=0$ in a potential $V$ to then be found at $\mathbf{x}_N$ at $t=T$?", the answer can be either expressed by a propagator, or by a path integral. Specifically:
$$
P(\mathbf{x}_N, T;\mathbf{x}_0, 0) = \langle\mathbf{x}_N|e^{-i\frac{H}{\hbar}T}|\mathbf{x_0}\rangle=A\int D[\mathbf{x}(t)]e^{i\frac{S}{\hbar}}
$$
Ok. Explanation time! The first formulation is your typical quantum bra-ket notation: you pick your initial state, you evolve it using a time propagator that is the exponential of the Hamiltonian (aka the energy of the system) divided by the Planck constant times $-iT$, then you project it against the desired final state. The overlap between the evolved state and the desired one is how likely it is for that specific process to happen. Still with me? Good.
The second formulation is the Feynman path integral notation. It says something slightly different: it says that the probability of finding the system in that final state is proportional (there's a normalization factor $A$ that we just won't worry about for now) to the sum of an imaginary exponential of the action, divided by the Planck constant, for all possible paths that connect the two states in that time. That's what $ D[\mathbf{x}(t)]$ means: it's not an integral over a variable, it's an integral over functions. And as you may guess, that's mighty hard to compute. More on that in a minute. First, let's consider what this path integral formalism means.
The action, like the energy, isn't an absolute value, it's defined up to a constant, being an integral of a very energy-like quantity, the Lagrangian, and all. So let's assume that there exists one path that connects our two events that has minimal action, and all other paths have bigger one. Note, this assumption isn't always true, and that means this reasoning and all its consequences can get messed up sometimes; but then again, so can the NEB if there are two equally possible saddle points, so bear with me. If we work under that assumption, then we can set $S=0$ for that path, and $S > 0$ for all the others. These other paths then will contribute with oscillatory terms to the overall integral, and the bigger $S$, the faster the oscillations, the more likely they will just cancel each other out. If we imagine shrinking $\hbar$ (which is a big no-no, it's not called a constant for nothing, but let's fashion ourselves gods and create our own versions of the universe for a moment), then these oscillations become wilder and wilder; and in the limit of $\hbar \rightarrow 0$, that is, in the limit of a perfectly classical, Newtonian, totally-not-quantum universe, they go completely out of control, and only one path remains to contribute: the one with the minimal action.
We have just retrieved the Principle of Least Action, which says that the (Newtonian) path between two points in space and time is always the one with the minimal possible action.
So, what does this have to do with NEB? Well, we need a few more steps, and a trick.
Let's suppose we have a classical system, and want to compute the least action path between two points in space and time. The thing is, all Newtonian trajectories are the least action path between where they start and where they arrive; but we don't know where they will arrive before we try them. Here, instead, we know both initial and final conditions, and we don't know anything about the path itself (including the initial velocity). So, how do we go about doing that, especially with a computer?
Well, I'd say, we discretize the integral to compute the action into $N$ steps, with a time step $dt = T/N$, separating it into a sum of intermediate steps $\mathbf{x}_1, \mathbf{x}_2, ...$. In this way, the action becomes:
$$
S = \sum_{i=0}^N\left[\frac{1}{2}m\dot{\mathbf{x_i}}^2-V(\mathbf{x_i})\right]dt
$$
This discretization by the way is also an excellent way of computing the integral above, and is often used for example in quantum field theory.
So how do we compute those velocities? Well, let's just assume that they are constant between each pair of steps, so
$$
\dot{\mathbf{x_i}} = \frac{\mathbf{x_i}-\mathbf{x_{i-1}}}{dt}
$$
and 
$$
S = \sum_{i=i}^N\frac{1}{2}m \frac{(\mathbf{x_i}-\mathbf{x_{i-1}})^2}{dt}-\sum_{i=0}^NV(\mathbf{x_i})dt
$$
Is this starting to look like your original NEB object function, if you take $k=m/dt$?
But see, there's still a problem, namely, that pernicious "minus" sign in front of the potential. In your original formula, it's a plus! That's an Hamiltonian, not a Lagrangian. So, what makes it disappear?
Last trick, I swear. Wick Rotation time. Another favourite of purveyors of QFT.
This one sounds a bit like magic, really. See that quantum propagator, above? I mean this:
$$
e^{-i\frac{H}{\hbar}t}
$$
Now, if you know the Hamiltonian is basically the system's energy, this looks a lot like a partition function. So let's make it look like one. Let's make a change of parameters: $T \rightarrow -i\tau$.
$$
e^{-i\frac{H}{\hbar}t} \rightarrow e^{-\frac{H}{\hbar}\tau}
$$
Ok, that's gotta be cheating, right? But in fact, it's perfectly okay, we're simply redefining one parameter in our maths, nothing's changed. We call $\tau$ the "imaginary time" and the only thing that's really important to remember is that it has nothing to do with real time and we should never relate the two as if they were the same, they're not. Now let's look at what that does for our action. We have to change its element of time, so $dt \rightarrow -id\tau$, but see what happens...
$$
S = i\sum_{i=i}^N\frac{1}{2}m \frac{(\mathbf{x_i}-\mathbf{x_{i-1}})^2}{d\tau}+i\sum_{i=0}^NV(\mathbf{x_i})d\tau
$$
Well, here we have it! There's also a lot of convenient consequences, like, if we go back to the path integral, now the non-minimal actions don't just oscillate, they vanish exponentially, and that makes the integral converge much better. But in the process we lost the original connection between paths and dynamics! These paths we get out of optimising this action aren't real paths, they're paths in "imaginary time", which frankly sounds something out of a bad Doctor Who episode. So does this time have to do with anything? Well, check the part where we first performed the Wick rotation. That looks a lot like a partition function, right? In fact, it would totally be a partition function if we set $\tau = \hbar\beta$ ($\beta$ here being the usual inverse of the temperature times the Boltzmann constant). So there you have: imaginary time is inverse temperature. When you compute that path, above, you're not looking for a specific path in time, you're looking for a path at a given temperature, and the higher $T$ (the final time), the lower that temperature... uh... $T$ I guess (ok, I realise here I actually used some slightly confusing notation. Sorry for that). Turns out, your $k$ in the NEB objective function is exactly proportional to the temperature. Set it high, and particles will cut corners: they've got enough kinetic energy to do that. Set it low, and particles will just slide back in their potential basins: they can't leave them.
And that's why the NEB uses that objective function, and what its physical sense is.
