Understanding derivation of discretized FIRE algorithm

Question

Introduction/Preamble

@SusiLehtola's answer to Basics of numerical energy minimization techniques used in molecular dynamics? mentions conjugate gradients, BFGS for energy minimization, Metropolis Monte Carlo, and the FIRE algorithm for dynamical simulations.

The FIRE algorithm was introduced in the 2006 Phys. Rev. Letter Structural Relaxation Made Simple (Bitzek, Koskinen, Gähler, Moseler and Gumbsch, PRL 97, 170201, also available 1, 2, 3)

Consider a blind skier searching for the fastest way to the bottom of a valley in an unknown mountain range described by the potential energy landscape $E(x)$ with $x = (x_1,x_2)$. Assuming that the skier is able to retard and steer we would recommend him to use the following equation of motion:

$$\mathbf{\dot{v}}(t) = \mathbf{F}(t)/m - \gamma(t) |\mathbf{v}(t)| \left(\mathbf{\hat{v}}(t) - \mathbf{\hat{F}}(t) \right) \tag{1}\label{eq1} $$

with the mass m, the velocity $\mathbf{v} = \mathbf{\dot{x}}$, the force $\mathbf{F} = \nabla E(\mathbf{x})$, and hat for a unit vector. We recommend as strategy that the skier introduces acceleration in a direction that is ‘steeper’ than the current direction of motion via the function $\gamma(t)$ if the power $P(t) = \mathbf{F}(t) \cdot \mathbf{v}(t)$is positive, and in order to avoid uphill motion he simply stops as soon as the power becomes negative. On the other hand, $\gamma(t)$ should not be too large, because the current velocities carry information about the reasonable ‘average’ descent direction and energy scale.

and later:

The numerical treatment of the algorithm is simple. Any common MD integrator can be used as the basis for the propagation due to the conservative forces. The MD trajectory is continuously re-adjusted by two kinds of velocity modifications: a) the above mentioned immediate stop upon uphill motion and b) a simple mixing of the global ($3N_{atoms}$ dimensional) velocity and force vectors $\mathbf{v} \rightarrow (1-\alpha) \mathbf{v} + \alpha \mathbf{\hat{F}} |\mathbf{v}|$ resulting from an Euler-discretization of the last term in eq. (1) with time step $\Delta t$ and $\alpha = \gamma \Delta t$. Both $\Delta t$ and $\alpha$ are treated as dynamically adaptive quantities.

Explicitly, the FIRE algorithm uses the following propagation rules (given initial values for $\Delta t, \alpha = \alpha_{start}$ and for the global vectors $\mathbf{x}$ and $\mathbf{v}=0$):

• MD: calculate $\mathbf{x}, \mathbf{F} = -\nabla E(\mathbf{x})$, and $\mathbf{v}$ using any common MD integrator; check for convergence
• F1: calculate $P = \mathbf{F} \cdot \mathbf{v}$
• F2: set $\mathbf{v} \rightarrow (1-\alpha) \mathbf{v} + \alpha \mathbf{\hat{F}} |\mathbf{v}|$
• F3: if $P > 0$ and the number of steps since $P$ was negative is larger than $N_{min}$, increase the time step $\Delta t \rightarrow \min(\Delta t f_{inc}, \Delta t_{max})$ and decrease $\alpha \rightarrow \alpha f_{\alpha}$
• F4: if $P \le 0$, decrease time step $\Delta t \rightarrow \Delta t f_{dec}$, freeze the system $\mathbf{v} \rightarrow 0$ and set $\alpha$ back to $\alpha_{start}$
• F5: return to MD

In relaxation an accurate calculation of the atomic trajectories is not necessary, and the adaptive time step allows FIRE to increase $\Delta t$ until either the largest stable time step $\Delta t_{max}$ is reached, or an energy minimum along the current direction of motion ($P < 0$) is encountered.

Actual Question

Despite the lengthy preamble, my question is fairly short.

Question: Understanding discretized FIRE; how did the authors get from Equation $\eqref{eq1}$ to the recommended algorithm?

I've bolded

...Euler-discretization of the last term in eq. (1)...

and the Euler method is an extremely simple if not particularly accurate ODE solver. The last bit of the quote points out that particularly high accuracy is not necessary in order to simply find our way to the bottom; we don't care if we take the exact same path as a the skier to get there, but only that we do.

Once could alternatively solve Equation $\eqref{eq1}$ with a fancy ODE solver as well, but then one would have to construct an appropriate $\gamma(t)$ and and manage step sizes accordingly and that would probably not end up being FIRE anymore.

But here I'm only asking how the authors get from Equation $\eqref{eq1}$ to the recommended algorithm.

Answer 1

Using the backward Euler method to generate velocities from Eq. \ref{eq1}, we would arrive at the expression

$$ \mathbf{v}_{i+1}= \mathbf{v}_i+\bigg(\frac{\mathbf{F}_{i+1}}{m} - \gamma_{i+1} |\mathbf{v}_{i+1}| \left(\mathbf{\hat{v}}_{i+1} - \mathbf{\hat{F}}_{i+1} \right)\bigg)\Delta t $$ where I replaced the time variable with an index. It is probably possible to solve this expression for $\mathbf{v}_{i+1}$, but due to the dependence on $\mathbf{v}_{i+1}$ on the right-hand side, it's more complicated then we would like.

A simple way around this that should get us to roughly the same result is to use a predictor-corrector like scheme. In the first step we generate a guess for the next velocity using an MD integrator $$ \tilde{\mathbf{v}}_{i+1}=\mathbf{v}_i+\bigg(\frac{\mathbf{F}_{i+1}}{m}\bigg)\Delta t $$ where I'm writing the guess with a tilde. To correct this guess and move it closer to the desired FIRE velocity, we add in the remaining terms from Eq. \ref{eq1}, using $\tilde{\mathbf{v}}_{i+1}$ in place of $\mathbf{v}_{i+1}$ on the right hand side. $$ \mathbf{v}_{i+1}=\tilde{\mathbf{v}}_{i+1}-\bigg(\gamma_{i+1} |\tilde{\mathbf{v}}_{i+1}| \left(\tilde{\mathbf{\hat{v}}}_{i+1} - \mathbf{\hat{F}}_{i+1} \right)\bigg)\Delta t $$

After simplifying, this gives step F2 of their proposed algorithm. I don't know if this is how they derived their approach, but it seems to fit their analogy well, taking a trajectory and incrementally adjusting it to follow a steeper path.