Performance: Avoid magnetic calculations being trapped in local minima

[WHUweiqingzhou]

Background

To end it, we need to clarify some issue:

1. how to tell SCF has been trapped in local minima?
2. how to help SCF jump over it?

Any discussion is welcomed.

Describe the solution you'd like

Design a routine to avoid magnetic calculations being trapped in local minima.

Task list only for developers

• [ ] Notice possible changes of behavior
• [ ] Explain the changes of codes in core modules of ESolver, HSolver, ElecState, Hamilt, Operator or Psi

Notice Possible Changes of Behavior (Reminder only for developers)

No response

Notice any changes of core modules (Reminder only for developers)

No response

Notice Possible Changes of Core Modules (Reminder only for developers)

No response

Additional Context

No response

Task list for Issue attackers (only for developers)

• [ ] Review and understand the proposed feature and its importance.
• [ ] Research on the existing solutions and relevant research articles/resources.
• [ ] Discuss with the team to evaluate the feasibility of implementing the feature.
• [ ] Create a design document outlining the proposed solution and implementation details.
• [ ] Get feedback from the team on the design document.
• [ ] Develop the feature following the agreed design.
• [ ] Write unit tests and integration tests for the feature.
• [ ] Update the documentation to include the new feature.
• [ ] Perform code review and address any issues.
• [ ] Merge the feature into the main branch.
• [ ] Monitor for any issues or bugs reported by users after the feature is released.
• [ ] Address any issues or bugs reported by users and continuously improve the feature.

[kirk0830]

I have some experience on handling magnetic system such as Fe and/or Co-contained systems. First I would like to share reviewer's comments on my previous work and my answer: Comment:
...Details regarding the computational modeling such as spin-states, alignment of spins at Fe and Co centers, terminations of the CoFe2O4 (110) facets need to be provided and properly treated... Answer:
... Because temperatures of reaction and calcination (500 oC) were close and far above Curie point respectively (~650 K, Figure R12a, J. Magn. Magn. Mater. 2009, 321(9): 1251-1255), at first supports in models we planned to construct were assumed to be paramagnetic, which means symmetry-breaking was not manually set in wavefunction optimization. We fully agree with you that spin-alignment is important for this kind of ferrimagnetic material, thus after re-consideration, we think because all our characterization have been performed at room temperature and our reaction temperature indeed has an non-negligible range below Curie point, symmetry broken system is a more appropriate choice for discussion. Therefore, we manually broke the symmetry in the revised manuscript, setting all magnetization of Co atoms as up and that of Fe atoms as down, according to published works (J. Phys.: Condens Matter. 2011, 23(42): 426004)...

The most important thing is, Actually, it is NOT appropriate to run a calculation with those possibly-to-be-magnetic atoms but without any priori knowledge. The first thing one should do is to check whether it is AFM, FM or with some other kinds of magnetic ordering, (and even Curie temperature) then rationally set starting magnetization, or it is very possible to be heavily criticized by reviewers. Therefore set starting magnetization precisely as much as possible is what user should do.

Then I try to directly answer questions:

1. try with other starting magnetization, then compare magnetic or charge population analysis results and energies, then one will get to know whether it is a system with several minima.
2. if we confine in context of the question "let SCF iteration jump out of the local minima", there are more subtle aspects, the first thing I would like to declare is, SCF can only converge to a point where the first derivative is zero, especially for algorithms like diagonalization, but for the second order derivatives which determines whether a solution is a local minimum or saddle point, it is actually a question of wavefunction stability. There are several methods dealing with unstable wavefunctions. @jinzx10 has notes for it. However, to really find a global minimum, it is a rather harder one compared with the stability and it is included in the field of "global optimization". A momentum factor would be of help but there are empirical parameters for it, @jinzx10 also find there is a technique called SCF metadynamics, which is analogous to the classical metadynamics during which more and more Gaussian hills are added to extended Largrangian. I also read some papers they set initial guess of local occupation matrix, like what ABACUS may support the omc keyword, as reported it would help UO2 system converge to state with lower energy, it then would be the content we used to discussed at the beginning of ABACUS new DFT+U implementation/development.

[WHUweiqingzhou]

@dyzheng and @hongriTianqi I will try so-called U-Ramping method to improve the convergence of DFTU. To this end, I will

1. add a new input parameter uramping
2. create a vector dftu.U0 to copy GlobalC::dftu.U in DFTU, and set GlobalC::dftu.U zero value.
3. modify GlobalC::dftu.U in ESolver_KS_LCAO::iter_init when SCF mixing restarts:

   
   if(mixing_restart)
   {
   for(int i = 0; i < GlobalC::dftu.U.size(); i++)
   {
   if (GlobalC::dftu.U[i] + udamping < GlobalC::dftu.U0[i] ) 
   {
       GlobalC::dftu.U[i] += udamping
   }
   else
   {
       GlobalC::dftu.U[i] = GlobalC::dftu.U0[i]
   }

} }