[lammps-users] About high accuracy Conjugate Gradient

To whom it may concern:

I'm trying to resolve a high accuracy CG relaxation.
I set

minimize 1.0e-12 100000 500000
min_style cg

In the input, and relaxation is done. Then I increase the accuracy and tries to continue from there:

minimize 1.0e-15 100000 500000
min_style cg

To my suprise, the relaxation runs two steps and stop. In the log file, I can see the total energies are identical and the xyz positions of the atoms are also identical.

I checked the maximum force and found it is 0.00044 eV/Ang. It is a mistery that the position doesn't change at all.

Is there something I did wrong?

Thanks for the help.

Sincerely,

To whom it may concern:

I'm trying to resolve a high accuracy CG relaxation.
I set

minimize 1.0e-12 100000 500000
min_style cg

In the input, and relaxation is done. Then I increase the accuracy and
tries to continue from there:

minimize 1.0e-15 100000 500000
min_style cg

To my suprise, the relaxation runs two steps and stop. In the log file, I
can see the total energies are identical and the xyz positions of the
atoms are also identical.

I checked the maximum force and found it is 0.00044 eV/Ang. It is a
mistery that the position doesn't change at all.

Is there something I did wrong?

yep. you forgot about the limitations of floating point numbers!

unlike real numbers, floating point numbers are not continuous
and have only about 14 digits accuracy. furthermore, 0.00044 ev/A
is a very small force. no wonder there is no "visible" movement.

if you also factor in the errors introduced by cutoff radii and alike,
there is very little point in trying to optimize beyond a certain numerical
limit. perhaps, if you explain why you need to push so far, people
here may have an alternate suggestion.

cheers,
   axel.

Thanks Axel for your inpput.

I tried to relax a system to its known ground state. In the numerical CG runs, the final state has 0.02eV higher than the actual ground state. And the final two steps in numerical runs has almost no difference in total energy(<1.e-10 eV). This is not what I expected.

Besides, the equilibirum bond length from the numerical CG run is 1.e-4 Angstrom larger than the ground state.

I need at least 1.e-6 Angstrom in bond length accuracy to get the physics right.

Is there any other way to get around?

Thanks,

There are several stopping criteria for the CG minimizer, involving
energy and force and change from step-to-step. One of those
must have been triggered by your problem. I don't often run
CG to those high tolerances so I haven't experiienced the issue.

Steve

Hi, Steve:

From the document I read:

The minimization stops if any of several criteria are met:

     * the change in energy between outer iterations is less than the tolerance
     * the number of outer iterations exceeds maxiter
     * the number of force evaluations exceeds maxeval
     * the 3N dimensional force vector goes (nearly) to zero

And in the input I specify:

minimize 1.0e-15 100000 500000
min_style cg

The program runs step 0 and step 1 and then stops.
Since the maximum force is 0.0004eV/Ang, I guess only the first criteria is achieved.

* the change in energy between outer iterations is less than the
tolerance

I can see from the output that the change in energy is exactly ZERO.
But this is not acceptable since the atoms didn't move at all. Of course the change in energy between outer iterations is ZERO.

There must be something wrong.

Thanks for the help.

Like others were saying, you could have hit machine precision on the energy, in which case, you can't do much more.

On the other hand, I'd check to see if forces are really consistent with finite differences of energy. This has nipped me in the bud many times, and it's probably something that's not carefully checked in all the LAMMPS potentials. If you started at a point where the forces were pointing in an uphill direction on the energy landscape, or you were using a discontinuous potential that didn't go to zero at cutoff, it would probably give the behavior you describe after you restart. Are there any pairs which are close to the cutoff?

Also, it's not clear from the docs whether the condition:

* the 3N dimensional force vector goes (nearly) to zero

means that it stops on max(force), average(mag(force)).
If you've got a big system, all of the other forces are much less than the max, and the latter interpretation of the stopping condition is correct, then you might be stopping because of this condition and not the energy condition.

--Craig