From 81b98f0d962c36fc7acd26acd874ecceaa4ded68 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Wed, 18 Sep 2024 16:48:24 +0800 Subject: [PATCH] center2 terms --- source/Makefile.Objects | 1 + .../hamilt_lcaodft/CMakeLists.txt | 1 + .../hamilt_lcaodft/FORCE_gamma.cpp | 14 +- .../hamilt_lcaodft/FORCE_k.cpp | 13 +- .../module_hamilt_lcao/hamilt_lcaodft/edm.cpp | 2 +- .../hamilt_lcaodft/force_stress_arrays.h | 4 + .../pulay_force_stress/pulay_force_stress.h | 57 ++++--- .../pulay_force_stress/pulay_force_stress.hpp | 140 ------------------ .../pulay_force_stress_center2.cpp | 103 +++++++++++++ .../pulay_force_stress_center2_template.hpp | 127 ++++++++++++++++ .../pulay_force_stress_gint.hpp | 39 +++++ 11 files changed, 337 insertions(+), 164 deletions(-) delete mode 100644 source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress.hpp create mode 100644 source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress_center2.cpp create mode 100644 source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress_center2_template.hpp create mode 100644 source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress_gint.hpp diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 66b0d9e261..51f4181913 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -570,6 +570,7 @@ OBJS_LCAO=evolve_elec.o\ FORCE_k.o\ stress_tools.o\ edm.o\ + pulay_force_stress_center2.o\ fedm_gamma.o\ fedm_k.o\ ftvnl_dphi_gamma.o\ diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt b/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt index 8ccb4c9f6e..0bf170672f 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt +++ b/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt @@ -15,6 +15,7 @@ if(ENABLE_LCAO) operator_lcao/td_nonlocal_lcao.cpp operator_lcao/sc_lambda_lcao.cpp operator_lcao/dftu_lcao.cpp + pulay_force_stress/pulay_force_stress_center2.cpp FORCE_STRESS.cpp FORCE_gamma.cpp FORCE_k.cpp diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp index 566ebbb6d0..2c4b45c312 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp @@ -216,12 +216,20 @@ void Force_LCAO::ftable(const bool isforce, // allocate DHloc_fixed_x, DHloc_fixed_y, DHloc_fixed_z this->allocate(pv, fsr, two_center_bundle, orb); + const double* dSx[3] = { fsr.DSloc_x, fsr.DSloc_y, fsr.DSloc_z }; + const double* dSxy[6] = { fsr.DSloc_11, fsr.DSloc_12, fsr.DSloc_13, fsr.DSloc_22, fsr.DSloc_23, fsr.DSloc_33 }; // calculate the force related to 'energy density matrix'. - this->cal_fedm(isforce, isstress, fsr, ucell, + // this->cal_fedm(isforce, isstress, fsr, ucell, + // this->cal_edm(pelec, *psi, *dm, *kv, pv, GlobalV::NSPIN, GlobalV::NBANDS, ucell, *ra), + // psi, pv, pelec, foverlap, soverlap); + PulayForceStress::cal_pulay_fs(foverlap, soverlap, this->cal_edm(pelec, *psi, *dm, *kv, pv, GlobalV::NSPIN, GlobalV::NBANDS, ucell, *ra), - psi, pv, pelec, foverlap, soverlap); + ucell, pv, dSx, dSxy, isforce, isstress); - this->cal_ftvnl_dphi(dm, pv, ucell, fsr, isforce, isstress, ftvnl_dphi, stvnl_dphi); + const double* dHx[3] = { fsr.DHloc_fixed_x, fsr.DHloc_fixed_y, fsr.DHloc_fixed_z }; + const double* dHxy[6] = { fsr.DHloc_fixed_11, fsr.DHloc_fixed_12, fsr.DHloc_fixed_13, fsr.DHloc_fixed_22, fsr.DHloc_fixed_23, fsr.DHloc_fixed_33 }; + // this->cal_ftvnl_dphi(dm, pv, ucell, fsr, isforce, isstress, ftvnl_dphi, stvnl_dphi); + PulayForceStress::cal_pulay_fs(ftvnl_dphi, stvnl_dphi, *dm, ucell, pv, dHx, dHxy, isforce, isstress); this->cal_fvnl_dbeta(dm, pv, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp index d22d6947c1..c96cbf708d 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp @@ -313,13 +313,20 @@ void Force_LCAO>::ftable(const bool isforce, kv->get_nks(), kv->kvec_d); + const double* dSx[3] = { fsr.DSloc_Rx, fsr.DSloc_Ry, fsr.DSloc_Rz }; // calculate the energy density matrix // and the force related to overlap matrix and energy density matrix. - this->cal_fedm(isforce, isstress, fsr, ucell, + // this->cal_fedm(isforce, isstress, fsr, ucell, + // this->cal_edm(pelec, *psi, *dm, *kv, pv, GlobalV::NSPIN, GlobalV::NBANDS, ucell, *ra), + // psi, pv, pelec, foverlap, soverlap, kv, ra); + PulayForceStress::cal_pulay_fs(foverlap, soverlap, this->cal_edm(pelec, *psi, *dm, *kv, pv, GlobalV::NSPIN, GlobalV::NBANDS, ucell, *ra), - psi, pv, pelec, foverlap, soverlap, kv, ra); + ucell, pv, dSx, fsr.DH_r, isforce, isstress, ra, -1.0, 1.0); - this->cal_ftvnl_dphi(dm, pv, ucell, fsr, isforce, isstress, ftvnl_dphi, stvnl_dphi, ra); + const double* dHx[3] = { fsr.DHloc_fixedR_x, fsr.DHloc_fixedR_y, fsr.DHloc_fixedR_z }; // T+Vnl + const double* dHxy[6] = { fsr.stvnl11, fsr.stvnl12, fsr.stvnl13, fsr.stvnl22, fsr.stvnl23, fsr.stvnl33 }; //T + // this->cal_ftvnl_dphi(dm, pv, ucell, fsr, isforce, isstress, ftvnl_dphi, stvnl_dphi, ra); + PulayForceStress::cal_pulay_fs(ftvnl_dphi, stvnl_dphi, *dm, ucell, pv, dHx, dHxy, isforce, isstress, ra, 1.0, -1.0); // doing on the real space grid. // this->cal_fvl_dphi(isforce, isstress, pelec->pot, gint, fvl_dphi, svl_dphi); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp index 8aecad21af..064e4fa52a 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp @@ -97,6 +97,6 @@ elecstate::DensityMatrix, double> Force_LCAO struct TGint { using type = Gint_Gamma; }; template <> struct TGint> { using type = Gint_k; }; #endif -/// calculate $Tr[D*dH/dx]$ and $1/V Tr[D*(dH/dx_a*x_b)] -/// where D can be either density matrix or energy density matrix +/// calculate the abstract formulas: +/// $Tr[D*dH/dx]$ (force) and $1/V Tr[D*(dH/dx_a*x_b)]$ (stress) +/// where D can be any (energy) density matrix +/// and H can be any operator namespace PulayForceStress { - /// for 2-center terms - // template - // void cal_pulay_fs( - // ModuleBase::matrix& f, ///< [out] force - // ModuleBase::matrix& s, ///< [out] stress - // const elecstate::DensityMatrix& dm, ///< [in] density matrix or energy density matrix - // const UnitCell& ucell, ///< [in] unit cell - // const ForceStressArrays& fsr, - // const bool& isstress - // ); - /// for grid terms + /// for 2-center-integration terms, provided force and stress derivatives + template + void cal_pulay_fs( + ModuleBase::matrix& f, ///< [out] force + ModuleBase::matrix& s, ///< [out] stress + const elecstate::DensityMatrix& dm, ///< [in] density matrix or energy density matrix + const UnitCell& ucell, ///< [in] unit cell + const Parallel_Orbitals& pv, ///< [in] parallel orbitals + const double* (&dHSx)[3], ///< [in] dHSx x, y, z, for force + const double* (&dHSxy)[6], ///< [in] dHSxy 11, 12, 13, 22, 23, 33, for stress + const bool& isforce, + const bool& isstress, + Record_adj* ra = nullptr, + const double& factor_force = 1.0, + const double& factor_stress = 1.0); + + /// for 2-center-integration terms, provided force derivatives and coordinate difference + template + void cal_pulay_fs( + ModuleBase::matrix& f, ///< [out] force + ModuleBase::matrix& s, ///< [out] stress + const elecstate::DensityMatrix& dm, ///< [in] density matrix or energy density matrix + const UnitCell& ucell, ///< [in] unit cell + const Parallel_Orbitals& pv, ///< [in] parallel orbitals + const double* (&dHSx)[3], ///< [in] dHSx x, y, z, for force and stress + const double* dtau, ///< [in] dr x, y, z, for stress + const bool& isforce, + const bool& isstress, + Record_adj* ra = nullptr, + const double& factor_force = 1.0, + const double& factor_stress = 1.0); + + /// for grid-integration terms template void cal_pulay_fs( ModuleBase::matrix& f, ///< [out] force @@ -39,8 +63,7 @@ namespace PulayForceStress typename TGint::type& gint, ///< [in] Gint object const bool& isforce, const bool& isstress, - const bool& set_dmr_gint = true - ); + const bool& set_dmr_gint = true); } - -#include "pulay_force_stress.hpp" \ No newline at end of file +#include "pulay_force_stress_center2_template.hpp" +#include "pulay_force_stress_gint.hpp" \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress.hpp deleted file mode 100644 index b51ec6552c..0000000000 --- a/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress.hpp +++ /dev/null @@ -1,140 +0,0 @@ -#pragma once -#include "pulay_force_stress.h" -#include "module_hamilt_lcao/hamilt_lcaodft/stress_tools.h" -namespace PulayForceStress -{ - // inline std::vector single_derivative( - // const TwoCenterBundle& c2, - // const char& dtype, - // const int& T1, const int& L1, const int& N1, const int& m1, - // const int& T2, const int& L2, const int& N2, const int& m2, - // const ModuleBase::Vector3& dtau_lat0) - // { - // std::vector dhs_xyz(3); - // // convert m (0,1,...2l) to M (-l, -l+1, ..., l-1, l) - // const int M1 = (m1 % 2 == 0) ? -m1 / 2 : (m1 + 1) / 2; - // const int M2 = (m2 % 2 == 0) ? -m2 / 2 : (m2 + 1) / 2; - // switch (dtype) - // { - // case 'S': - // two_center_bundle.overlap_orb->calculate(T1, L1, N1, M1, T2, L2, N2, M2, dtau_lat0, nullptr, dhs_xyz.data()); - // break; - // case 'T': - // two_center_bundle.kinetic_orb->calculate(T1, L1, N1, M1, T2, L2, N2, M2, dtau_lat0, nullptr, dhs_xyz.data()); - // break; - // // how about nonlocal? - // default: // not supposed to happen - // ModuleBase::WARNING_QUIT("LCAO_domain::build_ST_new", "dtype must be S or T"); - // } - // } - - // inline std::tuple>, std::vector>>> get_dHS - // ( - // const TwoCenterIntegrator& c2, - // const char& dtype, - // const bool& isstress - // ) // replace build_ST_new - // { - // std::vector> dHSx; - // std::vector>> dHSxy; - // return std::make_tuple(dHSx, dHSxy); - // } - // template - // void cal_pulay_fs( - // ModuleBase::matrix& f, ///< [out] force - // ModuleBase::matrix& s, ///< [out] stress - // const elecstate::DensityMatrix& dm, ///< [in] density matrix - // const UnitCell& ucell, ///< [in] unit cell - // const TwoCenterBundle& c2, - // const bool& isstress) - // { - // const int nspin = GlobalV::NSPIN; - // const int nbasis = GlobalV::NLOCAL; - // // any better container? - // const auto& dHSs = get_dHSx(c2); - // std::vector> dHSx = std::get<0>(dHSs); - // std::vector>> dHSxy = std::get<1>(dHSs); - // for (int i = 0; i < nbasis; i++) - // { - // const int iat = ucell.iwt2iat[i]; - // for (int j = 0; j < nbasis; j++) - // { - // const int mu = pv.global2local_row(j); - // const int nu = pv.global2local_col(i); - - // if (mu >= 0 && nu >= 0) - // { - // const int index = mu * pv.ncol + nu; - // double sum = 0.0; - // for (int is = 0; is < nspin; ++is) - // { - // sum += dm.get_DMK(is + 1, 0, nu, mu); - // } - // for (int i = 0;i < 3;++i) { f(iat, i) += 2 * sum * dHSx[i][index]; } - // if (isstress) - // { - // for (int i = 0;i < 3;++i) { for (int j = i;j < 3;++j) { s(i, j) += 2 * sum * dHSxy[i][j][index]; } } - // } - // } - // } - // } - - // if (isstress) - // { - // StressTools::stress_fill(ucell.lat0, ucell.omega, s); - // } - // } - template - void cal_pulay_fs( - ModuleBase::matrix& f, ///< [out] force - ModuleBase::matrix& s, ///< [out] stress - const elecstate::DensityMatrix& dm, ///< [in] density matrix - const UnitCell& ucell, ///< [in] unit cell - const elecstate::Potential* pot, ///< [in] potential on grid - typename TGint::type& gint, - const bool& isforce, - const bool& isstress, - const bool& set_dmr_gint) - { - if (set_dmr_gint) { gint.transfer_DM2DtoGrid(dm.get_DMR_vector()); } // 2d block to grid - - const int nspin = GlobalV::NSPIN; - for (int is = 0; is < nspin; ++is) - { - const double* vr_eff1 = pot->get_effective_v(is); - const double* vofk_eff1 = nullptr; - - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) - { - vofk_eff1 = pot->get_effective_vofk(is); - - Gint_inout inout( - is, - vr_eff1, - vofk_eff1, - isforce, - isstress, - &f, - &s, - Gint_Tools::job_type::force_meta); - - gint.cal_gint(&inout); - } - else - { - Gint_inout inout( - is, - vr_eff1, - isforce, - isstress, - &f, - &s, - Gint_Tools::job_type::force); - - gint.cal_gint(&inout); - } - } - - if (isstress) { StressTools::stress_fill(-1.0, ucell.omega, s); } - } -} \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress_center2.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress_center2.cpp new file mode 100644 index 0000000000..04c89a9163 --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress_center2.cpp @@ -0,0 +1,103 @@ +#include "pulay_force_stress.h" +namespace PulayForceStress +{ + template<> // gamma-only, provided xy + void cal_pulay_fs( + ModuleBase::matrix& f, + ModuleBase::matrix& s, + const elecstate::DensityMatrix& dm, + const UnitCell& ucell, + const Parallel_Orbitals& pv, + const double* (&dHSx)[3], + const double* (&dHSxy)[6], + const bool& isforce, + const bool& isstress, + Record_adj* ra, + const double& factor_force, + const double& factor_stress) + { + ModuleBase::TITLE("Force_LCAO", "cal_pulay_fs_center2"); + ModuleBase::timer::tick("Force_LCAO", "cal_pulay_fs_center2"); + + const int nspin = GlobalV::NSPIN; + const int nlocal = GlobalV::NLOCAL; + + for (int i = 0; i < nlocal; ++i) + { + const int iat = ucell.iwt2iat[i]; + for (int j = 0; j < nlocal; ++j) + { + const int mu = pv.global2local_row(j); + const int nu = pv.global2local_col(i); + + if (mu >= 0 && nu >= 0) + { + const int index = mu * pv.ncol + nu; + double sum = 0.0; + for (int is = 0; is < nspin; ++is) { sum += dm.get_DMK(is + 1, 0, nu, mu); } + if (isforce) + { + const double sumf = sum * factor_force; + for (int i = 0; i < 3; ++i) { f(iat, i) += sumf * 2.0 * dHSx[i][index]; } + } + if (isstress) + { + const double sums = sum * factor_stress; + int ij = 0; + for (int i = 0; i < 3;++i) { for (int j = i; j < 3; ++j) { s(i, j) += sums * dHSxy[ij++][index]; } } + } + } + } + } + + if (isstress) { StressTools::stress_fill(ucell.lat0, ucell.omega, s); } + + ModuleBase::timer::tick("Force_LCAO", "cal_pulay_fs_center2"); + } + + template<> //multi-k, provided xy + void cal_pulay_fs( + ModuleBase::matrix& f, + ModuleBase::matrix& s, + const elecstate::DensityMatrix, double>& dm, + const UnitCell& ucell, + const Parallel_Orbitals& pv, + const double* (&dHSx)[3], + const double* (&dHSxy)[6], + const bool& isforce, + const bool& isstress, + Record_adj* ra, + const double& factor_force, + const double& factor_stress) + { + auto stress_func = [](ModuleBase::matrix& local_s, const double& dm2d1_s, const double** dHSx, const double** dHSxy, const double* dtau, const int& irr) + { + int ij = 0; + for (int i = 0; i < 3; ++i) { for (int j = i; j < 3; ++j) { local_s(i, j) += dm2d1_s * dHSxy[ij++][irr]; } } + }; + cal_pulay_fs(f, s, dm, ucell, pv, dHSx, dHSxy, nullptr, isforce, isstress, ra, factor_force, factor_stress, stress_func); + } + + template<> // multi-k, provided x + void cal_pulay_fs( + ModuleBase::matrix& f, + ModuleBase::matrix& s, + const elecstate::DensityMatrix, double>& dm, + const UnitCell& ucell, + const Parallel_Orbitals& pv, + const double* (&dHSx)[3], + const double* dtau, + const bool& isforce, + const bool& isstress, + Record_adj* ra, + const double& factor_force, + const double& factor_stress) + { + auto stress_func = [](ModuleBase::matrix& local_s, const double& dm2d1_s, const double** dHSx, const double** dHSxy, const double* dtau, const int& irr) + { + for (int i = 0; i < 3; ++i) { for (int j = i; j < 3; ++j) { local_s(i, j) += dm2d1_s * dHSx[i][irr] * dtau[irr * 3 + j]; } } + }; + cal_pulay_fs(f, s, dm, ucell, pv, dHSx, nullptr, dtau, isforce, isstress, ra, factor_force, factor_stress, stress_func); + } + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress_center2_template.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress_center2_template.hpp new file mode 100644 index 0000000000..ba1345d90a --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress_center2_template.hpp @@ -0,0 +1,127 @@ +#pragma once +#include "pulay_force_stress.h" +#include "module_base/timer.h" +namespace PulayForceStress +{ + // common kernel + template + inline void cal_pulay_fs( + ModuleBase::matrix& f, + ModuleBase::matrix& s, + const elecstate::DensityMatrix& dm, + const UnitCell& ucell, + const Parallel_Orbitals& pv, + const double** dHSx, + const double** dHSxy, + const double* dtau, + const bool& isforce, + const bool& isstress, + Record_adj* ra, + const double& factor_force, + const double& factor_stress, + Tfunc& stress_func) + { + ModuleBase::TITLE("Force_LCAO", "cal_pulay_fs_center2"); + ModuleBase::timer::tick("Force_LCAO", "cal_pulay_fs_center2"); + + const int nspin_DMR = (GlobalV::NSPIN == 2) ? 2 : 1; + int total_irr = 0; +#ifdef _OPENMP +#pragma omp parallel + { + int num_threads = omp_get_num_threads(); + ModuleBase::matrix local_s(3, 3); + int local_total_irr = 0; +#else + ModuleBase::matrix& local_s = s; + int& local_total_irr = total_irr; +#endif + +#ifdef _OPENMP +#pragma omp for schedule(dynamic) +#endif + for (int iat = 0; iat < ucell.nat; iat++) + { + const int T1 = ucell.iat2it[iat]; + Atom* atom1 = &ucell.atoms[T1]; + const int I1 = ucell.iat2ia[iat]; + // get iat1 + int iat1 = ucell.itia2iat(T1, I1); + double* f_iat; + if (isforce) { f_iat = &f(iat, 0); } +#ifdef _OPENMP + // using local stack to avoid false sharing in multi-threaded case + double f_tmp[3] = { 0.0, 0.0, 0.0 }; + if (num_threads > 1) { f_iat = f_tmp; } +#endif + int irr = pv.nlocstart[iat]; + const int start1 = ucell.itiaiw2iwt(T1, I1, 0); + for (int cb = 0; cb < ra->na_each[iat]; ++cb) + { + const int T2 = ra->info[iat][cb][3]; + const int I2 = ra->info[iat][cb][4]; + const int start2 = ucell.itiaiw2iwt(T2, I2, 0); + Atom* atom2 = &ucell.atoms[T2]; + // get iat2 + int iat2 = ucell.itia2iat(T2, I2); + double Rx = ra->info[iat][cb][0]; + double Ry = ra->info[iat][cb][1]; + double Rz = ra->info[iat][cb][2]; + // get BaseMatrix + if (pv.get_row_size(iat1) <= 0 || pv.get_col_size(iat2) <= 0) { continue; } + std::vector*> tmp_matrix; + for (int is = 0; is < nspin_DMR; ++is) + { + tmp_matrix.push_back(dm.get_DMR_pointer(is + 1)->find_matrix(iat1, iat2, Rx, Ry, Rz)); + } + for (int mu = 0; mu < pv.get_row_size(iat1); ++mu) + { + for (int nu = 0; nu < pv.get_col_size(iat2); ++nu) + { + // the DMR should not be summed over spin, do the summation here + double dm2d1 = 0.0; + for (int is = 0; is < nspin_DMR; ++is) { dm2d1 += tmp_matrix[is]->get_value(mu, nu); } + double dm2d2 = 2.0 * dm2d1; + if (isforce) + { + const double dm2d2_f = dm2d2 * factor_force; + for (int i = 0; i < 3; ++i) { f_iat[i] += dm2d2_f * dHSx[i][irr]; } + } + if (isstress) + { + const double dm2d1_s = dm2d1 * factor_stress; + stress_func(local_s, dm2d1_s, dHSx, dHSxy, dtau, irr); + } + ++local_total_irr; + ++irr; + } + } + } +#ifdef _OPENMP + if (isforce && num_threads > 1) { for (int i = 0; i < 3; ++i) { f(iat, i) += f_iat[i]; } } +#endif + } // end iat +#ifdef _OPENMP +#pragma omp critical(cal_foverlap_k_reduce) + { + total_irr += local_total_irr; + if (isstress) + { + for (int i = 0; i < 3; ++i) { for (int j = i; j < 3; ++j) { s(i, j) += local_s(i, j); } } + } + } + } +#endif + + if (total_irr != pv.nnr) + { + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "wrong irr", total_irr); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "wrong pv.nnr", pv.nnr); + ModuleBase::WARNING_QUIT("Force_LCAO::cal_pulay_fs_center2", "irr!=pv.nnr"); + } + + if (isstress) { StressTools::stress_fill(ucell.lat0, ucell.omega, s); } + + ModuleBase::timer::tick("Force_LCAO", "cal_pulay_fs_center2"); + } +} diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress_gint.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress_gint.hpp new file mode 100644 index 0000000000..5b29d0e269 --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress_gint.hpp @@ -0,0 +1,39 @@ +#pragma once +#include "pulay_force_stress.h" +#include "module_hamilt_lcao/hamilt_lcaodft/stress_tools.h" +#include "module_hamilt_general/module_xc/xc_functional.h" +namespace PulayForceStress +{ + template + void cal_pulay_fs( + ModuleBase::matrix& f, ///< [out] force + ModuleBase::matrix& s, ///< [out] stress + const elecstate::DensityMatrix& dm, ///< [in] density matrix + const UnitCell& ucell, ///< [in] unit cell + const elecstate::Potential* pot, ///< [in] potential on grid + typename TGint::type& gint, + const bool& isforce, + const bool& isstress, + const bool& set_dmr_gint) + { + if (set_dmr_gint) { gint.transfer_DM2DtoGrid(dm.get_DMR_vector()); } // 2d block to grid + const int nspin = GlobalV::NSPIN; + for (int is = 0; is < nspin; ++is) + { + const double* vr_eff1 = pot->get_effective_v(is); + const double* vofk_eff1 = nullptr; + if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + { + vofk_eff1 = pot->get_effective_vofk(is); + Gint_inout inout(is, vr_eff1, vofk_eff1, isforce, isstress, &f, &s, Gint_Tools::job_type::force_meta); + gint.cal_gint(&inout); + } + else + { + Gint_inout inout(is, vr_eff1, isforce, isstress, &f, &s, Gint_Tools::job_type::force); + gint.cal_gint(&inout); + } + } + if (isstress) { StressTools::stress_fill(-1.0, ucell.omega, s); } + } +} \ No newline at end of file