Steepest descend try to optimize already optimized structures

[grhawk]

Hello,

@mahrossi @ceriottm

000_optimization contains the "real optimization". It took more than expected because I used a worse starting structure to be sure it works.

In 000_opt_optimized I used the previously optimized structure and I re-ran exactly the same optimization. It looks like the optimization never converges.

000_optimization.zip 000_opt_optimized_structure.zip

[mahrossi]

@litman90 since you did so many tests for the geometry optimizers, could you have a look at this issue?

[NirnaethArniedi]

It seems that in the SDOptimizer() class there is no update of the self.old_f object, so the forces differences evaluated in exitstep() never gets to zero so the code run until it gets to the maximum number of allowed steps.

[mahrossi]

@litman90 Could you look at this -- did you look at this?

[litman90]

I am having a look right now.

[mahrossi]

@NirnaethArniedi Could you provide inputs and outputs for what you are trying to do? In principle if everything was ok that condition shouldn't even be triggered.

[NirnaethArniedi]

I was running the example/lammps/paracetamol-geop files, using the sd optimisation. See input and output files attached.

paracetamol-geop-sd_output.zip paracetamol-geop-sd_input.zip

[mahrossi]

Okay, this may have to do with the tolerance of 1e-10 and perhaps some noise in lammps itself. The question is that the condition of making the force "flat" (i.e. that it doesn't change from the old to the new evaluation) is a bit strange. That usually doesn't happen if the numerical settings of the code you are evaluating the force in are tight. Could you confirm that if you use a tolerance of 1e-9 or 1e-8 the optimization goes well?

[ceriottm]

I agree, "force that does not change" is a weird requirement for convergence. I think we need to sanitize the bailout conditions throughout the GEOP routines.

On 8 May 2017 at 16:57, Mariana Rossi notifications@github.com wrote:

Okay, this may have to do with the tolerance of 1e-10 and perhaps some noise in lammps itself. The question is that the condition of making the force "flat" (i.e. that it doesn't change from the old to the new evaluation) is a bit strange. That usually doesn't happen if the numerical settings of the code you are evaluating the force in are tight. Could you confirm that if you use a tolerance of 1e-9 or 1e-8 the optimization goes well?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/cosmo-epfl/i-pi-dev/issues/174#issuecomment-299890981, or mute the thread https://github.com/notifications/unsubscribe-auth/ABESZwTLG2uH6pmyAUY0ipF2ImSxMLwiks5r3y12gaJpZM4Mj89D .

[NirnaethArniedi]

I ran the same geop with <tolerances> [1e-5, 1e-5, 1e-5] </tolerances> and <tolerances> [1e-8, 1e-8, 1e-8] </tolerances> in the i-pi input file without any changes in the final result (same final energy, simulation running until maximum number of steps is reached with system not changing at all after a few hundreds iterations).

In the output files it seems that those tolerances parameter are not taken into account and that the tolerances used are the ones defined by default in the GeopMotion class (that is 1e-06).

[mahrossi]

Alright, thanks a lot, we'll look into it further then. It is true that the forces just flatten though, right? (Meaning, they just stop changing). Does this happen with other optimizers too?

[NirnaethArniedi]

No for the other optimizers everything looks fine (you can look at the quick benchmarking report attached). From what I've seen the forces really flattens. Doing some debugging I have seen that at some point the init_step parameter passed to min_brent is null, and the method just return the initial point (what I understood is that with init_step=0 the bracketing return an interval reduced to the x0 initial point and then Brent's method just return this point). So the positions of the atoms are not moved any more and so the forces are not evolving either.

Benchmark.pdf

[mahrossi]

Perfect, great report and this may then be related to what Riccardo reported. What you report makes perfect sense (as in being a bug) and that should be the reason why this method is failing. Feel free to fix it if you see an immediate fix, otherwise I/Yair will change it in the next days.

I may add that the TRM method really benefits from an initial guess for the Hessian that is minimally better than identity (which is the current default). With that the TRM method should/could really win over all others. As an example, I quickly googled the documentation of molpro: http://www.molpro.net/info/2010.1/doc/manual/node568.html . We have a routine that we could add to the i-PI tools to produce Lindh hessians. They are not parametrized for all elements though.

[ceriottm]

I second the idea of having a preconditioner. Something easy I was thinking is using the VdW radii to build the Hessian

On 8 May 2017 at 18:14, Mariana Rossi notifications@github.com wrote:

Perfect, great report and this may then be related to what Riccardo reported. What you report makes perfect sense (as in being a bug) and that should be the reason why this method is failing. Feel free to fix it if you see an immediate fix, otherwise I/Yair will change it in the next days.

I may add that the TRM method really benefits from an initial guess for the Hessian that is minimally better than identity (which is the current default). With that the TRM method should/could really win over all others. As an example, I quickly googled the documentation of molpro: http://www.molpro.net/info/2010.1/doc/manual/node568.html . We have a routine that we could add to the i-PI tools to produce Lindh hessians. They are not parametrized for all elements though.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/cosmo-epfl/i-pi-dev/issues/174#issuecomment-299914186, or mute the thread https://github.com/notifications/unsubscribe-auth/ABESZ69xe-Bi0KRy0GuVR4ldPtZS0qBHks5r3z9-gaJpZM4Mj89D .

[litman90]

Thank you @NirnaethArniedi for finding that bug in min_brent, but I think it is a signal for a bug in another place. I join to @mahrossi comment if you find it, feel free to fix it, otherwise, we fix it soon. Regarding tolerances, that part works properly. You are just not using the correct format ;) The correct one is for example:

1e-3 1e-5 1e-4

I am not completly sure why i-pi doesn't complain about the wrong format.

[litman90]

< tolerances > < position > 1e-3 < /position > < energy > 1e-5 < /energy > < force > 1e-4 < /force > < /tolerances >

[NirnaethArniedi]

Thank you for the tolerances, I will try that.

I agree with the fact that the bug doesn't come from min_brent, since it is used in the CGOptimizer were everything works fine. I just managed to track the origin of flattened forces up to this point, but it could be normal (if we are at an energy minimum it is logical that nothing moves afterwards) and only a problem in the bailout conditions.

[litman90]

Note that in the CG test the "force that does not change" is triggered for the exit. The problem is that the SD and CG use the energy difference to find an approximate minimum inside each step (the 1-D line search or to be more specific inside "min_brent" loop). After a certain point, the displacement is so small and the energy difference reaches the noise level (this is system/ force evaluation code dependent). In this point, the 1D-line search doesn't work anymore. This doesn't happen with the other methods like bfgs, lbfgs or bfgstrm since they use gradient information in the approximate minimization in each step. Maybe we have to change that algorithm. On top of that, the "zeps" value in the min_brent function wasn't small enough for the tolerance needed. I tested with the H2O qtip4pf and pswater and I reached "| max component force |" < 10-8. (There is a little bug in qtip4pf that I will report in another issue). I can't reach that value with the paracetamol crystal. After getting an energy difference of ~1e-11 and a displacement of ~1e-07 the algorithm breaks. However, the maximum force component is ~10^-4 (5e-3ev/a) which is quite high. I have just opened a new branch for this issue.