An Algorithm Predicting the Optimal Mechanical Response of Electronic Energy Difference

The use of mechanical forces at the molecular level has been shown to be an interesting tool for modulating different chemical and physical molecular properties. The so-called covalent mechanochemistry deals with the application of precise mechanical forces that induce specific changes in the structure, stability, reactivity, and other physical properties. The use of this kind of force to modulate photophysical properties and photochemical reactivity has also been studied. Nevertheless, the general problem of mechanical modulation of the energy gap between two electronic states has been addressed only with the development of simple theoretical models. Here, we develop and implement an algorithm providing the Largest energy Gap variation with Minimal mechanical Force (LGMF) that allows the determination of the optimal mechanical forces tuning the electronic energy gap, as well as to identify the maximum mechanical response of a molecular system to the application of any mechanical stimulus. The algorithm has been implemented for diverse molecular systems showing different degrees of flexibility. The phyton code of the algorithm is available in a public repository.

The algorithm developed in this study deals with the general problem of finding critical points for a differentiable scalar function, which is the case of a potential energy surfaces for a given electronic state or any combination of them.The most commonly used methods for finding these characteristic points are based on Newton's methods.The aim of these methods is to search for the zeros of a given function.For the stationary points, it is assumed that the gradient at that point is 0 within a convergence threshold.
Let us assume a scalar function () where  ∈ ℝ  .For determining the zeros of this function, we use the secant method.To do this, the function is expanded at the point   using the Taylor polynomial of first order.

𝐉(𝐱 𝟎 ) = ( 𝜕𝑔(𝐱 𝟎 )
1 ⋮ (  )   ) Known the () function, we look for the point that satisfies () = 0. 0 = (  ) + (  )( −   ) (S2) By solving, we get: In this way we obtain the value of the zero of the function.However, since we do not know the exact function (we have approximated it using a Taylor polynomial of order 1, see Eq. S1), it is necessary to check if this point is really a zero of the function.To do this, it is necessary to evaluate the function at the new point obtained using Eq.S1 and re-compute the point where () = 0 is fulfilled.Therefore, by iterating according to this procedure, a sequential approach to the real solution is fulfilled.Iteratively, the zero of the function is obtained: + =   − [(  )] −1 (  ) (S4) A zero of the function is obtained when ‖ + −   ‖ < , where  is a critical value that ensures the convergence of the procedure.
The search for critical points is based on the same principle.Let's assume a scalar function () as above.It is assumed that the function can be locally approximated to a quadratic function.
Therefore, by deriving Eq.S5 (  + ∆  ) = (  ) + (  )∆  (S6) Equating Eq. 6 give rise to: As can be seen the equation is similar to Eq. S4.The difference is that in this case we are looking for the zeros of a vector function : ℝ  → ℝ  instead of a scalar function : ℝ  → ℝ.So instead of having the matrix (  ) of (1xn) dimension, we have the matrix (  ) of (nxn) dimension.So, iteratively using Eq.S7, we will obtain the critical point of the function.
Newton's method, however, can be included in a group of optimisation methods known as linear search methods.In this group, the search for the minimum is carried out by performing a displacement known as the descent direction, which minimises the value of the function, converging to the minimum.
where (  ) is the search vector or direction of descent and   ∈ ℝ + is a constant known as the step length, which scales the displacement vector in order to ensure certain conditions and accelerate convergence.Depending on the descent vector and the step length constant,   , we can obtain different optimisation methods.
As shown in Fig. 2, for each of the iteration steps it is necessary to calculate both the gradient and the Hessian matrix.In cases where the calculation of the Hessian matrix is too timeconsuming, there are certain methods for updating the Hessian matrix using the displacement and gradient difference information.∆ −1 and gradient difference information (  ) − ( −1 ).
These methods are known as Quasi-Newton methods.One of the most widely used methods is the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method where   = ( + ) − (  ) y   = ∆  .According to Eq. S9, for the update to be possible and the process to yield a positive definite matrix  +1 defined positive matrix, the curvature condition must be      > .If necessary, to ensure the correct update of the Hessian, we can force the displacement to meet Wolfe's conditions [S1].To ensure that these conditions are met, the optimisation is performed using the equation where   allows us to scale the vector to ensure that the conditions are optimal for the Hessian update.
It can be shown that these conditions are satisfied for convex or locally convex functions.In the method proposed in this article, although we do not do an unconstrained optimisation leading to the minimum of the function, throughout the constrained optimisation, the points are located always around a minimum of a scalar function (PES).Near such minima the second-order approximations are adequate: such second-order surfaces fulfil the condition of convex functions.
Therefore, the BFGS method is a good method for updating the Hessian for displacements along local minima of PESs.In order to check how good the Hessian update is compared to the exact Hessian, we performed the computation of optimal forces using both.These results can be seen in figure S1.

Algorithm convergence criteria
The problem we address in this article is the conditional stationary point search problem.
Solving these problems can be done with the Lagrange multiplier method, which states that a constrained problem of n variables gives rise to an unconstrained problem of n+k variables where k is the number of constraints.For a scalar function () where { ∈ ℝ  } and a constraint () = , it can be shown that it is possible to construct a function (, ) =  −  which will be a linear combination of the two above.The restrained stationary points fulfil  = .It is possible to use the above method for searching for zeros and use Eq.S6 and S7 to iteratively find the solutions.However, the above condition implies that the solutions must satisfy where  is a scalar known as the Lagrange multiplier.Then, we will resort to linear search methods, making use of Eq.S8 to find those points that satisfy condition given by S11.In particular, as we discuss in the main text, the direction of displacement will be determined by the expression The linear search methods following Eq.8 are based on the fact that the direction of the displacement must minimize the function to be optimised.Therefore, they are methods used for the search of minima.But the actual problem is about finding conditional extrema, both minima and maxima.So to find local maxima, we must impose the opposite condition on the displacement vector, that is, maximizing the function.This means that, for the search for the local maximum, the displacement vector must be Once we have defined the displacement vector, the convergence have to be checked.By means of the iteratively linear search processes, different points that are sequentially closer to the solution are found, reaching a point that is close to it within a certain threshold.Therefore, the solution has been reached when ‖ + −   ‖ < , where  is a critical value that ensures meaningful convergence criteria.The solutions satisfy the condition of conditional extremum  =  which implies that ∆  ⊥ = .Therefore, it is possible to use as a convergence criteria: where  is again a threshold value.∆  ⊥ provides a measure of the angle between the vectors ∆  and ‖  ‖ 2 .It tends to 1, when ‖∆  ⊥ ( + )‖ becomes smaller enough.The convergence is reached when the cosine of the angle is close to 1.
In linear search methods, the choice of initial search direction is important, but also the modulus of the displacement.If the displacement is too small, convergence to the point will be slow.On the other hand, if it is too large, it does not converge to be solution, but rather provokes an oscillation along the desired point.This displacement is scaled by the constant   .There is a wide variety of methods for the calculation of the   which are based on the minimization of the function to be optimised at each step, in this case   (  − ∆  ⊥ (  ))[S1].However, in this case, we implement a constant parameter  = 0.0035 obtained by optimizing the efficiency of the algorithm in different tested molecular systems.The displacement then has been taken as:

Figure S1 .
Figure S1.Predictions made with the LGMF algorithm using updated hessians (red points)