[lammps-users] Slow convergence of energy minimization

Dear Users,

I have here a question about system relaxation by means of energy minimization, i.e. regarding the Lammps minimize command.

My system contains several millions of Si atoms (I simulate with Stillinger-Weber potential), in which a specified concentration of defects in form of Frenkel-pairs are created. I intend to relax such a system to its lowest energy before studying further physical processes we are interested in. However, when I used minimize command together with style cg (or sd) to relax it, I found such a relaxation proceeds very very slowly if I set energy tolerance as 1e-11 or 1e-12; even when TOL = 1e-12, the actual dE/|E| is still around 5e-12 after 400,000 timesteps are performed. As we want to relax dozens of samples, such a slow speed is not that we can accept.

And something I found interesting is that, during most time of the relaxation, energy falls very very slowly; but sometimes, within dozens or hundreds of timesteps, energy can fall rather significantly--but as this quick energy descending is not dominating the process, it doesn't help much to speed up the relaxation.

So I want to ask scientists who are experienced in this lammps minimize: is the code relating to energy minimization in Lammps efficiently optimized or not? What is the physical reason of the quick energy descending sometimes happens within short time (perhaps because the code is written not in the most optimized way)? And any good suggestions for speeding up the relaxation process by means simulation setup optimization (I mean avoiding modifying the lammps code)?

Thanks and with best regards!

Xiang Gu

Dear Users,

dear xiang gu,

I have here a question about system relaxation by means of energy
minimization, i.e. regarding the Lammps minimize command.

My system contains several millions of Si atoms (I simulate with
Stillinger-Weber potential), in which a specified concentration of defects
in form of Frenkel-pairs are created. I intend to relax such a system to
its lowest energy before studying further physical processes we are
interested in. However, when I used minimize command together with style
cg (or sd) to relax it, I found such a relaxation proceeds very very

if you'd make a quick survey of the convergence characteristics
of a conjugate gradient and even more so of a steepest descend
algorithm, you'll see that none of them is well suited for your
project. from what i've seen, you best shot would be a linear
scaling BFGS variant. of course, one important question that you
should consider first is, how exhaustive isyour minimization
supposed should be, i.e. how tight you want to converge. any
minimization scheme will only lead you to the nearest local minimum
and depending how rugger you potential hypersurface is and how
energetically close various minimal are, you may be setting your
convergence parameter needlessly too tightly.

slowly if I set energy tolerance as 1e-11 or 1e-12; even when TOL = 1e-12,
the actual dE/|E| is still around 5e-12 after 400,000 timesteps are

...and you also should consider the limitations of double precision
floating point numbers. you have only 14 digits accuracy in your
mantissa. so going below a relative convergence of less than 1e-11 or
less for millions of particles is asking for a _lot_.

performed. As we want to relax dozens of samples, such a slow speed is not
that we can accept.

And something I found interesting is that, during most time of the
relaxation, energy falls very very slowly; but sometimes, within dozens or
hundreds of timesteps, energy can fall rather significantly--but as this
quick energy descending is not dominating the process, it doesn't help
much to speed up the relaxation.

that is just a sign of the characteristics of your potential
hypersurface in combination with the algorithm.

So I want to ask scientists who are experienced in this lammps
minimize: is the code relating to energy minimization in Lammps
efficiently optimized or not? What is the physical reason of the quick
energy descending sometimes happens within short time (perhaps because the
code is written not in the most optimized way)? And any good suggestions

optimization _cannot_ have anything to do with that, since
with or without optimization you should get the same result
(within the limitations of floating point accuracy).

for speeding up the relaxation process by means simulation setup
optimization (I mean avoiding modifying the lammps code)?

i have three suggestions:

a) don't go for an accuracy that is by far exeeding
the accuracy of your model and the quality of information
that the method you choose can provide.

b) you may want to try an MD with annealing and
a friction term as an alternative.

c) any "optimization" in that sense that you are looking for
would mean to implement a better algorithm and that requires
that you (or somebody else) get your hands "dirty",
i.e. modify lammps.

you have to keep in mind, that the purpose of minimizer algorithms
in most MD codes is simply to remove excessive potential
energy after assembling a system from fragments and there
CG and even SD are efficient enough, so that MD runs don't
crash or have atoms/molecules being "shot around" because of
excessively high forces.

cheers,
   axel.

I'm going to release some new code for minimization this
week, which might improve the situation. Does the
SW potential go smoothly to zero at the cutoff both
for energy and force? If not, this causes inherent
problems with minimization that affect convergence,
just like it does for a cutoff LJ or Coulomb potential.

Steve

Seems like discontinuity in the PES (or its curvature) would confuse the CG algorithm and cause it to terminate early (PE goes up when moving in the supposed downhill direction) or get stuck. My guess is that this is a real feature of his model since the algorithm doesn't terminate or get stuck.

I've never been able to quench binary LJ glasses of more than a few tens of thousands of particles to their Inherent Structures owing to effects just like Xiang Gu describes. In the binary LJ glasses, there are bursts of energy (and stress) relaxation events even when the PES gets very very (but not absolutely) flat. These relaxation events happen with CG or simple continuous gradient dynamics (with longer running time between relaxation events in the latter case). They are "physical". Changing the minimizer routine would not "help".

Doing real inertial continuous MD with fix_NVE and fix_viscous could help with converging to *some* energy minimum, but could change the minimum to which one would converge. This is probably no less artificial than tweaking parameters in the CG algorithm, anyway, since CG is not guaranteed to converge uniquely to anything like gradient dynamics is. One probably needs to think about the *physics* of the initial relaxation process for the problem of interest (temperature (if any), viscosity, etc...). CG is not a physical representation of the dynamics of anything.

I think Xiang Gu's relaxation events might be interesting in and of themselves.

My $0.02.

--Craig