Some Mathematical Reasoning on the Artificial Force Induced Reaction Method

There are works of the Maeda–Morokuma group, which propose the artificial force induced reaction (AFIR) method (Maeda et al., J. Comput. Chem. 2014, 35, 166 and 2018, 39, 233). We study this important method from a theoretical point of view. The understanding of the proposers does not use the barrier breakdown point of the AFIR parameter, which usually is half of the reaction path between the minimum and the transition state which is searched for. Based on a comparison with the theory of Newton trajectories, we could better understand the method. It allows us to follow along some reaction pathways from minimum to saddle point, or vice versa. We discuss some well‐known two‐dimensional test surfaces where we calculate full AFIR pathways. If one has special AFIR curves at hand, one can also study the behavior of the ansatz. © 2019 The Authors. Journal of Computational Chemistry published by Wiley Periodicals, Inc.


Introduction
Considerable interest is attached to the search of reaction pathways in chemistry, especially the points which govern these ways: minimums and saddle points of index one (SP 1 ) on the potential energy surface (PES) of a reaction system. The reaction pathway is defined as a one-dimensional description of a chemical reaction through a sequence of molecular geometries in an M-dimensional configuration space.
The AFIR method is an ansatz, which disturbs the given PES by an external force. [1,2] It is a generalized case of the treatment in mechanochemistry. [3][4][5] It has some similarity with the SEGO method (standard and enforced geometry optimization). [6] By the disturbance, one moves the stationary points of the former PES to new locations. By following the successive forcedisplaced stationary points, one gets a curve which can, in good cases, connect a minimum and a SP 1 by a kind of reaction path. The AFIR path has analogous properties. This paper has the following sections: next, we refer to the AFIR method, and we calculate a reaction path by pieces of a curve by consecutive AFIR points. A more theoretical tool is obtained by a variational formula for full AFIR curves. Further special properties like dependence of the AFIR curve on the coordinates, and avoided crossing (AC), are discussed separately by examples. At the end, we add a discussion and some conclusions.

The AFIR Method
The proposal of the Maeda-Morokuma group is to use an effective PES [1,2,[7][8][9][10] with internal coordinates r Ri + R j rij p r ij P N k < l R k + R l r kl p : Here V (r) is the original PES, α is a factor which plays the role of a numerical parameter which drives the calculation, R i and R j are covalent radii of atoms i and j. The vector r with the components r ij contains the distances between the corresponding atoms i and j of the studied chemical system. The dimension of all r ij is maximally M = N N−1 ð Þ 2 for a molecule with N atoms. It is possible to include a lower number of distances only. [1,11] Of course, all r ij > 0. To imagine the external force, f, directly, we write the components with numbers i, j and the effective PES is thus, the force, f, acts on the current point, r. If the force is zero, α = 0, then we have the original PES. Interesting are the cases with increasing amount of α. Equation (2) controls the relation between different bonds. It means that the larger distances nearly disappear in the extra force for p > 2 but only the smallest distances make a contribution to the resulting direction of the force. So, eventually, the small distances of H-atom bonds, which do not react in the system of question, should not be used in f. [7] Of course, if the extra force moves all stationary points of the PES out of their former places then a minimum and an SP can coalesce, and a former barrier can disappear. Such a situation occurs in a point labeled barrier breakdown point (BBP) with α = α max , and it is instructive to compare it for Newton trajectories (NTs). [12,13] So a new valley opens for a contact between former distant minimums. Thus, one can use the ansatz to detect reaction valleys. [1,2,9,11,[14][15][16][17][18][19][20][21] To the purpose, one has to choose α ≥ α max . The important relation is not discussed in the AFIR papers. If it holds then the former initial minimum disappears and a new minimum exists on the corresponding effective PES, which may be near a searched minimum of the original PES. Many examples are drawn recently for NTs. [4,5,13,[22][23][24][25] Here, in contrast, we propose to calculate the full "reaction path" between the initial minimum and the next SP of index 1 by AFIR. The aim will allow us to better understand the behavior of the AFIR method and its possible improvement. One starts at a known minimum with α = 0. Like for NTs, [4,5,26] a continuous increase of the strength of the force, α, will move the stationary points of the effective new PES. In the first papers to AFIR, Meada et al. only searched for an increase of the force. [1,10,11,14] In the last great review, [2] they propose to use only one fixed value of α. We first remark that in the case the value should be larger than α max because then the optimization on the effective PES does not go back near to the original minimum.
Here, we propose to improve the method by two alternating pieces of the curve of new stationary AFIR points. We propose to use an increase of the parameter, α, up to the BBP at α max of the AFIR curve, and a decrease of the parameter, α, after the BBP. Then the obtained curve points could fully describe the curve between two original stationary points over a BBP, like in the case of NTs. [4] The maximal α determines the BBP. Note that the BBP is not an approximation of the original SP of the PES. The BBP is usually anywhere between the initial minimum and the next SP, see some instructive discussions for NTs. [4] At the next stationary point the parameter α has again to converge to zero.
Since in the AFIR method only one test-α is used [2] this has to be greater than α max . Because then the optimization can jump along a new valley to a minimum near to the searched one. If it is test-α < α max then the optimization of the effective PES will get a minimum before the BBP, near to the original initial minimum. Figure 1 shows the result of calculations for changing values of α for the Rhee-Pande test surface. [4,22,27] We use (x, y) for only two abstract coordinates, thus dimension M = 2, and the exponent p = 6 in eq. (1) and put formally R 1 = 1/2 and R 2 = 1/2. Thus, the equation becomes The calculation goes step by step. We put a value of the parameter α near zero, and search the next stationary point of the effective PES, V eff . This point is used for the next α-step as the initial value for the optimization. We proceed in the same way after the convergence is reached. The AFIR points start at the right minimum, R at (4.03, 0.97) at zero energy, and they go with negative values of α up to α max = −9.98 at the BBP (pink) at point (4.7, 2.35). Then the parameter is again decreased (in absolute values). With the decreasing amounts of the parameter after the BBP, the curve at least correctly finds the rightmost SP at (4.25, 2.96). The same can continued in direction of the intermediate minimum; however, at the corresponding next BBP, the curve ends. There is a small gap to the intermediate where no AFIR points exist. An explanation of the fact comes below in the next section.
The other half of the reaction path from the product minimum, P at point (1.0, 4.0) and energy 3.64, is analogous. Up to a values of α = −6.5 at the leftmost BBP (pink), at (2.8, 5.15), we increase the parameter, but the way to the SP before the intermediate, at (3.0, 4.56), is again get by a decrease. And the piece between the leftmost SP of index 1 and the intermediate minimum at (3.76, 4.03) at energy 5.96 is a next ordinary part of an AFIR curve.
The procedure to go up by small steps for α and optimize the corresponding stationary point, up to a BBP, and then go back to zero for α to find an SP 1 , this procedure will work in every dimension M.

AFIR Curves by a Variational Formula
We use a first variational structure [28] of the AFIR model where g(r) is the gradient of the PES, α is the Lagrange multiplier and φ(r) is the derivation of the extra term, r r (f(r) T Á r). If we assume that φ(r) 6 ¼ 0, we can write the variation ansatz in another form where U is the unit matrix. Here, the task would be to derive the implicit tangent from the given term. However, the problem in the model ansatz, eq. (1), is the quite complicated expression of φ(r). To make an attempt to achieve an expression, we execute the following derivation. First we rewrite eq. (1) in the Now we define H(r): = r r g T (r) and G(r): = r r φ T (r). After some mathematical manipulations we obtain where r r α(r) has the form We note that if φ(r) is independent of r then eq. (6) reduces to the tangent of the reduced gradient following (RGF) model, [26,29] which is another version of the NT model. From this point of view, we can say that the AFIR method is a generalization of the NT model.
In the case of M = 2, we can use eq. (5) for a numeric search of a solution of the AFIR curves. We employ a Mathematica contour plot in Figure 2 for the zero contour of the square of the norm of the left hand side of eq. (5). The point by point optimized AFIR points of Figure 1 fit well in the resulting curves. The gap between the rightmost SP of index 1 and the intermediate is an avoided crossing of two AFIR curves.
A next problematic property is, at the other side of the minimums, for positive α values, that the AFIR curves escape into the mountains. They do not converge to the upper SP 1 at (1.59, 1.45). This SP is very higher in energy, it is 24 units, than the pathways through the intermediate where the left SP 1 is at 12.84 energy units. Thus, a reaction will proceed over the lower reaction path, and not over the upper one; however, from a theoretical point of view, one would like to know also the higher energy pathway. But the AFIR points from the two global minimums do not converge to this SP. An AFIR curve can start in this SP, but it connects the SP to the summit of the surface, an SP of index two. In contrast, steepest descent from the SP at (1.59, 1.45) will find the two global minimums.
Note that here the origin (0, 0) of the plane of the coordinates plays an exceptional role which is an artifact of this twodimensional test surface. (The origin was also excluded for the application of eq. (5).) It seems that not all AFIR curves through the zero point have a geometric sense. In the AFIR, ansatz of eq. (1) is used only atomic distances of the chemical system which cannot be zero. Thus, in real chemical calculations the zero of the distance coordinates does not appear.

Coordinate Dependence of AFIR Curves
Because of the nonlinearity of the AFIR ansatz, eq. (1), the resulting curves not only depend on the PES, but also depend on the used coordinates. We demonstrate it with a very simple test surface with one minimum and two SP of index 1, the Konda-Avdoshenko-Makarov surface. [12,13] It is a surface with two reaction pathways between a reactant and the exit. Across to the exit it has a very flat ridge. Unfortunately, the minimum is the zero of the coordinates system. Nevertheless, an AFIR curve in Figure 3a connects the minimum and the lower SP region, but near the lower SP it suffers from a small avoided crossing. Another AFIR curve follows nicely the flat ridge. Near point (0.9, 0.8) this part of an AFIR curve has a turning point (TP). In this example, all AFIR curve directions Figure 2. Test surface of Rhee and Pande [27] with AFIR curves (blue). The minimums and saddle points of index one are shown in black, the BBPs in red. [Color figure can be viewed at wileyonlinelibrary.com] Figure 3. (a) KAM test surface [12] with AFIR curves (blue). Green lines are the Det(H eff ) = 0 points of the surface. The crossing with an AFIR curve is a BBP of this curve. (b) The same surface but all coordinates moved by (2,2). By this transformation, we avoid the point (0, 0) in eq. (1). The new set of AFIR curves is different with respect to the previous set of AFIR curves. [Color figure can be viewed at wileyonlinelibrary.com] through the minimum are different from the eigenvalues of the Hessian. It seems that not all AFIR curves through the zero point have a geometric sense.
In Figure 3b, we use the same surface but moved symmetrically by a coordinate transformation by a linearly moved origin with a distance (2,2). The AFIR curves have a quite other form! But they reflect well the ridge of the surface. A former direct pathway from the minimum to the lower SP of index 1 is missing. The AFIR curve now connects indirectly both SP with the minimum over turning points.
The green lines in Figure 3 (and in the following figures) are the Det(H eff ) = 0 points of the surface where H eff is the Hessian matrix of the second derivatives of the effective PES with respect of the coordinates. There the gradient norm of the effective PES has a maximum or a minimum if one goes along the AFIR curve. The crossing of a green line with an AFIR curve is a BBP of this curve.
The next example has another set of scaled coordinates. Figure 4a shows the AFIR curves for a modified threeminimums surface [13,30] obtained by the variational formula, eq. (5). The modified surface is defined by The three minimums may mean one reactant minimum, Min R , and two different product minimums, Min P1 and Min P2 . The corresponding saddles are also so depicted. The example is chosen because the range of the coordinates is extended by a factor of 10. Nevertheless, here all AFIR curves suffer from avoided crossings (AC). No two stationary points are truly connected. Especially the Min R is far away from the origin, but it has a large AC to the SP 2 . AC means that the AFIR method can fail because the stationary points cannot coalesce on such separated AFIR curves. Figure 4b shows the AFIR curves for the same surface but changed coordinates by the symmetric distance (20,20). Now the picture again changes totally. The movement of the coordinates origin out of the global bowl improves the situation. AFIR curves connect one each minimum with one next SP, but not with the corresponding other SP. The Min R is connected to SP 1 , the Min P1 is connected to SP P , and the Min P2 is connected to SP 2 . Between SP P and SP 2 emerges an AFIR arc over the maximum. Again, some avoided crossings exist. No two minimums are truly connected. Figure 5a shows the AFIR curves for the Eckhardt test surface [31] obtained by the variational formula. The example is used because here the "forbidden" origin of the coordinates is the maximum of the surface reflecting to a certain degree the case of distance coordinates, r of eq. (1), where the case r = 0 is a really forbidden singularity. The Eckhardt surface is again a surface with two different reaction pathways between reactant and product minimums like the Rhee-Pande case. Here again Figure 4. (a) Three minimums [30] BQC test surface with blue AFIR curves. Green lines are the BBP curves with Det(H eff ) = 0. (b) The same surface moved by the coordinate pair (20,20). By this transformation, we avoid the point (0, 0) in eq. (1). Again, AFIR curves have strongly changed. [Color figure can be viewed at wileyonlinelibrary.com] Figure 5. (a) Eckhardt test surface [31] with blue AFIR curves. Green lines are the BBP curves Det(H eff ) = 0. (b) The same surface where the origin is moved by (3,3). By this transformation, we avoid the point (0, 0) in eq. (1). The AFIR curves qualitatively change. [Color figure can be viewed at wileyonlinelibrary.com] all AFIR curves suffer from avoided crossings. No two stationary points are truly connected. But again, the upper SP u seems to be more isolated than the lower one, SP l . Figure 5b shows the AFIR curves for a moved by (3, 3) Eckhardt surface. This situation is not better than under panel (a): the product minimum has no connection to other stationary points by an AFIR curve, the reactant minimum is connected to the lower SP l by a strange AFIR curve, but the upper SP u again is isolated.

Further Examples
Many AFIR curves show an AC. We could not assign any useful property of the PES to such ACs. It is in contrast to NTs. There the ACs indicate the neighborhood of a valley-ridge inflection (VRI) point which is crossed by a bifurcating, a singular NT. Singular NTs divide the "regions of influence" of the different stationary points. However, here, so to say, "singular" AFIR curves with a bifurcation are very seldom because these curves do not form a dense family of curves. They are unique curves. One cannot try to change the "search-direction" of the AFIR curve to get a nearby "singular" AFIR curve like a singular NT. The bifurcation of NTs is quite easier to calculate [33] and depends directly on the Hessian of the PES. Because of the nonlinearity of the ansatz of eq. (1), the connection to the effective Hessian will be quite more complicated. Figure 6 shows the AFIR curves for the well-known Müller-Braun surface [32] [MB] obtained by the variational formula, eq. (5). Here, only the main SP 1 and the intermediate minimum are connected by an AFIR curve. All other AFIR curves suffer from avoided crossings. It is in contrast to the case of NTs. [23] At every stationary point we detect one AFIR curve which coarsely follows an eigenvector direction of the Hessian. At two minimums, they follow the smaller eigenvalue direction, thus the "reaction valley," but at the third minimum, the corresponding AFIR curve follows the larger eigenvalue direction. The rule for this pattern is not clear. At the main SP 1 near point (−0.8, 0.65) the AFIR curve here crosses along the ridge, not along the reaction path direction.

Discussion
The examples demonstrate that the AFIR method can follow a valley from a minimum to an SP 1 , or vice versa, at least in good cases.
There are some specialties: 1. There are often gaps by an avoided crossing of the AFIR curves. The hypothetical bifurcation points inside an avoided crossing seem to have no geometrical meaning. In contrast, regular NTs connect the minimums with the SP 1 of the PES. Bifurcation points of NTs are valley-ridge inflection points. Additionally, an AC can destroy the planed action of the AFIR method. 2. The AFIR curve can have a turning point. This means that the curve touches a level line. Such behavior is also known from NTs. If a turning point emerges then the corresponding curve should not serve for a model of a reaction path because the TP has usually a higher energy than the next SP 1 . 3. A problematic property of the AFIR method, at least in the example of Figure 2, as bad as in others, is that here an unsatisfactory behavior emerges into the inverse directions of the two global minimums. The corresponding AFIR curves for positive values of α escape into the left and right mountains; however, they do not find the SP 1 at point (1.59, 1.45). Thus, not every SP 1 which is connected to a minimum by a steepest descent is also connected with this minimum by an AFIR curve. 4. Usually one AFIR curve leads through the stationary points, correspondingly, to positive or negative values of the parameter, α. It is again like for NTs, but there we can choose any direction which then is the leading direction of the NT. The NTs have a quite greater variability because around a stationary point all search directions are possible. The NTs form a dense net of curves on the PES. And the NTs are a linear ansatz, thus very easier to handle than the AFIR method. 5. A search for optimal BBPs [34] is not possible with AFIR curves because they have their fixed direction at every point. To determine an optimal direction, the search direction must be continuously changeable to determine the optimal NT.

Conclusions
In former applications, the AFIR method is handled as a "black box." It is not discussed that the α max of the BBP plays the decisive role for the planed action. The use of only a fixed value of the parameter α for a test calculation [2] where then it is hoped to find a next minimum more or less accidentally, gives away possibilities of this ansatz. In good cases, a consecutive use of small α-steps can follow a reaction path up to the searched SP 1 directly. But one has to be careful: α has to increase up to α max at a "barrier breakdown point" BBP and then to decrease back to zero at the next stationary point. However, the emergence of "avoided crossings" of AFIR curves can destroy their exploit-ability for a full reaction pathway.
The dependence on the coordinates of the external force makes the method somewhat tricky! The last conclusion is that one should better use the simpler Newton trajectories, thus a fixed, constant force, f, in eq. (3). NTs are better adapted to the task of the AFIR method. Figure 6. MB test surface [32] with blue AFIR curves. Green lines are the curves Det(H eff ) = 0. The crossing point between a green line and an AFIR curve is the BBP of this curve. [Color figure can be viewed at wileyonlinelibrary.com]