From 722d3bbf1e062839749f15c720fc4efdc8782558 Mon Sep 17 00:00:00 2001 From: Qianruipku Date: Thu, 23 Nov 2023 13:24:43 +0800 Subject: [PATCH 1/7] Refactor: charge_mixing Before: only one mixing_beta can be used for this mixing After: We can define any form of mixing_beta function. --- source/Makefile.Objects | 4 + source/module_base/CMakeLists.txt | 4 + .../module_mixing/broyden_mixing.cpp | 200 ++++++++++++ .../module_mixing/broyden_mixing.h | 234 ++------------ source/module_base/module_mixing/mixing.cpp | 113 +++++++ source/module_base/module_mixing/mixing.h | 91 ++---- .../module_base/module_mixing/mixing_data.cpp | 21 +- .../module_mixing/plain_mixing.cpp | 99 ++++++ .../module_base/module_mixing/plain_mixing.h | 130 ++++---- .../module_mixing/pulay_mixing.cpp | 182 +++++++++++ .../module_base/module_mixing/pulay_mixing.h | 219 ++----------- .../module_mixing/test/mixing_test.cpp | 88 +++-- .../module_elecstate/module_charge/charge.h | 2 +- .../module_charge/charge_mixing.cpp | 303 +++++++++++------- .../module_charge/charge_mixing.h | 50 ++- .../test/charge_mixing_test.cpp | 169 ++++++++-- source/module_esolver/esolver_ks.cpp | 35 +- .../module_io/test/for_testing_input_conv.h | 9 +- 18 files changed, 1179 insertions(+), 774 deletions(-) create mode 100644 source/module_base/module_mixing/broyden_mixing.cpp create mode 100644 source/module_base/module_mixing/mixing.cpp create mode 100644 source/module_base/module_mixing/plain_mixing.cpp create mode 100644 source/module_base/module_mixing/pulay_mixing.cpp diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 84b4f6fe17..67ada634d5 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -147,6 +147,10 @@ OBJS_BASE=abfs-vector3_order.o\ cubic_spline.o\ spherical_bessel_transformer.o\ mixing_data.o\ + mixing.o\ + plain_mixing.o\ + pulay_mixing.o\ + broyden_mixing.o\ OBJS_CELL=atom_pseudo.o\ diff --git a/source/module_base/CMakeLists.txt b/source/module_base/CMakeLists.txt index 5fb8bd0dd5..b6d362d5ec 100644 --- a/source/module_base/CMakeLists.txt +++ b/source/module_base/CMakeLists.txt @@ -54,6 +54,10 @@ add_library( formatter_table.cpp formatter_contextfmt.cpp module_mixing/mixing_data.cpp + module_mixing/mixing.cpp + module_mixing/plain_mixing.cpp + module_mixing/pulay_mixing.cpp + module_mixing/broyden_mixing.cpp ${LIBM_SRC} ) diff --git a/source/module_base/module_mixing/broyden_mixing.cpp b/source/module_base/module_mixing/broyden_mixing.cpp new file mode 100644 index 0000000000..95caba2b5b --- /dev/null +++ b/source/module_base/module_mixing/broyden_mixing.cpp @@ -0,0 +1,200 @@ +#include "broyden_mixing.h" + +#include "module_base/lapack_connector.h" +#include "module_base/memory.h" +#include "module_base/module_container/base/third_party/blas.h" +#include "module_base/timer.h" +#include "module_base/tool_title.h" +namespace Base_Mixing +{ +template void Broyden_Mixing::tem_push_data(Mixing_Data& mdata, + const double* data_in, + const double* data_out, + std::function screen, + std::function mix, + const bool& need_calcoef); +template void Broyden_Mixing::tem_push_data( + Mixing_Data& mdata, + const std::complex* data_in, + const std::complex* data_out, + std::function*)> screen, + std::function*, const std::complex*, const std::complex*)> mix, + const bool& need_calcoef); + +template +void Broyden_Mixing::tem_push_data(Mixing_Data& mdata, + const FPTYPE* data_in, + const FPTYPE* data_out, + std::function screen, + std::function mix, + const bool& need_calcoef) +{ + const size_t length = mdata.length; + std::vector F_tmp(length); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) +#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 data(length); + mix(data.data(), data_in, F_tmp.data()); + + 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(dF); + FPTYPE* FP_F = static_cast(F); + if (mdata.ndim_use == 1) + { + address = &mdata; + // allocate + if (F != nullptr) + free(F); + F = malloc(sizeof(FPTYPE) * length); + FP_F = static_cast(F); + if (dF != nullptr) + free(dF); + dF = malloc(sizeof(FPTYPE) * length * mixing_ndim); + FP_dF = static_cast(dF); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) +#endif + for (int i = 0; i < length; ++i) + { + FP_F[i] = F_tmp[i]; + } + } + else + { + this->ndim_cal_dF = std::min(this->ndim_cal_dF + 1, this->mixing_ndim); + start_dF = (this->start_dF + 1) % this->mixing_ndim; +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) +#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[start_dF * length + i] -= FP_F[i]; + } + } +}; + +template void Broyden_Mixing::tem_cal_coef(const Mixing_Data& mdata, + std::function inner_product); +template void Broyden_Mixing::tem_cal_coef( + const Mixing_Data& mdata, + std::function*, std::complex*)> inner_product); + +template +void Broyden_Mixing::tem_cal_coef(const Mixing_Data& mdata, std::function inner_product) +{ + ModuleBase::TITLE("Charge_Mixing", "Simplified_Broyden_mixing"); + ModuleBase::timer::tick("Charge", "Broyden_mixing"); + if (address != &mdata && address != nullptr) + ModuleBase::WARNING_QUIT( + "Broyden_mixing", + "One Broyden_Mixing object can only bind one Mixing_Data object to calculate coefficients"); + const int length = mdata.length; + FPTYPE* FP_dF = static_cast(dF); + FPTYPE* FP_F = static_cast(F); + if (ndim_cal_dF > 0) + { + ModuleBase::matrix beta_tmp(ndim_cal_dF, ndim_cal_dF); + // beta(i, j) = + for (int i = 0; i < ndim_cal_dF; ++i) + { + FPTYPE* dFi = FP_dF + i * length; + for (int j = i; j < ndim_cal_dF; ++j) + { + if (i != start_dF && j != start_dF) + { + beta_tmp(i, j) = beta(i, j); + } + else + { + FPTYPE* dFj = FP_dF + j * length; + beta(i, j) = beta_tmp(i, j) = inner_product(dFi, dFj); + } + if (j != i) + { + beta_tmp(j, i) = beta_tmp(i, j); + } + } + } + double* work = new double[ndim_cal_dF]; + int* iwork = new int[ndim_cal_dF]; + char uu = 'U'; + int info; + dsytrf_(&uu, &ndim_cal_dF, beta_tmp.c, &ndim_cal_dF, iwork, work, &ndim_cal_dF, &info); + if (info != 0) + ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when factorizing beta."); + dsytri_(&uu, &ndim_cal_dF, beta_tmp.c, &ndim_cal_dF, iwork, work, &info); + if (info != 0) + ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when DSYTRI beta."); + for (int i = 0; i < ndim_cal_dF; ++i) + { + for (int j = i + 1; j < ndim_cal_dF; ++j) + { + beta_tmp(i, j) = beta_tmp(j, i); + } + } + for (int i = 0; i < ndim_cal_dF; ++i) + { + FPTYPE* dFi = FP_dF + i * length; + work[i] = inner_product(dFi, FP_F); + } + // gamma[i] = \sum_j beta_tmp(i,j) * work[j] + std::vector gamma(ndim_cal_dF); + container::BlasConnector::gemv('N', + ndim_cal_dF, + ndim_cal_dF, + 1.0, + beta_tmp.c, + ndim_cal_dF, + work, + 1, + 0.0, + gamma.data(), + 1); + coef[mdata.start] = 1 + gamma[dFindex_move(0)]; + for (int i = 1; i < ndim_cal_dF; ++i) + { + coef[mdata.index_move(-i)] = gamma[dFindex_move(-i)] - gamma[dFindex_move(-i + 1)]; + } + coef[mdata.index_move(-ndim_cal_dF)] = -gamma[dFindex_move(-ndim_cal_dF + 1)]; + + delete[] work; + delete[] iwork; + } + else + { + coef[0] = 1.0; + } + + FPTYPE* dFnext = FP_dF + dFindex_move(1) * length; +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) +#endif + for (int i = 0; i < length; ++i) + { + dFnext[i] = FP_F[i]; + } + ModuleBase::timer::tick("Charge", "Broyden_mixing"); +}; +} // namespace Base_Mixing \ No newline at end of file diff --git a/source/module_base/module_mixing/broyden_mixing.h b/source/module_base/module_mixing/broyden_mixing.h index 1b3ee66ff9..d3fbbd82f3 100755 --- a/source/module_base/module_mixing/broyden_mixing.h +++ b/source/module_base/module_mixing/broyden_mixing.h @@ -1,12 +1,7 @@ #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" -#include "module_base/global_variable.h" namespace Base_Mixing { @@ -30,21 +25,16 @@ namespace Base_Mixing class Broyden_Mixing : public Mixing { public: - Broyden_Mixing(const int& mixing_ndim, const double& mixing_beta) + Broyden_Mixing(const int& mixing_ndim) { this->mixing_ndim = mixing_ndim; this->data_ndim = mixing_ndim + 1; - this->mixing_beta = mixing_beta; this->coef = std::vector(mixing_ndim + 1); this->beta = ModuleBase::matrix(mixing_ndim, mixing_ndim, true); - if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4) - { - this->two_beta = 0; - } - else if (GlobalV::NSPIN == 2) - { - this->two_beta = 1; - } + } + Broyden_Mixing(const int& mixing_ndim, const double& mixing_beta) : Broyden_Mixing(mixing_ndim) + { + this->mixing_beta = mixing_beta; } virtual ~Broyden_Mixing() override { @@ -71,6 +61,7 @@ class Broyden_Mixing : public Mixing * @param data_in x_in * @param data_out x_out = f(x_in) * @param screen pointer to the screen function for Ker-Ker mixing + * @param mix (double* out, const double* in, const double* residual) calculate 'out' with 'in' and residual * @param need_calcoef whether need to calculate the coef * */ @@ -78,17 +69,20 @@ class Broyden_Mixing : public Mixing const double* data_in, const double* data_out, std::function screen, + std::function mix, const bool& need_calcoef) override { - this->tem_push_data(mdata, data_in, data_out, screen, need_calcoef); + this->tem_push_data(mdata, data_in, data_out, screen, mix, need_calcoef); }; - virtual void push_data(Mixing_Data& mdata, - const std::complex* data_in, - const std::complex* data_out, - std::function*)> screen, - const bool& need_calcoef) override + virtual void push_data( + Mixing_Data& mdata, + const std::complex* data_in, + const std::complex* data_out, + std::function*)> screen, + std::function*, const std::complex*, const std::complex*)> mix, + const bool& need_calcoef) override { - this->tem_push_data(mdata, data_in, data_out, screen, need_calcoef); + this->tem_push_data(mdata, data_in, data_out, screen, mix, need_calcoef); }; /** @@ -115,6 +109,7 @@ class Broyden_Mixing : public Mixing * @param data_in x_in * @param data_out x_out = f(x_in) * @param screen pointer to the screen function for Ker-Ker mixing + * @param mix (double* out, const double* in, const double* residual) calculate 'out' with 'in' and residual * @param need_calcoef whether need to calculate the coef * */ @@ -123,101 +118,8 @@ class Broyden_Mixing : public Mixing const FPTYPE* data_in, const FPTYPE* data_out, std::function screen, - const bool& need_calcoef) - { - const size_t length = mdata.length; - std::vector F_tmp(length); -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) -#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 data(length); - // mix density and magnetic density sperately - if (this->two_beta == 0) - { - // rho_tot -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) -#endif - for (int i = 0; i < length; ++i) - { - data[i] = data_in[i] + this->mixing_beta * F_tmp[i]; - } - } - else if (this->two_beta == 1) - { - // rho_tot -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) -#endif - for (int i = 0; i < length / 2; ++i) - { - data[i] = data_in[i] + this->mixing_beta * F_tmp[i]; - } - // magnetism -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) -#endif - for (int i = length / 2; i < length; ++i) - { - data[i] = data_in[i] + GlobalV::MIXING_BETA_MAG * 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(dF); - FPTYPE* FP_F = static_cast(F); - if (mdata.ndim_use == 1) - { - address = &mdata; - // allocate - if (F != nullptr) - free(F); - F = malloc(sizeof(FPTYPE) * length); - FP_F = static_cast(F); - if (dF != nullptr) - free(dF); - dF = malloc(sizeof(FPTYPE) * length * mixing_ndim); - FP_dF = static_cast(dF); -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) -#endif - for (int i = 0; i < length; ++i) - { - FP_F[i] = F_tmp[i]; - } - } - else - { - this->ndim_cal_dF = std::min(this->ndim_cal_dF + 1, this->mixing_ndim); - start_dF = (this->start_dF + 1) % this->mixing_ndim; -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) -#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[start_dF * length + i] -= FP_F[i]; - } - } - }; + std::function mix, + const bool& need_calcoef); /** * @brief calculate coeficients for mixing @@ -226,101 +128,7 @@ class Broyden_Mixing : public Mixing * @param inner_product pointer to the inner dot function */ template - void tem_cal_coef(const Mixing_Data& mdata, std::function inner_product) - { - ModuleBase::TITLE("Charge_Mixing", "Simplified_Broyden_mixing"); - ModuleBase::timer::tick("Charge", "Broyden_mixing"); - if (address != &mdata && address != nullptr) - ModuleBase::WARNING_QUIT( - "Broyden_mixing", - "One Broyden_Mixing object can only bind one Mixing_Data object to calculate coefficients"); - const int length = mdata.length; - FPTYPE* FP_dF = static_cast(dF); - FPTYPE* FP_F = static_cast(F); - if (ndim_cal_dF > 0) - { - ModuleBase::matrix beta_tmp(ndim_cal_dF, ndim_cal_dF); - // beta(i, j) = - for (int i = 0; i < ndim_cal_dF; ++i) - { - FPTYPE* dFi = FP_dF + i * length; - for (int j = i; j < ndim_cal_dF; ++j) - { - if (i != start_dF && j != start_dF) - { - beta_tmp(i, j) = beta(i, j); - } - else - { - FPTYPE* dFj = FP_dF + j * length; - beta(i, j) = beta_tmp(i, j) = inner_product(dFi, dFj); - } - if (j != i) - { - beta_tmp(j, i) = beta_tmp(i, j); - } - } - } - double* work = new double[ndim_cal_dF]; - int* iwork = new int[ndim_cal_dF]; - char uu = 'U'; - int info; - dsytrf_(&uu, &ndim_cal_dF, beta_tmp.c, &ndim_cal_dF, iwork, work, &ndim_cal_dF, &info); - if (info != 0) - ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when factorizing beta."); - dsytri_(&uu, &ndim_cal_dF, beta_tmp.c, &ndim_cal_dF, iwork, work, &info); - if (info != 0) - ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when DSYTRI beta."); - for (int i = 0; i < ndim_cal_dF; ++i) - { - for (int j = i + 1; j < ndim_cal_dF; ++j) - { - beta_tmp(i, j) = beta_tmp(j, i); - } - } - for (int i = 0; i < ndim_cal_dF; ++i) - { - FPTYPE* dFi = FP_dF + i * length; - work[i] = inner_product(dFi, FP_F); - } - // gamma[i] = \sum_j beta_tmp(i,j) * work[j] - std::vector gamma(ndim_cal_dF); - container::BlasConnector::gemv('N', - ndim_cal_dF, - ndim_cal_dF, - 1.0, - beta_tmp.c, - ndim_cal_dF, - work, - 1, - 0.0, - gamma.data(), - 1); - coef[mdata.start] = 1 + gamma[dFindex_move(0)]; - for (int i = 1; i < ndim_cal_dF; ++i) - { - coef[mdata.index_move(-i)] = gamma[dFindex_move(-i)] - gamma[dFindex_move(-i + 1)]; - } - coef[mdata.index_move(-ndim_cal_dF)] = -gamma[dFindex_move(-ndim_cal_dF + 1)]; - - delete[] work; - delete[] iwork; - } - else - { - coef[0] = 1.0; - } - - FPTYPE* dFnext = FP_dF + dFindex_move(1) * length; -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) -#endif - for (int i = 0; i < length; ++i) - { - dFnext[i] = FP_F[i]; - } - ModuleBase::timer::tick("Charge", "Broyden_mixing"); - }; + void tem_cal_coef(const Mixing_Data& mdata, std::function inner_product); private: // F = data_out - data_in @@ -343,8 +151,6 @@ class Broyden_Mixing : public Mixing } // the number of calculated dF int ndim_cal_dF = 0; - // if we should considere two beta - int two_beta = 0; }; } // namespace Base_Mixing #endif \ No newline at end of file diff --git a/source/module_base/module_mixing/mixing.cpp b/source/module_base/module_mixing/mixing.cpp new file mode 100644 index 0000000000..bd522bc60b --- /dev/null +++ b/source/module_base/module_mixing/mixing.cpp @@ -0,0 +1,113 @@ +#include "mixing.h" + +#include "module_base/module_container/base/third_party/blas.h" +namespace Base_Mixing +{ + +void Mixing::push_data(Mixing_Data& mdata, + const double* data_in, + const double* data_out, + std::function screen, + const bool& need_calcoef) +{ + const size_t length = mdata.length; + this->push_data( + mdata, + data_in, + data_out, + screen, + [this, length](double* out, const double* in, const double* sres) { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 512) +#endif + for (int i = 0; i < length; ++i) + { + out[i] = in[i] + this->mixing_beta * sres[i]; + } + }, + need_calcoef); + return; +} + +void Mixing::push_data(Mixing_Data& mdata, + const std::complex* data_in, + const std::complex* data_out, + std::function*)> screen, + const bool& need_calcoef) +{ + const size_t length = mdata.length; + this->push_data( + mdata, + data_in, + data_out, + screen, + [this, length](std::complex* out, const std::complex* in, const std::complex* sres) { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 256) +#endif + for (int i = 0; i < length; ++i) + { + out[i] = in[i] + this->mixing_beta * sres[i]; + } + }, + need_calcoef); + return; +} + +void Mixing::mix_data(const Mixing_Data& mdata, double* data_mix) +{ + if (mdata.length <= 0) + return; + double* FP_data = static_cast(mdata.data); + if (mdata.ndim_use == 1) + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 512) +#endif + for (int i = 0; i < mdata.length; ++i) + data_mix[i] = FP_data[i]; + return; + } + container::BlasConnector::gemv('N', + mdata.length, + mdata.ndim_use, + 1.0, + FP_data, + mdata.length, + coef.data(), + 1, + 0.0, + data_mix, + 1); +} +void Mixing::mix_data(const Mixing_Data& mdata, std::complex* data_mix) +{ + if (mdata.length <= 0) + return; + std::complex* FP_data = static_cast*>(mdata.data); + if (mdata.ndim_use == 1) + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 256) +#endif + for (int i = 0; i < mdata.length; ++i) + data_mix[i] = FP_data[i]; + return; + } + // conver coef to complex + std::vector> coef_complex(coef.size()); + for (int i = 0; i < coef.size(); ++i) + coef_complex[i] = coef[i]; + container::BlasConnector::gemv('N', + mdata.length, + mdata.ndim_use, + 1.0, + FP_data, + mdata.length, + coef_complex.data(), + 1, + 0.0, + data_mix, + 1); +} +} // namespace Base_Mixing \ No newline at end of file diff --git a/source/module_base/module_mixing/mixing.h b/source/module_base/module_mixing/mixing.h index ff2e9d76eb..ba0226ee79 100644 --- a/source/module_base/module_mixing/mixing.h +++ b/source/module_base/module_mixing/mixing.h @@ -3,7 +3,7 @@ #include #include "mixing_data.h" -#include "module_base/module_container/base/third_party/blas.h" + namespace Base_Mixing { @@ -41,6 +41,7 @@ class Mixing * @param data_in x_in * @param data_out x_out = f(x_in) * @param screen pointer to the screen function for Ker-Ker mixing + * @param mix (double* out, const double* in, const double* residual) calculate 'out' with 'in' and residual * @param need_calcoef whether need to calculate the coef * */ @@ -48,15 +49,39 @@ class Mixing const double* data_in, const double* data_out, std::function screen, + std::function mix, const bool& need_calcoef) = 0; + virtual void push_data( + Mixing_Data& mdata, + const std::complex* data_in, + const std::complex* data_out, + std::function*)> screen, + std::function*, const std::complex*, const std::complex*)> mix, + const bool& need_calcoef) + = 0; + + /** + * @brief + * + * @param mdata store information of this iterative step + * @param data_in x_in + * @param data_out x_out = f(x_in) + * @param screen pointer to the screen function for Ker-Ker mixing + * @param need_calcoef whether need to calculate the coef + * + */ + virtual void push_data(Mixing_Data& mdata, + const double* data_in, + const double* data_out, + std::function screen, + const bool& need_calcoef); virtual void push_data(Mixing_Data& mdata, const std::complex* data_in, const std::complex* data_out, std::function*)> screen, - const bool& need_calcoef) - = 0; - + const bool& need_calcoef); + /** * @brief calculate coeficients for mixing * @@ -74,62 +99,8 @@ class Mixing * @param mdata Mixing_Data * @param data_mix output data */ - void mix_data(const Mixing_Data& mdata, double* data_mix) - { - if (mdata.length <= 0) - return; - double* FP_data = static_cast(mdata.data); - if (mdata.ndim_use == 1) - { -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 512) -#endif - for (int i = 0; i < mdata.length; ++i) - data_mix[i] = FP_data[i]; - return; - } - container::BlasConnector::gemv('N', - mdata.length, - mdata.ndim_use, - 1.0, - FP_data, - mdata.length, - coef.data(), - 1, - 0.0, - data_mix, - 1); - } - void mix_data(const Mixing_Data& mdata, std::complex* data_mix) - { - if (mdata.length <= 0) - return; - std::complex* FP_data = static_cast*>(mdata.data); - if (mdata.ndim_use == 1) - { -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 256) -#endif - for (int i = 0; i < mdata.length; ++i) - data_mix[i] = FP_data[i]; - return; - } - // conver coef to complex - std::vector> coef_complex(coef.size()); - for (int i = 0; i < coef.size(); ++i) - coef_complex[i] = coef[i]; - container::BlasConnector::gemv('N', - mdata.length, - mdata.ndim_use, - 1.0, - FP_data, - mdata.length, - coef_complex.data(), - 1, - 0.0, - data_mix, - 1); - } + void mix_data(const Mixing_Data& mdata, double* data_mix); + void mix_data(const Mixing_Data& mdata, std::complex* data_mix); /** * @brief reset mixing diff --git a/source/module_base/module_mixing/mixing_data.cpp b/source/module_base/module_mixing/mixing_data.cpp index a257f4e788..6c41757ee3 100644 --- a/source/module_base/module_mixing/mixing_data.cpp +++ b/source/module_base/module_mixing/mixing_data.cpp @@ -3,30 +3,35 @@ namespace Base_Mixing { -Mixing_Data::Mixing_Data(const int& ndim, const int& length, const size_t &type_size) +Mixing_Data::Mixing_Data(const int& ndim, const int& length, const size_t& type_size) { this->ndim_tot = ndim; this->length = length; - this->data = malloc(ndim * length * type_size); + if (ndim * length > 0) + { + this->data = malloc(ndim * length * type_size); + } } Mixing_Data::~Mixing_Data() { - if(this->data != nullptr) + if (this->data != nullptr) free(this->data); } -void Mixing_Data::resize(const int& ndim, const int& length, const size_t &type_size) +void Mixing_Data::resize(const int& ndim, const int& length, const size_t& type_size) { this->ndim_tot = ndim; this->length = length; - if(this->data != nullptr) + if (this->data != nullptr) free(this->data); - this->data = malloc(ndim * length * type_size); + if (ndim * length > 0) + { + this->data = malloc(ndim * length * type_size); + } this->start = -1; this->ndim_use = 0; this->ndim_history = 0; } - -} +} // namespace Base_Mixing diff --git a/source/module_base/module_mixing/plain_mixing.cpp b/source/module_base/module_mixing/plain_mixing.cpp new file mode 100644 index 0000000000..025cc67b1e --- /dev/null +++ b/source/module_base/module_mixing/plain_mixing.cpp @@ -0,0 +1,99 @@ +#include "plain_mixing.h" + +#include "module_base/memory.h" +#include "module_base/timer.h" +#include "module_base/tool_title.h" +namespace Base_Mixing +{ +template void Plain_Mixing::tem_push_data(Mixing_Data& mdata, + const double* data_in, + const double* data_out, + std::function screen, + std::function mix, + const bool& need_calcoef); +template void Plain_Mixing::tem_push_data( + Mixing_Data& mdata, + const std::complex* data_in, + const std::complex* data_out, + std::function*)> screen, + std::function*, const std::complex*, const std::complex*)> mix, + const bool& need_calcoef); + +template +void Plain_Mixing::tem_push_data(Mixing_Data& mdata, + const FPTYPE* data_in, + const FPTYPE* data_out, + std::function screen, + std::function mix, + const bool& need_calcoef) +{ + const size_t length = mdata.length; + std::vector F_tmp(length); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) +#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 data(length); + mix(data.data(), data_in, F_tmp.data()); + + mdata.push(data.data()); +}; + +template void Plain_Mixing::simple_mix(double* data_new, + const double* data_in, + const double* data_out, + const int& length, + std::function screen); +template void Plain_Mixing::simple_mix(std::complex* data_new, + const std::complex* data_in, + const std::complex* data_out, + const int& length, + std::function*)> screen); +template +void Plain_Mixing::simple_mix(FPTYPE* data_new, + const FPTYPE* data_in, + const FPTYPE* data_out, + const int& length, + std::function screen) +{ + if (screen == nullptr) + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) +#endif + for (int ig = 0; ig < length; ig++) + { + data_new[ig] = data_in[ig] + this->mixing_beta * (data_out[ig] - data_in[ig]); + } + } + else + { + std::vector F_tmp(length); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) +#endif + for (int i = 0; i < length; ++i) + { + F_tmp[i] = data_out[i] - data_in[i]; + } + screen(F_tmp.data()); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) +#endif + for (int i = 0; i < length; ++i) + { + data_new[i] = data_in[i] + this->mixing_beta * F_tmp[i]; + } + } +} + +} // namespace Base_Mixing \ No newline at end of file diff --git a/source/module_base/module_mixing/plain_mixing.h b/source/module_base/module_mixing/plain_mixing.h index c87e2ee4a2..e7e3c87def 100644 --- a/source/module_base/module_mixing/plain_mixing.h +++ b/source/module_base/module_mixing/plain_mixing.h @@ -1,10 +1,6 @@ #ifndef PLAIN_MIXING_H_ #define PLAIN_MIXING_H_ #include "mixing.h" -#include "module_base/memory.h" -#include "module_base/timer.h" -#include "module_base/tool_title.h" -#include "module_base/global_variable.h" namespace Base_Mixing { @@ -15,19 +11,13 @@ namespace Base_Mixing class Plain_Mixing : public Mixing { public: - Plain_Mixing(const double& mixing_beta) + Plain_Mixing() { - this->mixing_beta = mixing_beta; - this->data_ndim = 1; this->coef = std::vector(1, 1.0); - if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4) - { - this->two_beta = 0; - } - else if (GlobalV::NSPIN == 2) - { - this->two_beta = 1; - } + } + Plain_Mixing(const double& mixing_beta) : Plain_Mixing() + { + this->mixing_beta = mixing_beta; } virtual ~Plain_Mixing() override{}; @@ -47,6 +37,7 @@ class Plain_Mixing : public Mixing * @param data_in x_in * @param data_out x_out = f(x_in) * @param screen pointer to the screen function for Ker-Ker mixing + * @param mix (double* out, const double* in, const double* residual) calculate 'out' with 'in' and residual * @param need_calcoef whether need to calculate the coef * */ @@ -54,17 +45,20 @@ class Plain_Mixing : public Mixing const double* data_in, const double* data_out, std::function screen, + std::function mix, const bool& need_calcoef) override { - this->tem_push_data(mdata, data_in, data_out, screen, need_calcoef); + this->tem_push_data(mdata, data_in, data_out, screen, mix, need_calcoef); }; - virtual void push_data(Mixing_Data& mdata, - const std::complex* data_in, - const std::complex* data_out, - std::function*)> screen, - const bool& need_calcoef) override + virtual void push_data( + Mixing_Data& mdata, + const std::complex* data_in, + const std::complex* data_out, + std::function*)> screen, + std::function*, const std::complex*, const std::complex*)> mix, + const bool& need_calcoef) override { - this->tem_push_data(mdata, data_in, data_out, screen, need_calcoef); + this->tem_push_data(mdata, data_in, data_out, screen, mix, need_calcoef); }; /** @@ -83,6 +77,29 @@ class Plain_Mixing : public Mixing return; } + /** + * @brief Simple plain mixing + * data_new = data_in + mixing_beta * (data_out - data_in) + * + * @param data_new can be the same as data_in or data_out + */ + void plain_mix(double* data_new, + const double* data_in, + const double* data_out, + const int& length, + std::function screen) + { + this->simple_mix(data_new, data_in, data_out, length, screen); + } + void plain_mix(std::complex* data_new, + const std::complex* data_in, + const std::complex* data_out, + const int& length, + std::function*)> screen) + { + this->simple_mix(data_new, data_in, data_out, length, screen); + } + private: /** * @brief @@ -91,6 +108,7 @@ class Plain_Mixing : public Mixing * @param data_in x_in * @param data_out x_out = f(x_in) * @param screen pointer to the screen function for Ker-Ker mixing + * @param mix (double* out, const double* in, const double* residual) calculate 'out' with 'in' and residual * @param need_calcoef whether need to calculate the coef * */ @@ -99,61 +117,21 @@ class Plain_Mixing : public Mixing const FPTYPE* data_in, const FPTYPE* data_out, std::function screen, - const bool& need_calcoef) - { - const size_t length = mdata.length; - std::vector F_tmp(length); -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) -#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 data(length); - if (this->two_beta == 0) - { - // rho_tot -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) -#endif - for (int i = 0; i < length; ++i) - { - data[i] = data_in[i] + this->mixing_beta * F_tmp[i]; - } - } - else if (this->two_beta == 1) - { - // rho_tot -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) -#endif - for (int i = 0; i < length / 2; ++i) - { - data[i] = data_in[i] + this->mixing_beta * F_tmp[i]; - } - // magnetism -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) -#endif - for (int i = length / 2; i < length; ++i) - { - data[i] = data_in[i] + GlobalV::MIXING_BETA_MAG * F_tmp[i]; - } - - } + std::function mix, + const bool& need_calcoef); - mdata.push(data.data()); - }; - - // if we should considere two beta - int two_beta = 0; + /** + * @brief Simple plain mixing + * data_new = data_in + mixing_beta * (data_out - data_in) + * + * @param data_new can be the same as data_in or data_out + */ + template + void simple_mix(FPTYPE* data_new, + const FPTYPE* data_in, + const FPTYPE* data_out, + const int& length, + std::function screen); }; } // namespace Base_Mixing #endif \ No newline at end of file diff --git a/source/module_base/module_mixing/pulay_mixing.cpp b/source/module_base/module_mixing/pulay_mixing.cpp new file mode 100644 index 0000000000..0c85e8e071 --- /dev/null +++ b/source/module_base/module_mixing/pulay_mixing.cpp @@ -0,0 +1,182 @@ +#include "pulay_mixing.h" + +#include "module_base/lapack_connector.h" +#include "module_base/memory.h" +#include "module_base/timer.h" +#include "module_base/tool_title.h" +namespace Base_Mixing +{ +template void Pulay_Mixing::tem_push_data(Mixing_Data& mdata, + const double* data_in, + const double* data_out, + std::function screen, + std::function mix, + const bool& need_calcoef); +template void Pulay_Mixing::tem_push_data( + Mixing_Data& mdata, + const std::complex* data_in, + const std::complex* data_out, + std::function*)> screen, + std::function*, const std::complex*, const std::complex*)> mix, + const bool& need_calcoef); + +template +void Pulay_Mixing::tem_push_data(Mixing_Data& mdata, + const FPTYPE* data_in, + const FPTYPE* data_out, + std::function screen, + std::function mix, + const bool& need_calcoef) +{ + const size_t length = mdata.length; + std::vector F_tmp(length); + +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) +#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 data(length); + mix(data.data(), data_in, F_tmp.data()); + + mdata.push(data.data()); + + if (!need_calcoef) + return; + + if (address != &mdata && address != nullptr) + ModuleBase::WARNING_QUIT( + "Pulay_Mixing", + "One Pulay_Mixing object can only bind one Mixing_Data object to calculate coefficients"); + + FPTYPE* FP_F = static_cast(F); + if (mdata.ndim_use == 1) + { + address = &mdata; + // allocate + if (F != nullptr) + free(F); + F = malloc(sizeof(FPTYPE) * length * mixing_ndim); + FP_F = static_cast(F); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) +#endif + for (int i = 0; i < length; ++i) + { + FP_F[i] = F_tmp[i]; + } + } + else + { + start_F = (this->start_F + 1) % this->mixing_ndim; + FPTYPE* FP_startF = FP_F + start_F * length; +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) +#endif + for (int i = 0; i < length; ++i) + { + FP_startF[i] = F_tmp[i]; + } + } +}; + +template void Pulay_Mixing::tem_cal_coef(const Mixing_Data& mdata, + std::function inner_product); +template void Pulay_Mixing::tem_cal_coef( + const Mixing_Data& mdata, + std::function*, std::complex*)> inner_product); + +template +void Pulay_Mixing::tem_cal_coef(const Mixing_Data& mdata, std::function inner_product) +{ + ModuleBase::TITLE("Charge_Mixing", "Pulay_mixing"); + ModuleBase::timer::tick("Charge", "Pulay_mixing"); + if (address != &mdata && address != nullptr) + ModuleBase::WARNING_QUIT( + "Pulay_mixing", + "One Pulay_Mixing object can only bind one Mixing_Data object to calculate coefficients"); + const int length = mdata.length; + FPTYPE* FP_F = static_cast(F); + + if (mdata.ndim_use > 1) + { + const int ndim_use = mdata.ndim_use; + ModuleBase::matrix beta_tmp(ndim_use, ndim_use); + // beta(i, j) = + for (int i = 0; i < ndim_use; ++i) + { + FPTYPE* Fi = FP_F + i * length; + for (int j = i; j < ndim_use; ++j) + { + if (i != start_F && j != start_F) + { + beta_tmp(i, j) = beta(i, j); + } + else + { + FPTYPE* Fj = FP_F + j * length; + beta(i, j) = beta_tmp(i, j) = inner_product(Fi, Fj); + } + if (j != i) + { + beta(j, i) = beta_tmp(j, i) = beta_tmp(i, j); + } + } + } + + double* work = new double[ndim_use]; + int* iwork = new int[ndim_use]; + char uu = 'U'; + int info; + dsytrf_(&uu, &ndim_use, beta_tmp.c, &ndim_use, iwork, work, &ndim_use, &info); + if (info != 0) + ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when factorizing beta."); + dsytri_(&uu, &ndim_use, beta_tmp.c, &ndim_use, iwork, work, &info); + if (info != 0) + ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when DSYTRI beta."); + for (int i = 0; i < ndim_use; ++i) + { + for (int j = i + 1; j < ndim_use; ++j) + { + beta_tmp(i, j) = beta_tmp(j, i); + } + } + + // coef{i} = \sum_j beta{ij} / \sum_k \sum_j beta{kj} + double sum_beta = 0.; + for (int i = 0; i < ndim_use; ++i) + { + for (int j = 0; j < ndim_use; ++j) + { + sum_beta += beta_tmp(j, i); + } + } + for (int i = 0; i < ndim_use; ++i) + { + coef[i] = 0.; + for (int j = 0; j < ndim_use; ++j) + { + coef[i] += beta_tmp(i, j); + } + coef[i] /= sum_beta; + } + delete[] work; + delete[] iwork; + } + else + { + beta(0, 0) = inner_product(FP_F, FP_F); + coef[0] = 1.0; + } + + ModuleBase::timer::tick("Charge", "Pulay_mixing"); +}; +} // namespace Base_Mixing \ No newline at end of file diff --git a/source/module_base/module_mixing/pulay_mixing.h b/source/module_base/module_mixing/pulay_mixing.h index 04bb0e4713..0e6c6de595 100644 --- a/source/module_base/module_mixing/pulay_mixing.h +++ b/source/module_base/module_mixing/pulay_mixing.h @@ -1,12 +1,7 @@ #ifndef PULAY_MIXING_H_ #define PULAY_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" -#include "module_base/global_variable.h" namespace Base_Mixing { @@ -24,21 +19,16 @@ namespace Base_Mixing class Pulay_Mixing : public Mixing { public: - Pulay_Mixing(const int& mixing_ndim, const double& mixing_beta) + Pulay_Mixing(const int& mixing_ndim) { this->mixing_ndim = mixing_ndim; this->data_ndim = mixing_ndim; - this->mixing_beta = mixing_beta; this->coef = std::vector(mixing_ndim); this->beta = ModuleBase::matrix(mixing_ndim, mixing_ndim, true); - if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4) - { - this->two_beta = 0; - } - else if (GlobalV::NSPIN == 2) - { - this->two_beta = 1; - } + } + Pulay_Mixing(const int& mixing_ndim, const double& mixing_beta) : Pulay_Mixing(mixing_ndim) + { + this->mixing_beta = mixing_beta; } virtual ~Pulay_Mixing() override { @@ -62,6 +52,7 @@ class Pulay_Mixing : public Mixing * @param data_in x_in * @param data_out x_out = f(x_in) * @param screen pointer to the screen function for Ker-Ker mixing + * @param mix (double* out, const double* in, const double* residual) calculate 'out' with 'in' and residual * @param need_calcoef whether need to calculate the coef * */ @@ -69,17 +60,20 @@ class Pulay_Mixing : public Mixing const double* data_in, const double* data_out, std::function screen, + std::function mix, const bool& need_calcoef) override { - this->tem_push_data(mdata, data_in, data_out, screen, need_calcoef); + this->tem_push_data(mdata, data_in, data_out, screen, mix, need_calcoef); }; - virtual void push_data(Mixing_Data& mdata, - const std::complex* data_in, - const std::complex* data_out, - std::function*)> screen, - const bool& need_calcoef) override + virtual void push_data( + Mixing_Data& mdata, + const std::complex* data_in, + const std::complex* data_out, + std::function*)> screen, + std::function*, const std::complex*, const std::complex*)> mix, + const bool& need_calcoef) override { - this->tem_push_data(mdata, data_in, data_out, screen, need_calcoef); + this->tem_push_data(mdata, data_in, data_out, screen, mix, need_calcoef); }; /** @@ -106,6 +100,7 @@ class Pulay_Mixing : public Mixing * @param data_in x_in * @param data_out x_out = f(x_in) * @param screen pointer to the screen function for Ker-Ker mixing + * @param mix (double* out, const double* in, const double* residual) calculate 'out' with 'in' and residual * @param need_calcoef whether need to calculate the coef * */ @@ -114,97 +109,8 @@ class Pulay_Mixing : public Mixing const FPTYPE* data_in, const FPTYPE* data_out, std::function screen, - const bool& need_calcoef) - { - const size_t length = mdata.length; - std::vector F_tmp(length); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) -#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 data(length); - if (this->two_beta == 0) - { - // rho_tot -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) -#endif - for (int i = 0; i < length; ++i) - { - data[i] = data_in[i] + this->mixing_beta * F_tmp[i]; - } - } - else if (this->two_beta == 1) - { - // rho_tot -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) -#endif - for (int i = 0; i < length / 2; ++i) - { - data[i] = data_in[i] + this->mixing_beta * F_tmp[i]; - } - // magnetism -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) -#endif - for (int i = length / 2; i < length; ++i) - { - data[i] = data_in[i] + GlobalV::MIXING_BETA_MAG * F_tmp[i]; - } - - } - - mdata.push(data.data()); - - if (!need_calcoef) - return; - - if (address != &mdata && address != nullptr) - ModuleBase::WARNING_QUIT( - "Pulay_Mixing", - "One Pulay_Mixing object can only bind one Mixing_Data object to calculate coefficients"); - - FPTYPE* FP_F = static_cast(F); - if (mdata.ndim_use == 1) - { - address = &mdata; - // allocate - if (F != nullptr) - free(F); - F = malloc(sizeof(FPTYPE) * length * mixing_ndim); - FP_F = static_cast(F); -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) -#endif - for (int i = 0; i < length; ++i) - { - FP_F[i] = F_tmp[i]; - } - } - else - { - start_F = (this->start_F + 1) % this->mixing_ndim; - FPTYPE* FP_startF = FP_F + start_F * length; -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE)) -#endif - for (int i = 0; i < length; ++i) - { - FP_startF[i] = F_tmp[i]; - } - } - }; + std::function mix, + const bool& need_calcoef); /** * @brief calculate coeficients for mixing @@ -213,90 +119,7 @@ class Pulay_Mixing : public Mixing * @param inner_product pointer to the inner dot function */ template - void tem_cal_coef(const Mixing_Data& mdata, std::function inner_product) - { - ModuleBase::TITLE("Charge_Mixing", "Pulay_mixing"); - ModuleBase::timer::tick("Charge", "Pulay_mixing"); - if (address != &mdata && address != nullptr) - ModuleBase::WARNING_QUIT( - "Pulay_mixing", - "One Pulay_Mixing object can only bind one Mixing_Data object to calculate coefficients"); - const int length = mdata.length; - FPTYPE* FP_F = static_cast(F); - - if (mdata.ndim_use > 1) - { - const int ndim_use = mdata.ndim_use; - ModuleBase::matrix beta_tmp(ndim_use, ndim_use); - // beta(i, j) = - for (int i = 0; i < ndim_use; ++i) - { - FPTYPE* Fi = FP_F + i * length; - for (int j = i; j < ndim_use; ++j) - { - if (i != start_F && j != start_F) - { - beta_tmp(i, j) = beta(i, j); - } - else - { - FPTYPE* Fj = FP_F + j * length; - beta(i, j) = beta_tmp(i, j) = inner_product(Fi, Fj); - } - if (j != i) - { - beta(j, i) = beta_tmp(j, i) = beta_tmp(i, j); - } - } - } - - double* work = new double[ndim_use]; - int* iwork = new int[ndim_use]; - char uu = 'U'; - int info; - dsytrf_(&uu, &ndim_use, beta_tmp.c, &ndim_use, iwork, work, &ndim_use, &info); - if (info != 0) - ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when factorizing beta."); - dsytri_(&uu, &ndim_use, beta_tmp.c, &ndim_use, iwork, work, &info); - if (info != 0) - ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when DSYTRI beta."); - for (int i = 0; i < ndim_use; ++i) - { - for (int j = i + 1; j < ndim_use; ++j) - { - beta_tmp(i, j) = beta_tmp(j, i); - } - } - - // coef{i} = \sum_j beta{ij} / \sum_k \sum_j beta{kj} - double sum_beta = 0.; - for (int i = 0; i < ndim_use; ++i) - { - for (int j = 0; j < ndim_use; ++j) - { - sum_beta += beta_tmp(j, i); - } - } - for (int i = 0; i < ndim_use; ++i) - { - coef[i] = 0.; - for (int j = 0; j < ndim_use; ++j) - { - coef[i] += beta_tmp(i, j); - } - coef[i] /= sum_beta; - } - delete[] work; - delete[] iwork; - } - else - { - beta(0, 0) = inner_product(FP_F, FP_F); - coef[0] = 1.0; - } - - ModuleBase::timer::tick("Charge", "Pulay_mixing"); - }; + void tem_cal_coef(const Mixing_Data& mdata, std::function inner_product); // F = data_out - data_in void* F = nullptr; @@ -308,8 +131,6 @@ class Pulay_Mixing : public Mixing int mixing_ndim = -1; // start index for F int start_F = 0; - // if we should considere two beta - int two_beta = 0; }; } // namespace Base_Mixing #endif \ No newline at end of file diff --git a/source/module_base/module_mixing/test/mixing_test.cpp b/source/module_base/module_mixing/test/mixing_test.cpp index 972c12ece9..ea0e2415d5 100755 --- a/source/module_base/module_mixing/test/mixing_test.cpp +++ b/source/module_base/module_mixing/test/mixing_test.cpp @@ -1,9 +1,10 @@ +#include + #include "../broyden_mixing.h" #include "../plain_mixing.h" #include "../pulay_mixing.h" #include "gmock/gmock.h" #include "gtest/gtest.h" -#include #define DOUBLETHRESHOLD 1e-8 double ext_inner_product_mock(double* x1, double* x2) @@ -65,7 +66,7 @@ class Mixing_Test : public testing::Test * [x3] [-6/12 -3/12 36/12][x3] */ template - void solve_linear_eq(FPTYPE* x_in, FPTYPE* x_out) + void solve_linear_eq(FPTYPE* x_in, FPTYPE* x_out, bool diff_beta = false) { this->mixing->init_mixing_data(xdata, 3, sizeof(FPTYPE)); std::vector delta_x(3); @@ -96,8 +97,25 @@ class Mixing_Test : public testing::Test { break; } - - this->mixing->push_data(this->xdata, x_in, x_out, screen, true); + if (diff_beta) + { + this->mixing->push_data( + this->xdata, + x_in, + x_out, + screen, + // mixing can use different mixing_beta for one vector + [](FPTYPE* out, const FPTYPE* in, const FPTYPE* sres) { + out[0] = in[0] + 0.5 * sres[0]; + out[1] = in[1] + 0.6 * sres[1]; + out[2] = in[2] + 0.5 * sres[2]; + }, + true); + } + else + { + this->mixing->push_data(this->xdata, x_in, x_out, screen, true); + } this->mixing->cal_coef(this->xdata, inner_product); @@ -136,10 +154,10 @@ TEST_F(Mixing_Test, BroydenSolveLinearEq) init_method("broyden"); std::vector x_in = xd_ref; std::vector x_out(3); - solve_linear_eq(x_in.data(), x_out.data()); - EXPECT_NEAR(x_out[0], 2.9999999999999996, DOUBLETHRESHOLD); - EXPECT_NEAR(x_out[1], 2.0000000000000004, DOUBLETHRESHOLD); - EXPECT_NEAR(x_out[2], 1.0000000000000000, DOUBLETHRESHOLD); + solve_linear_eq(x_in.data(), x_out.data(), true); + EXPECT_NEAR(x_out[0], 3.0, DOUBLETHRESHOLD); + EXPECT_NEAR(x_out[1], 2.0, DOUBLETHRESHOLD); + EXPECT_NEAR(x_out[2], 1.0, DOUBLETHRESHOLD); ASSERT_EQ(niter, 5); this->mixing->reset(); @@ -147,10 +165,10 @@ TEST_F(Mixing_Test, BroydenSolveLinearEq) std::vector> xc_in = xc_ref; std::vector> xc_out(3); - solve_linear_eq>(xc_in.data(), xc_out.data()); - EXPECT_NEAR(xc_out[0].real(), 3.0000000000000009, DOUBLETHRESHOLD); - EXPECT_NEAR(xc_out[1].real(), 1.9999999999999998, DOUBLETHRESHOLD); - EXPECT_NEAR(xc_out[2].real(), 0.99999999999999944, DOUBLETHRESHOLD); + solve_linear_eq>(xc_in.data(), xc_out.data(), true); + EXPECT_NEAR(xc_out[0].real(), 3.0, DOUBLETHRESHOLD); + EXPECT_NEAR(xc_out[1].real(), 2.0, DOUBLETHRESHOLD); + EXPECT_NEAR(xc_out[2].real(), 1.0, DOUBLETHRESHOLD); ASSERT_EQ(niter, 5); std::string output; Base_Mixing::Mixing_Data testdata; @@ -166,9 +184,7 @@ TEST_F(Mixing_Test, BroydenSolveLinearEq) testing::HasSubstr("One Broyden_Mixing object can only bind one Mixing_Data object to calculate coefficients")); testing::internal::CaptureStdout(); - EXPECT_EXIT(this->mixing->cal_coef(testdata, ext_inner_product_mock), - ::testing::ExitedWithCode(0), - ""); + EXPECT_EXIT(this->mixing->cal_coef(testdata, ext_inner_product_mock), ::testing::ExitedWithCode(0), ""); output = testing::internal::GetCapturedStdout(); EXPECT_THAT( output, @@ -214,9 +230,7 @@ TEST_F(Mixing_Test, PulaySolveLinearEq) testing::HasSubstr("One Pulay_Mixing object can only bind one Mixing_Data object to calculate coefficients")); testing::internal::CaptureStdout(); - EXPECT_EXIT(this->mixing->cal_coef(testdata, ext_inner_product_mock), - ::testing::ExitedWithCode(0), - ""); + EXPECT_EXIT(this->mixing->cal_coef(testdata, ext_inner_product_mock), ::testing::ExitedWithCode(0), ""); output = testing::internal::GetCapturedStdout(); EXPECT_THAT( output, @@ -248,6 +262,24 @@ TEST_F(Mixing_Test, PlainSolveLinearEq) EXPECT_NEAR(xc_out[2].real(), 1.0000211250842703, DOUBLETHRESHOLD); ASSERT_EQ(niter, 10); + // test mix_data of plain_mixing + std::vector x_tmp(3); + this->mixing->push_data(this->xdata, x_in.data(), x_out.data(), nullptr, true); + this->mixing->mix_data(this->xdata, x_tmp.data()); + Base_Mixing::Plain_Mixing plain_mix(mixing_beta); + plain_mix.plain_mix(x_in.data(), x_in.data(), x_out.data(), 3, [](double* x) {}); + EXPECT_NEAR(x_tmp[0], x_in[0], DOUBLETHRESHOLD); + EXPECT_NEAR(x_tmp[1], x_in[1], DOUBLETHRESHOLD); + EXPECT_NEAR(x_tmp[2], x_in[2], DOUBLETHRESHOLD); + + std::vector> xc_tmp(3); + this->mixing->push_data(this->xdata, xc_in.data(), xc_out.data(), nullptr, true); + this->mixing->mix_data(this->xdata, xc_tmp.data()); + plain_mix.plain_mix(xc_in.data(), xc_in.data(), xc_out.data(), 3, nullptr); + EXPECT_NEAR(xc_tmp[0].real(), xc_in[0].real(), DOUBLETHRESHOLD); + EXPECT_NEAR(xc_tmp[1].real(), xc_in[1].real(), DOUBLETHRESHOLD); + EXPECT_NEAR(xc_tmp[2].real(), xc_in[2].real(), DOUBLETHRESHOLD); + this->mixing->reset(); clear(); @@ -255,14 +287,14 @@ TEST_F(Mixing_Test, PlainSolveLinearEq) TEST_F(Mixing_Test, OtherCover) { - this->mixing = new Base_Mixing::Broyden_Mixing(2, 0.7); - Base_Mixing::Mixing_Data nodata; - this->mixing->init_mixing_data(nodata, 0, sizeof(double)); - this->mixing->push_data(nodata, (double*)nullptr, (double*)nullptr, nullptr, false); - this->mixing->push_data(nodata, (double*)nullptr, (double*)nullptr, nullptr, false); - this->mixing->mix_data(nodata, (double*)nullptr); - this->mixing->mix_data(nodata, (std::complex*)nullptr); - EXPECT_EQ(nodata.length, 0); - - clear(); + this->mixing = new Base_Mixing::Broyden_Mixing(2, 0.7); + Base_Mixing::Mixing_Data nodata; + this->mixing->init_mixing_data(nodata, 0, sizeof(double)); + this->mixing->push_data(nodata, (double*)nullptr, (double*)nullptr, nullptr, false); + this->mixing->push_data(nodata, (double*)nullptr, (double*)nullptr, nullptr, false); + this->mixing->mix_data(nodata, (double*)nullptr); + this->mixing->mix_data(nodata, (std::complex*)nullptr); + EXPECT_EQ(nodata.length, 0); + + clear(); } \ No newline at end of file diff --git a/source/module_elecstate/module_charge/charge.h b/source/module_elecstate/module_charge/charge.h index b0ca0b0cdd..55bf2ee369 100644 --- a/source/module_elecstate/module_charge/charge.h +++ b/source/module_elecstate/module_charge/charge.h @@ -178,7 +178,7 @@ class Charge int nxyz; // total nuber of r vectors int ngmc; // number of g vectors in this processor int nspin; // number of spins - ModulePW::PW_Basis* rhopw = nullptr; + ModulePW::PW_Basis* rhopw = nullptr;// When double_grid is used, rhopw = rhodpw (dense grid) private: void destroy(); // free arrays liuyu 2023-03-12 diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index ff175c85f3..44f1b55ed2 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -16,16 +16,19 @@ Charge_Mixing::Charge_Mixing() Charge_Mixing::~Charge_Mixing() { delete this->mixing; + delete this->mixing_highf; } void Charge_Mixing::set_mixing(const std::string& mixing_mode_in, const double& mixing_beta_in, const int& mixing_ndim_in, const double& mixing_gg0_in, - const bool& mixing_tau_in) + const bool& mixing_tau_in, + const double& mixing_beta_mag_in) { this->mixing_mode = mixing_mode_in; this->mixing_beta = mixing_beta_in; + this->mixing_beta_mag = mixing_beta_mag_in; this->mixing_ndim = mixing_ndim_in; this->mixing_gg0 = mixing_gg0_in; this->mixing_tau = mixing_tau_in; @@ -61,6 +64,13 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in, { ModuleBase::WARNING_QUIT("Charge_Mixing", "This Mixing mode is not implemended yet,coming soon."); } + if (GlobalV::double_grid) + { + // ONLY smooth part of charge density is mixed by specific mixing method + // The high_frequency part is mixed by plain mixing method. + delete this->mixing_highf; + this->mixing_highf = new Base_Mixing::Plain_Mixing(this->mixing_beta); + } if (GlobalV::SCF_THR_TYPE == 1) { @@ -82,54 +92,55 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in, return; } -void Charge_Mixing::set_rhopw(ModulePW::PW_Basis* rhopw_in) +void Charge_Mixing::set_rhopw(ModulePW::PW_Basis* rhopw_in, ModulePW::PW_Basis* rhodpw_in) { this->rhopw = rhopw_in; + this->rhodpw = rhodpw_in; } -void Charge_Mixing::need_auto_set() -{ - this->autoset = true; -} - -void Charge_Mixing::auto_set(const double& bandgap_in, const UnitCell& ucell_) -{ - // auto set parameters once - if (!this->autoset) - { - return; - } - else - { - this->autoset = false; - } - GlobalV::ofs_running << "--------------AUTO-SET---------------" << std::endl; - // 0.8 for nspin=1 and 0.4 for others - if (GlobalV::NSPIN == 1) - { - this->mixing->mixing_beta = this->mixing_beta = 0.8; - GlobalV::ofs_running << " Autoset mixing_beta to " << this->mixing_beta << std::endl; - } - else - { - this->mixing->mixing_beta = this->mixing_beta = 0.4; - GlobalV::MIXING_BETA_MAG = 4 * this->mixing_beta; - GlobalV::ofs_running << " Autoset mixing_beta to " << this->mixing_beta << std::endl; - GlobalV::ofs_running << " Autoset mixing_beta_mag to " << GlobalV::MIXING_BETA_MAG << std::endl; - } - GlobalV::MIXING_BETA = mixing_beta; - - // auto set kerker mixing_gg0 = 1.0 as default - this->mixing_gg0 = 1.0; - GlobalV::ofs_running << " Autoset mixing_gg0 to " << this->mixing_gg0 << std::endl; - if (GlobalV::NSPIN == 1) - { - GlobalV::MIXING_GG0_MAG = 0.0; - GlobalV::ofs_running << " Autoset mixing_gg0_mag to " << GlobalV::MIXING_GG0_MAG << std::endl; - } - - GlobalV::ofs_running << "-------------------------------------" << std::endl; -} +// void Charge_Mixing::need_auto_set() +// { +// this->autoset = true; +// } + +// void Charge_Mixing::auto_set(const double& bandgap_in, const UnitCell& ucell_) +// { +// // auto set parameters once +// if (!this->autoset) +// { +// return; +// } +// else +// { +// this->autoset = false; +// } +// GlobalV::ofs_running << "--------------AUTO-SET---------------" << std::endl; +// // 0.8 for nspin=1 and 0.4 for others +// if (GlobalV::NSPIN == 1) +// { +// this->mixing->mixing_beta = this->mixing_beta = 0.8; +// GlobalV::ofs_running << " Autoset mixing_beta to " << this->mixing_beta << std::endl; +// } +// else +// { +// this->mixing->mixing_beta = this->mixing_beta = 0.4; +// GlobalV::MIXING_BETA_MAG = 4 * this->mixing_beta; +// GlobalV::ofs_running << " Autoset mixing_beta to " << this->mixing_beta << std::endl; +// GlobalV::ofs_running << " Autoset mixing_beta_mag to " << GlobalV::MIXING_BETA_MAG << std::endl; +// } +// GlobalV::MIXING_BETA = mixing_beta; + +// // auto set kerker mixing_gg0 = 1.0 as default +// this->mixing_gg0 = 1.0; +// GlobalV::ofs_running << " Autoset mixing_gg0 to " << this->mixing_gg0 << std::endl; +// if (GlobalV::NSPIN == 1) +// { +// GlobalV::MIXING_GG0_MAG = 0.0; +// GlobalV::ofs_running << " Autoset mixing_gg0_mag to " << GlobalV::MIXING_GG0_MAG << std::endl; +// } + +// GlobalV::ofs_running << "-------------------------------------" << std::endl; +// } double Charge_Mixing::get_drho(Charge* chr, const double nelec) { @@ -204,116 +215,83 @@ void Charge_Mixing::mix_rho_recip(Charge* chr) // ONLY smooth part of charge density is mixed by specific mixing method // The high_frequency part is mixed by plain mixing method. // NOTE: chr->rhopw is dense, while this->rhopw is smooth - std::complex* rhog_in = nullptr; - std::complex* rhog_out = nullptr; + std::complex*rhogs_in = chr->rhog_save[0], *rhogs_out = chr->rhog[0]; + std::complex*rhoghf_in = nullptr, *rhoghf_out = nullptr; if (GlobalV::double_grid) { - rhog_in = new std::complex[GlobalV::NSPIN * this->rhopw->npw]; - rhog_out = new std::complex[GlobalV::NSPIN * this->rhopw->npw]; - for (int is = 0; is < GlobalV::NSPIN; is++) - { - std::memcpy(&rhog_in[is * rhopw->npw], chr->rhog_save[is], rhopw->npw * sizeof(std::complex)); - std::memcpy(&rhog_out[is * rhopw->npw], chr->rhog[is], rhopw->npw * sizeof(std::complex)); - } - } - else - { - rhog_in = chr->rhog_save[0]; - rhog_out = chr->rhog[0]; + // divide into smooth part and high_frequency part + divide_data(chr->rhog_save[0], rhogs_in, rhoghf_in); + divide_data(chr->rhog[0], rhogs_out, rhoghf_out); } auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip, this, std::placeholders::_1); - this->mixing->push_data(this->rho_mdata, rhog_in, rhog_out, screen, true); + this->mixing->push_data(this->rho_mdata, rhogs_in, rhogs_out, screen, true); auto inner_product = std::bind(&Charge_Mixing::inner_product_recip, this, std::placeholders::_1, std::placeholders::_2); this->mixing->cal_coef(this->rho_mdata, inner_product); - this->mixing->mix_data(this->rho_mdata, rhog_out); + this->mixing->mix_data(this->rho_mdata, rhogs_out); if (GlobalV::double_grid) { - // simple mixing for high_frequencies - this->high_freq_mix(chr->rhog[0], chr->rhog_save[0], GlobalV::NSPIN * chr->rhopw->npw); + // plain mixing for high_frequencies + const int ndimhf = (this->rhodpw->npw - this->rhopw->npw) * GlobalV::NSPIN; + this->mixing_highf->plain_mix(rhoghf_out, rhoghf_in, rhoghf_out, ndimhf, nullptr); // combine smooth part and high_frequency part - for (int is = 0; is < GlobalV::NSPIN; is++) - { - std::memcpy(chr->rhog[is], &rhog_out[is * rhopw->npw], rhopw->npw * sizeof(std::complex)); - } - - delete[] rhog_in; - delete[] rhog_out; + combine_data(chr->rhog[0], rhogs_out, rhoghf_out); + clean_data(rhogs_in, rhoghf_in); } // rhog to rho for (int is = 0; is < GlobalV::NSPIN; is++) { - chr->rhopw->recip2real(chr->rhog[is], chr->rho[is]); + rhodpw->recip2real(chr->rhog[is], chr->rho[is]); } // For kinetic energy density if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) { - std::vector> kin_g(GlobalV::NSPIN * chr->rhopw->npw); - std::vector> kin_g_save(GlobalV::NSPIN * chr->rhopw->npw); - std::complex* taug_in = nullptr; - std::complex* taug_out = nullptr; - if (GlobalV::double_grid) + std::vector> kin_g(GlobalV::NSPIN * rhodpw->npw); + std::vector> kin_g_save(GlobalV::NSPIN * rhodpw->npw); + + for (int is = 0; is < GlobalV::NSPIN; ++is) { - taug_in = new std::complex[GlobalV::NSPIN * this->rhopw->npw]; - taug_out = new std::complex[GlobalV::NSPIN * this->rhopw->npw]; + rhodpw->real2recip(chr->kin_r[is], &kin_g[is * rhodpw->npw]); + rhodpw->real2recip(chr->kin_r_save[is], &kin_g_save[is * rhodpw->npw]); } - else + std::complex*taugs_in = kin_g_save.data(), *taugs_out = kin_g.data(); + std::complex*taughf_in = nullptr, *taughf_out = nullptr; + if (GlobalV::double_grid) { - taug_in = kin_g_save.data(); - taug_out = kin_g.data(); + // divide into smooth part and high_frequency part + divide_data(kin_g_save.data(), taugs_in, taughf_in); + divide_data(kin_g.data(), taugs_out, taughf_out); } - for (int is = 0; is < GlobalV::NSPIN; ++is) - { - chr->rhopw->real2recip(chr->kin_r[is], &kin_g[is * chr->rhopw->npw]); - chr->rhopw->real2recip(chr->kin_r_save[is], &kin_g_save[is * chr->rhopw->npw]); - - // similar to rhog and rhog_save - if (GlobalV::double_grid) - { - std::memcpy(&taug_in[is * rhopw->npw], - &kin_g_save[is * chr->rhopw->npw], - rhopw->npw * sizeof(std::complex)); - std::memcpy(&taug_out[is * rhopw->npw], - &kin_g[is * chr->rhopw->npw], - rhopw->npw * sizeof(std::complex)); - } - } // Note: there is no kerker modification for tau because I'm not sure // if we should have it. If necessary we can try it in the future. - this->mixing->push_data(this->tau_mdata, taug_in, taug_out, nullptr, false); + this->mixing->push_data(this->tau_mdata, taugs_in, taugs_out, nullptr, false); - this->mixing->mix_data(this->tau_mdata, taug_out); + this->mixing->mix_data(this->tau_mdata, taugs_out); if (GlobalV::double_grid) { // simple mixing for high_frequencies - this->high_freq_mix(kin_g.data(), kin_g_save.data(), GlobalV::NSPIN * chr->rhopw->npw); + const int ndimhf = (this->rhodpw->npw - this->rhopw->npw) * GlobalV::NSPIN; + this->mixing_highf->plain_mix(taughf_out, taughf_in, taughf_out, ndimhf, nullptr); // combine smooth part and high_frequency part - for (int is = 0; is < GlobalV::NSPIN; is++) - { - std::memcpy(&kin_g[is * chr->rhopw->npw], - &taug_out[is * rhopw->npw], - rhopw->npw * sizeof(std::complex)); - } - - delete[] taug_in; - delete[] taug_out; + combine_data(kin_g.data(), taugs_out, taughf_out); + clean_data(taugs_in, taughf_in); } // kin_g to kin_r for (int is = 0; is < GlobalV::NSPIN; is++) { - chr->rhopw->recip2real(&kin_g[is * chr->rhopw->npw], chr->kin_r[is]); + rhodpw->recip2real(&kin_g[is * rhodpw->npw], chr->kin_r[is]); } } @@ -349,6 +327,8 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr) { rhog_in = chr->rhog_save[0]; rhog_out = chr->rhog[0]; + auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip_new, this, std::placeholders::_1); + this->mixing->push_data(this->rho_mdata, rhog_in, rhog_out, screen, true); } else if (GlobalV::NSPIN == 2) { @@ -356,10 +336,29 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr) chr->get_rhog_mag(); rhog_in = chr->rhog_mag_save; rhog_out = chr->rhog_mag; + const int npw = this->rhopw->npw; + auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip_new, this, std::placeholders::_1); + auto twobeta_mix + = [this, npw](std::complex* out, const std::complex* in, const std::complex* sres) { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 256) +#endif + for (int i = 0; i < npw; ++i) + { + out[i] = in[i] + this->mixing_beta * sres[i]; + } + // magnetism +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 256) +#endif + for (int i = npw; i < 2 * npw; ++i) + { + out[i] = in[i] + this->mixing_beta_mag * sres[i]; + } + }; + this->mixing->push_data(this->rho_mdata, rhog_in, rhog_out, screen, twobeta_mix, true); } - auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip_new, this, std::placeholders::_1); - this->mixing->push_data(this->rho_mdata, rhog_in, rhog_out, screen, true); // can choose inner_product_recip_new1 or inner_product_recip_new2 auto inner_product_new = std::bind(&Charge_Mixing::inner_product_recip_new1, this, std::placeholders::_1, std::placeholders::_2); @@ -557,11 +556,11 @@ void Charge_Mixing::mix_rho(Charge* chr) // --------------------Mixing Body-------------------- if (GlobalV::SCF_THR_TYPE == 1) { - if (GlobalV::double_grid) + if (GlobalV::double_grid || GlobalV::use_paw) { mix_rho_recip(chr); } - else // new mixing method do not support double_grid yet + else // new mixing method do not support double_grid or paw yet { mix_rho_recip_new(chr); } @@ -1058,20 +1057,76 @@ double Charge_Mixing::rhog_dot_product(const std::complex* const* const return sum; } -void Charge_Mixing::high_freq_mix(std::complex* data, - const std::complex* data_save, - const int& number) const +void Charge_Mixing::divide_data(std::complex* data_d, + std::complex*& data_s, + std::complex*& data_hf) { - ModuleBase::TITLE("Charge_Mixing", "high_freq_mix"); - ModuleBase::timer::tick("Charge_Mixing", "high_freq_mix"); - - if (GlobalV::double_grid) + ModuleBase::TITLE("Charge_Mixing", "divide_data"); + if (GlobalV::NSPIN == 1) { - for (int ig = 0; ig < number; ig++) + data_s = data_d; + data_hf = data_d + this->rhopw->npw; + } + else + { + const int ndimd = this->rhodpw->npw; + const int ndims = this->rhopw->npw; + const int ndimhf = ndimd - ndims; + data_s = new std::complex[GlobalV::NSPIN * ndims]; + data_hf = nullptr; + if (ndimhf > 0) { - data[ig] = data_save[ig] + mixing_beta * (data[ig] - data_save[ig]); + data_hf = new std::complex[GlobalV::NSPIN * ndimhf]; } + for (int is = 0; is < GlobalV::NSPIN; ++is) + { + std::memcpy(data_s + is * ndims, data_d + is * ndimd, ndims * sizeof(std::complex)); + std::memcpy(data_hf + is * ndimhf, data_d + is * ndimd + ndims, ndimhf * sizeof(std::complex)); + } + } +} +void Charge_Mixing::combine_data(std::complex* data_d, + std::complex*& data_s, + std::complex*& data_hf) +{ + ModuleBase::TITLE("Charge_Mixing", "combine_data"); + if (GlobalV::NSPIN == 1) + { + data_s = nullptr; + data_hf = nullptr; + return; } + else + { + const int ndimd = this->rhodpw->npw; + const int ndims = this->rhopw->npw; + const int ndimhf = ndimd - ndims; + for (int is = 0; is < GlobalV::NSPIN; ++is) + { + std::memcpy(data_d + is * ndimd, data_s + is * ndims, ndims * sizeof(std::complex)); + std::memcpy(data_d + is * ndimd + ndims, data_hf + is * ndimhf, ndimhf * sizeof(std::complex)); + } + delete[] data_s; + delete[] data_hf; + data_s = nullptr; + data_hf = nullptr; + } +} - ModuleBase::timer::tick("Charge_Mixing", "high_freq_mix"); +void Charge_Mixing::clean_data(std::complex*& data_s, std::complex*& data_hf) +{ + ModuleBase::TITLE("Charge_Mixing", "clean_data"); + if (GlobalV::NSPIN == 1) + { + data_s = nullptr; + data_hf = nullptr; + return; + } + else + { + delete[] data_s; + delete[] data_hf; + data_s = nullptr; + data_hf = nullptr; + } } \ No newline at end of file diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index f101c916b1..f941a7fe5e 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -5,6 +5,7 @@ #include "module_base/global_function.h" #include "module_base/global_variable.h" #include "module_base/module_mixing/mixing.h" +#include "module_base/module_mixing/plain_mixing.h" #include "module_cell/unitcell.h" class Charge_Mixing { @@ -16,6 +17,8 @@ class Charge_Mixing Base_Mixing::Mixing_Data tau_mdata; Base_Mixing::Mixing_Data nhat_mdata; + Base_Mixing::Plain_Mixing* mixing_highf = nullptr; ///< The high_frequency part is mixed by plain mixing method. + /** * @brief reset mixing * @@ -82,19 +85,20 @@ class Charge_Mixing const double& mixing_beta_in, const int& mixing_ndim_in, const double& mixing_gg0_in, - const bool& mixing_tau_in); // mohan add mixing_gg0_in 2014-09-27 + const bool& mixing_tau_in, + const double& mixing_beta_mag_in = 1.6); // mohan add mixing_gg0_in 2014-09-27 - /** - * @brief use auto set - * - */ - void need_auto_set(); + // /** + // * @brief use auto set + // * + // */ + // void need_auto_set(); - /** - * @brief auto set mixing gg0 and mixing_beta - * - */ - void auto_set(const double& bandgap_in, const UnitCell& ucell_); + // /** + // * @brief auto set mixing gg0 and mixing_beta + // * + // */ + // void auto_set(const double& bandgap_in, const UnitCell& ucell_); /** * @brief Get the drho @@ -102,8 +106,8 @@ class Charge_Mixing */ double get_drho(Charge* chr, const double nelec); - // init pwrho, sunliang add 2023-05-08 - void set_rhopw(ModulePW::PW_Basis* rhopw_in); + // init pwrho and rhodpw + void set_rhopw(ModulePW::PW_Basis* rhopw_in, ModulePW::PW_Basis* rhodpw_in); // extracting parameters // normally these parameters will not be used @@ -132,20 +136,36 @@ class Charge_Mixing //====================================== std::string mixing_mode = "broyden"; double mixing_beta = 0.8; + double mixing_beta_mag = 1.6; int mixing_ndim = 8; double mixing_gg0 = 0.0; // mohan add 2014-09-27 bool mixing_tau = false; bool new_e_iteration = true; - ModulePW::PW_Basis* rhopw = nullptr; + ModulePW::PW_Basis* rhopw = nullptr; ///< smooth grid + ModulePW::PW_Basis* rhodpw = nullptr; ///< dense grid, same as rhopw for ncpp. bool autoset = false; private: double rhog_dot_product(const std::complex* const* const rhog1, const std::complex* const* const rhog2) const; - void high_freq_mix(std::complex* data, const std::complex* data_save, const int& number) const; + /** + * @brief divide rho/tau to smooth and high frequency parts + * + */ + void divide_data(std::complex* data_d, std::complex*& data_s, std::complex*& data_hf); + /** + * @brief gather smooth and high frequency parts to rho/tau + * + */ + void combine_data(std::complex* data_d, std::complex*& data_s, std::complex*& data_hf); + /** + * @brief clean smooth and high frequency parts + * + */ + void clean_data(std::complex*& data_s, std::complex*& data_hf); }; #endif diff --git a/source/module_elecstate/test/charge_mixing_test.cpp b/source/module_elecstate/test/charge_mixing_test.cpp index 78e2f65edd..794b4c7345 100644 --- a/source/module_elecstate/test/charge_mixing_test.cpp +++ b/source/module_elecstate/test/charge_mixing_test.cpp @@ -241,6 +241,10 @@ UnitCell ucell; * Charge_Mixing::mix_rho_recip(chr) * Charge_Mixing::mix_rho_real(chr) * - mix rho with different methods + * - MixDivCombTest: Charge_Mixing::divide_data + * Charge_Mixing::combine_data + * Charge_Mixing::clean_data + * - divide and combine data * */ @@ -254,8 +258,13 @@ class ChargeMixingTest : public ::testing::Test pw_basis.initparameters(false, 20); pw_basis.setuptransform(); pw_basis.collect_local_pw(); + pw_dbasis.initgrids(4, ModuleBase::Matrix3(1, 0, 0, 0, 1, 0, 0, 0, 1), 40); + pw_dbasis.initparameters(false, 40); + pw_dbasis.setuptransform(&pw_basis); + pw_dbasis.collect_local_pw(); } ModulePW::PW_Basis pw_basis; + ModulePW::PW_Basis_Sup pw_dbasis; Charge charge; }; @@ -264,7 +273,7 @@ TEST_F(ChargeMixingTest, SetMixingTest) omp_set_num_threads(1); GlobalV::NSPIN = 1; Charge_Mixing CMtest; - CMtest.set_rhopw(&pw_basis); + CMtest.set_rhopw(&pw_basis, &pw_basis); double beta = 1.0; int dim = 1; double gg0 = 1; @@ -310,7 +319,7 @@ TEST_F(ChargeMixingTest, SetMixingTest) TEST_F(ChargeMixingTest, AutoSetTest) { Charge_Mixing CMtest; - CMtest.set_rhopw(&pw_basis); + CMtest.set_rhopw(&pw_basis, &pw_basis); GlobalV::SCF_THR_TYPE = 1; GlobalV::NSPIN = 1; @@ -343,7 +352,7 @@ TEST_F(ChargeMixingTest, AutoSetTest) TEST_F(ChargeMixingTest, KerkerScreenRecipTest) { Charge_Mixing CMtest; - CMtest.set_rhopw(&pw_basis); + CMtest.set_rhopw(&pw_basis, &pw_basis); GlobalV::NSPIN = 1; GlobalC::ucell.tpiba = 1.0; @@ -399,7 +408,7 @@ TEST_F(ChargeMixingTest, KerkerScreenRecipTest) TEST_F(ChargeMixingTest, KerkerScreenRecipNewTest) { Charge_Mixing CMtest; - CMtest.set_rhopw(&pw_basis); + CMtest.set_rhopw(&pw_basis, &pw_basis); GlobalC::ucell.tpiba = 1.0; // for new kerker @@ -460,7 +469,7 @@ TEST_F(ChargeMixingTest, InnerDotTest) { // REAL Charge_Mixing CMtest; - CMtest.set_rhopw(&pw_basis); + CMtest.set_rhopw(&pw_basis, &pw_basis); GlobalV::NSPIN = 1; std::vector drhor1(pw_basis.nrxx); std::vector drhor2(pw_basis.nrxx); @@ -528,7 +537,7 @@ TEST_F(ChargeMixingTest, InnerDotTest) TEST_F(ChargeMixingTest, InnerDotNewTest) { Charge_Mixing CMtest; - CMtest.set_rhopw(&pw_basis); + CMtest.set_rhopw(&pw_basis, &pw_basis); GlobalV::NSPIN = 1; // inner_product_recip_new1 @@ -621,7 +630,7 @@ TEST_F(ChargeMixingTest, MixRhoTest) //--------------------------------MAIN BODY-------------------------------- // RECIPROCAL Charge_Mixing CMtest_recip; - CMtest_recip.set_rhopw(&pw_basis); + CMtest_recip.set_rhopw(&pw_basis, &pw_basis); GlobalV::SCF_THR_TYPE = 1; CMtest_recip.set_mixing(mode, beta, dim, gg0, mixingtau); CMtest_recip.mix_reset(); @@ -652,7 +661,7 @@ TEST_F(ChargeMixingTest, MixRhoTest) // REAL Charge_Mixing CMtest_real; GlobalV::SCF_THR_TYPE = 2; - CMtest_real.set_rhopw(&pw_basis); + CMtest_real.set_rhopw(&pw_basis, &pw_basis); CMtest_real.set_mixing(mode, beta, dim, gg0, mixingtau); CMtest_real.mix_reset(); for(int i = 0 ; i < nspin * nrxx; ++i) @@ -685,33 +694,139 @@ TEST_F(ChargeMixingTest, MixRhoTest) delete[] charge.kin_r_save; } -TEST_F(ChargeMixingTest, HighFreqMixTest) +TEST_F(ChargeMixingTest, MixDoubleGridRhoTest) { - const int number = 10; - std::complex* data = new std::complex[number]; - std::complex* data_save = new std::complex[number]; - - // initialize - for (int i = 0; i < number; i++) + GlobalV::double_grid = true; + charge.set_rhopw(&pw_dbasis); + const int nspin = GlobalV::NSPIN = 1; + GlobalV::DOMAG_Z = false; + FUNC_TYPE = 3; + double beta = 0.7; + int dim = 1; + double gg0 = 0.0; + bool mixingtau = true; + std::string mode = "plain"; + const int nrxx = pw_dbasis.nrxx; + const int npw = pw_dbasis.npw; + charge._space_rho = new double[nspin * nrxx]; + charge._space_rho_save = new double[nspin * nrxx]; + charge._space_rhog = new std::complex[nspin * npw]; + charge._space_rhog_save = new std::complex[nspin * npw]; + charge._space_kin_r = new double[nspin * nrxx]; + charge._space_kin_r_save = new double[nspin * nrxx]; + charge.rho = new double*[nspin]; + charge.rhog = new std::complex*[nspin]; + charge.rho_save = new double*[nspin]; + charge.rhog_save = new std::complex*[nspin]; + charge.kin_r = new double*[nspin]; + charge.kin_r_save = new double*[nspin]; + for (int is = 0; is < nspin; is++) { - data[i] = {1.0, 1.0}; - data_save[i] = {2.0, 2.0}; + charge.rho[is] = charge._space_rho + is * nrxx; + charge.rhog[is] = charge._space_rhog + is * npw; + charge.rho_save[is] = charge._space_rho_save + is * nrxx; + charge.rhog_save[is] = charge._space_rhog_save + is * npw; + charge.kin_r[is] = charge._space_kin_r + is * nrxx; + charge.kin_r_save[is] = charge._space_kin_r_save + is * nrxx; + } + std::vector real_ref(nspin * nrxx); + std::vector real_save_ref(nspin * nrxx); + std::vector> recip_ref(nspin * npw); + std::vector> recip_save_ref(nspin * npw); + for (int i = 0; i < nspin * npw; ++i) + { + recip_ref[i] = std::complex(double(i), 1.0); + recip_save_ref[i] = std::complex(double(i), 0.0); + } + for (int i = 0; i < nspin; ++i) + { + pw_dbasis.recip2real(recip_ref.data() + i * npw, real_ref.data() + i * nrxx); + pw_dbasis.recip2real(recip_save_ref.data() + i * npw, real_save_ref.data() + i * nrxx); + } + //--------------------------------MAIN BODY-------------------------------- + // RECIPROCAL + Charge_Mixing CMtest_recip; + CMtest_recip.set_rhopw(&pw_basis, &pw_dbasis); + GlobalV::SCF_THR_TYPE = 1; + CMtest_recip.set_mixing(mode, beta, dim, gg0, mixingtau); + CMtest_recip.mix_reset(); + for (int i = 0; i < nspin * npw; ++i) + { + charge._space_rhog[i] = recip_ref[i]; + charge._space_rhog_save[i] = recip_save_ref[i]; + } + for (int i = 0; i < nspin * nrxx; ++i) + { + charge._space_rho[i] = real_ref[i]; + charge._space_rho_save[i] = real_save_ref[i]; + } + CMtest_recip.mix_rho(&charge); + for (int is = 0; is < nspin; ++is) + { + for (int ir = 0; ir < nrxx; ++ir) + { + EXPECT_NEAR(charge.rho_save[is][ir], real_ref[is * nrxx + ir], 1e-8); + } + for (int ig = 0; ig < npw; ++ig) + { + EXPECT_NEAR(charge.rhog[is][ig].real(), recip_save_ref[is * npw + ig].real(), 1e-8); + EXPECT_NEAR(charge.rhog[is][ig].imag(), recip_save_ref[is * npw + ig].imag() + 0.7, 1e-8); + } } - Charge_Mixing mixer; - mixer.set_rhopw(&pw_basis); - mixer.set_mixing("broyden", 1.0, 1, 0.2, false); - mixer.high_freq_mix(data, data_save, number); + //------------------------------------------------------------------------- + delete[] charge._space_rho; + delete[] charge._space_rho_save; + delete[] charge._space_rhog; + delete[] charge._space_rhog_save; + delete[] charge._space_kin_r; + delete[] charge._space_kin_r_save; + delete[] charge.rho; + delete[] charge.rhog; + delete[] charge.rho_save; + delete[] charge.rhog_save; + delete[] charge.kin_r; + delete[] charge.kin_r_save; +} - for (int i = 0; i < number; i++) +TEST_F(ChargeMixingTest, MixDivCombTest) +{ + // NSPIN = 1 + GlobalV::NSPIN = 1; + Charge_Mixing CMtest; + CMtest.set_rhopw(&pw_basis, &pw_dbasis); + std::vector> data(pw_dbasis.npw, 1.0); + std::complex*datas, *datahf; + std::complex*datas2, *datahf2; + CMtest.divide_data(data.data(), datas, datahf); + EXPECT_EQ(datas, data.data()); + EXPECT_EQ(datahf, data.data() + pw_basis.npw); + CMtest.combine_data(data.data(), datas, datahf); + EXPECT_EQ(datas, nullptr); + EXPECT_EQ(datahf, nullptr); + + CMtest.divide_data(data.data(), datas2, datahf2); + CMtest.clean_data(datas2, datahf2); + EXPECT_EQ(datas2, nullptr); + EXPECT_EQ(datahf2, nullptr); + + // NSPIN = 2 + GlobalV::NSPIN = 2; + data.resize(pw_dbasis.npw * 2, 1.0); + std::vector> dataout(pw_dbasis.npw * 2, 1.0); + CMtest.divide_data(data.data(), datas, datahf); + CMtest.combine_data(dataout.data(), datas, datahf); + EXPECT_EQ(datas, nullptr); + EXPECT_EQ(datahf, nullptr); + for (int i = 0; i < pw_dbasis.npw * 2; ++i) { - std::complex expected = data_save[i] + mixer.mixing_beta * (data[i] - data_save[i]); - EXPECT_DOUBLE_EQ(data[i].real(), expected.real()); - EXPECT_DOUBLE_EQ(data[i].imag(), expected.imag()); + EXPECT_EQ(dataout[i], data[i]); } - delete[] data; - delete[] data_save; + CMtest.divide_data(data.data(), datas2, datahf2); + CMtest.clean_data(datas2, datahf2); + EXPECT_EQ(datas2, nullptr); + EXPECT_EQ(datahf2, nullptr); } #undef private diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index a16b24065c..e532492cd0 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -47,7 +47,7 @@ namespace ModuleESolver /// charge mixing ///---------------------------------------------------------- p_chgmix = new Charge_Mixing(); - p_chgmix->set_rhopw(this->pw_rho); + p_chgmix->set_rhopw(this->pw_rho, this->pw_rhod); ///---------------------------------------------------------- /// wavefunc @@ -78,7 +78,8 @@ namespace ModuleESolver GlobalV::MIXING_BETA, GlobalV::MIXING_NDIM, GlobalV::MIXING_GG0, - GlobalV::MIXING_TAU); + GlobalV::MIXING_TAU, + GlobalV::MIXING_BETA_MAG); // I use default value to replace autoset // using bandgap to auto set mixing_beta // if (std::abs(GlobalV::MIXING_BETA + 10.0) < 1e-6) @@ -420,21 +421,21 @@ namespace ModuleESolver { //----------charge mixing--------------- //before first calling mix_rho(), bandgap and cell structure can be analyzed to get better default parameters - if(iter == 1) - { - double bandgap_for_autoset = 0.0; - if (!GlobalV::TWO_EFERMI) - { - this->pelec->cal_bandgap(); - bandgap_for_autoset = this->pelec->bandgap; - } - else - { - this->pelec->cal_bandgap_updw(); - bandgap_for_autoset = std::min(this->pelec->bandgap_up, this->pelec->bandgap_dw); - } - p_chgmix->auto_set(bandgap_for_autoset, GlobalC::ucell); - } + // if(iter == 1) + // { + // double bandgap_for_autoset = 0.0; + // if (!GlobalV::TWO_EFERMI) + // { + // this->pelec->cal_bandgap(); + // bandgap_for_autoset = this->pelec->bandgap; + // } + // else + // { + // this->pelec->cal_bandgap_updw(); + // bandgap_for_autoset = std::min(this->pelec->bandgap_up, this->pelec->bandgap_dw); + // } + // p_chgmix->auto_set(bandgap_for_autoset, GlobalC::ucell); + // } p_chgmix->mix_rho(pelec->charge); //----------charge mixing done----------- diff --git a/source/module_io/test/for_testing_input_conv.h b/source/module_io/test/for_testing_input_conv.h index 0004b3c89e..0cdaf8b82d 100644 --- a/source/module_io/test/for_testing_input_conv.h +++ b/source/module_io/test/for_testing_input_conv.h @@ -286,11 +286,10 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in, { return; } -// void Charge_Mixing::need_auto_set(){} -void Charge_Mixing::need_auto_set() -{ - this->autoset = true; -} +// void Charge_Mixing::need_auto_set() +// { +// this->autoset = true; +// } void Occupy::decision(const std::string& name, const std::string& smearing_method, const double& smearing_sigma) { return; From 25381a9c35dda52dfc4dcc9769fba6bfad17d69a Mon Sep 17 00:00:00 2001 From: Qianruipku Date: Thu, 23 Nov 2023 07:29:50 +0000 Subject: [PATCH 2/7] update results --- .../101_PW_15_pseudopots_LibxcLDA/result.ref | 8 ++++---- tests/integrate/103_PW_15_CS_CF/result.ref | 10 +++++----- .../103_PW_15_CS_CF_bspline/result.ref | 10 +++++----- tests/integrate/103_PW_OU_CS_CF/result.ref | 10 +++++----- tests/integrate/107_PW_outWfcPw/result.ref | 18 +++++++++--------- tests/integrate/107_PW_outWfcR/result.ref | 16 ++++++++-------- tests/integrate/108_PW_RE/result.ref | 10 +++++----- tests/integrate/108_PW_RE_MG/result.ref | 8 ++++---- tests/integrate/109_PW_CR_fix_ac/result.ref | 10 +++++----- tests/integrate/109_PW_CR_fix_b/result.ref | 10 +++++----- tests/integrate/109_PW_CR_fix_bc/result.ref | 10 +++++----- tests/integrate/109_PW_CR_fix_c/result.ref | 10 +++++----- tests/integrate/109_PW_CR_moveatoms/result.ref | 10 +++++----- tests/integrate/110_PW_SY/result.ref | 10 +++++----- tests/integrate/112_PW_dipole/result.ref | 10 +++++----- tests/integrate/112_PW_efield/result.ref | 10 +++++----- tests/integrate/113_PW_gatefield/result.ref | 10 +++++----- tests/integrate/115_PW_sol_H2O/result.ref | 12 ++++++------ tests/integrate/120_PW_KP_MD_MSST/result.ref | 10 +++++----- .../120_PW_KP_MD_MSST_level2/result.ref | 10 +++++----- tests/integrate/120_PW_KP_MD_NPT/result.ref | 10 +++++----- tests/integrate/120_PW_KP_MD_NVT/result.ref | 10 +++++----- tests/integrate/170_PW_MD_1O/result.ref | 10 +++++----- tests/integrate/170_PW_MD_2O/result.ref | 10 +++++----- .../201_NO_KP_DJ_CF_CS_GaAs/result.ref | 10 +++++----- tests/integrate/204_NO_KP_AFM/result.ref | 10 +++++----- tests/integrate/270_NO_MD_1O/result.ref | 10 +++++----- 27 files changed, 141 insertions(+), 141 deletions(-) diff --git a/tests/integrate/101_PW_15_pseudopots_LibxcLDA/result.ref b/tests/integrate/101_PW_15_pseudopots_LibxcLDA/result.ref index d91f58b519..8b7c06f8d7 100644 --- a/tests/integrate/101_PW_15_pseudopots_LibxcLDA/result.ref +++ b/tests/integrate/101_PW_15_pseudopots_LibxcLDA/result.ref @@ -1,8 +1,8 @@ -etotref -196.4803276789654660 -etotperatomref -98.2401638395 +etotref -196.4803276808400483 +etotperatomref -98.2401638404 totalforceref 0.000000 -totalstressref 1409.543919 +totalstressref 1409.559594 pointgroupref T_d spacegroupref O_h nksibzref 1 -totaltimeref +0.11872 +totaltimeref 0.38 diff --git a/tests/integrate/103_PW_15_CS_CF/result.ref b/tests/integrate/103_PW_15_CS_CF/result.ref index b4a60eb90c..a404e1378c 100644 --- a/tests/integrate/103_PW_15_CS_CF/result.ref +++ b/tests/integrate/103_PW_15_CS_CF/result.ref @@ -1,8 +1,8 @@ -etotref -4206.8778432386488930 -etotperatomref -2103.4389216193 -totalforceref 8.104034 -totalstressref 81776.142537 +etotref -4206.8778432376702767 +etotperatomref -2103.4389216188 +totalforceref 8.101896 +totalstressref 81776.176460 pointgroupref C_1 spacegroupref C_1 nksibzref 8 -totaltimeref +1.32321 +totaltimeref 6.35 diff --git a/tests/integrate/103_PW_15_CS_CF_bspline/result.ref b/tests/integrate/103_PW_15_CS_CF_bspline/result.ref index 4f6f705ab2..d5836d68ac 100644 --- a/tests/integrate/103_PW_15_CS_CF_bspline/result.ref +++ b/tests/integrate/103_PW_15_CS_CF_bspline/result.ref @@ -1,8 +1,8 @@ -etotref -4200.2532789945953482 -etotperatomref -2100.1266394973 -totalforceref 10.586354 -totalstressref 82981.350483 +etotref -4200.2532789747046991 +etotperatomref -2100.1266394874 +totalforceref 10.537940 +totalstressref 82981.598052 pointgroupref C_1 spacegroupref C_1 nksibzref 2 -totaltimeref +0.45091 +totaltimeref 2.35 diff --git a/tests/integrate/103_PW_OU_CS_CF/result.ref b/tests/integrate/103_PW_OU_CS_CF/result.ref index 870e6ed958..ae273f866f 100644 --- a/tests/integrate/103_PW_OU_CS_CF/result.ref +++ b/tests/integrate/103_PW_OU_CS_CF/result.ref @@ -1,8 +1,8 @@ -etotref -1687.3372942255016369 -etotperatomref -843.6686471128 -totalforceref 4.988094 -totalstressref 34076.542350 +etotref -1687.3372941874322350 +etotperatomref -843.6686470937 +totalforceref 4.989240 +totalstressref 34076.373392 pointgroupref C_1 spacegroupref C_1 nksibzref 8 -totaltimeref +0.98146 +totaltimeref 4.71 diff --git a/tests/integrate/107_PW_outWfcPw/result.ref b/tests/integrate/107_PW_outWfcPw/result.ref index 5e24c222da..75e833385b 100644 --- a/tests/integrate/107_PW_outWfcPw/result.ref +++ b/tests/integrate/107_PW_outWfcPw/result.ref @@ -1,12 +1,12 @@ -etotref -197.1405644408571 -etotperatomref -98.5702822204 -Max_wfc_1 0.7386 -Max_wfc_2 0.3996 -Max_wfc_3 0.3811 -Max_wfc_4 0.3619 -Max_wfc_5 0.4156 -Max_wfc_6 0.3917 +etotref -197.1405644417894 +etotperatomref -98.5702822209 +Max_wfc_1 0.8113 +Max_wfc_2 0.4031 +Max_wfc_3 0.4073 +Max_wfc_4 0.3141 +Max_wfc_5 0.4158 +Max_wfc_6 0.4405 pointgroupref T_d spacegroupref O_h nksibzref 1 -totaltimeref +totaltimeref 0.37 diff --git a/tests/integrate/107_PW_outWfcR/result.ref b/tests/integrate/107_PW_outWfcR/result.ref index 00601d5cda..d8c86e8633 100644 --- a/tests/integrate/107_PW_outWfcR/result.ref +++ b/tests/integrate/107_PW_outWfcR/result.ref @@ -1,12 +1,12 @@ -etotref -197.1405644408571 -etotperatomref -98.5702822204 +etotref -197.1405644417894 +etotperatomref -98.5702822209 variance_wfc_r_0_0 0.31340 -variance_wfc_r_0_1 2.15785 -variance_wfc_r_0_2 2.17777 -variance_wfc_r_0_3 2.60806 -variance_wfc_r_0_4 0.78311 -variance_wfc_r_0_5 0.84843 +variance_wfc_r_0_1 1.71056 +variance_wfc_r_0_2 2.39604 +variance_wfc_r_0_3 1.66607 +variance_wfc_r_0_4 1.05190 +variance_wfc_r_0_5 1.29386 pointgroupref T_d spacegroupref O_h nksibzref 1 -totaltimeref 0.11908 +totaltimeref 0.42 diff --git a/tests/integrate/108_PW_RE/result.ref b/tests/integrate/108_PW_RE/result.ref index 6798abbece..ef666be794 100644 --- a/tests/integrate/108_PW_RE/result.ref +++ b/tests/integrate/108_PW_RE/result.ref @@ -1,5 +1,5 @@ -etotref -211.6098046542826410 -etotperatomref -105.8049023271 -totalforceref 11.823856 -totalstressref 941.983060 -totaltimeref +0.87209 +etotref -211.6098046531405146 +etotperatomref -105.8049023266 +totalforceref 11.824018 +totalstressref 941.977060 +totaltimeref 1.40 diff --git a/tests/integrate/108_PW_RE_MG/result.ref b/tests/integrate/108_PW_RE_MG/result.ref index 1042672674..badc240369 100644 --- a/tests/integrate/108_PW_RE_MG/result.ref +++ b/tests/integrate/108_PW_RE_MG/result.ref @@ -1,7 +1,7 @@ -etotref -210.7981143003123634 -etotperatomref -105.3990571502 -totalforceref 25.503030 +etotref -210.7981142976093167 +etotperatomref -105.3990571488 +totalforceref 25.502040 pointgroupref C_3v spacegroupref D_3d nksibzref 4 -totaltimeref 0.70155 +totaltimeref 1.10 diff --git a/tests/integrate/109_PW_CR_fix_ac/result.ref b/tests/integrate/109_PW_CR_fix_ac/result.ref index 0a66bc3f9b..bfbb04c2b9 100644 --- a/tests/integrate/109_PW_CR_fix_ac/result.ref +++ b/tests/integrate/109_PW_CR_fix_ac/result.ref @@ -1,5 +1,5 @@ -etotref -211.8190977232769967 -etotperatomref -105.9095488616 -totalforceref 0.065800 -totalstressref 360.021672 -totaltimeref +0.37137 +etotref -211.8190977198273686 +etotperatomref -105.9095488599 +totalforceref 0.066072 +totalstressref 360.004787 +totaltimeref 1.65 diff --git a/tests/integrate/109_PW_CR_fix_b/result.ref b/tests/integrate/109_PW_CR_fix_b/result.ref index 622d4c8e4c..46474ecddb 100644 --- a/tests/integrate/109_PW_CR_fix_b/result.ref +++ b/tests/integrate/109_PW_CR_fix_b/result.ref @@ -1,5 +1,5 @@ -etotref -211.8207812193545010 -etotperatomref -105.9103906097 -totalforceref 0.046652 -totalstressref 355.694939 -totaltimeref +0.35793 +etotref -211.8207812490062452 +etotperatomref -105.9103906245 +totalforceref 0.046818 +totalstressref 355.692789 +totaltimeref 1.69 diff --git a/tests/integrate/109_PW_CR_fix_bc/result.ref b/tests/integrate/109_PW_CR_fix_bc/result.ref index 0d2cc131d7..b3885d766b 100644 --- a/tests/integrate/109_PW_CR_fix_bc/result.ref +++ b/tests/integrate/109_PW_CR_fix_bc/result.ref @@ -1,5 +1,5 @@ -etotref -211.8190977077987611 -etotperatomref -105.9095488539 -totalforceref 0.065846 -totalstressref 360.020805 -totaltimeref +0.36126 +etotref -211.8190977227684471 +etotperatomref -105.9095488614 +totalforceref 0.066048 +totalstressref 360.009033 +totaltimeref 1.68 diff --git a/tests/integrate/109_PW_CR_fix_c/result.ref b/tests/integrate/109_PW_CR_fix_c/result.ref index bc2328fb19..101075b45c 100644 --- a/tests/integrate/109_PW_CR_fix_c/result.ref +++ b/tests/integrate/109_PW_CR_fix_c/result.ref @@ -1,5 +1,5 @@ -etotref -211.8207812361457343 -etotperatomref -105.9103906181 -totalforceref 0.046674 -totalstressref 355.690008 -totaltimeref +0.36365 +etotref -211.8207812442472857 +etotperatomref -105.9103906221 +totalforceref 0.046844 +totalstressref 355.678950 +totaltimeref 1.66 diff --git a/tests/integrate/109_PW_CR_moveatoms/result.ref b/tests/integrate/109_PW_CR_moveatoms/result.ref index 85ff3664bd..b5586a7d34 100644 --- a/tests/integrate/109_PW_CR_moveatoms/result.ref +++ b/tests/integrate/109_PW_CR_moveatoms/result.ref @@ -1,5 +1,5 @@ -etotref -211.8220727412902988 -etotperatomref -105.9110363706 -totalforceref 0.000018 -totalstressref 348.867117 -totaltimeref +0.31793 +etotref -211.8220727631340026 +etotperatomref -105.9110363816 +totalforceref 0.000000 +totalstressref 348.869952 +totaltimeref 0.93 diff --git a/tests/integrate/110_PW_SY/result.ref b/tests/integrate/110_PW_SY/result.ref index a5d2583afd..fe004bbaaa 100644 --- a/tests/integrate/110_PW_SY/result.ref +++ b/tests/integrate/110_PW_SY/result.ref @@ -1,8 +1,8 @@ -etotref -1713.7222698260661673 -etotperatomref -107.1076418641 -totalforceref 0.000324 -totalstressref 140.793438 +etotref -1713.7222698220132315 +etotperatomref -107.1076418639 +totalforceref 0.001974 +totalstressref 140.790489 pointgroupref T_d spacegroupref O_h nksibzref 3 -totaltimeref +4.93845 +totaltimeref 8.08 diff --git a/tests/integrate/112_PW_dipole/result.ref b/tests/integrate/112_PW_dipole/result.ref index d4b5184cb6..7f764fb0fe 100644 --- a/tests/integrate/112_PW_dipole/result.ref +++ b/tests/integrate/112_PW_dipole/result.ref @@ -1,5 +1,5 @@ -etotref -441.6269887886720085 -etotperatomref -147.2089962629 -totalforceref 7.601276 -totalstressref 2040.163610 -totaltimeref +0.33318 +etotref -441.6269887882587000 +etotperatomref -147.2089962628 +totalforceref 7.601019 +totalstressref 2040.168902 +totaltimeref 0.63 diff --git a/tests/integrate/112_PW_efield/result.ref b/tests/integrate/112_PW_efield/result.ref index 58a1f39874..e05a307a13 100644 --- a/tests/integrate/112_PW_efield/result.ref +++ b/tests/integrate/112_PW_efield/result.ref @@ -1,5 +1,5 @@ -etotref -441.7143153004764713 -etotperatomref -147.2381051002 -totalforceref 7.550200 -totalstressref 2052.354473 -totaltimeref +0.30156 +etotref -441.7143153067928552 +etotperatomref -147.2381051023 +totalforceref 7.548332 +totalstressref 2052.331678 +totaltimeref 0.62 diff --git a/tests/integrate/113_PW_gatefield/result.ref b/tests/integrate/113_PW_gatefield/result.ref index 5aed1ee411..8138df826b 100644 --- a/tests/integrate/113_PW_gatefield/result.ref +++ b/tests/integrate/113_PW_gatefield/result.ref @@ -1,5 +1,5 @@ -etotref -432.8463288659988280 -etotperatomref -144.2821096220 -totalforceref 20.870422 -totalstressref 2257.718361 -totaltimeref +0.49608 +etotref -432.8463288663492108 +etotperatomref -144.2821096221 +totalforceref 20.870256 +totalstressref 2257.704502 +totaltimeref 0.58 diff --git a/tests/integrate/115_PW_sol_H2O/result.ref b/tests/integrate/115_PW_sol_H2O/result.ref index ee7466b218..d11bd84ca6 100644 --- a/tests/integrate/115_PW_sol_H2O/result.ref +++ b/tests/integrate/115_PW_sol_H2O/result.ref @@ -1,6 +1,6 @@ -etotref -443.1237068197164604 -etotperatomref -147.7079022732 -totalforceref 16.865518 -esolelref -1.02249508526 -esolcavref +0.0325034001928 -totaltimeref 44.58741 +etotref -443.1237199062131253 +etotperatomref -147.7079066354 +totalforceref 16.858746 +esolelref -1.0225228008 +esolcavref 0.0325034862 +totaltimeref 53.61 diff --git a/tests/integrate/120_PW_KP_MD_MSST/result.ref b/tests/integrate/120_PW_KP_MD_MSST/result.ref index 42bc85bfff..13752bd260 100644 --- a/tests/integrate/120_PW_KP_MD_MSST/result.ref +++ b/tests/integrate/120_PW_KP_MD_MSST/result.ref @@ -1,5 +1,5 @@ -etotref -250.4498851597248 -etotperatomref -125.2249425799 -totalforceref 7.137234 -totalstressref 322.547564 -totaltimeref 5.60 +etotref -250.4498856072767 +etotperatomref -125.2249428036 +totalforceref 7.137474 +totalstressref 322.546737 +totaltimeref 10.48 diff --git a/tests/integrate/120_PW_KP_MD_MSST_level2/result.ref b/tests/integrate/120_PW_KP_MD_MSST_level2/result.ref index 1db0d6f115..d65862306e 100644 --- a/tests/integrate/120_PW_KP_MD_MSST_level2/result.ref +++ b/tests/integrate/120_PW_KP_MD_MSST_level2/result.ref @@ -1,5 +1,5 @@ -etotref -209.0846269525608 -etotperatomref -104.5423134763 -totalforceref 0.015086 -totalstressref 60.102536 -totaltimeref 5.07 +etotref -209.0846269675987 +etotperatomref -104.5423134838 +totalforceref 0.015008 +totalstressref 60.102361 +totaltimeref 11.67 diff --git a/tests/integrate/120_PW_KP_MD_NPT/result.ref b/tests/integrate/120_PW_KP_MD_NPT/result.ref index 5369dbdda3..98ff9fd55f 100644 --- a/tests/integrate/120_PW_KP_MD_NPT/result.ref +++ b/tests/integrate/120_PW_KP_MD_NPT/result.ref @@ -1,5 +1,5 @@ -etotref -211.8470586921492 -etotperatomref -105.9235293461 -totalforceref 1.994280 -totalstressref 446.726145 -totaltimeref 1.41 +etotref -211.8470631132377 +etotperatomref -105.9235315566 +totalforceref 1.994252 +totalstressref 446.654926 +totaltimeref 3.41 diff --git a/tests/integrate/120_PW_KP_MD_NVT/result.ref b/tests/integrate/120_PW_KP_MD_NVT/result.ref index 16161ef5cb..e433978802 100644 --- a/tests/integrate/120_PW_KP_MD_NVT/result.ref +++ b/tests/integrate/120_PW_KP_MD_NVT/result.ref @@ -1,5 +1,5 @@ -etotref -211.789138597789 -etotperatomref -105.8945692989 -totalforceref 1.818226 -totalstressref 419.846974 -totaltimeref 1.16 +etotref -211.7891383209277 +etotperatomref -105.8945691605 +totalforceref 1.818148 +totalstressref 419.850434 +totaltimeref 2.44 diff --git a/tests/integrate/170_PW_MD_1O/result.ref b/tests/integrate/170_PW_MD_1O/result.ref index dfe1770831..3dbc283de2 100644 --- a/tests/integrate/170_PW_MD_1O/result.ref +++ b/tests/integrate/170_PW_MD_1O/result.ref @@ -1,5 +1,5 @@ -etotref -211.7989766829926 -etotperatomref -105.8994883415 -totalforceref 0.673682 -totalstressref 387.295780 -totaltimeref +1.233 +etotref -211.7989760460549 +etotperatomref -105.8994880230 +totalforceref 0.674278 +totalstressref 387.284907 +totaltimeref 2.05 diff --git a/tests/integrate/170_PW_MD_2O/result.ref b/tests/integrate/170_PW_MD_2O/result.ref index be1016e1da..e561ba9a48 100644 --- a/tests/integrate/170_PW_MD_2O/result.ref +++ b/tests/integrate/170_PW_MD_2O/result.ref @@ -1,5 +1,5 @@ -etotref -211.7989767295324 -etotperatomref -105.8994883648 -totalforceref 0.673568 -totalstressref 387.288868 -totaltimeref +1.2656 +etotref -211.7989761048663 +etotperatomref -105.8994880524 +totalforceref 0.673900 +totalstressref 387.307779 +totaltimeref 1.94 diff --git a/tests/integrate/201_NO_KP_DJ_CF_CS_GaAs/result.ref b/tests/integrate/201_NO_KP_DJ_CF_CS_GaAs/result.ref index 012575f46e..f535918521 100644 --- a/tests/integrate/201_NO_KP_DJ_CF_CS_GaAs/result.ref +++ b/tests/integrate/201_NO_KP_DJ_CF_CS_GaAs/result.ref @@ -1,5 +1,5 @@ -etotref -4884.3328502754275178 -etotperatomref -2442.1664251377 -totalforceref 771.987060 -totalstressref 59946.597702 -totaltimeref +0.69368 +etotref -4884.332850270051 +etotperatomref -2442.1664251350 +totalforceref 771.983676 +totalstressref 59946.600686 +totaltimeref 1.72 diff --git a/tests/integrate/204_NO_KP_AFM/result.ref b/tests/integrate/204_NO_KP_AFM/result.ref index 2d8ddcf572..397570e39c 100644 --- a/tests/integrate/204_NO_KP_AFM/result.ref +++ b/tests/integrate/204_NO_KP_AFM/result.ref @@ -1,9 +1,9 @@ -etotref -6437.8376615092847715 -etotperatomref -3218.9188307546 +etotref -6437.837660182953 +etotperatomref -3218.9188300915 totalforceref 0.000000 -totalstressref 559.239585 -Compare_mulliken_pass 0 +totalstressref 559.242540 +Compare_mulliken_pass 1 pointgroupref O_h spacegroupref O_h nksibzref 4 -totaltimeref +11.96992 +totaltimeref 26.97 diff --git a/tests/integrate/270_NO_MD_1O/result.ref b/tests/integrate/270_NO_MD_1O/result.ref index 7e743f0d8d..8f398579b6 100644 --- a/tests/integrate/270_NO_MD_1O/result.ref +++ b/tests/integrate/270_NO_MD_1O/result.ref @@ -1,5 +1,5 @@ -etotref -211.2974370287481 -etotperatomref -105.6487185144 -totalforceref 3.467318 -totalstressref 555.120337 -totaltimeref +2.3953 +etotref -211.2974366541717 +etotperatomref -105.6487183271 +totalforceref 3.467360 +totalstressref 555.123731 +totaltimeref 5.00 From e75d70ed98baf36704850c0c7584d371fd7360ad Mon Sep 17 00:00:00 2001 From: Qianruipku Date: Thu, 23 Nov 2023 23:21:14 +0800 Subject: [PATCH 3/7] fix UTs compile --- source/module_io/test/for_testing_input_conv.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/source/module_io/test/for_testing_input_conv.h b/source/module_io/test/for_testing_input_conv.h index 0cdaf8b82d..f39953f7a1 100644 --- a/source/module_io/test/for_testing_input_conv.h +++ b/source/module_io/test/for_testing_input_conv.h @@ -282,7 +282,8 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in, const double& mixing_beta_in, const int& mixing_ndim_in, const double& mixing_gg0_in, - const bool& mixing_tau_in) + const bool& mixing_tau_in, + const double& mixing_beta_mag_in) { return; } From faf2d0918df0b84f8ff50afbc935d2ee295720ff Mon Sep 17 00:00:00 2001 From: Qianruipku Date: Fri, 24 Nov 2023 11:25:26 +0800 Subject: [PATCH 4/7] fix: mix_rho_real --- .../module_charge/charge_mixing.cpp | 28 ++- .../integrate/204_NO_KP_AFM/mulliken.txt.ref | 165 +++++++++--------- tests/integrate/204_NO_KP_AFM/result.ref | 2 +- 3 files changed, 108 insertions(+), 87 deletions(-) mode change 100644 => 100755 tests/integrate/204_NO_KP_AFM/result.ref diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index 44f1b55ed2..edcbc622fd 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -445,7 +445,9 @@ void Charge_Mixing::mix_rho_real(Charge* chr) if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4) { rhor_in = chr->rho_save[0]; - rhor_out = chr->rho[0]; + rhor_out = chr->rho[0]; + auto screen = std::bind(&Charge_Mixing::Kerker_screen_real, this, std::placeholders::_1); + this->mixing->push_data(this->rho_mdata, rhor_in, rhor_out, screen, true); } else if (GlobalV::NSPIN == 2) { @@ -453,9 +455,29 @@ void Charge_Mixing::mix_rho_real(Charge* chr) chr->get_rho_mag(); rhor_in = chr->rho_mag_save; rhor_out = chr->rho_mag; + const int nrxx = this->rhopw->nrxx; + auto screen = std::bind(&Charge_Mixing::Kerker_screen_real, this, std::placeholders::_1); + auto twobeta_mix + = [this, nrxx](double* out, const double* in, const double* sres) { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 256) +#endif + for (int i = 0; i < nrxx; ++i) + { + out[i] = in[i] + this->mixing_beta * sres[i]; + } + // magnetism +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 256) +#endif + for (int i = nrxx; i < 2 * nrxx; ++i) + { + out[i] = in[i] + this->mixing_beta_mag * sres[i]; + } + }; + this->mixing->push_data(this->rho_mdata, rhor_in, rhor_out, screen, twobeta_mix, true); } - auto screen = std::bind(&Charge_Mixing::Kerker_screen_real, this, std::placeholders::_1); - this->mixing->push_data(this->rho_mdata, rhor_in, rhor_out, screen, true); + auto inner_product = std::bind(&Charge_Mixing::inner_product_real, this, std::placeholders::_1, std::placeholders::_2); this->mixing->cal_coef(this->rho_mdata, inner_product); diff --git a/tests/integrate/204_NO_KP_AFM/mulliken.txt.ref b/tests/integrate/204_NO_KP_AFM/mulliken.txt.ref index af52d175a4..a9e10a5032 100644 --- a/tests/integrate/204_NO_KP_AFM/mulliken.txt.ref +++ b/tests/integrate/204_NO_KP_AFM/mulliken.txt.ref @@ -5,93 +5,92 @@ CALCULATE THE MULLIkEN ANALYSIS FOR EACH ATOM Total charge: 32 Decomposed Mulliken populations 0 Zeta of Fe Spin 1 Spin 2 Sum Diff -s 0 0.5144 0.3474 0.8618 0.167 - sum over m 0.5144 0.3474 0.8618 0.167 -s 1 0.927 0.9006 1.828 0.0264 - sum over m 0.927 0.9006 1.828 0.0264 -s 2 0.03713 0.01672 0.05384 0.02041 - sum over m 0.03713 0.01672 0.05384 0.02041 -s 3 -0.0003171 8.473e-06 -0.0003086 -0.0003256 - sum over m -0.0003171 8.473e-06 -0.0003086 -0.0003256 - sum over m+zeta 1.478 1.265 2.743 0.2135 -pz 0 0.9969 0.997 1.994 -3.411e-05 -px 0 0.9969 0.997 1.994 -8.355e-05 -py 0 0.9971 0.9971 1.994 -2.298e-05 - sum over m 2.991 2.991 5.982 -0.0001406 -pz 1 0.1918 0.1466 0.3384 0.04523 -px 1 0.0226 0.01029 0.03289 0.01231 -py 1 0.08327 0.06054 0.1438 0.02272 - sum over m 0.2977 0.2174 0.5151 0.08026 - sum over m+zeta 3.289 3.209 6.497 0.08011 -dz^2 0 0.9872 0.1572 1.144 0.83 -dxz 0 0.9322 0.4583 1.391 0.4739 -dyz 0 0.9206 0.6848 1.605 0.2358 -dx^2-y^2 0 0.988 0.4926 1.481 0.4953 -dxy 0 0.9347 0.1697 1.104 0.765 - sum over m 4.763 1.963 6.725 2.8 -dz^2 1 0.007172 0.008624 0.0158 -0.001452 -dxz 1 0.001558 0.002921 0.00448 -0.001363 -dyz 1 -0.002762 -0.001656 -0.004419 -0.001106 -dx^2-y^2 1 0.001388 0.0006564 0.002045 0.0007318 -dxy 1 0.001891 0.004545 0.006437 -0.002654 - sum over m 0.009248 0.01509 0.02434 -0.005843 - sum over m+zeta 4.772 1.978 6.75 2.794 -fz^3 0 0.0001958 6.926e-05 0.0002651 0.0001265 -fxz^2 0 5.282e-05 3.967e-05 9.249e-05 1.315e-05 -fyz^2 0 0.0003679 0.0002939 0.0006618 7.395e-05 -fzx^2-zy^2 0 0.003579 0.001847 0.005426 0.001732 -fxyz 0 0.0001303 0.0004518 0.000582 -0.0003215 -fx^3-3*xy^2 0 0.001085 0.000658 0.001743 0.0004271 -f3yx^2-y^3 0 0.001208 0.0007619 0.00197 0.0004459 - sum over m 0.006619 0.004121 0.01074 0.002497 - sum over m+zeta 0.006619 0.004121 0.01074 0.002497 +s 0 0.51 0.3574 0.8674 0.1526 + sum over m 0.51 0.3574 0.8674 0.1526 +s 1 0.9204 0.897 1.817 0.02333 + sum over m 0.9204 0.897 1.817 0.02333 +s 2 0.04825 0.02184 0.07009 0.0264 + sum over m 0.04825 0.02184 0.07009 0.0264 +s 3 -0.0004055 -5.52e-05 -0.0004607 -0.0003503 + sum over m -0.0004055 -5.52e-05 -0.0004607 -0.0003503 + sum over m+zeta 1.478 1.276 2.754 0.202 +pz 0 0.997 0.9971 1.994 -3.896e-05 +px 0 0.997 0.9971 1.994 -8.588e-05 +py 0 0.9972 0.9972 1.994 -2.241e-05 + sum over m 2.991 2.991 5.983 -0.0001473 +pz 1 0.1919 0.1459 0.3379 0.04601 +px 1 0.02289 0.01048 0.03337 0.01242 +py 1 0.08382 0.06127 0.1451 0.02255 + sum over m 0.2987 0.2177 0.5163 0.08097 + sum over m+zeta 3.29 3.209 6.499 0.08083 +dz^2 0 0.9859 0.1591 1.145 0.8268 +dxz 0 0.931 0.4532 1.384 0.4778 +dyz 0 0.9202 0.678 1.598 0.2422 +dx^2-y^2 0 0.9877 0.4979 1.486 0.4898 +dxy 0 0.9328 0.168 1.101 0.7648 + sum over m 4.758 1.956 6.714 2.801 +dz^2 1 0.008353 0.006887 0.01524 0.001466 +dxz 1 0.002451 0.00305 0.005501 -0.0005986 +dyz 1 -0.004139 -0.004893 -0.009032 0.0007534 +dx^2-y^2 1 0.00147 0.0006257 0.002096 0.0008448 +dxy 1 0.003618 0.004275 0.007893 -0.0006573 + sum over m 0.01175 0.009945 0.0217 0.001808 + sum over m+zeta 4.769 1.966 6.736 2.803 +fz^3 0 0.0001724 0.00011 0.0002824 6.232e-05 +fxz^2 0 5.322e-05 3.895e-05 9.217e-05 1.427e-05 +fyz^2 0 0.0003704 0.0002977 0.0006681 7.269e-05 +fzx^2-zy^2 0 0.003643 0.00188 0.005524 0.001763 +fxyz 0 0.0001469 0.0005171 0.000664 -0.0003702 +fx^3-3*xy^2 0 0.00111 0.0006851 0.001796 0.0004253 +f3yx^2-y^3 0 0.001232 0.000789 0.002021 0.0004427 + sum over m 0.006728 0.004318 0.01105 0.00241 + sum over m+zeta 0.006728 0.004318 0.01105 0.00241 Total Charge on atom: Fe 16 -Total Magnetism on atom: Fe 3.089 +Total Magnetism on atom: Fe 3.088 1 Zeta of Fe Spin 1 Spin 2 Sum Diff -s 0 0.3472 0.5142 0.8614 -0.167 - sum over m 0.3472 0.5142 0.8614 -0.167 -s 1 0.9006 0.927 1.828 -0.0264 - sum over m 0.9006 0.927 1.828 -0.0264 -s 2 0.01669 0.03713 0.05382 -0.02044 - sum over m 0.01669 0.03713 0.05382 -0.02044 -s 3 7.194e-06 -0.0003187 -0.0003115 0.0003259 - sum over m 7.194e-06 -0.0003187 -0.0003115 0.0003259 - sum over m+zeta 1.264 1.478 2.742 -0.2136 -pz 0 0.997 0.9969 1.994 3.345e-05 -px 0 0.997 0.9969 1.994 8.279e-05 -py 0 0.9971 0.9971 1.994 2.236e-05 - sum over m 2.991 2.991 5.982 0.0001386 -pz 1 0.1466 0.1918 0.3384 -0.04524 -px 1 0.01029 0.0226 0.03289 -0.01231 -py 1 0.06054 0.08327 0.1438 -0.02273 - sum over m 0.2174 0.2977 0.5151 -0.08028 - sum over m+zeta 3.209 3.289 6.497 -0.08014 -dz^2 0 0.157 0.9872 1.144 -0.8302 -dxz 0 0.4584 0.9322 1.391 -0.4738 -dyz 0 0.6849 0.9206 1.605 -0.2357 -dx^2-y^2 0 0.4925 0.988 1.48 -0.4955 -dxy 0 0.1697 0.9347 1.104 -0.765 - sum over m 1.962 4.763 6.725 -2.8 -dz^2 1 0.008577 0.007111 0.01569 0.001466 -dxz 1 0.002919 0.001562 0.004481 0.001358 -dyz 1 -0.001651 -0.002764 -0.004415 0.001113 -dx^2-y^2 1 0.0006373 0.001364 0.002001 -0.0007266 -dxy 1 0.004545 0.001894 0.006438 0.002651 - sum over m 0.01503 0.009166 0.02419 0.005861 - sum over m+zeta 1.977 4.772 6.749 -2.794 -fz^3 0 6.958e-05 0.0001962 0.0002657 -0.0001266 -fxz^2 0 3.967e-05 5.281e-05 9.248e-05 -1.313e-05 -fyz^2 0 0.0002939 0.0003679 0.0006617 -7.403e-05 -fzx^2-zy^2 0 0.001846 0.003578 0.005424 -0.001732 -fxyz 0 0.0004519 0.0001302 0.000582 0.0003217 -fx^3-3*xy^2 0 0.0006582 0.001085 0.001743 -0.0004268 -f3yx^2-y^3 0 0.0007621 0.001208 0.00197 -0.0004456 - sum over m 0.004121 0.006618 0.01074 -0.002496 - sum over m+zeta 0.004121 0.006618 0.01074 -0.002496 - +s 0 0.3573 0.5099 0.8672 -0.1526 + sum over m 0.3573 0.5099 0.8672 -0.1526 +s 1 0.897 0.9203 1.817 -0.02334 + sum over m 0.897 0.9203 1.817 -0.02334 +s 2 0.02181 0.04828 0.07009 -0.02646 + sum over m 0.02181 0.04828 0.07009 -0.02646 +s 3 -5.646e-05 -0.0004059 -0.0004624 0.0003494 + sum over m -5.646e-05 -0.0004059 -0.0004624 0.0003494 + sum over m+zeta 1.276 1.478 2.754 -0.202 +pz 0 0.9971 0.997 1.994 3.917e-05 +px 0 0.9971 0.997 1.994 8.514e-05 +py 0 0.9972 0.9972 1.994 2.265e-05 + sum over m 2.991 2.991 5.983 0.000147 +pz 1 0.1457 0.1923 0.338 -0.04664 +px 1 0.01048 0.0229 0.03337 -0.01242 +py 1 0.06098 0.08418 0.1452 -0.0232 + sum over m 0.2171 0.2994 0.5165 -0.08226 + sum over m+zeta 3.209 3.291 6.499 -0.08211 +dz^2 0 0.159 0.9859 1.145 -0.827 +dxz 0 0.4533 0.931 1.384 -0.4777 +dyz 0 0.6781 0.9202 1.598 -0.2421 +dx^2-y^2 0 0.4977 0.9877 1.485 -0.49 +dxy 0 0.1681 0.9328 1.101 -0.7648 + sum over m 1.956 4.758 6.714 -2.802 +dz^2 1 0.006832 0.008304 0.01514 -0.001472 +dxz 1 0.003046 0.002453 0.0055 0.0005928 +dyz 1 -0.004903 -0.004143 -0.009046 -0.0007601 +dx^2-y^2 1 0.0006069 0.001452 0.002058 -0.0008446 +dxy 1 0.004271 0.003618 0.007889 0.000653 + sum over m 0.009853 0.01168 0.02154 -0.001831 + sum over m+zeta 1.966 4.769 6.735 -2.803 +fz^3 0 0.0001101 0.0001728 0.0002829 -6.268e-05 +fxz^2 0 3.895e-05 5.321e-05 9.216e-05 -1.427e-05 +fyz^2 0 0.0002961 0.0003722 0.0006683 -7.611e-05 +fzx^2-zy^2 0 0.001878 0.003645 0.005523 -0.001767 +fxyz 0 0.0005171 0.0001469 0.000664 0.0003702 +fx^3-3*xy^2 0 0.0006852 0.00111 0.001796 -0.0004253 +f3yx^2-y^3 0 0.0007887 0.001232 0.002021 -0.0004433 + sum over m 0.004314 0.006732 0.01105 -0.002418 + sum over m+zeta 0.004314 0.006732 0.01105 -0.002418 Total Charge on atom: Fe 16 -Total Magnetism on atom: Fe -3.089 +Total Magnetism on atom: Fe -3.09 diff --git a/tests/integrate/204_NO_KP_AFM/result.ref b/tests/integrate/204_NO_KP_AFM/result.ref old mode 100644 new mode 100755 index 397570e39c..98b8096c53 --- a/tests/integrate/204_NO_KP_AFM/result.ref +++ b/tests/integrate/204_NO_KP_AFM/result.ref @@ -2,7 +2,7 @@ etotref -6437.837660182953 etotperatomref -3218.9188300915 totalforceref 0.000000 totalstressref 559.242540 -Compare_mulliken_pass 1 +Compare_mulliken_pass 0 pointgroupref O_h spacegroupref O_h nksibzref 4 From 5059bcded7f3d8440455c8a9d7e3fd8cf934ec98 Mon Sep 17 00:00:00 2001 From: Qianruipku Date: Fri, 24 Nov 2023 23:27:15 +0800 Subject: [PATCH 5/7] add some time for Ctest since Autotest nearly cost 1500s --- .github/workflows/test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 58bd6ea54f..419b080ff1 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -25,4 +25,4 @@ jobs: GTEST_COLOR: 'yes' OMP_NUM_THREADS: '2' run: | - cmake --build build --target test ARGS="-V" + cmake --build build --target test ARGS="-V --timeout 1600" From 4785ee3b23efd6271c7ea55f35dd27691c28cce6 Mon Sep 17 00:00:00 2001 From: Qianruipku Date: Sat, 25 Nov 2023 11:11:12 +0800 Subject: [PATCH 6/7] add explanation for parameters --- .../module_charge/charge_mixing.h | 41 +++++++++++++------ 1 file changed, 29 insertions(+), 12 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index f941a7fe5e..f6492e7d29 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -12,10 +12,10 @@ class Charge_Mixing public: Charge_Mixing(); ~Charge_Mixing(); - Base_Mixing::Mixing* mixing = nullptr; - Base_Mixing::Mixing_Data rho_mdata; - Base_Mixing::Mixing_Data tau_mdata; - Base_Mixing::Mixing_Data nhat_mdata; + Base_Mixing::Mixing* mixing = nullptr; ///< Mixing object to mix charge density, kinetic energy density and compensation density + Base_Mixing::Mixing_Data rho_mdata; ///< Mixing data for charge density + Base_Mixing::Mixing_Data tau_mdata; ///< Mixing data for kinetic energy density + Base_Mixing::Mixing_Data nhat_mdata; ///< Mixing data for compensation density Base_Mixing::Plain_Mixing* mixing_highf = nullptr; ///< The high_frequency part is mixed by plain mixing method. @@ -80,6 +80,7 @@ class Charge_Mixing * @param mixing_ndim_in mixing ndim * @param mixing_gg0_in mixing gg0 for Kerker screen * @param mixing_tau_in whether to use tau mixing + * @param mixing_beta_mag_in mixing beta for magnetism */ void set_mixing(const std::string& mixing_mode_in, const double& mixing_beta_in, @@ -107,6 +108,13 @@ class Charge_Mixing double get_drho(Charge* chr, const double nelec); // init pwrho and rhodpw + + /** + * @brief Set the smooth and dense grids + * + * @param rhopw_in smooth grid + * @param rhodpw_in dense grid when double grid is used, otherwise same as rhopw + */ void set_rhopw(ModulePW::PW_Basis* rhopw_in, ModulePW::PW_Basis* rhodpw_in); // extracting parameters @@ -134,18 +142,18 @@ class Charge_Mixing //====================================== // General parameters //====================================== - std::string mixing_mode = "broyden"; - double mixing_beta = 0.8; - double mixing_beta_mag = 1.6; - int mixing_ndim = 8; - double mixing_gg0 = 0.0; // mohan add 2014-09-27 - bool mixing_tau = false; + std::string mixing_mode = "broyden"; ///< mixing mode: "plain", "broyden", "pulay" + double mixing_beta = 0.8; ///< mixing beta for density + double mixing_beta_mag = 1.6; ///< mixing beta for magnetism + int mixing_ndim = 8; ///< mixing ndim for broyden and pulay + double mixing_gg0 = 0.0; ///< mixing gg0 for Kerker screen + bool mixing_tau = false; ///< whether to use tau mixing bool new_e_iteration = true; ModulePW::PW_Basis* rhopw = nullptr; ///< smooth grid ModulePW::PW_Basis* rhodpw = nullptr; ///< dense grid, same as rhopw for ncpp. - bool autoset = false; + // bool autoset = false; private: double rhog_dot_product(const std::complex* const* const rhog1, @@ -153,16 +161,25 @@ class Charge_Mixing /** * @brief divide rho/tau to smooth and high frequency parts + * @param data_d dense data + * @param data_s smooth data + * @param data_hf high frequency data = dense data - smooth data * */ void divide_data(std::complex* data_d, std::complex*& data_s, std::complex*& data_hf); /** * @brief gather smooth and high frequency parts to rho/tau - * + * @param data_d dense data + * @param data_s smooth data + * @param data_hf high frequency data = dense data - smooth data + * */ void combine_data(std::complex* data_d, std::complex*& data_s, std::complex*& data_hf); /** * @brief clean smooth and high frequency parts + * @param data_d dense data + * @param data_s smooth data + * @param data_hf high frequency data = dense data - smooth data * */ void clean_data(std::complex*& data_s, std::complex*& data_hf); From 2c356e27d3dec461b00a33465bd5658a421b7d21 Mon Sep 17 00:00:00 2001 From: Qianruipku Date: Sat, 25 Nov 2023 12:03:39 +0800 Subject: [PATCH 7/7] delete default beta_mag in set_mixing --- .../module_charge/charge_mixing.h | 2 +- .../module_elecstate/test/charge_mixing_test.cpp | 16 ++++++++-------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index f6492e7d29..993c229057 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -87,7 +87,7 @@ class Charge_Mixing const int& mixing_ndim_in, const double& mixing_gg0_in, const bool& mixing_tau_in, - const double& mixing_beta_mag_in = 1.6); // mohan add mixing_gg0_in 2014-09-27 + const double& mixing_beta_mag_in); // mohan add mixing_gg0_in 2014-09-27 // /** // * @brief use auto set diff --git a/source/module_elecstate/test/charge_mixing_test.cpp b/source/module_elecstate/test/charge_mixing_test.cpp index 794b4c7345..fb20b2519a 100644 --- a/source/module_elecstate/test/charge_mixing_test.cpp +++ b/source/module_elecstate/test/charge_mixing_test.cpp @@ -282,12 +282,12 @@ TEST_F(ChargeMixingTest, SetMixingTest) bool mixingtau = false; GlobalV::SCF_THR_TYPE = 1; std::string mode = "broyden"; - CMtest.set_mixing(mode, beta, dim, gg0, mixingtau); + CMtest.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6); EXPECT_EQ(CMtest.rho_mdata.length, pw_basis.npw); GlobalV::SCF_THR_TYPE = 2; mode = "broyden"; - CMtest.set_mixing(mode, beta, dim, gg0, mixingtau); + CMtest.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6); EXPECT_EQ(CMtest.rho_mdata.length, pw_basis.nrxx); EXPECT_EQ(CMtest.get_mixing_mode(), "broyden"); EXPECT_EQ(CMtest.get_mixing_beta(), 1.0); @@ -298,19 +298,19 @@ TEST_F(ChargeMixingTest, SetMixingTest) mixingtau = true; mode = "plain"; GlobalV::SCF_THR_TYPE = 1; - CMtest.set_mixing(mode, beta, dim, gg0, mixingtau); + CMtest.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6); CMtest.mix_reset(); EXPECT_EQ(CMtest.tau_mdata.length, pw_basis.npw); GlobalV::SCF_THR_TYPE = 2; - CMtest.set_mixing(mode, beta, dim, gg0, mixingtau); + CMtest.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6); CMtest.mix_reset(); EXPECT_EQ(CMtest.tau_mdata.length, pw_basis.nrxx); mode = "nothing"; std::string output; testing::internal::CaptureStdout(); - EXPECT_EXIT(CMtest.set_mixing(mode, beta, dim, gg0, mixingtau);, ::testing::ExitedWithCode(0), ""); + EXPECT_EXIT(CMtest.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6);, ::testing::ExitedWithCode(0), ""); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("This Mixing mode is not implemended yet,coming soon.")); } @@ -632,7 +632,7 @@ TEST_F(ChargeMixingTest, MixRhoTest) Charge_Mixing CMtest_recip; CMtest_recip.set_rhopw(&pw_basis, &pw_basis); GlobalV::SCF_THR_TYPE = 1; - CMtest_recip.set_mixing(mode, beta, dim, gg0, mixingtau); + CMtest_recip.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6); CMtest_recip.mix_reset(); for(int i = 0 ; i < nspin * npw; ++i) { @@ -662,7 +662,7 @@ TEST_F(ChargeMixingTest, MixRhoTest) Charge_Mixing CMtest_real; GlobalV::SCF_THR_TYPE = 2; CMtest_real.set_rhopw(&pw_basis, &pw_basis); - CMtest_real.set_mixing(mode, beta, dim, gg0, mixingtau); + CMtest_real.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6); CMtest_real.mix_reset(); for(int i = 0 ; i < nspin * nrxx; ++i) { @@ -748,7 +748,7 @@ TEST_F(ChargeMixingTest, MixDoubleGridRhoTest) Charge_Mixing CMtest_recip; CMtest_recip.set_rhopw(&pw_basis, &pw_dbasis); GlobalV::SCF_THR_TYPE = 1; - CMtest_recip.set_mixing(mode, beta, dim, gg0, mixingtau); + CMtest_recip.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6); CMtest_recip.mix_reset(); for (int i = 0; i < nspin * npw; ++i) {