Dear all,

I want to know how “quadratic line search” reduce forces on each atom to zero. I read cpp file and saw this:

relerr = fabs(1.0+(0.5*alpha*(alpha-alphaprev)*(fh+fhprev)-ecurrent)/engprev);

But I really don’t understand what is calculated by this equation and how it bring forces to zero.

Thank you very much.

Weina

I'll let Aidan respond to this Q.

Steve

The formula for relerr is a sanity check, to provide really bad choices

for alpha0. Often, it has not effect. I have changed the formula a little,

and added some comments explaining where it comes from. Steve will post a

patch shortly.

--- min_linesearch.cpp (revision 3487)

+++ min_linesearch.cpp (working copy)

+ // linemin: quadratic line search (adapted from Dennis and Schnabel)

+ // The objective function is approximated by a quadratic

+ // function in alpha, for sufficiently small alpha.

+ // This idea is the same as that used in the well-known secant

+ // method. However, since the change in the objective function

+ // (difference of two finite numbers) is not known as accurately

+ // as the gradient (which is close to zero), all the expressions

+ // are written in terms of gradients. In this way, we can converge

+ // the LAMMPS forces much closer to zero.

+ //

+ // We know E,Eprev,fh,fhprev. The Taylor series about alpha_prev

+ // truncated at the quadratic term is:

+ //

+ // E = Eprev - del_alpha*fhprev + (1/2)del_alpha^2*Hprev

+ //

+ // and

+ //

+ // fh = fhprev - del_alpha*Hprev

+ //

+ // where del_alpha = alpha-alpha_prev

+ //

+ // We solve these two equations for Hprev and E=Esolve, giving:

+ //

+ // Esolve = Eprev - del_alpha*(f+fprev)/2

+ //

+ // We define relerr to be:

+ //

+ // relerr = |(Esolve-E)/Eprev|

+ // = |1.0 - (0.5*del_alpha*(f+fprev)+E)/Eprev|

+ //

+ // If this is accurate to within a reasonable tolerance, then

+ // we go ahead and use a secant step to fh = 0:

+ //

+ // alpha0 = alpha - (alpha-alphaprev)*fh/delfh;

+ //

- relerr = fabs(1.0+(0.5*alpha*(alpha-alphaprev)*

- (fh+fhprev)-ecurrent)/engprev);

+ relerr = fabs(1.0-(0.5*(alpha-alphaprev)*

+ (fh+fhprev)+ecurrent)/engprev);

+ alpha0 = alpha - (alpha-alphaprev)*fh/delfh;

Aidan's update was part of the 5Dec09 patch.

Steve

I think I meant prevent, not provide. How strange.