From 54fb9cfad394ec9e0949583376c597a0cfb9b460 Mon Sep 17 00:00:00 2001 From: Qianruipku Date: Sat, 16 Sep 2023 11:25:26 +0800 Subject: [PATCH] Refactor: add Mixing module temporarily; abacus_pw can be compiled now --- source/Makefile.Objects | 5 +- source/module_base/CMakeLists.txt | 1 + source/module_base/lapack_connector.h | 2 +- .../module_mixing/broyden_mixing.h | 260 ++++++ source/module_base/module_mixing/mixing.h | 140 ++++ .../module_base/module_mixing/mixing_data.cpp | 31 + .../module_base/module_mixing/mixing_data.h | 96 +++ .../module_base/module_mixing/plain_mixing.h | 87 ++ source/module_elecstate/CMakeLists.txt | 2 - .../module_elecstate/module_charge/charge.cpp | 36 +- .../module_elecstate/module_charge/charge.h | 6 + .../module_charge/charge_broyden.cpp | 181 ---- .../module_charge/charge_mixing.cpp | 778 +++++++++--------- .../module_charge/charge_mixing.h | 156 +--- .../module_charge/charge_pulay.cpp | 696 ---------------- source/module_esolver/esolver_ks_lcao.cpp | 6 - source/module_esolver/esolver_ks_pw.cpp | 4 - 17 files changed, 1073 insertions(+), 1414 deletions(-) create mode 100644 source/module_base/module_mixing/broyden_mixing.h create mode 100644 source/module_base/module_mixing/mixing.h create mode 100644 source/module_base/module_mixing/mixing_data.cpp create mode 100644 source/module_base/module_mixing/mixing_data.h create mode 100644 source/module_base/module_mixing/plain_mixing.h delete mode 100644 source/module_elecstate/module_charge/charge_broyden.cpp delete mode 100644 source/module_elecstate/module_charge/charge_pulay.cpp diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 8e11febe5e..5054394b70 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -22,6 +22,7 @@ VPATH=./src_global:\ ./module_base/module_container/base/core:\ ./module_base/module_container/ATen/core:\ ./module_base/module_container/ATen/kernels:\ +./module_base/module_mixing:\ ./module_md:\ ./module_basis/module_pw:\ ./module_esolver:\ @@ -137,6 +138,8 @@ OBJS_BASE=abfs-vector3_order.o\ formatter_contextfmt.o\ cubic_spline.o\ spherical_bessel_transformer.o\ + mixing_data.o\ + OBJS_CELL=atom_pseudo.o\ atom_spec.o\ @@ -472,10 +475,8 @@ OBJS_SRCPW=H_Ewald_pw.o\ charge.o\ charge_init.o\ charge_mpi.o\ - charge_broyden.o\ charge_extra.o\ charge_mixing.o\ - charge_pulay.o\ fp_energy.o\ forces.o\ force_op.o\ diff --git a/source/module_base/CMakeLists.txt b/source/module_base/CMakeLists.txt index f3aa16ec2e..795ff0e771 100644 --- a/source/module_base/CMakeLists.txt +++ b/source/module_base/CMakeLists.txt @@ -52,6 +52,7 @@ add_library( formatter_physfmt.cpp formatter_table.cpp formatter_contextfmt.cpp + module_mixing/mixing_data.cpp ${LIBM_SRC} ) diff --git a/source/module_base/lapack_connector.h b/source/module_base/lapack_connector.h index 29d70e8e34..5865e178be 100644 --- a/source/module_base/lapack_connector.h +++ b/source/module_base/lapack_connector.h @@ -94,7 +94,7 @@ extern "C" // dsytrf_ computes the Bunch-Kaufman factorization of a double precision // symmetric matrix, while dsytri takes its output to perform martrix inversion void dsytrf_(const char* uplo, const int* n, double * a, const int* lda, - int *ipiv,double *work, int* lwork ,int *info); + int *ipiv,double *work, const int* lwork ,int *info); void dsytri_(const char* uplo,const int* n,double *a, const int *lda, int *ipiv, double * work,int *info); // Peize Lin add dsptrf and dsptri 2016-06-21, to compute inverse real symmetry indefinit matrix. diff --git a/source/module_base/module_mixing/broyden_mixing.h b/source/module_base/module_mixing/broyden_mixing.h new file mode 100644 index 0000000000..46d6979993 --- /dev/null +++ b/source/module_base/module_mixing/broyden_mixing.h @@ -0,0 +1,260 @@ +#ifndef BROYDEN_MIXING_H_ +#define BROYDEN_MIXING_H_ +#include "mixing.h" +#include "module_base/lapack_connector.h" +#include "module_base/matrix.h" +#include "module_base/memory.h" +#include "module_base/timer.h" +#include "module_base/tool_title.h" + +namespace Base_Mixing +{ +/** + * @brief Simplified modified broyden_mixing method. + * Ref: D.D. Johnson PRB 38, 12807 (1988) + * Here the weight w0 of the error of the inverse Jacobian is set to 0 and the weight wn of + * the error of each previous iteration is set to same. + * @note Formula: + * F = n_out - n_in + * dF{i} = F_{i-1} - F{i} //different from Ref + * dn_in{i} = n_in_{i-1} - n_in{i} //different from Ref + * alpha{ij} = + * beta{ij} = inv(alpha){ij} + * c{mk} = + * gamma{mn} = \sum_k c{mk} * beta{kn} + * n{m+1} = n_in{m} + mixing_beta*F{m} - \sum_n gamma{mn} * (dn_in{n} + mixing_beta*dF{n}) + * mixing_data{i} = n_in{i} + mixing_beta*F{i} + * n{m+1} = \sum_i coef{i} * mixing_data{i} + */ +class Broyden_Mixing : public Mixing +{ + public: + Broyden_Mixing(const int& mixing_ndim) + { + this->mixing_ndim = mixing_ndim; + this->coef = std::vector(mixing_ndim); + this->beta = ModuleBase::matrix(mixing_ndim, mixing_ndim, true); + } + virtual ~Broyden_Mixing() override + { + if (F != nullptr) + free(F); + if (dF != nullptr) + free(dF); + }; + virtual void push_data(Mixing_Data& mdata, + const double* data_in, + const double* data_out, + std::function screen, + const bool& need_calcoef) override + { + this->tem_push_data(mdata, data_in, data_out, screen, need_calcoef); + }; + virtual void push_data(Mixing_Data& mdata, + const std::complex* data_in, + const std::complex* data_out, + std::function*)> screen, + const bool& need_calcoef) override + { + this->tem_push_data(mdata, data_in, data_out, screen, need_calcoef); + }; + virtual void cal_coef(const Mixing_Data& mdata, std::function inner_dot) override + { + tem_cal_coef(mdata, inner_dot); + } + virtual void cal_coef(const Mixing_Data& mdata, + std::function*, std::complex*)> inner_dot) override + { + tem_cal_coef(mdata, inner_dot); + } + + private: + template + void tem_push_data(Mixing_Data& mdata, + 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, 128) +#endif + for (int i = 0; i < length; ++i) + { + F_tmp[i] = data_out[i] - data_in[i]; + } + + // get screened F + if (screen != nullptr) + screen(F_tmp.data()); + + // container::Tensor data = data_in + mixing_beta * F; + std::vector data(length); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 128) +#endif + for (int i = 0; i < length; ++i) + { + data[i] = data_in[i] + this->mixing_beta * F_tmp[i]; + } + + mdata.push(data.data()); + + if (!need_calcoef) + return; + + if (address != &mdata && address != nullptr) + ModuleBase::WARNING_QUIT( + "Broyden_Mixing", + "One Broyden_Mixing object can only bind one Mixing_Data object to calculate coefficients"); + + FPTYPE* FP_dF = static_cast(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, 128) +#endif + for (int i = 0; i < length; ++i) + { + FP_F[i] = F_tmp[i]; + } + } + else + { + const int previous = mdata.index_move(-1); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 128) +#endif + for (int i = 0; i < length; ++i) + { + FP_F[i] = F_tmp[i]; + // dF{n} = F{n-1} - F{n} = -(F{n} - F{n-1}) + FP_dF[previous * length + i] -= FP_F[i]; + } + } + }; + + template + void tem_cal_coef(const Mixing_Data& mdata, std::function inner_dot) + { + ModuleBase::TITLE("Charge_Mixing", "Simplified_Broyden_mixing"); + ModuleBase::timer::tick("Charge", "Broyden_mixing"); + if (address != &mdata && address != nullptr) + ModuleBase::WARNING_QUIT( + "Mixing", + "One Mixing object can only bind one Mixing_Data object to calculate coefficients"); + const int length = mdata.length; + const int start = mdata.start; + FPTYPE* FP_dF = static_cast(dF); + FPTYPE* FP_F = static_cast(F); + if (mdata.ndim_previous > 0) + { + const int ndim_previous = mdata.ndim_previous; + ModuleBase::matrix beta_tmp(ndim_previous, ndim_previous); + //beta(i, j) = + for (int i = 0; i < ndim_previous; ++i) + { + FPTYPE* dFi = FP_dF + i * length; + for (int j = i; j < ndim_previous; ++j) + { + if (i < ndim_previous - 1 && j < ndim_previous - 1) + { + beta_tmp(i, j) = beta(i, j); + } + FPTYPE* dFj = FP_dF + j * length; + beta(i, j) = beta_tmp(i, j) = inner_dot(dFi, dFj); + if (j != i) + { + beta(j, i) = beta_tmp(j, i) = beta_tmp(i, j); + } + } + } + double* work = new double[ndim_previous]; + int* iwork = new int[ndim_previous]; + char uu = 'U'; + int info; + dsytrf_(&uu, &ndim_previous, beta_tmp.c, &ndim_previous, iwork, work, &ndim_previous, &info); + if (info != 0) + ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when factorizing beta."); + dsytri_(&uu, &ndim_previous, beta_tmp.c, &ndim_previous, iwork, work, &info); + if (info != 0) + ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when DSYTRI beta."); + for (int i = 0; i < ndim_previous; ++i) + { + for (int j = i + 1; j < ndim_previous; ++j) + { + beta_tmp(i, j) = beta_tmp(j, i); + } + } + for (int i = 0; i < ndim_previous; ++i) + { + FPTYPE* dFi = FP_dF + i * length; + work[i] = inner_dot(dFi, FP_F); + } + //gamma[i] = \sum_j beta_tmp(i,j) * work[j] + std::vector gamma(ndim_previous); + container::BlasConnector::gemv('N', + ndim_previous, + ndim_previous, + 1.0, + beta_tmp.c, + ndim_previous, + work, + 1, + 0.0, + gamma.data(), + 1); + + coef[start] = 1 + gamma[mdata.index_move(-1)]; + for (int i = 1; i < ndim_previous - 1; ++i) + { + coef[mdata.index_move(-i)] = gamma[mdata.index_move(-i - 1)] - gamma[mdata.index_move(-i)]; + } + coef[mdata.index_move(-ndim_previous)] = -gamma[mdata.index_move(-ndim_previous)]; + + delete[] work; + delete[] iwork; + } + else + { + coef[0] = 1.0; + } + + FPTYPE* dFstart = FP_dF + start * length; +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 128) +#endif + for (int i = 0; i < length; ++i) + { + dFstart[i] = FP_F[i]; + } + ModuleBase::timer::tick("Charge", "Broyden_mixing"); + }; + + private: + // F = data_out - data_in + void* F = nullptr; + // dF = F_{n+1} - F_n + void* dF = nullptr; + // binded mixing_data + Mixing_Data* address = nullptr; + // beta_ij = + ModuleBase::matrix beta; +}; +} // namespace Base_Mixing +#endif \ No newline at end of file diff --git a/source/module_base/module_mixing/mixing.h b/source/module_base/module_mixing/mixing.h new file mode 100644 index 0000000000..ed94b40dfe --- /dev/null +++ b/source/module_base/module_mixing/mixing.h @@ -0,0 +1,140 @@ +#ifndef MIXING_H_ +#define MIXING_H_ +#include + +#include "mixing_data.h" +#include "module_base/module_container/base/third_party/blas.h" +namespace Base_Mixing +{ + +/** + * @brief Mixing class can mixing different steps of data to solver the iteration problem. + * For equation x = f(x), we can use iterative process to solve it: + * x_{n+1} = f(x_n), n = 0, 1, 2, ... + * To acceralte the convergence, we can use mixing method, we need information including + * x_in, x_out=f(x_in) in each iterative step. + * + */ +class Mixing +{ + public: + Mixing(){}; + virtual ~Mixing(){}; + + /** + * @brief init mixing data + * + * @param mdata mixing data + * @param length the length of each vector + * @param type_size size of type + * + */ + virtual void init_mixing_data(Mixing_Data& mdata, const int& length, const size_t& type_size) + { + mdata.resize(mixing_ndim, length, type_size); + } + + /** + * @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) + = 0; + 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; + /** + * @brief calculate + * + * @tparam T class type + * @param inner_dot pointer to the inner dot function + */ + virtual void cal_coef(const Mixing_Data& mdata, std::function inner_dot) = 0; + virtual void cal_coef(const Mixing_Data& mdata, + std::function*, std::complex*)> inner_dot) + = 0; + + /** + * @brief calculate the mixing data + * + * @param mdata Mixing_Data + * @param data_mix output data + */ + void mix_data(const Mixing_Data& mdata, double* data_mix) + { + double* FP_data = static_cast(mdata.data); + if(mdata.ndim_use == 1) + { + 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) + { + std::complex* FP_data = static_cast*>(mdata.data); + if(mdata.ndim_use == 1) + { + 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); + } + + /** + * @brief reset mixing + * + */ + void reset(); + + public: + // mixing_beta from INPUT + double mixing_beta = 0.7; + // coeficients for mixing + std::vector coef; + // ndim for mixing + int mixing_ndim = 1; +}; + +} // namespace Base_Mixing + +#endif \ No newline at end of file diff --git a/source/module_base/module_mixing/mixing_data.cpp b/source/module_base/module_mixing/mixing_data.cpp new file mode 100644 index 0000000000..6e329b9782 --- /dev/null +++ b/source/module_base/module_mixing/mixing_data.cpp @@ -0,0 +1,31 @@ +#include "mixing_data.h" + +namespace Base_Mixing +{ + +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); +} + +Mixing_Data::~Mixing_Data() +{ + if(this->data != nullptr) + free(this->data); +} + +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) + free(this->data); + this->data = malloc(ndim * length * type_size); + this->start = -1; + this->ndim_use = 0; +} + + +} diff --git a/source/module_base/module_mixing/mixing_data.h b/source/module_base/module_mixing/mixing_data.h new file mode 100644 index 0000000000..9b2435ac89 --- /dev/null +++ b/source/module_base/module_mixing/mixing_data.h @@ -0,0 +1,96 @@ +#ifndef MIXING_DATA_H_ +#define MIXING_DATA_H_ +#include + +#include "module_base/module_container/ATen/tensor.h" +namespace Base_Mixing +{ + +/** + * @brief data for Mixing class + * + */ +class Mixing_Data +{ + public: + Mixing_Data() = default; + /** + * @brief Construct a new Mixing_Data object + * + * @param ndim store ndim vectors for mixing + * @param length the length of each vector + * @param type_size size of type + * + */ + Mixing_Data(const int& ndim, const int& length, const size_t& type_size); + + /** + * @brief Destroy the Mixing_Data object + * + */ + ~Mixing_Data(); + + /** + * @brief resize the data + * + * @param ndim store ndim vectors for mixing + * @param length the length of each vector + * @param type_size size of type + * + */ + void resize(const int& ndim, const int& length, const size_t& type_size); + + /** + * @brief push data to the tensor + * + */ + template + void push(const FPTYPE* data_in) + { + this->start = (this->start + 1) % this->ndim_tot; + this->ndim_use = std::min(this->ndim_use + 1, this->ndim_tot); + this->ndim_previous = std::min(this->ndim_previous + 1, this->ndim_tot); + FPTYPE* FP_data = static_cast(this->data); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 128) +#endif + for (int i = 0; i < this->length; ++i) + FP_data[this->start * this->length + i] = data_in[i]; + } + + /** + * @brief reset mixing + * + */ + void reset_mixing() + { + this->ndim_use = 0; + this->start = -1; + } + + /** + * @brief get the index of i-th vector + * + */ + int index_move(const int& n) const + { + return (n + this->start + ndim_tot) % ndim_tot; + } + + public: + // Tensor pointer to store the data + void* data = nullptr; + // the number of vectors for mixing + int ndim_tot = 0; + // the length of each vector + int length = 0; + // the start index for vector: start = this->index_move(0) + int start = -1; + // the number of used vectors for mixing + int ndim_use = 0; + // the number of previous vectors for mixing + int ndim_previous = -1; +}; + +} // namespace Base_Mixing +#endif \ 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 new file mode 100644 index 0000000000..6104b4d78b --- /dev/null +++ b/source/module_base/module_mixing/plain_mixing.h @@ -0,0 +1,87 @@ +#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" + +namespace Base_Mixing +{ +/** + * @brief Plain mixing : rho_new = rho_in + mixing_beta * (rho_out - rho_in) + * + */ +class Plain_Mixing : public Mixing +{ + public: + Plain_Mixing() + { + this->mixing_ndim = 1; + this->coef = std::vector(1, 1.0); + + } + virtual ~Plain_Mixing() override{}; + virtual void push_data(Mixing_Data& mdata, + const double* data_in, + const double* data_out, + std::function screen, + const bool& need_calcoef) override + { + this->tem_push_data(mdata, data_in, data_out, screen, need_calcoef); + }; + virtual void push_data(Mixing_Data& mdata, + const std::complex* data_in, + const std::complex* data_out, + std::function*)> screen, + const bool& need_calcoef) override + { + this->tem_push_data(mdata, data_in, data_out, screen, need_calcoef); + }; + virtual void cal_coef(const Mixing_Data& mdata, std::function inner_dot) override + { + return; + } + virtual void cal_coef(const Mixing_Data& mdata, + std::function*, std::complex*)> inner_dot) override + { + return; + } + + private: + template + void tem_push_data(Mixing_Data& mdata, + 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, 128) +#endif + for (int i = 0; i < length; ++i) + { + F_tmp[i] = data_out[i] - data_in[i]; + } + + // get screened F + if (screen != nullptr) + screen(F_tmp.data()); + + // container::Tensor data = data_in + mixing_beta * F; + std::vector data(length); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 128) +#endif + for (int i = 0; i < length; ++i) + { + data[i] = data_in[i] + this->mixing_beta * F_tmp[i]; + } + + mdata.push(data.data()); + }; +}; +} // namespace Base_Mixing +#endif \ No newline at end of file diff --git a/source/module_elecstate/CMakeLists.txt b/source/module_elecstate/CMakeLists.txt index 82f96b581b..3738f99f8d 100644 --- a/source/module_elecstate/CMakeLists.txt +++ b/source/module_elecstate/CMakeLists.txt @@ -18,10 +18,8 @@ list(APPEND objects module_charge/charge.cpp module_charge/charge_init.cpp module_charge/charge_mpi.cpp - module_charge/charge_broyden.cpp module_charge/charge_extra.cpp module_charge/charge_mixing.cpp - module_charge/charge_pulay.cpp module_charge/symmetry_rho.cpp module_charge/symmetry_rhog.cpp fp_energy.cpp diff --git a/source/module_elecstate/module_charge/charge.cpp b/source/module_elecstate/module_charge/charge.cpp index ec4ffeb546..6db6ac8f81 100644 --- a/source/module_elecstate/module_charge/charge.cpp +++ b/source/module_elecstate/module_charge/charge.cpp @@ -62,15 +62,6 @@ void Charge::destroy() { for (int i = 0; i < GlobalV::NSPIN; i++) { - delete[] rho[i]; - delete[] rhog[i]; - delete[] rho_save[i]; - delete[] rhog_save[i]; - if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5) - { - delete[] kin_r[i]; - delete[] kin_r_save[i]; - } if(GlobalV::use_paw) { delete[] nhat[i]; @@ -83,6 +74,12 @@ void Charge::destroy() delete[] rhog_save; delete[] rho_core; delete[] rhog_core; + delete[] _space_rho; + delete[] _space_rho_save; + delete[] _space_rhog; + delete[] _space_rhog_save; + delete[] _space_kin_r; + delete[] _space_kin_r_save; if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5) { delete[] kin_r; @@ -120,6 +117,15 @@ void Charge::allocate(const int& nspin_in) } // allocate memory + _space_rho = new double[nspin * nrxx]; + _space_rho_save = new double[nspin * nrxx]; + _space_rhog = new std::complex[nspin * ngmc]; + _space_rhog_save = new std::complex[nspin * ngmc]; + if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5) + { + _space_kin_r = new double[nspin * nrxx]; + _space_kin_r_save = new double[nspin * nrxx]; + } rho = new double*[nspin]; rhog = new std::complex*[nspin]; rho_save = new double*[nspin]; @@ -137,19 +143,19 @@ void Charge::allocate(const int& nspin_in) for (int is = 0; is < nspin; is++) { - rho[is] = new double[nrxx]; - rhog[is] = new std::complex[ngmc]; - rho_save[is] = new double[nrxx]; - rhog_save[is] = new std::complex[ngmc]; + rho[is] = _space_rho + is * nrxx; + rhog[is] = _space_rhog + is * ngmc; + rho_save[is] = _space_rho_save + is * nrxx; + rhog_save[is] = _space_rhog_save + is * ngmc; ModuleBase::GlobalFunc::ZEROS(rho[is], nrxx); ModuleBase::GlobalFunc::ZEROS(rhog[is], ngmc); ModuleBase::GlobalFunc::ZEROS(rho_save[is], nrxx); ModuleBase::GlobalFunc::ZEROS(rhog_save[is], ngmc); if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5) { - kin_r[is] = new double[nrxx]; + kin_r[is] = _space_kin_r + is * nrxx; ModuleBase::GlobalFunc::ZEROS(kin_r[is], nrxx); - kin_r_save[is] = new double[nrxx]; + kin_r_save[is] = _space_kin_r_save + is * nrxx; ModuleBase::GlobalFunc::ZEROS(kin_r_save[is], nrxx); } if(GlobalV::use_paw) diff --git a/source/module_elecstate/module_charge/charge.h b/source/module_elecstate/module_charge/charge.h index 829e2ff777..5abcb64530 100644 --- a/source/module_elecstate/module_charge/charge.h +++ b/source/module_elecstate/module_charge/charge.h @@ -42,7 +42,13 @@ class Charge double **kin_r = nullptr; // kinetic energy density in real space, for meta-GGA double **kin_r_save = nullptr; // kinetic energy density in real space, for meta-GGA // wenfei 2021-07-28 + private: + //temporary + double *_space_rho = nullptr, *_space_rho_save = nullptr; + std::complex *_space_rhog = nullptr, *_space_rhog_save = nullptr; + double *_space_kin_r = nullptr, *_space_kin_r_save = nullptr; + public: double **nhat = nullptr; //compensation charge for PAW double **nhat_save = nullptr; //compensation charge for PAW // wenfei 2023-09-05 diff --git a/source/module_elecstate/module_charge/charge_broyden.cpp b/source/module_elecstate/module_charge/charge_broyden.cpp deleted file mode 100644 index 728c770441..0000000000 --- a/source/module_elecstate/module_charge/charge_broyden.cpp +++ /dev/null @@ -1,181 +0,0 @@ -#include "charge_mixing.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" -#include "module_base/global_variable.h" -#include "module_base/inverse_matrix.h" -#include "module_base/lapack_connector.h" -#include "module_base/parallel_reduce.h" -#include "module_base/memory.h" -#include "module_base/timer.h" - -void Charge_Mixing::Simplified_Broyden_mixing(const int &iter, - Charge* chr) -{ - ModuleBase::TITLE("Charge_Mixing","Simplified_Broyden_mixing"); - ModuleBase::timer::tick("Charge", "Broyden_mixing"); - //It is a simplified modified broyden_mixing method. - //Ref: D.D. Johnson PRB 38, 12807 (1988) - //Here the weight w0 of the error of the inverse Jacobian is set to 0 and the weight wn of - //the error of each previous iteration is set to same. - - // (1) - this->allocate_Broyden(); - - int iter_used = std::min(iter-1, mixing_ndim); - int ipos = iter-2 - int((iter-2)/mixing_ndim) * mixing_ndim; - if(iter > 1) - { -#ifdef _OPENMP -#pragma omp parallel for collapse(2) schedule(static, 128) -#endif - for(int is=0; isrhopw->npw; ++ig) - { - dF[ipos][is][ig] -= chr->rhog[is][ig]; - dn[ipos][is][ig] -= chr->rhog_save[is][ig]; - } - } - } -#ifdef _OPENMP -#pragma omp parallel for collapse(2) schedule(static, 128) -#endif - for(int is=0; isrhopw->npw; ++ig) - { - dF[mixing_ndim][is][ig] = chr->rhog[is][ig]; - dn[mixing_ndim][is][ig] = chr->rhog_save[is][ig]; - } - } - - if(iter_used > 0) - { - this->beta.create(iter_used, iter_used,false); - for(int i = 0; i < iter_used; ++i) - { - for(int j = i; j < iter_used; ++j) - { - beta(i,j) = rhog_dot_product( this->dF[i], this->dF[j] ); - if(j != i) - { - beta(j,i)=beta(i,j); - } - } - } - double * work = new double [iter_used]; - int * iwork = new int [iter_used]; - char uu='U'; - int info; - dsytrf_(&uu,&iter_used,beta.c,&iter_used,iwork,work,&iter_used,&info); - if(info != 0) ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when factorizing beta."); - dsytri_(&uu,&iter_used,beta.c,&iter_used,iwork,work,&info); - if(info != 0) ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when DSYTRI beta."); - for(int i = 0; i < iter_used; ++i) - { - for(int j = i + 1; j < iter_used; ++j) - { - beta(i,j) = beta(j,i); - } - } - for(int i = 0 ; i < iter_used ; ++i) - { - work[i] = rhog_dot_product( this->dF[i], chr->rhog ); - } - for(int i = 0 ; i < iter_used ; ++i) - { - double gamma0 = 0; - for(int j = 0; j < iter_used ; ++j) - { - gamma0 += beta(i,j) * work[j]; - } -#ifdef _OPENMP -#pragma omp parallel for collapse(2) schedule(static, 512) -#endif - for(int is=0; isrhopw->npw; ++ig) - { - chr->rhog[is][ig] -= gamma0 * dF[i][is][ig]; - chr->rhog_save[is][ig] -= gamma0 * dn[i][is][ig]; - } - } - - } - delete[] work; - delete[] iwork; - } - int inext = iter-1 - int((iter-1)/mixing_ndim) * mixing_ndim; - -#ifdef _OPENMP -#pragma omp parallel for collapse(2) schedule(static, 128) -#endif - for(int is=0; isrhopw->npw; ++ig) - { - dF[inext][is][ig] = dF[mixing_ndim][is][ig]; - dn[inext][is][ig] = dn[mixing_ndim][is][ig]; - } - } - - - for(int is=0; isrhopw->npw; ++ig) - { - chr->rhog_save[is][ig] += mixing_beta * chr->rhog[is][ig]; - } - this->rhopw->recip2real( chr->rhog_save[is], chr->rho[is]); - } - ModuleBase::timer::tick("Charge", "Broyden_mixing"); - return; -} - -void Charge_Mixing::allocate_Broyden() -{ - if(!initb) - { - int npdim = mixing_ndim + 1; // another array is used for temporarily store - this->dF = new std::complex**[npdim]; - this->dn = new std::complex**[npdim]; - - for (int i=0; i*[GlobalV::NSPIN]; - dn[i] = new std::complex*[GlobalV::NSPIN]; - for (int is=0; is[this->rhopw->npw]; - dn[i][is] = new std::complex[this->rhopw->npw]; - } - } - ModuleBase::Memory::record("ChgMix::dF", sizeof(std::complex) * GlobalV::NSPIN*npdim*this->rhopw->npw); - ModuleBase::Memory::record("ChgMix::dn", sizeof(std::complex) * GlobalV::NSPIN*npdim*this->rhopw->npw); - this->initb = true; - } - - return; -} - -void Charge_Mixing::deallocate_Broyden() -{ - if (initb) - { - for (int i=0; iinitb = false; - } -} \ No newline at end of file diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index 05799ad9d5..398e49716a 100644 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -1,489 +1,483 @@ #include "charge_mixing.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" + +#include "module_base/element_elec_config.h" #include "module_base/inverse_matrix.h" +#include "module_base/module_mixing/broyden_mixing.h" +#include "module_base/module_mixing/plain_mixing.h" #include "module_base/parallel_reduce.h" #include "module_base/timer.h" -#include "module_base/element_elec_config.h" - +#include "module_hamilt_pw/hamilt_pwdft/global.h" Charge_Mixing::Charge_Mixing() { - rstep = 0; - dstep = rstep - 1;//alway like this. - initp = false; - initb = false; } Charge_Mixing::~Charge_Mixing() { - if (initp) - { - deallocate_Pulay(); - } - if (initb) - { - deallocate_Broyden(); - } + delete this->mixing; } -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 -) +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) { this->mixing_mode = mixing_mode_in; this->mixing_beta = mixing_beta_in; this->mixing_ndim = mixing_ndim_in; - this->mixing_gg0 = mixing_gg0_in; //mohan add 2014-09-27 - this->mixing_tau = mixing_tau_in; + this->mixing_gg0 = mixing_gg0_in; // mohan add 2014-09-27 + this->mixing_tau = mixing_tau_in; + if (GlobalV::SCF_THR_TYPE == 1) + { + this->rho_mdata.resize(GlobalV::NSPIN, this->rhopw->npw, sizeof(std::complex)); + this->tau_mdata.resize(GlobalV::NSPIN, this->rhopw->npw, sizeof(std::complex)); + } + else + { + this->rho_mdata.resize(GlobalV::NSPIN, this->rhopw->nrxx, sizeof(double)); + this->tau_mdata.resize(GlobalV::NSPIN, this->rhopw->nrxx, sizeof(double)); + } - if(mixing_tau && mixing_mode == "broyden") - { - GlobalV::ofs_running << "Note : mixing_tau has only been implemented for plain and pulay mixing" << std::endl; - } + if (this->mixing_mode == "broyden") + { + this->mixing = new Base_Mixing::Broyden_Mixing(this->mixing_ndim); + } + else if (this->mixing_mode == "plain") + { + this->mixing = new Base_Mixing::Plain_Mixing(); + } + // else if(this->mixing_mode == "pulay") + // { + // this->mixing = new Base_Mixing::Pulay_Mixing(this->mixing_ndim); + // } + else + { + ModuleBase::WARNING_QUIT("Charge_Mixing", "This Mixing mode is not implemended yet,coming soon."); + } + if (GlobalV::SCF_THR_TYPE == 1) + { + this->mixing->init_mixing_data(this->rho_mdata, + this->rhopw->npw * GlobalV::NSPIN, + sizeof(std::complex)); + if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) + this->mixing->init_mixing_data(this->tau_mdata, + this->rhopw->npw * GlobalV::NSPIN, + sizeof(std::complex)); + } + else + { + this->mixing->init_mixing_data(this->rho_mdata, this->rhopw->nrxx * GlobalV::NSPIN, sizeof(double)); + if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) + this->mixing->init_mixing_data(this->tau_mdata, this->rhopw->nrxx * GlobalV::NSPIN, sizeof(double)); + } return; } +void Charge_Mixing::set_rhopw(ModulePW::PW_Basis* rhopw_in) +{ + this->rhopw = rhopw_in; +} + void Charge_Mixing::need_auto_set() { - this->autoset = true; + 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---------------"<mixing_beta = 0.2; - } - else - { - this->mixing_beta = 0.7; - } - GlobalV::ofs_running<<" Autoset mixing_beta to "<mixing_beta<mixing_gg0 = 1.5; - } - else - { - this->mixing_gg0 = 0.0; - } - GlobalV::ofs_running<<" Autoset mixing_gg0 to "<mixing_gg0<autoset) + { + return; + } + else + { + this->autoset = false; + } + GlobalV::ofs_running << "--------------AUTO-SET---------------" << std::endl; + // 0.2 for metal and 0.7 for others + if (bandgap_in * ModuleBase::Ry_to_eV < 1.0) + { + this->mixing_beta = 0.2; + } + else + { + this->mixing_beta = 0.7; + } + GlobalV::ofs_running << " Autoset mixing_beta to " << this->mixing_beta << std::endl; + bool has_trans_metal = false; + // find elements of cell + for (int it = 0; it < ucell_.ntype; it++) + { + if (ModuleBase::IsTransMetal.find(ucell_.atoms[it].ncpp.psd) != ModuleBase::IsTransMetal.end()) + { + if (ModuleBase::IsTransMetal.at(ucell_.atoms[it].ncpp.psd)) + { + has_trans_metal = true; + } + } + } + // auto set kerker mixing for trans metal system + if (has_trans_metal) + { + this->mixing_gg0 = 1.5; + } + else + { + this->mixing_gg0 = 0.0; + } + GlobalV::ofs_running << " Autoset mixing_gg0 to " << this->mixing_gg0 << std::endl; + GlobalV::ofs_running << "-------------------------------------" << std::endl; + // auto set for inhomogeneous system } double Charge_Mixing::get_drho(Charge* chr, const double nelec) { - for (int is=0; isrhopw->real2recip(chr->rho[is], chr->rhog[is]); + ModuleBase::TITLE("Charge_Mixing", "get_drho"); + ModuleBase::timer::tick("Charge_Mixing", "get_drho"); + double drho = 0.0; - ModuleBase::GlobalFunc::NOTE("Perform FFT on rho_save(r) to obtain rho_save(G)."); - this->rhopw->real2recip(chr->rho_save[is], chr->rhog_save[is]); + if (GlobalV::SCF_THR_TYPE == 1) + { + for (int is = 0; is < GlobalV::NSPIN; ++is) + { + ModuleBase::GlobalFunc::NOTE("Perform FFT on rho(r) to obtain rho(G)."); + this->rhopw->real2recip(chr->rho[is], chr->rhog[is]); + ModuleBase::GlobalFunc::NOTE("Perform FFT on rho_save(r) to obtain rho_save(G)."); + this->rhopw->real2recip(chr->rho_save[is], chr->rhog_save[is]); + } - ModuleBase::GlobalFunc::NOTE("Calculate the charge difference between rho(G) and rho_save(G)"); + ModuleBase::GlobalFunc::NOTE("Calculate the charge difference between rho(G) and rho_save(G)"); + std::vector> drhog(GlobalV::NSPIN * this->rhopw->npw); #ifdef _OPENMP -#pragma omp parallel for schedule(static, 512) +#pragma omp parallel for collapse(2) schedule(static, 512) #endif - for (int ig=0; igrhopw->npw; ig++) + for (int is = 0; is < GlobalV::NSPIN; ++is) { - chr->rhog[is][ig] -= chr->rhog_save[is][ig]; + for (int ig = 0; ig < this->rhopw->npw; ig++) + { + drhog[is * rhopw->npw + ig] = chr->rhog[is][ig] - chr->rhog_save[is][ig]; + } } - } - ModuleBase::GlobalFunc::NOTE("Calculate the norm of the Residual std::vector: < R[rho] | R[rho_save] >"); - double scf_thr = this->rhog_dot_product( chr->rhog, chr->rhog); - - if(GlobalV::test_charge)GlobalV::ofs_running << " scf_thr from rhog_dot_product is " << scf_thr << std::endl; - - // This is a temporary method to choose the type of scf_thr, in order to test different scf_thr_type. - // Charge mixing will be refactored later. - if(GlobalV::SCF_THR_TYPE == 1 && !GlobalV::test_charge) - { - return scf_thr; - } - - // scf_thr calculated from real space. - double scf_thr2 = 0.0; - for(int is=0; is"); + drho = this->inner_dot_recip(drhog.data(), drhog.data()); + } + else + { + // Note: Maybe it is wrong. + // The inner_dot_real function (L1-norm) is different from that (L2-norm) in mixing. + for (int is = 0; is < GlobalV::NSPIN; is++) + { + if (is != 0 && is != 3 && GlobalV::DOMAG_Z) + { + continue; + } #ifdef _OPENMP -#pragma omp parallel for reduction(+:scf_thr2) +#pragma omp parallel for reduction(+ : drho) #endif - for(int ir=0; irrhopw->nrxx; ir++) - { - scf_thr2 += std::abs( chr->rho[is][ir] - chr->rho_save[is][ir] ); - } - } - - Parallel_Reduce::reduce_double_pool( scf_thr2 ); - assert( nelec != 0); - assert( GlobalC::ucell.omega > 0); - assert( this->rhopw->nxyz > 0); - scf_thr2 *= GlobalC::ucell.omega / static_cast( this->rhopw->nxyz ); - scf_thr2 /= nelec; - if(GlobalV::test_charge)GlobalV::ofs_running << " scf_thr from real space grid is " << scf_thr2 << std::endl; - - // mohan add 2011-01-22 - //if(LINEAR_SCALING && LOCAL_BASIS) xiaohui modify 2013-09-01 - if(GlobalV::SCF_THR_TYPE == 2) - { - scf_thr = scf_thr2; - } - return scf_thr; + for (int ir = 0; ir < this->rhopw->nrxx; ir++) + { + drho += std::abs(chr->rho[is][ir] - chr->rho_save[is][ir]); + } + } + Parallel_Reduce::reduce_double_pool(drho); + assert(nelec != 0); + assert(GlobalC::ucell.omega > 0); + assert(this->rhopw->nxyz > 0); + drho *= GlobalC::ucell.omega / static_cast(this->rhopw->nxyz); + drho /= nelec; + } + + ModuleBase::timer::tick("Charge_Mixing", "get_drho"); + return drho; } -void Charge_Mixing::mix_rho(const int &iter, Charge* chr) +void Charge_Mixing::mix_rho_recip(const int& iter, Charge* chr) { - ModuleBase::TITLE("Charge_Mixing","mix_rho"); - ModuleBase::timer::tick("Charge", "mix_rho"); - - // the charge before mixing. - double **rho123 = new double*[GlobalV::NSPIN]; - for(int is=0; isrhopw->nrxx]; - if(is!=0 && is!=3 && GlobalV::DOMAG_Z) - { - ModuleBase::GlobalFunc::ZEROS(rho123[is], this->rhopw->nrxx); - continue; - } -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 512) -#endif - for(int ir=0; irrhopw->nrxx; ++ir) - { - rho123[is][ir] = chr->rho[is][ir]; - } - } - - double **kin_r123; - if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) + // electronic density + // rhog and rhog_save are calculated in get_drho() function + std::complex* rhog_in = chr->rhog_save[0]; + std::complex* rhog_out = chr->rhog[0]; + + 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); + + auto inner_dot = std::bind(&Charge_Mixing::inner_dot_recip, this, std::placeholders::_1, std::placeholders::_2); + this->mixing->cal_coef(this->rho_mdata, inner_dot); + + this->mixing->mix_data(this->rho_mdata, rhog_out); + + for (int is = 0; is < GlobalV::NSPIN; is++) { - kin_r123 = new double*[GlobalV::NSPIN]; - for(int is=0; isrhopw->nrxx]; -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 512) -#endif - for(int ir=0; irrhopw->nrxx; ++ir) - { - kin_r123[is][ir] = chr->kin_r[is][ir]; - } - } - } - - if ( this->mixing_mode == "plain") + this->rhopw->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) { - this->plain_mixing(chr); + std::vector> taug_in(GlobalV::NSPIN * this->rhopw->npw); + std::vector> taug_out(GlobalV::NSPIN * this->rhopw->npw); + for (int is = 0; is < GlobalV::NSPIN; ++is) + { + this->rhopw->real2recip(chr->kin_r[is], &taug_out[is * rhopw->npw]); + this->rhopw->real2recip(chr->kin_r_save[is], &taug_in[is * rhopw->npw]); + } + // 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.data(), taug_out.data(), nullptr, false); + + this->mixing->mix_data(this->tau_mdata, taug_out.data()); + + for (int is = 0; is < GlobalV::NSPIN; is++) + { + this->rhopw->recip2real(&taug_out[is * rhopw->npw], chr->kin_r[is]); + } } - else if ( this->mixing_mode == "pulay") + + return; +} + +void Charge_Mixing::mix_rho_real(const int& iter, Charge* chr) +{ + // electronic density + double* rhor_in = chr->rho_save[0]; + double* 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); + + auto inner_dot = std::bind(&Charge_Mixing::inner_dot_real, this, std::placeholders::_1, std::placeholders::_2); + this->mixing->cal_coef(this->rho_mdata, inner_dot); + + this->mixing->mix_data(this->rho_mdata, rhor_out); + + double *taur_out, *taur_in; + if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) { - this->Pulay_mixing(chr); + taur_in = chr->kin_r_save[0]; + taur_out = chr->kin_r[0]; + // 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, taur_in, taur_out, nullptr, false); + + this->mixing->mix_data(this->tau_mdata, taur_out); } - else if ( this->mixing_mode == "broyden") +} + +void Charge_Mixing::mix_rho(const int& iter, Charge* chr) +{ + ModuleBase::TITLE("Charge_Mixing", "mix_rho"); + ModuleBase::timer::tick("Charge", "mix_rho"); + + if (GlobalV::SCF_THR_TYPE == 1) { - this->Simplified_Broyden_mixing(iter, chr); + mix_rho_recip(iter, chr); } else { - ModuleBase::WARNING_QUIT("Charge_Mixing","Not implemended yet,coming soon."); + mix_rho_real(iter, chr); } - // mohan add 2012-06-05 - // rho_save is the charge before mixing - for(int is=0; isrhopw->nrxx; ++ir) - { - chr->rho_save[is][ir] = rho123[is][ir]; - } - } - - for(int is=0; is* drhog) +{ + if (this->mixing_gg0 <= 0.0) + return; + const double fac = this->mixing_gg0; + const double gg0 = std::pow(fac * 0.529177 / GlobalC::ucell.tpiba, 2); #ifdef _OPENMP #pragma omp parallel for collapse(2) schedule(static, 512) #endif - for(int is=0; isrhopw->nrxx; ++ir) - { - chr->kin_r_save[is][ir] = kin_r123[is][ir]; - } - } - - for(int is=0; isrhopw->npw; ++ig) + { + double gg = this->rhopw->gg[ig]; + double filter_g = std::max(gg / (gg + gg0), 0.1); + drhog[is * this->rhopw->npw + ig] *= filter_g; + } + } + return; +} - if(new_e_iteration) new_e_iteration = false; +void Charge_Mixing::Kerker_screen_real(double* drhor) +{ + std::vector> drhog(this->rhopw->npw * GlobalV::NSPIN); + for (int is = 0; is < GlobalV::NSPIN; ++is) + { + this->rhopw->real2recip(drhor + is * this->rhopw->nrxx, drhog.data() + is * this->rhopw->npw); + } + Kerker_screen_recip(drhog.data()); + for (int is = 0; is < GlobalV::NSPIN; ++is) + { + this->rhopw->recip2real(drhog.data() + is * this->rhopw->npw, drhor + is * this->rhopw->nrxx); + } +} - ModuleBase::timer::tick("Charge","mix_rho"); - return; +double Charge_Mixing::inner_dot_recip(std::complex* rho1, std::complex* rho2) +{ + std::complex** rho1_2d = new std::complex*[GlobalV::NSPIN]; + std::complex** rho2_2d = new std::complex*[GlobalV::NSPIN]; + for (int is = 0; is < GlobalV::NSPIN; is++) + { + rho1_2d[is] = rho1 + is * this->rhopw->npw; + rho2_2d[is] = rho2 + is * this->rhopw->npw; + } + double result = this->rhog_dot_product(rho1_2d, rho2_2d); + delete[] rho1_2d; + delete[] rho2_2d; + return result; } -void Charge_Mixing::plain_mixing(Charge* chr) const +double Charge_Mixing::inner_dot_real(double* rho1, double* rho2) { - ModuleBase::timer::tick("Charge", "plain_mixing"); - // if mixing_beta == 1, each electron iteration, - // use all new charge density, - // on the contrary, if mixing_beta == 0, - // no new charge will be generated! - const double mix_old = 1 - mixing_beta; - -//xiaohui add 2014-12-09 - for (int is=0; ismixing_gg0 > 0.0) - { - double* Rrho = new double[this->rhopw->nrxx]; - std::complex *kerpulay = new std::complex[this->rhopw->npw]; - double* kerpulayR = new double[this->rhopw->nrxx]; -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 512) -#endif - for(int ir=0; irrhopw->nrxx; ir++) - { - Rrho[ir] = chr->rho[is][ir] - chr->rho_save[is][ir]; - } - this->rhopw->real2recip(Rrho, kerpulay); - - const double fac = this->mixing_gg0; - const double gg0 = std::pow(fac * 0.529177 / GlobalC::ucell.tpiba, 2); - double* filter_g = new double[this->rhopw->npw]; -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 256) -#endif - for(int ig=0; igrhopw->npw; ig++) - { - double gg = this->rhopw->gg[ig]; - filter_g[ig] = std::max(gg / (gg + gg0), 0.1); - - kerpulay[ig] = (1 - filter_g[ig]) * kerpulay[ig]; - } - this->rhopw->recip2real(kerpulay, kerpulayR); -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 128) -#endif - for(int ir=0; irrhopw->nrxx; ir++) - { - Rrho[ir] = Rrho[ir] - kerpulayR[ir]; - chr->rho[is][ir] = Rrho[ir] * mixing_beta + chr->rho_save[is][ir]; - chr->rho_save[is][ir] = chr->rho[is][ir]; - } - - delete[] Rrho; - delete[] kerpulay; - delete[] kerpulayR; - delete[] filter_g; - } - else - { -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 256) -#endif - for (int ir=0; irrhopw->nrxx; ir++) - { - chr->rho[is][ir] = chr->rho[is][ir]*mixing_beta + mix_old*chr->rho_save[is][ir]; - chr->rho_save[is][ir] = chr->rho[is][ir]; - } - } - - if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) - { + double rnorm = 0.0; #ifdef _OPENMP -#pragma omp parallel for schedule(static, 256) +#pragma omp parallel for reduction(+ : rnorm) #endif - for (int ir=0; irrhopw->nrxx; ir++) - { - chr->kin_r[is][ir] = chr->kin_r[is][ir]*mixing_beta + mix_old*chr->kin_r_save[is][ir]; - chr->kin_r_save[is][ir] = chr->kin_r[is][ir]; - } - } - - } - ModuleBase::timer::tick("Charge", "plain_mixing"); - return; + for (int ir = 0; ir < this->rhopw->nrxx * GlobalV::NSPIN; ++ir) + { + rnorm += rho1[ir] * rho2[ir]; + } + Parallel_Reduce::reduce_double_pool(rnorm); + return rnorm; } -double Charge_Mixing::rhog_dot_product( - const std::complex*const*const rhog1, - const std::complex*const*const rhog2 -) const +double Charge_Mixing::rhog_dot_product(const std::complex* const* const rhog1, + const std::complex* const* const rhog2) const { - ModuleBase::TITLE("Charge_Mixing","rhog_dot_product"); - ModuleBase::timer::tick("Charge_Mixing","rhog_dot_product"); + ModuleBase::TITLE("Charge_Mixing", "rhog_dot_product"); + ModuleBase::timer::tick("Charge_Mixing", "rhog_dot_product"); static const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / GlobalC::ucell.tpiba2; static const double fac2 = ModuleBase::e2 * ModuleBase::FOUR_PI / (ModuleBase::TWO_PI * ModuleBase::TWO_PI); double sum = 0.0; - - auto part_of_noncolin = [&]() // Peize Lin change goto to function at 2020.01.31 - { - double sum = 0.0; + + auto part_of_noncolin = [&]() // Peize Lin change goto to function at 2020.01.31 + { + double sum = 0.0; #ifdef _OPENMP -#pragma omp parallel for reduction(+:sum) +#pragma omp parallel for reduction(+ : sum) #endif - for (int ig=0; igrhopw->npw; ++ig) - { - if(this->rhopw->gg[ig]<1e-8) continue; - sum += ( conj( rhog1[0][ig] )* rhog2[0][ig] ).real() / this->rhopw->gg[ig]; - } - sum *= fac; - return sum; - }; - - switch ( GlobalV::NSPIN ) + for (int ig = 0; ig < this->rhopw->npw; ++ig) + { + if (this->rhopw->gg[ig] < 1e-8) + continue; + sum += (conj(rhog1[0][ig]) * rhog2[0][ig]).real() / this->rhopw->gg[ig]; + } + sum *= fac; + return sum; + }; + + switch (GlobalV::NSPIN) { - case 1: - sum += part_of_noncolin(); - break; + case 1: + sum += part_of_noncolin(); + break; - case 2: - { - // (1) First part of density error. + case 2: { + // (1) First part of density error. #ifdef _OPENMP -#pragma omp parallel for reduction(+:sum) +#pragma omp parallel for reduction(+ : sum) #endif - for (int ig=0; igrhopw->npw; ++ig) - { - if(this->rhopw->gg[ig]<1e-8) continue; - sum += ( conj( rhog1[0][ig]+rhog1[1][ig] ) * (rhog2[0][ig]+rhog2[1][ig]) ).real() / this->rhopw->gg[ig]; - } - sum *= fac; + for (int ig = 0; ig < this->rhopw->npw; ++ig) + { + if (this->rhopw->gg[ig] < 1e-8) + continue; + sum += (conj(rhog1[0][ig] + rhog1[1][ig]) * (rhog2[0][ig] + rhog2[1][ig])).real() / this->rhopw->gg[ig]; + } + sum *= fac; - if(GlobalV::GAMMA_ONLY_PW) - { - sum *= 2.0; - } + if (GlobalV::GAMMA_ONLY_PW) + { + sum *= 2.0; + } - // (2) Second part of density error. - // including |G|=0 term. - double sum2 = 0.0; + // (2) Second part of density error. + // including |G|=0 term. + double sum2 = 0.0; - sum2 += fac2 * ( conj( rhog1[0][0]-rhog1[1][0] ) * ( rhog2[0][0]-rhog2[1][0] ) ).real(); + sum2 += fac2 * (conj(rhog1[0][0] - rhog1[1][0]) * (rhog2[0][0] - rhog2[1][0])).real(); - double mag = 0.0; + double mag = 0.0; #ifdef _OPENMP -#pragma omp parallel for reduction(+:mag) +#pragma omp parallel for reduction(+ : mag) #endif - for (int ig=0; igrhopw->npw; ig++) - { - mag += ( conj( rhog1[0][ig]-rhog1[1][ig] ) * ( rhog2[0][ig]-rhog2[1][ig] ) ).real(); - } - mag *= fac2; - - //if(GlobalV::GAMMA_ONLY_PW); - if(GlobalV::GAMMA_ONLY_PW) // Peize Lin delete ; 2020.01.31 - { - mag *= 2.0; - } - - //std::cout << " sum=" << sum << " mag=" << mag << std::endl; - sum2 += mag; - sum += sum2; - break; - } - case 4: - // non-collinear spin, added by zhengdy - if(!GlobalV::DOMAG&&!GlobalV::DOMAG_Z) - sum += part_of_noncolin(); - else - { - //another part with magnetization + for (int ig = 0; ig < this->rhopw->npw; ig++) + { + mag += (conj(rhog1[0][ig] - rhog1[1][ig]) * (rhog2[0][ig] - rhog2[1][ig])).real(); + } + mag *= fac2; + + // if(GlobalV::GAMMA_ONLY_PW); + if (GlobalV::GAMMA_ONLY_PW) // Peize Lin delete ; 2020.01.31 + { + mag *= 2.0; + } + + // std::cout << " sum=" << sum << " mag=" << mag << std::endl; + sum2 += mag; + sum += sum2; + break; + } + case 4: + // non-collinear spin, added by zhengdy + if (!GlobalV::DOMAG && !GlobalV::DOMAG_Z) + sum += part_of_noncolin(); + else + { + // another part with magnetization #ifdef _OPENMP -#pragma omp parallel for reduction(+:sum) +#pragma omp parallel for reduction(+ : sum) #endif - for (int ig=0; igrhopw->npw; ig++) - { - if(ig==this->rhopw->ig_gge0) continue; - sum += ( conj( rhog1[0][ig] )* rhog2[0][ig] ).real() / this->rhopw->gg[ig]; - } - sum *= fac; - const int ig0 = this->rhopw->ig_gge0; - if(ig0 > 0) - { - sum += fac2 * ((conj( rhog1[1][ig0])*rhog2[1][ig0]).real() + - (conj( rhog1[2][ig0])*rhog2[2][ig0]).real() + - (conj( rhog1[3][ig0])*rhog2[3][ig0]).real()); - } - double fac3 = fac2; - if(GlobalV::GAMMA_ONLY_PW) - { - fac3 *= 2.0; - } + for (int ig = 0; ig < this->rhopw->npw; ig++) + { + if (ig == this->rhopw->ig_gge0) + continue; + sum += (conj(rhog1[0][ig]) * rhog2[0][ig]).real() / this->rhopw->gg[ig]; + } + sum *= fac; + const int ig0 = this->rhopw->ig_gge0; + if (ig0 > 0) + { + sum += fac2 + * ((conj(rhog1[1][ig0]) * rhog2[1][ig0]).real() + (conj(rhog1[2][ig0]) * rhog2[2][ig0]).real() + + (conj(rhog1[3][ig0]) * rhog2[3][ig0]).real()); + } + double fac3 = fac2; + if (GlobalV::GAMMA_ONLY_PW) + { + fac3 *= 2.0; + } #ifdef _OPENMP -#pragma omp parallel for reduction(+:sum) +#pragma omp parallel for reduction(+ : sum) #endif - for (int ig=0; igrhopw->npw; ig++) - { - if(ig == ig0) continue; - sum += fac3 * ((conj( rhog1[1][ig])*rhog2[1][ig]).real() + - (conj( rhog1[2][ig])*rhog2[2][ig]).real() + - (conj( rhog1[3][ig])*rhog2[3][ig]).real()); - } - } - break; + for (int ig = 0; ig < this->rhopw->npw; ig++) + { + if (ig == ig0) + continue; + sum += fac3 + * ((conj(rhog1[1][ig]) * rhog2[1][ig]).real() + (conj(rhog1[2][ig]) * rhog2[2][ig]).real() + + (conj(rhog1[3][ig]) * rhog2[3][ig]).real()); + } + } + break; } - Parallel_Reduce::reduce_double_pool( sum ); + Parallel_Reduce::reduce_double_pool(sum); - ModuleBase::timer::tick("Charge_Mixing","rhog_dot_product"); + ModuleBase::timer::tick("Charge_Mixing", "rhog_dot_product"); - sum *= GlobalC::ucell.omega * 0.5; + sum *= GlobalC::ucell.omega * 0.5; return sum; } - diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index 5cfb6bf492..994081e1f0 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -1,51 +1,46 @@ -//================================= -// Mohan add 2010-02-01 -//================================= + #ifndef CHARGE_MIXING_H #define CHARGE_MIXING_H - -//================================== -// (1) Plain Mixing -// (2) KerKer Mixing -// (3) Pulay Mixing -// (4) Modified Broden Mixing -//=================================== +#include "charge.h" #include "module_base/global_function.h" #include "module_base/global_variable.h" -#include "module_base/matrix.h" +#include "module_base/module_mixing/mixing.h" #include "module_cell/unitcell.h" -#include "charge.h" class Charge_Mixing { - public: - Charge_Mixing(); - ~Charge_Mixing(); + public: + Charge_Mixing(); + ~Charge_Mixing(); + Base_Mixing::Mixing* mixing = nullptr; + Base_Mixing::Mixing_Data rho_mdata; + Base_Mixing::Mixing_Data tau_mdata; + + void mix_rho(const int& iter, Charge* chr); + void mix_rho_recip(const int& iter, Charge* chr); + void mix_rho_real(const int& iter, Charge* chr); -//====================================== -// General interfaces, in charge_mixing.cpp -//====================================== - void 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 - );//mohan add mixing_gg0_in 2014-09-27 + void Kerker_screen_recip(std::complex* rhog); + void Kerker_screen_real(double* rho); - void need_auto_set(); - void auto_set(const double& bandgap_in, const UnitCell& ucell_); + double inner_dot_recip(std::complex* rho1, std::complex* rho2); + double inner_dot_real(double* rho1, double* rho2); - double get_drho(Charge* chr, const double nelec); + //====================================== + // General interfaces, in charge_mixing.cpp + //====================================== + void 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); // mohan add mixing_gg0_in 2014-09-27 - void mix_rho(const int &iter, Charge* chr);// mix rho + void need_auto_set(); + void auto_set(const double& bandgap_in, const UnitCell& ucell_); - //for Pulay method - //if first electronic step, then reset charge mixing - void reset(); + 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, sunliang add 2023-05-08 + void set_rhopw(ModulePW::PW_Basis* rhopw_in); // extracting parameters // normally these parameters will not be used @@ -56,93 +51,24 @@ class Charge_Mixing int get_mixing_ndim() const {return mixing_ndim;} double get_mixing_gg0() const {return mixing_gg0;} - int get_totstep() const {return totstep;} - int get_rstep() const {return rstep;} - int get_dstep() const {return dstep;} - int get_idstep() const {return idstep;} - double* get_alpha() const {return alpha;} - - private: - -//====================================== -// General parameters -//====================================== + private: + //====================================== + // General parameters + //====================================== std::string mixing_mode; double mixing_beta; int mixing_ndim; - double mixing_gg0; //mohan add 2014-09-27 - bool mixing_tau; + double mixing_gg0; // mohan add 2014-09-27 + bool mixing_tau; bool new_e_iteration; - ModulePW::PW_Basis* rhopw = nullptr; - -//====================================== -// simple plain mixing method, in charge_mixing.cpp -//====================================== - void plain_mixing(Charge* chr) const; - - double rhog_dot_product(const std::complex*const*const rhog1, const std::complex*const*const rhog2) const; - -//====================================== -// Pulay mixing method, in charge_pulay.cpp -//====================================== - - void Pulay_mixing(Charge* chr); - - bool initp; // p stands for pulay algorithms - void allocate_Pulay(); - void deallocate_Pulay(); - - // auxiliary variables / subroutines - int irstep; //mohan add 2012-02-10 - int idstep; - int totstep; - int rstep; // the record step; - int dstep; // Delta step " dstep = rstep-1 ". - double* alpha; // - sum (Abar * dRR) - - double*** Rrho;// Rrho(i) = rho(i) - rho_save(i), (GlobalV::NSPIN, rstep, pw.nrxx) - double*** dRrho;// dRrho(i) = Rrho(i+1) - Rrho(i), (GlobalV::NSPIN, dstep, pw.nrxx) - double*** drho;// drho(i)= rho_save(i+1) - rho_save2(i), (GlobalV::NSPIN, dstep, pw.nrxx) - double** rho_save2;//rho_save: rho_in, rho_save2: rho_in(last step) - - double*** Rtau;//same things, but for kinetic energy density - double*** dRtau; - double*** dtau; - double** tau_save2; - - ModuleBase::matrix Abar; // ^{-1} - double* dRR; // - - void generate_datas(const int &irstep, const int &idstep, const int &totstep, Charge* chr); - void generate_Abar(ModuleBase::matrix &A)const; - void inverse_preA(const int &dim, ModuleBase::matrix &preA)const; - void inverse_real_symmetry_matrix(ModuleBase::matrix &A)const; // indicate the spin. - void generate_dRR(const int &m); - void generate_alpha(); - void generate_new_rho(const int &is,const int &m, Charge* chr); - - void generate_residual_vector(double *residual, const double* rho_out, const double* rho_in)const; - -//====================================== -// Broyden mixing method, in charge_broyden.cpp -//====================================== - - void Simplified_Broyden_mixing(const int &iter, - Charge* chr); //qianrui created 2021-5-15 - - bool initb; // b stands for Broyden algorithms. - void allocate_Broyden(); - void deallocate_Broyden(); - - ModuleBase::matrix beta; // (dstep, dstep) - - std::complex*** dF; // dF(i) = rhog(i) - rhog_save(i), (GlobalV::NSPIN, rstep, rhopw->npw) - std::complex*** dn; // dn(i) = rhog(i+1) - rhog(i), (GlobalV::NSPIN, rstep, rhopw->npw) + ModulePW::PW_Basis* rhopw = nullptr; + bool autoset = false; - private: - bool autoset = false; + private: + double rhog_dot_product(const std::complex* const* const rhog1, + const std::complex* const* const rhog2) const; }; #endif diff --git a/source/module_elecstate/module_charge/charge_pulay.cpp b/source/module_elecstate/module_charge/charge_pulay.cpp deleted file mode 100644 index 899a055e1d..0000000000 --- a/source/module_elecstate/module_charge/charge_pulay.cpp +++ /dev/null @@ -1,696 +0,0 @@ -#include "charge_mixing.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" -#include "module_base/inverse_matrix.h" -#include "module_base/parallel_reduce.h" -#include "module_base/memory.h" -#include "module_base/tool_threading.h" -#include "module_base/timer.h" - -static inline double calculate_residual_norm(int nrxx, double *residual1, double* residual2) -{ - // calculate the norm of the residual std::vector: - // (the target to minimize in Pulay's algorithm) - double rnorm = 0.0; -#ifdef _OPENMP -#pragma omp parallel for reduction(+:rnorm) -#endif - for(int ir=0; irmixing_ndim; - dstep = this->mixing_ndim - 1; - assert(dstep>0); - - // (1) allocate - this->allocate_Pulay(); - - if (irstep==rstep) irstep=0; - if (idstep==dstep) idstep=0; - - if(GlobalV::test_charge)ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"irstep",irstep); - if(GlobalV::test_charge)ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"idstep",idstep); - if(GlobalV::test_charge)ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"totstep",totstep); - - //---------------------------------------------- - // calculate "dR^{i} = R^{i+1} - R^{i}" - // calculate "drho^{i} = rho^{i+1} - rho^{i}" - //---------------------------------------------- - this->generate_datas(irstep, idstep, totstep, chr); - - // not enough steps, not full matrix. - if(totstep < dstep) - { - int premix = 1; - - // mohan update 2011-06-14 - if(totstep==0) premix = 2; - - if(premix == 1) - { - this->generate_Abar(Abar); - - //----------------------------- - // inverse part of the matrix. - //----------------------------- - const int predim = idstep+1; - ModuleBase::matrix preA(predim, predim); - for(int i=0; iinverse_preA(predim,preA); - - - //---------------------------------------------- - // get the information from part of the matrix. - //---------------------------------------------- - Abar.zero_out(); - for(int i=0; igenerate_dRR(irstep); - - this->generate_alpha(); - - for(int is=0; isgenerate_new_rho(is,irstep,chr); - } - } - - if(premix == 2) - { - // if not enough step, take kerker mixing method. - this->plain_mixing(chr); - } - - if(totstep>0) - { - ++idstep; - } - ++irstep; - ++totstep; - } - else - { - // generate matrix A = - this->generate_Abar(Abar); - - // inverse A matrix to become Abar. - this->inverse_real_symmetry_matrix(Abar); - - this->generate_dRR(irstep); - - this->generate_alpha(); - - for(int is=0; isgenerate_new_rho(is,irstep,chr); - } - - ++irstep; - ++idstep; - ++totstep; - } - - ModuleBase::timer::tick("Charge", "Pulay_mixing"); - return; -} - -void Charge_Mixing::reset() // Peize Lin add 2018-11-01 -{ - this->new_e_iteration = true; - - irstep = 0; - idstep = 0; - totstep = 0; - - // liuyu add 2023-03-29 - // if md_prec_level == 2, charge mixing should re-allocate - // due to the change of FFT grids - if (GlobalV::md_prec_level == 2) - { - if (this->mixing_mode == "pulay") - { - this->deallocate_Pulay(); - } - else if (this->mixing_mode == "broyden") - { - this->deallocate_Broyden(); - } - } -} - -void Charge_Mixing::set_rhopw(ModulePW::PW_Basis* rhopw_in) -{ - this->rhopw = rhopw_in; -} - -void Charge_Mixing::allocate_Pulay() -{ - auto zeros_kernel = [&](int num_threads, int thread_id) - { - int beg, len; - ModuleBase::TASK_DIST_1D(num_threads, thread_id, this->rhopw->nrxx, beg, len); - for(int i=0; iinitp) - { - ModuleBase::TITLE("Charge_Mixing","allocate_pulay"); - ModuleBase::GlobalFunc::NOTE("rstep is used to record Rrho"); - if(GlobalV::test_charge)ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"rstep",rstep); - ModuleBase::GlobalFunc::NOTE("dstep is used to record dRrho, drho"); - if(GlobalV::test_charge)ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"dstep",dstep); - assert(rstep>1); - dstep = rstep - 1; - - // (1) allocate Rrho[i]: rho_out[i] - rho_in[i] - this->Rrho = new double**[GlobalV::NSPIN]; - for (int is=0; isrhopw->nrxx]; - } - } - ModuleBase::Memory::record("ChgMix::Rrho", sizeof(double) * GlobalV::NSPIN*rstep*this->rhopw->nrxx); - - // (2) allocate "dRrho[i] = Rrho[i+1] - Rrho[i]" of the last few steps. - // allocate "drho[i] = rho[i+1] - rho[i]" of the last few steps. - // allocate rho_save2: rho[i] of the last few steps. - this->dRrho = new double**[GlobalV::NSPIN]; - this->drho = new double**[GlobalV::NSPIN]; - this->rho_save2 = new double*[GlobalV::NSPIN]; - - for (int is=0; isrhopw->nrxx]; - - for (int i=0; irhopw->nrxx]; - drho[is][i] = new double[this->rhopw->nrxx]; - } - } - ModuleBase::Memory::record("ChgMix::dRrho", sizeof(double) * GlobalV::NSPIN*dstep*this->rhopw->nrxx); - ModuleBase::Memory::record("ChgMix::drho", sizeof(double) * GlobalV::NSPIN*dstep*this->rhopw->nrxx); - ModuleBase::Memory::record("ChgMix::rho_save2", sizeof(double) * GlobalV::NSPIN*this->rhopw->nrxx); - - if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) - { - this->Rtau = new double**[GlobalV::NSPIN]; - for (int is=0; isrhopw->nrxx]; - } - } - ModuleBase::Memory::record("ChgMix::Rtau", sizeof(double) * GlobalV::NSPIN*rstep*this->rhopw->nrxx); - - this->dRtau = new double**[GlobalV::NSPIN]; - this->dtau = new double**[GlobalV::NSPIN]; - this->tau_save2 = new double*[GlobalV::NSPIN]; - - for (int is=0; isrhopw->nrxx]; - - for (int i=0; irhopw->nrxx]; - dtau[is][i] = new double[this->rhopw->nrxx]; - } - } - ModuleBase::Memory::record("ChgMix::dRtau", sizeof(double) * GlobalV::NSPIN*dstep*this->rhopw->nrxx); - ModuleBase::Memory::record("ChgMix::dtau", sizeof(double) * GlobalV::NSPIN*dstep*this->rhopw->nrxx); - ModuleBase::Memory::record("ChgMix::tau_save2", sizeof(double) * GlobalV::NSPIN*this->rhopw->nrxx); - } - - ModuleBase::GlobalFunc::NOTE("Allocate Abar = , dimension = dstep."); - this->Abar.create(dstep, dstep); - ModuleBase::Memory::record("ChgMix::Abar", sizeof(double) * dstep*dstep); - - // (4) allocate dRR = - ModuleBase::GlobalFunc::NOTE("Allocate dRR = < dR | R >, dimension = dstep"); - this->dRR = new double[dstep]; - - // (5) allocate alpha - ModuleBase::GlobalFunc::NOTE("Allocate alpha, dimension = dstep"); - this->alpha = new double[dstep]; - - // (6) zeros all arrays - ModuleBase::OMP_PARALLEL(zeros_kernel); - - this->initp = true; - } - - // mohan add 2010-07-16 - if(this->new_e_iteration) - { - ModuleBase::OMP_PARALLEL(zeros_kernel); - } - return; -} - -void Charge_Mixing::deallocate_Pulay() -{ - if (!this->initp) return; - // delete: Rrho[i] = rho_out[i] - rho_in[i]; - for (int is=0; isrhopw->nrxx) - for (int is=0; isinitp = false; -} - -// calculate < dR | dR > -// if spin is considered, double the size. -// < dR1,dR2 | dR1,dR2 > = < dR1 | dR1 > + < dR2 | dR2 > -void Charge_Mixing::generate_Abar(ModuleBase::matrix &A)const -{ - int step = 0; - - step=dstep; - - A.zero_out(); - - for(int is=0; isrhopw->nrxx, this->dRrho[is][j], this->dRrho[is][i] ); - if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) - { - A(i,j) += calculate_residual_norm(this->rhopw->nrxx, this->dRtau[is][j], this->dRtau[is][i] ); - } - A(j,i) = A(i,j); - } - } - } - Parallel_Reduce::reduce_double_pool(A.c, A.nr * A.nc); - return; -} - -#include "module_base/complexmatrix.h" -void Charge_Mixing::inverse_preA(const int &dim, ModuleBase::matrix &preA)const -{ - ModuleBase::ComplexMatrix B(dim, dim); - ModuleBase::ComplexMatrix C(dim, dim); - for(int i=0; i (preA(i,j), 0.0); - } - } - ModuleBase::Inverse_Matrix_Complex IMC; - IMC.init(dim); - IMC.using_zheev(B,C); - - for(int i=0; i (A(i,j),0.0); - } - } - - ModuleBase::Inverse_Matrix_Complex IMC; - IMC.init(step); - IMC.using_zheev(B,C); - - for(int i=0; idRR[i] += calculate_residual_norm(this->rhopw->nrxx, this->dRrho[is][i], this->Rrho[is][m]); - if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) - { - this->dRR[i] += calculate_residual_norm(this->rhopw->nrxx, this->dRtau[is][i], this->Rtau[is][m]); - } - } - } - Parallel_Reduce::reduce_double_pool(dRR, dstep); - - return; -} - -// use dstep to genearte Abar(dstep, dstep) -void Charge_Mixing::generate_alpha() -{ - for(int i=0; ialpha[i] = 0; - for(int j=0; jalpha[i] -= this->Abar(j,i) * this->dRR[j]; - } - } - - return; -} - -void Charge_Mixing::generate_new_rho(const int &is, const int &m, Charge* chr) -{ - double mixp = this->mixing_beta; -#ifdef _OPENMP -#pragma omp parallel for -#endif - for(int ir=0; irrhopw->nrxx; ir++) - { - double rhonew = chr->rho_save[is][ir] + mixp * this->Rrho[is][m][ir]; - for(int i=0; ialpha[i] * ( this->drho[is][i][ir] + mixp * this->dRrho[is][i][ir] ); - } - chr->rho[is][ir] = rhonew; - } - - if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) - { -#ifdef _OPENMP -#pragma omp parallel for -#endif - for(int ir=0; irrhopw->nrxx; ir++) - { - double rhonew = chr->kin_r_save[is][ir] + mixp * this->Rtau[is][m][ir]; - for(int i=0; ialpha[i] * ( this->dtau[is][i][ir] + mixp * this->dRtau[is][i][ir] ); - } - chr->kin_r[is][ir] = rhonew; - } - } - return; -} - -void Charge_Mixing::generate_residual_vector(double *residual, const double* rho_out, const double* rho_in)const -{ -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 512) -#endif - for (int ir=0; irrhopw->nrxx; ir++) - { - residual[ir]= rho_out[ir] - rho_in[ir]; - } - return; -} - - -// calculate "dR^{i} = R^{i+1} - R^{i}" -// calculate "drho^{i} = rho^{i+1} - rho^{i}" -void Charge_Mixing::generate_datas(const int &irstep, const int &idstep, const int &totstep, Charge* chr) -{ - //=============================================== - // calculate the important "Rrho". - // init the "Rrho = rho - rho_save" - // which Rrho to be update now? answer: irstep - //=============================================== - - ModuleBase::GlobalFunc::NOTE("Generate Residual std::vector from rho and rho_save."); - for (int is=0; isgenerate_residual_vector( this->Rrho[is][irstep], chr->rho[is], chr->rho_save[is]); - if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) - { - this->generate_residual_vector( this->Rtau[is][irstep], chr->kin_r[is], chr->kin_r_save[is]); - } - - // 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. - - if(this->mixing_gg0 > 0.0) - { - std::complex *kerpulay = new std::complex[this->rhopw->npw]; - double* kerpulayR = new double[this->rhopw->nrxx]; - - this->rhopw->real2recip(Rrho[is][irstep], kerpulay); - - const double fac = this->mixing_gg0; - const double gg0 = std::pow(fac * 0.529177 /GlobalC::ucell.tpiba, 2); -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 128) -#endif - for(int ig=0; igrhopw->npw; ig++) - { - double gg = this->rhopw->gg[ig]; - double filter_g = std::max(gg / (gg + gg0), 0.1); - kerpulay[ig] = (1 - filter_g) * kerpulay[ig]; - } - - this->rhopw->recip2real(kerpulay, kerpulayR); -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 256) -#endif - for(int ir=0; irrhopw->nrxx; ir++) - { - Rrho[is][irstep][ir] = Rrho[is][irstep][ir] - kerpulayR[ir]; - } - - delete[] kerpulay; - delete[] kerpulayR; - } - } - - ModuleBase::OMP_PARALLEL([&](int num_threads, int thread_id) - { - int irbeg, irend; - ModuleBase::TASK_DIST_1D(num_threads, thread_id, this->rhopw->nrxx, irbeg, irend); - irend = irbeg + irend; - if(totstep==0) - { - // don't need to calculate 'dRrho' and 'drho'. - - // mohan fix the bug 2010/03/26. - // save 'rho_save2' in order to calculate drho in - // the next iteration. - for(int is=0; isrho_save2[is][ir] = chr->rho_save[is][ir]; - } - } - - if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) - { - for(int is=0; istau_save2[is][ir] = chr->kin_r_save[is][ir]; - } - } - } - } - else if(totstep>0) - { - // which dRrho to be update now? answer: dstep=istep-1; - // irstep range: [0, rstep) - const int nowR = irstep; - int lastR = irstep - 1; - - if (thread_id == 0) - { - if(GlobalV::test_charge)ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "now irstep", nowR); - if(GlobalV::test_charge)ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "last irstep",lastR); - } - - if(lastR < 0) lastR += rstep; - for (int is=0; isdRrho[is][idstep][ir] = this->Rrho[is][nowR][ir] - this->Rrho[is][lastR][ir]; - this->drho[is][idstep][ir] = chr->rho_save[is][ir] - this->rho_save2[is][ir]; - // save 'rho_save2' in order to calculate drho in - this->rho_save2[is][ir] = chr->rho_save[is][ir]; - } - } - - if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) - { - for (int is=0; isdRtau[is][idstep][ir] = this->Rtau[is][nowR][ir] - this->Rtau[is][lastR][ir]; - this->dtau[is][idstep][ir] = chr->kin_r_save[is][ir] - this->tau_save2[is][ir]; - // save 'tau_save2' in order to calculate drho in - this->tau_save2[is][ir] = chr->kin_r_save[is][ir]; - } - } - } - } - }); - - ModuleBase::GlobalFunc::NOTE("Calculate drho = rho_{in}^{i+1} - rho_{in}^{i}"); - return; -} - diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 72e07e56af..ee1c17f67a 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -445,12 +445,6 @@ void ESolver_KS_LCAO::Init_Basis_lcao(ORB_control& orb_con, Input& inp, UnitCell void ESolver_KS_LCAO::eachiterinit(const int istep, const int iter) { - - // mohan add 2010-07-16 - // used for pulay mixing. - if (iter == 1) - this->p_chgmix->reset(); - // mohan update 2012-06-05 this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband(); diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 1e046cc24a..89c6dec4c2 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -388,10 +388,6 @@ void ESolver_KS_PW::othercalculation(const int istep) template void ESolver_KS_PW::eachiterinit(const int istep, const int iter) { - // mohan add 2010-07-16 - if (iter == 1) - this->p_chgmix->reset(); - // mohan move harris functional to here, 2012-06-05 // use 'rho(in)' and 'v_h and v_xc'(in) this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband();