I am attempting to implement an algorithm for transition state finding (in Python, but that's not important). The algorithm is called BITSS (Binary Image-pair Transition State Search). It defines a composite potential from reactant and product structures: $$E_\mathrm{BITSS} = E_1 + E_2 + k_e(E_1-E_2)^2 + k_d(d-d_i)^2$$

Where 1 and 2 are reactant and product respectively. $E$ defines the molecular energy and d is the distance between them (defined as RMSD or something similar). The $d_i$ is the target distance criteria. Then the composite potential is minimized across all Cartesian coordinates of both images, while the distance criteria is tightened in steps. This essentially means that there is an artificial force that pushes the reactant and product towards each other (from the distance term), while not allowing one to jump over the barrier (due to the energy constraint that ensures the two are at similar energies). At convergence, both structures (images) reach the transition state.

Now, in the paper, they only provide the equation for BITSS energy. I have been able to derive an analytic expression for the BITSS gradient myself, from the individual molecular gradients. However, for a BFGS or RFO like optimisation, I also need the Hessian of the BITSS energy, which has a very complicated expression.

I remember reading somewhere, that if the optimisation is started with an identity matrix (1's on diagonal and 0 elsewhere) as Hessian guess, and then BFGS hessian updates are performed, then the Hessian will automatically get closer to the real value due to getting the gradient information in BFGS update.

Is this true? Is it fine to start an optimisation with an identity matrix Hessian guess and then rely on the Hessian update procedure? I am not optimising a molecule per se, but a potential defined from two molecules. I am using Cartesian coordinates throughout.

(Also, I have tried gradient descent steps, but they are too slow to converge)

Reference - Samuel J. Avis, Jack R. Panter, and Halim Kusumaatmaja , "A robust and memory-efficient transition state search method for complex energy landscapes", J. Chem. Phys. 157, 124107 (2022) https://doi.org/10.1063/5.0102145

  • $\begingroup$ Have you tried conjugate gradients? $\endgroup$ Jan 18, 2023 at 20:48
  • $\begingroup$ @ShernRenTee I think it will be too slow. This is not using forcefields, but using QM calculations, so every step is very costly. I have used steepest descent as I wrote and it was quite slow. I don't think conjugate gradient will be much faster than that. $\endgroup$
    – S R Maiti
    Jan 18, 2023 at 21:20
  • 1
    $\begingroup$ Algorithmically conjugate gradients can be appreciably faster than steepest descent. Why not just try it if? CG is just a few lines of code, all but trivial to implement, if it is not fast enough for your test cases you can then look into something more complicated. $\endgroup$
    – Ian Bush
    Jan 19, 2023 at 7:21
  • $\begingroup$ @IanBush I already have the code for RFO optimisation, the only remaining problem was handling the initial Hessian. I guess I could try CG, but I have high doubts it can surpass RFO. Edit: I am working on a pre-existing codebase which implements RFO, BFGS, Steepest Descent etc. I am using those algorithms on my potential, so it is easier to use those implementations than writing something from scratch myself. $\endgroup$
    – S R Maiti
    Jan 19, 2023 at 10:11
  • 1
    $\begingroup$ If all you are interested in is the minimum, it does not matter how you get there. You can use BFGS with an initial Hessian initialized to the identity matrix; however, if the spectrum of the Hessian is very broad, the minimization may take forever to converge. I would suggest looking into using automatic differentiation for the gradient and Hessian; it will make it way easier to get a working algorithm in place. You may even find that you don't need to determine the Hessian in closed form... $\endgroup$ Jan 22, 2023 at 18:18

1 Answer 1


This is a difficult question to answer, because it really depends what you mean by "OK". There is certainly no fundamental problem with an initial guess Hessian which is the identity matrix, and any quasi-Newton method will attempt to improve this iteratively. The improvements are done by approximating the second derivative matrix (Hessian) based on information from two consecutive force evaluations, by a term which is like the change in force divided by the change in configuration. I say "like" because both force and configuration are $3N$-dimensional vectors, so the division operation is not uniquely defined. In most schemes the Hessian update is explicit, but some schemes (e.g. conjugate gradients) use an implicit Hessian update.

Although there is no fundamental technical problem to this choice of initial Hessian, the performance and even stability of the method is likely to depend quite strongly on how reasonable your initial Hessian is. Quasi-Newton methods assume your system is "near" to the minimum, in the sense that the energy landscape is well approximated by a quadratic (quasi-Newton methods ignore cubic and higher-order terms), and it's important that your system stays within this quadratic region. If you stray into areas of configuration space where the energy surface is highly non-quadratic and/or nearer to a maximum, then the Hessian updates can become dominated by the higher-order terms (effectively noise) which leads to poorer and poorer Hessian approximations.

The efficiency of an optimisation method is generally related to what's called the condition number, $\theta$, which is the ratio between the largest and smallest eigenvalues ($\lambda_i$) of the Hessian,

$$ \theta = \left\vert\frac{\lambda_\mathrm{max}}{\lambda_\mathrm{min}}\right\vert. \tag{1} $$ If $\theta=1$, then even a steepest descent method should converge to the minimum in a single iteration. More generally, the number of iterations steepest descent methods take scales as $\theta$, whereas the number of iterations quasi-Newton methods take scales as $\sqrt{\theta}$ (although these scaling arguments assume the initial configuration is in the quadratic well).

Choosing a good initial Hessian is equivalent to preconditioning, a technique to improve the condition number by changing the coordinate system to one in which the Hessian is "more spherical" (i.e. the eigenvalues are closer to each other). The identity matrix doesn't change the coordinates at all, and so it doesn't help (but it doesn't hurt either); however, even a reasonable guess of the extreme eigenvalues might improve the condition number by a factor of 10 or so, resulting in the method only needing approximately a third ($\frac{1}{\sqrt{10}}$) of the iterations and also being more stable.

In the case of molecules, a large condition number can arise because stretching or compressing a bond is typically much more energetic than bond-bending or twisting; i.e. the diagonal elements of the Hessian corresponding to displacement a bond are large, but those relating to tangential displacements are small - a classic case of a poorly conditioned problem.

There are lots of ways you could think of preconditioning it, you could estimate the bond directions and the ratio of the energy curvature along them to a tangential direction. One approach which can be very effective is to use a simple proxy model forcefield, like Lennard-Jones, and calculate its inverse Hessian. It won't be correct, but it has some of the right behaviour and might improve the condition number.

You might find it helpful to look at the work of Christoph Ortner & James Kermode in this area, which includes work on Hessians for transition state searches:

A "universal preconditioner" for materials' Hessians, Packwood et al, https://doi.org/10.1063/1.4947024

Preconditioning NEB transition state searches, Makri et al, https://doi.org/10.1063/1.5064465 .

  • $\begingroup$ Thanks for the detailed answer! Is there a Hessian update method that is stable even if the surface is non-quadratic? I am doing transition state search from the reactant and product, so there is no doubt that it will go through regions where PES is not quadratic. $\endgroup$
    – S R Maiti
    Jan 31, 2023 at 10:06
  • 2
    $\begingroup$ @SRMaiti you could investigate the Stabilised Quasi-Newton Method (SQNM) arxiv.org/pdf/1412.0935.pdf or there are also a range of methods called Trust Region Methods. In principle, you could actually combine the two, but I'm not aware of anyone having done that. $\endgroup$ Feb 2, 2023 at 2:00
  • 2
    $\begingroup$ +1 I think many molecular packages use the Lindh model Hessian which is essentially a very simple force field, but is clearly better than the identity matrix. I've been surprised there aren't more examples considering a better model Hessian would help in TS searches for the reasons you mention. $\endgroup$ Feb 2, 2023 at 14:56
  • 1
    $\begingroup$ @GeoffHutchison For materials, we've found the method of Pfrommer et al. to be pretty good: doi.org/10.1006/jcph.1996.5612. There's also the "universal forcefield" of Soler, and the more recent work of Kermode et al, including TS searches doi.org/10.1063/1.5064465 $\endgroup$ Feb 4, 2023 at 0:28

You must log in to answer this question.

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