@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()?

  • $\begingroup$ companion question: Understanding derivation of discretized FIRE algorithm $\endgroup$
    – uhoh
    Commented Sep 24, 2022 at 6:47
  • 1
    $\begingroup$ +1. I'm guessing that F is 2xN because there's N atoms and 2 spatial dimensions? $\endgroup$ Commented Sep 24, 2022 at 8:45
  • $\begingroup$ @NikeDattani yep! "In my pure python implementation of FIRE... in two dimensions..." It might be very interesting to do an N-dimensional analysis of 2 atoms in hypermolecular dynamics, but not today. $\endgroup$
    – uhoh
    Commented Sep 24, 2022 at 9:23

2 Answers 2


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.

  • 2
    $\begingroup$ @uhoh why are your arguments shape (2,N)? F and v should both be vectors with 3Natoms elements. $\endgroup$ Commented Sep 24, 2022 at 10:37
  • $\begingroup$ Thanks for the edit! Also thanks for being a more careful reader than I am it seems. Yes I understand now. $\endgroup$
    – uhoh
    Commented Sep 25, 2022 at 0:12
  • 1
    $\begingroup$ They're shape (2,N) because @uhoh is using a 2D model system, rather than a 3D one; although I agree that a vector of length 2Natoms would be simpler in this part of the code. $\endgroup$ Commented Oct 4, 2022 at 0:07
  • $\begingroup$ @PhilHasnip yes, it was later clarified that uhoh's code is two-dimensional even though this critical piece of information is missing from the unnecessarily long question; uhoh also stated in the comment above above that they understood that the real problem was that the algorithm was implemented wrong. P is supposed to be a scalar, like in the paper, and the vector dot product was implemented wrong. $\endgroup$ Commented Oct 5, 2022 at 8:45

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.

  • $\begingroup$ Thanks for your help, and for .reshape(-1). I have a simple, hexagonal lattice for the adlayer so each atom (except on the edges) has six bonds. I'm using numpy's advanced indexing to calculate all bond energies (simple Hooke's law force); accessing all the positions individually is needed every time force or energy has to be evaluated. Maybe some day I will try some compiled code (or just "upgrade" to Julia), but it's still a fun challenged to see just how fast it can be using everything numpy has to offer. $\endgroup$
    – uhoh
    Commented Sep 25, 2022 at 0:19
  • $\begingroup$ I'm slowly gaining ground learning LAMMPS as well, and comparing their results will also be instructional. $\endgroup$
    – uhoh
    Commented Sep 25, 2022 at 0:20

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .