How to apply FIRE to many atoms where P=F·v seems to be a vector rather than a scalar?

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.

In my pure python implementation of FIRE (for educational purposes) in two dimensions, force and velocity vectors $\mathbf{F}$ and $\mathbf{v}$ are 2 x N arrays. P = (F * v).sum(axis=0) returns a one dimensional array of length N.

Do I go atom by atom and test its P(i) to decide which branch of the condition to use for it, or should I use just a single scalar P = (F * v).sum()?

Answer 1

Your code is wrong. Your excerpt of the FIRE paper states "a simple mixing of the global (3Natoms dimensional) velocity and force vectors", meaning $\mathbf{F}$ and $\mathbf{v}$ are both vectors with dimension $3N_\text{atoms}$.

$P(t) = \mathbf{F}(t) \cdot \mathbf{v}(t)$ is a scalar function. The correct implementation in Python is simply

P = numpy.dot(F, v)

when you use the proper array types.

Answer 2

As Susi noted, for your 2D case (as a subset of the more general 3D case) this dot product is assumed to be using the velocity/force written as $2N$ dimensional vectors rather than $(2,N)$ arrays.

However, it can certainly be convenient to have them as arrays when you want to more easily index into the array and grab the, say, the $y$ dimension of the 5th atom. If you want to keep the $(2,N)$ arrays while still being able to do the appropriate dot product, the quickest way is just to temporarily reshape your array as needed, e.g. P=np.dot(F.flatten(),v.flatten())

If performance is a concern, you should be able to use F.reshape(-1) (same for v) to obtain a view of the array with the correct dimensions rather than constructing a copy.