Refactor: add Mixing module temporarily; abacus_pw can be compiled now

source/Makefile.Objects

@@ -22,6 +22,7 @@ VPATH=./src_global:\
 ./module_base/module_container/base/core:\
 ./module_base/module_container/ATen/core:\
 ./module_base/module_container/ATen/kernels:\
+./module_base/module_mixing:\
 ./module_md:\
 ./module_basis/module_pw:\
 ./module_esolver:\
@@ -137,6 +138,8 @@ OBJS_BASE=abfs-vector3_order.o\
     formatter_contextfmt.o\
 	cubic_spline.o\
 	spherical_bessel_transformer.o\
+    mixing_data.o\
+
 
 OBJS_CELL=atom_pseudo.o\
     atom_spec.o\
@@ -472,10 +475,8 @@ OBJS_SRCPW=H_Ewald_pw.o\
     charge.o\
     charge_init.o\
     charge_mpi.o\
-    charge_broyden.o\
     charge_extra.o\
     charge_mixing.o\
-    charge_pulay.o\
     fp_energy.o\
     forces.o\
     force_op.o\

source/module_base/CMakeLists.txt

@@ -52,6 +52,7 @@ add_library(
     formatter_physfmt.cpp
     formatter_table.cpp
     formatter_contextfmt.cpp
+    module_mixing/mixing_data.cpp
     ${LIBM_SRC}
 )
 

source/module_base/lapack_connector.h

@@ -94,7 +94,7 @@ extern "C"
     // dsytrf_ computes the Bunch-Kaufman factorization of a double precision
     // symmetric matrix, while dsytri takes its output to perform martrix inversion
     void dsytrf_(const char* uplo, const int* n, double * a, const int* lda,
-                 int *ipiv,double *work, int* lwork ,int *info);
+                 int *ipiv,double *work, const int* lwork ,int *info);
     void dsytri_(const char* uplo,const int* n,double *a, const int *lda,
                  int *ipiv, double * work,int *info);
     // Peize Lin add dsptrf and dsptri 2016-06-21, to compute inverse real symmetry indefinit matrix.

source/module_base/module_mixing/broyden_mixing.h

@@ -0,0 +1,260 @@
+#ifndef BROYDEN_MIXING_H_
+#define BROYDEN_MIXING_H_
+#include "mixing.h"
+#include "module_base/lapack_connector.h"
+#include "module_base/matrix.h"
+#include "module_base/memory.h"
+#include "module_base/timer.h"
+#include "module_base/tool_title.h"
+
+namespace Base_Mixing
+{
+/**
+ * @brief Simplified modified broyden_mixing method.
+ *        Ref: D.D. Johnson PRB 38, 12807 (1988)
+ *        Here the weight w0 of the error of the inverse Jacobian is set to 0 and the weight wn of
+ *        the error of each previous iteration is set to same.
+ * @note  Formula:
+ *        F = n_out - n_in
+ *        dF{i} = F_{i-1} - F{i}  //different from Ref
+ *        dn_in{i} = n_in_{i-1} - n_in{i}  //different from Ref
+ *        alpha{ij} = <dF{i}, dF{j}>
+ *        beta{ij} = inv(alpha){ij}
+ *        c{mk} = <dF{k}, F{m}>
+ *        gamma{mn} = \sum_k c{mk} * beta{kn} 
+ *        n{m+1} = n_in{m} + mixing_beta*F{m} - \sum_n gamma{mn} * (dn_in{n} + mixing_beta*dF{n})
+ *        mixing_data{i} = n_in{i} + mixing_beta*F{i}
+ *        n{m+1} = \sum_i coef{i} * mixing_data{i} 
+ */
+class Broyden_Mixing : public Mixing
+{
+  public:
+    Broyden_Mixing(const int& mixing_ndim)
+    {
+        this->mixing_ndim = mixing_ndim;
+        this->coef = std::vector<double>(mixing_ndim);
+        this->beta = ModuleBase::matrix(mixing_ndim, mixing_ndim, true);
+    }
+    virtual ~Broyden_Mixing() override
+    {
+        if (F != nullptr)
+            free(F);
+        if (dF != nullptr)
+            free(dF);
+    };
+    virtual void push_data(Mixing_Data& mdata,
+                           const double* data_in,
+                           const double* data_out,
+                           std::function<void(double*)> screen,
+                           const bool& need_calcoef) override
+    {
+        this->tem_push_data(mdata, data_in, data_out, screen, need_calcoef);
+    };
+    virtual void push_data(Mixing_Data& mdata,
+                           const std::complex<double>* data_in,
+                           const std::complex<double>* data_out,
+                           std::function<void(std::complex<double>*)> screen,
+                           const bool& need_calcoef) override
+    {
+        this->tem_push_data(mdata, data_in, data_out, screen, need_calcoef);
+    };
+    virtual void cal_coef(const Mixing_Data& mdata, std::function<double(double*, double*)> inner_dot) override
+    {
+        tem_cal_coef(mdata, inner_dot);
+    }
+    virtual void cal_coef(const Mixing_Data& mdata,
+                          std::function<double(std::complex<double>*, std::complex<double>*)> inner_dot) override
+    {
+        tem_cal_coef(mdata, inner_dot);
+    }
+
+  private:
+    template <class FPTYPE>
+    void tem_push_data(Mixing_Data& mdata,
+                       const FPTYPE* data_in,
+                       const FPTYPE* data_out,
+                       std::function<void(FPTYPE*)> screen,
+                       const bool& need_calcoef)
+    {
+        const size_t length = mdata.length;
+        std::vector<FPTYPE> F_tmp(length);
+
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 128)
+#endif
+        for (int i = 0; i < length; ++i)
+        {
+            F_tmp[i] = data_out[i] - data_in[i];
+        }
+
+        // get screened F
+        if (screen != nullptr)
+            screen(F_tmp.data());
+
+        // container::Tensor data = data_in + mixing_beta * F;
+        std::vector<FPTYPE> data(length);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 128)
+#endif
+        for (int i = 0; i < length; ++i)
+        {
+            data[i] = data_in[i] + this->mixing_beta * F_tmp[i];
+        }
+
+        mdata.push(data.data());
+
+        if (!need_calcoef)
+            return;
+
+        if (address != &mdata && address != nullptr)
+            ModuleBase::WARNING_QUIT(
+                "Broyden_Mixing",
+                "One Broyden_Mixing object can only bind one Mixing_Data object to calculate coefficients");
+
+        FPTYPE* FP_dF = static_cast<FPTYPE*>(dF);
+        FPTYPE* FP_F = static_cast<FPTYPE*>(F);
+        if (mdata.ndim_use == 1)
+        {
+            address = &mdata;
+            // allocate
+            if (F != nullptr)
+                free(F);
+            F = malloc(sizeof(FPTYPE) * length);
+            FP_F = static_cast<FPTYPE*>(F);
+            if (dF != nullptr)
+                free(dF);
+            dF = malloc(sizeof(FPTYPE) * length * mixing_ndim);
+            FP_dF = static_cast<FPTYPE*>(dF);
+
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 128)
+#endif
+            for (int i = 0; i < length; ++i)
+            {
+                FP_F[i] = F_tmp[i];
+            }
+        }
+        else
+        {
+            const int previous = mdata.index_move(-1);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 128)
+#endif
+            for (int i = 0; i < length; ++i)
+            {
+                FP_F[i] = F_tmp[i];
+                // dF{n} = F{n-1} - F{n} = -(F{n} - F{n-1})
+                FP_dF[previous * length + i] -= FP_F[i];
+            }
+        }
+    };
+
+    template <class FPTYPE>
+    void tem_cal_coef(const Mixing_Data& mdata, std::function<double(FPTYPE*, FPTYPE*)> inner_dot)
+    {
+        ModuleBase::TITLE("Charge_Mixing", "Simplified_Broyden_mixing");
+        ModuleBase::timer::tick("Charge", "Broyden_mixing");
+        if (address != &mdata && address != nullptr)
+            ModuleBase::WARNING_QUIT(
+                "Mixing",
+                "One Mixing object can only bind one Mixing_Data object to calculate coefficients");
+        const int length = mdata.length;
+        const int start = mdata.start;
+        FPTYPE* FP_dF = static_cast<FPTYPE*>(dF);
+        FPTYPE* FP_F = static_cast<FPTYPE*>(F);
+        if (mdata.ndim_previous > 0)
+        {
+            const int ndim_previous = mdata.ndim_previous;
+            ModuleBase::matrix beta_tmp(ndim_previous, ndim_previous);
+            //beta(i, j) = <dF_i, dF_j>
+            for (int i = 0; i < ndim_previous; ++i)
+            {
+                FPTYPE* dFi = FP_dF + i * length;
+                for (int j = i; j < ndim_previous; ++j)
+                {
+                    if (i < ndim_previous - 1 && j < ndim_previous - 1)
+                    {
+                        beta_tmp(i, j) = beta(i, j);
+                    }
+                    FPTYPE* dFj = FP_dF + j * length;
+                    beta(i, j) = beta_tmp(i, j) = inner_dot(dFi, dFj);
+                    if (j != i)
+                    {
+                        beta(j, i) = beta_tmp(j, i) = beta_tmp(i, j);
+                    }
+                }
+            }
+            double* work = new double[ndim_previous];
+            int* iwork = new int[ndim_previous];
+            char uu = 'U';
+            int info;
+            dsytrf_(&uu, &ndim_previous, beta_tmp.c, &ndim_previous, iwork, work, &ndim_previous, &info);
+            if (info != 0)
+                ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when factorizing beta.");
+            dsytri_(&uu, &ndim_previous, beta_tmp.c, &ndim_previous, iwork, work, &info);
+            if (info != 0)
+                ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when DSYTRI beta.");
+            for (int i = 0; i < ndim_previous; ++i)
+            {
+                for (int j = i + 1; j < ndim_previous; ++j)
+                {
+                    beta_tmp(i, j) = beta_tmp(j, i);
+                }
+            }
+            for (int i = 0; i < ndim_previous; ++i)
+            {
+                FPTYPE* dFi = FP_dF + i * length;
+                work[i] = inner_dot(dFi, FP_F);
+            }
+            //gamma[i] = \sum_j beta_tmp(i,j) * work[j]
+            std::vector<double> gamma(ndim_previous);
+            container::BlasConnector::gemv('N',
+                                           ndim_previous,
+                                           ndim_previous,
+                                           1.0,
+                                           beta_tmp.c,
+                                           ndim_previous,
+                                           work,
+                                           1,
+                                           0.0,
+                                           gamma.data(),
+                                           1);
+
+            coef[start] = 1 + gamma[mdata.index_move(-1)];
+            for (int i = 1; i < ndim_previous - 1; ++i)
+            {
+                coef[mdata.index_move(-i)] = gamma[mdata.index_move(-i - 1)] - gamma[mdata.index_move(-i)];
+            }
+            coef[mdata.index_move(-ndim_previous)] =  -gamma[mdata.index_move(-ndim_previous)];
+
+            delete[] work;
+            delete[] iwork;
+        }
+        else
+        {
+            coef[0] = 1.0;
+        }
+
+        FPTYPE* dFstart = FP_dF + start * length;
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 128)
+#endif
+        for (int i = 0; i < length; ++i)
+        {
+            dFstart[i] = FP_F[i];
+        }
+        ModuleBase::timer::tick("Charge", "Broyden_mixing");
+    };
+
+  private:
+    // F = data_out - data_in
+    void* F = nullptr;
+    // dF = F_{n+1} - F_n
+    void* dF = nullptr;
+    // binded mixing_data
+    Mixing_Data* address = nullptr;
+    // beta_ij = <dF_i, dF_j>
+    ModuleBase::matrix beta;
+};
+} // namespace Base_Mixing
+#endif