DM converges to wrong state in MD simulations

[seungchul]

Dear Developers,

In my MD simulation, I found that energy and atomic forces jump suddenly (e.g. few tens of eV between adjacent step).
At certain step, starting DM is far different from correct one, and it converges to some wrong state.
I have checked this by comparing with energies from "Single Point Calculation" with the atomic structures at that MD step.
You can check in 'compare_eng.dat' (attached).

The attachment contains, fdf and stdout of MD simulations, and stdout files of several Single Point Calculation with structures from MD.

I've used
-- SIESTA 4.0 stable version
-- Linux cluster
-- Intel Composer XE 2013.5.192
-- Open MPI

Thank you very much,

[Nick Papior (nickpapior)]

Could you try:

1) increase mesh-cutoff, the forces are sensitive to the mesh-cutoff for various reasons. I would suggest you try around 250-300 Ry.
2) Your DM.tolerance is relatively high. It would perhaps be ideal to use the 1e-4
3) Why don't you have the slab-dipole correction turned on?

[seungchul (k-seungchul)]

Dear Nick Papior,

Thank you very much for your comment.
I've done several tests. The results are...

1. As you suggested, I've used more robust parameters for SCF convergence (e.g. larger mesh cutoff, small DM.tolerence). It gives better results, but I found sudden jump in E_KS. And of course, it becomes very slow.

2. Mesh-cutoff is less important than tolerance (DM, H, FreeE) for this problem.

3. In fact, in MD simulations, we usually use smaller mesh-cutoff and larger tolerance, due to the speed. So, it is not good idea to use similar parameters for usual structural relaxation and electronic structure calculation.

4. The solution I found is changing program slightly.
We can check if the E_KS jumps suddenly (i.e. unphysically). For example, it is nonsense that E_KS changes few electron volts between successive MD step.
I re-initialized DM when E_KS changes too much. (This is based on the observation that 'Single Point Calculation' gives reasonable energies.)
In this way we can do fast MD simulations without strange behaviour in energy.

Thank you very much.

SKim

[Nick Papior (nickpapior)]

Thanks for your feedback.

1) Slower, yes, but also more accurate. It is of course a trade-off. But an important aspect of MD simulations when performing MD calculations.

3) Yes, I see why you do it, but if the forces/energies are "incorrect" I think stabilising the simulation via increased mesh-cutoff would be beneficial.

Regarding 4), did you change the source to restart DM from atomic fillings if the change E_KS was too high? Or did you do this manually? If you changed the code, could you attach a patch here for us to review, perhaps we might consider adding it to the code base. It sounds like a possibly good idea.

[seungchul (k-seungchul)]

Dear Nick Papior,

After more test, I found that such a sudden change is due to 'subroutine post_scf_work'.
If you check MD.out file in the attachment of the first posting, you will see that E_KS at the last scf step and printed E_KS (which is 'corrected' by subroutine post_scf_work) are very different when the E_KS jumps.
I am attaching data from other case; same structure from the my first post, but different convergence parameters.

It seems that such behaviour happens less when better convergence parameters used and more GridCellSampling used.

I will attach my modification of SIESTA after few more test.

Thank you very much.

[Nick Papior (nickpapior)]

Could you try with the latest 4.0 commit:

https://bazaar.launchpad.net/~siesta-maint/siesta/rel-4.0/tarball/559

I think some of these issues may be fixed through some other fixes.

Please return with your findings.