I am currently trying to implement a GS2 (Gonzalez-Schlegel second order) IRC algorithm in a python code. I am following the original paper ref(1).

The main problem is in the constrained optimisation step of the GS2 algorithm. In this step, the geometry optimisation is done by holding the distance of the coordinates (euclidean distance) constant from the coordinates of a pivot point (fixed throughout optimisation).

I will reproduce the main equations for the GS2 constrained optimisation step here: $$p' = x' - x^* \tag{1}$$ $$\Delta x = x - x'\tag{2}$$ $$p = p' + \Delta x = x - x^* \tag{3}$$

where $x'$ and $x$ are the old and new coordinates for the constrained optimisation; $x^*$ is the pivot point that is constant. The $p'$ and $p$ similarly are the vectors from the old and new coordinates to the pivot point, respectively.

The energy is minimised under the constraint $p^T\cdot p = k$ (where k is some fixed scalar number). This is equivalent to optimization on the surface of a hypersphere. They express the optimisation problem with a Lagrange multiplier form: $$L = E' + g^T \cdot \Delta x + \frac{1}{2} \Delta x^T\cdot H \cdot \Delta x - \frac{1}{2} \lambda [p^T\cdot p - k]\tag{4}$$

Where the potential energy surface is expressed as a truncated Taylor series (as normal for QM stuff) and the $g$ and $H$ are the gradient vector and Hessian matrix as usual.

What I don't understand here is the half sign in front of the Lagrange muliplier $\lambda$.

Also, after this, they mention that at convergence, that $$\Delta x = - (H - \lambda I)^{-1} \cdot (g- \lambda p')\tag{5}$$

They then put that expression into the expression for $p$ (eqn. 3) and find the $\lambda$ (by doing a bracketing 1D scalar root search) for which the $\Delta x$ satisfies the constraint.

How does this expression follow from the above expression? I don't quite understand. Could anyone explain how this derivation works mathematically?

Also, another problem with these equations is that the step size is too high in most iterations so there is unpredictable behaviour. Is there any way to restrict the step size for the constrained optimisation step? (In the paper, they have a linear interpolation formula to stabilize the quasi-NR step, but it does not help too much)

References: (1) C. Gonzalez, H. B. Schlegel, J Chem. Phys., 1989, 90(4), 2155

  • 2
    $\begingroup$ Perhaps you forgot to add the reference paper (1) $\endgroup$ Commented Mar 18, 2023 at 19:27
  • $\begingroup$ @VandanRevanur Thanks for catching that! fixed it. $\endgroup$
    – S R Maiti
    Commented Mar 19, 2023 at 9:54

1 Answer 1


I'll propose a slightly different Lagrangian expression:

$$ L = E' + g^T \Delta x + \frac{1}{2} \Delta x^T H \Delta x + \frac{1}{2} \gamma \Delta x^T \Delta x - \frac{1}{2} \lambda (p^T p - k). \tag{1}$$

Here $\gamma$ is really a spring constant (with the corresponding units) that penalizes large norms for $\Delta x$. The sign of $\lambda$ doesn't really matter since it's enforcing a constraint, so we keep it as per the expression given above in equation 4.

Now, we would like to minimize the Lagrangian with respect to $\Delta x$, and we will do so by completing squares. First we substitute in $p = p' + \Delta x$ (equation 3), and collect terms in like powers of $\Delta x$:

$$ L = \frac{1}{2}(\Delta x^T[H + \gamma I - \lambda I]\Delta x) - (\lambda p' - g)^T \Delta x + \cdots \tag{2}$$

where the dotted terms are constant with respect to $\Delta x$. The term linear in $\Delta x$ combines a contribution from the forces $g$ and a contribution from the $p^Tp$ term of equation 4.

Now hark back to your high school algebra: a quadratic expression $ax^2 - bx + c$ is minimized at $x^* = (2a)^{-1} b$ (by completing the square), and this is no different:

$$ \Delta x^* \equiv \mathrm{arg min}\,L(\Delta x) = (H + \gamma I - \lambda I)^{-1}(\lambda p' - g) \tag{3}$$

Other than my proposed modification (equivalently, when $\gamma = 0$) this expression is exactly equation 5.

So hopefully this answers your three questions:

  1. Why half $\lambda$? Scaling your quadratic term by half leaves you with just $\lambda$ in the final answer, rather than $2 \lambda$. It's purely aesthetic -- but it also does leave you with one less factor of two to muck up when you're programming.

  2. How to get to equation 5? Expand in terms of $\Delta x$ and complete the square.

  3. How to constrain the step size? To me the most natural way is to put in an explicit norm restraint. From my final equation you can see that when $\gamma \gg \lambda$ and $\gamma \gg ||H||$, $\Delta x^*$ will just be $(\lambda p' - g)/\gamma$ -- that is, you'll simply be following the forces downwards, in ever decreasing step size as $\gamma$ increases.

  • $\begingroup$ Thanks! This is very useful to me. Could you please explain how you get from eqn.(1) to (2) in your answer? I am afraid I don't have much experience with vector algebra. Also, is there a way to solve for lambda and gamma directly? And how would you deal with cases where the Hessian is not positive definite? $\endgroup$
    – S R Maiti
    Commented Mar 20, 2023 at 10:00
  • $\begingroup$ @SRMaiti (1) have you tried deriving it yourself? You need to use a combination of expanding squares (which works just like expanding $(a+b)^2$ in normal algebra), recognising that for any vector $v$, $v^T v = v^T I v$, and grouping like matrix multiplications. Gamma is a constant you choose. If it is large enough the effective Hessian $H + \gamma I$ is always positive definite (answering your last question). It looks like I reinvented Hessian modification (nmr-relax.com/manual/Hessian_modifications.html). $\endgroup$ Commented Mar 21, 2023 at 9:32
  • $\begingroup$ Lambda is calculated by re-imposing the constraint $p^T p = k$ on the solution for $ \Delta x^*$ and finding the value which permits that. You will usually need some sort of iteration to accomplish this. $\endgroup$ Commented Mar 21, 2023 at 9:35
  • $\begingroup$ Right, I know that lambda is calculated by reimposing the constraint. That's the original method written in GS2 paper. They do it by a Newton root search, by starting from a value that's lower than the first eigenvalue of the Hessian, which ensures the level-shifted Hessian is positive definite. But how do I do root search with two variables and two constraints? Do they both need to be less than the first eigenvalue, or is it ok if the sum is less? Because lambda also appears in the gradient term. $\endgroup$
    – S R Maiti
    Commented Mar 21, 2023 at 10:27
  • 1
    $\begingroup$ Once again, $\gamma$ is a constant for each minimisation step. It is a number you specify. So if you use equation (3) of my answer you have one (modified) constraint in one variable, $\lambda$, not two constraints in two variables. You may find it easier to think of this as a Hessian modification $H \rightarrow H+\gamma I$, and if $\gamma$ is larger than the most negative (or zero) eigenvalue of $H$, then $H+\gamma I$ is guaranteed positive definite (and clearly still symmetric). $\endgroup$ Commented Mar 21, 2023 at 11:30

You must log in to answer this question.

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