Hi Steve & Aidan,
I'm getting a lot of failed linesearches in trying to minimize my
Lennard-Jones crystal. Going through the latest version of
linemin_quadratic() in min_linesearch.cpp, I couldn't understand the
'relerr' based criterion used for switching to quadratic mode.
Specifically, I can't understand the role of 'relerr' in following code
block (line 417):
I'll let Aidan respond to these Qs.
Steve
Hasan,
I am delighted that someone else is willing to take a close look at these
expressions. They originate with the idea that in the linesearch, 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.
So much for the math. In practice, many of the details were arrived in a
highly Edisonian manner. The final code does seem to work in most practical
cases, although it may be less robust that the backtrack linesearch in
situations where the forces do not exactly represent dE/dx, which can occur
for some potentials with discontinuities
Let's see if I can answer your detailed questions:
You are right. The expression for relerr is wrong. It should be:
de_ideal = -BACKTRACK_SLOPE*alpha*fdothall;
if (de_ideal >= -IDEAL_TOL) {
relerr = 0.0;
} else {
relerr = fabs(1.0+0.5*(alpha-alphaprev)*(fh+fhprev)/(ecurrent-engprev));
}
The derivation of this is given below. For whatever reason, both version
seem to work equally well. We are going to make further improvements to the
minimizer soon and will include that change.
The motivation behind IDEAL_TOL is that if the ideal energy change is too
small, then we do not expect to be able to measure the real energy reliably,
due to machine precision. IDEAL_TOL defines what we think is too small.
Finally, if you are trying to relax a crystal with fix box/relax aniso, with
different target pressures in x, y, and z, it is not unusual to encounter
difficulties. There are several reason for this, but it relates to the fact
that this case, unlike the isotropic case, can not be posed as a
minimization problem, unless the gradients include second-derivatives of
energy. We have some ideas for improvement, but have not tried them out yet.
Aidan
*** relerr derivation ***
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
The only unknown is Hprev, but we have two equations. So we can choose one
of E,Eprev,fh,fhprev as the second unknown, and compare the solution with
the known value. E seems like the best choice, giving:
Esolve = Eprev - del_alpha*(f+fprev)/2
relerr follows from that.
relerr = (Esolve-E)/(E-Eprev)
= -0.5*del_alpha*(f+fprev)/(E-Eprev) - 1