Feature: enable ABACUS can finish SCF if charge density oscillation is found

[WHUweiqingzhou]

Fix #5406

Motivation

For systems that are difficult to converge, the SCF process may exhibit oscillations in charge density, preventing further progress toward the specified convergence criteria and resulting in continuous oscillation until the maximum number of steps is reached; this greatly wastes computational resources. To address this issue, I have designed a new feature that allows ABACUS to terminate the SCF process early upon detecting oscillations, thus reducing subsequent meaningless calculations. The detection of oscillations is based on the slope of the logarithm of historical drho values.

What's Changed?

1. add a new parameter scf_os_stop. If true, SCF will stop if charge density oscillation is found. The default is false
2. add a new parameter to control the slope threshold scf_os_thr. The default value is -0.01
3. add a new parameter scf_os_ndim to control the number of past iterations used in slope calculation. Generally speaking, if the number of oscillation steps persists beyond the steps considered for charge density mixing mixing_ndim, the SCF is unlikely to break free from oscillations. So, I set the default value as the same as mixing_ndim, you can change it.
4. add a new function Charge_Mixing::if_scf_oscillate to determine if oscillation happens.
5. add a new member oscillate_esolver, similar to conv_esolver.
6. add some docs and comments for users.

[WHUweiqingzhou]

If anyone is unsure how to choose scf_thr_os, run this script and check your slope of your drho calculated. I think -0.01 is a good choice.

#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <sstream>
#include <cmath>

bool if_scf_oscillate(std::vector<double>& drho_history, const int iteration, const double drho, const int iternum_used, const double threshold)
{

    if(threshold >= 0) // close the function
    {
        return false;
    }

    // add drho into history
    drho_history[iteration - 1] = drho;

    // check if the history is long enough
    if(iteration < iternum_used)
    {
        return false;
    }

    // calculate the slope of the last iternum_used iterations' drho
    double slope = 0.0;

    // Least Squares Method
    double sumX = 0, sumY = 0, sumXY = 0, sumXX = 0;
    for (int i = iteration - iternum_used; i < iteration; i++)
    {
        sumX += i;
        sumY += std::log10(drho_history[i]);
        sumXY += i * std::log10(drho_history[i]);
        sumXX += i * i;
    }
    double numerator = iternum_used * sumXY - sumX * sumY;
    double denominator = iternum_used * sumXX - sumX * sumX;
    if (denominator == 0) {
        return false;
    }
    slope =  numerator / denominator;

    std::cout << iteration << "-th slope: " << slope << std::endl;

    // if the slope is less than the threshold, return true
    if(slope > threshold)
    {
        return true;
    }

    return false; 
}

int main() {
    std::string filename = "data.txt"; // 替换为你的文件名
    std::ifstream file(filename);
    std::vector<double> numbers;
    std::string line;

    if (file.is_open()) {
        while (getline(file, line)) {
            std::istringstream iss(line);
            double number;
            if (!(iss >> number)) { // 如果转换失败，跳过这个数字
                continue;
            }
            numbers.push_back(number);
        }
        file.close();
    } else {
        std::cerr << "Unable to open file" << std::endl;
        return 1;
    }

    int scf_nmax = 10;
    int scf_os_ndim = 3;
    double scf_thr_os = -0.05;

    std::vector<double> drho_history;
    drho_history.resize(scf_nmax);

    bool oscillate = false;

    for (int i = 1; i <= numbers.size(); i++) {
        oscillate = if_scf_oscillate(drho_history, i, numbers[i-1], scf_os_ndim, scf_thr_os);
        if (oscillate) {
            std::cout << "Oscillation detected at iteration " << i << std::endl;
        }
    }

    return 0;
}