diff --git a/source/Makefile.Objects b/source/Makefile.Objects index d99287243d..f6ed91aac4 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -568,12 +568,9 @@ OBJS_LCAO=evolve_elec.o\ FORCE_STRESS.o\ FORCE_gamma.o\ FORCE_k.o\ - fvl_dphi_gamma.o\ - fvl_dphi_k.o\ - fedm_gamma.o\ - fedm_k.o\ - ftvnl_dphi_gamma.o\ - ftvnl_dphi_k.o\ + stress_tools.o\ + edm.o\ + pulay_force_stress_center2.o\ fvnl_dbeta_gamma.o\ fvnl_dbeta_k.o\ grid_init.o\ diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt b/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt index bc72dcca3b..72209f066e 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt +++ b/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt @@ -15,15 +15,12 @@ if(ENABLE_LCAO) operator_lcao/td_nonlocal_lcao.cpp operator_lcao/sc_lambda_lcao.cpp operator_lcao/dftu_lcao.cpp + pulay_force_stress_center2.cpp FORCE_STRESS.cpp FORCE_gamma.cpp FORCE_k.cpp - fvl_dphi_gamma.cpp - fvl_dphi_k.cpp - fedm_gamma.cpp - fedm_k.cpp - ftvnl_dphi_gamma.cpp - ftvnl_dphi_k.cpp + stress_tools.cpp + edm.cpp fvnl_dbeta_gamma.cpp fvnl_dbeta_k.cpp grid_init.cpp diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h index 8644bf804b..51d840089d 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h @@ -94,7 +94,7 @@ class Force_LCAO const bool isstress, ForceStressArrays& fsr, const UnitCell& ucell, - const elecstate::DensityMatrix* dm, + const elecstate::DensityMatrix& dm, const psi::Psi* psi, const Parallel_Orbitals& pv, const elecstate::ElecState* pelec, @@ -136,12 +136,16 @@ class Force_LCAO typename TGint::type& gint, ModuleBase::matrix& fvl_dphi, ModuleBase::matrix& svl_dphi); + + elecstate::DensityMatrix cal_edm(const elecstate::ElecState* pelec, + const psi::Psi& psi, + const elecstate::DensityMatrix& dm, + const K_Vectors& kv, + const Parallel_Orbitals& pv, + const int& nspin, + const int& nbands, + const UnitCell& ucell, + Record_adj& ra) const; }; -// this namespace used to store global function for some stress operation -namespace StressTools -{ -// set upper matrix to whole matrix -void stress_fill(const double& lat0_, const double& omega_, ModuleBase::matrix& stress_matrix); -} // namespace StressTools #endif diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp index 9f72bf544c..db631841cd 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp @@ -13,6 +13,7 @@ #include "module_elecstate/elecstate_lcao.h" #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" #include "module_io/write_HS.h" +#include "module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress.h" template <> void Force_LCAO::allocate(const Parallel_Orbitals& pv, @@ -215,10 +216,17 @@ 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, dm, psi, pv, pelec, foverlap, soverlap); + PulayForceStress::cal_pulay_fs(foverlap, soverlap, + this->cal_edm(pelec, *psi, *dm, *kv, pv, PARAM.inp.nspin, PARAM.inp.nbands, ucell, *ra), + 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 }; + //tvnl_dphi + PulayForceStress::cal_pulay_fs(ftvnl_dphi, stvnl_dphi, *dm, ucell, pv, dHx, dHxy, isforce, isstress); this->cal_fvnl_dbeta(dm, pv, @@ -231,9 +239,9 @@ void Force_LCAO::ftable(const bool isforce, fvnl_dbeta, svnl_dbeta); - this->cal_fvl_dphi(isforce, isstress, pelec->pot, gint, fvl_dphi, svl_dphi); + // vl_dphi + PulayForceStress::cal_pulay_fs(fvl_dphi, svl_dphi, *dm, ucell, pelec->pot, gint, isforce, isstress, false/*reset dm to gint*/); - // caoyu add for DeePKS #ifdef __DEEPKS if (PARAM.inp.deepks_scf) { @@ -314,23 +322,3 @@ void Force_LCAO::ftable(const bool isforce, ModuleBase::timer::tick("Force_LCAO", "ftable"); return; } - -namespace StressTools -{ -void stress_fill(const double& lat0_, const double& omega_, ModuleBase::matrix& stress_matrix) -{ - assert(omega_ > 0.0); - double weight = lat0_ / omega_; - for (int i = 0; i < 3; ++i) - { - for (int j = 0; j < 3; ++j) - { - if (j > i) - { - stress_matrix(j, i) = stress_matrix(i, j); - } - stress_matrix(i, j) *= weight; - } - } -} -} // namespace StressTools diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp index 81e10d3bbf..1c4b10f251 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp @@ -12,6 +12,7 @@ #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/write_HS.h" +#include "module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress.h" #include #include @@ -312,14 +313,21 @@ 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, dm, psi, pv, pelec, foverlap, soverlap, kv, ra); + PulayForceStress::cal_pulay_fs(foverlap, soverlap, + this->cal_edm(pelec, *psi, *dm, *kv, pv, PARAM.inp.nspin, PARAM.inp.nbands, ucell, *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 + // tvnl_dphi + 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); + // vl_dphi + PulayForceStress::cal_pulay_fs(fvl_dphi, svl_dphi, *dm, ucell, pelec->pot, gint, isforce, isstress, false/*reset dm to gint*/); this->cal_fvnl_dbeta(dm, pv, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp new file mode 100644 index 0000000000..aacaa742dc --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp @@ -0,0 +1,103 @@ +#include "FORCE.h" +#include "module_elecstate/module_dm/cal_dm_psi.h" +#include "module_base/memory.h" +#include "module_parameter/parameter.h" +template<> +elecstate::DensityMatrix Force_LCAO::cal_edm(const elecstate::ElecState* pelec, + const psi::Psi& psi, + const elecstate::DensityMatrix& dm, + const K_Vectors& kv, + const Parallel_Orbitals& pv, + const int& nspin, + const int& nbands, + const UnitCell& ucell, + Record_adj& ra) const +{ + ModuleBase::matrix wg_ekb; + wg_ekb.create(nspin, nbands); + + for(int is=0; iswg(is,ib) * pelec->ekb(is, ib); + } + } + + // construct a DensityMatrix for Gamma-Only + elecstate::DensityMatrix edm(&pv, nspin); + +#ifdef __PEXSI + if (PARAM.inp.ks_solver == "pexsi") + { + auto pes = dynamic_cast*>(pelec); + for (int ik = 0; ik < nspin; ik++) + { + edm.set_DMK_pointer(ik, pes->get_DM()->pexsi_EDM[ik]); + } + + } + else +#endif + { + elecstate::cal_dm_psi(edm.get_paraV_pointer(), wg_ekb, psi, edm); + } + return edm; +} + +template<> +elecstate::DensityMatrix, double> Force_LCAO>::cal_edm(const elecstate::ElecState* pelec, + const psi::Psi>& psi, + const elecstate::DensityMatrix, double>& dm, + const K_Vectors& kv, + const Parallel_Orbitals& pv, + const int& nspin, + const int& nbands, + const UnitCell& ucell, + Record_adj& ra) const +{ + + // construct a DensityMatrix object + elecstate::DensityMatrix, double> edm(&kv, &pv, nspin); + + //-------------------------------------------- + // calculate the energy density matrix here. + //-------------------------------------------- + + ModuleBase::matrix wg_ekb; + wg_ekb.create(kv.get_nks(), nbands); + ModuleBase::Memory::record("Force::wg_ekb", sizeof(double) * kv.get_nks() * nbands); +#ifdef _OPENMP +#pragma omp parallel for collapse(2) schedule(static, 1024) +#endif + for (int ik = 0; ik < kv.get_nks(); ik++) + { + for (int ib = 0; ib < nbands; ib++) + { + wg_ekb(ik, ib) = pelec->wg(ik, ib) * pelec->ekb(ik, ib); + } + } + + // use the original formula (Hamiltonian matrix) to calculate energy density matrix + if (dm.EDMK.size()) + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (int ik = 0; ik < kv.get_nks(); ++ik) + { + edm.set_DMK_pointer(ik, dm.EDMK[ik].c); + } + } + else + { + // cal_dm_psi + elecstate::cal_dm_psi(edm.get_paraV_pointer(), wg_ekb, psi, edm); + } + + // cal_dm_2d + edm.init_DMR(ra, &ucell); + edm.cal_DMR(); + // edm.sum_DMR_spin(); + return edm; +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/fedm_gamma.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/fedm_gamma.cpp deleted file mode 100644 index 6fe6420d6d..0000000000 --- a/source/module_hamilt_lcao/hamilt_lcaodft/fedm_gamma.cpp +++ /dev/null @@ -1,112 +0,0 @@ -#include "FORCE.h" -#include "module_elecstate/elecstate_lcao.h" -#include "module_parameter/parameter.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" -#include "module_base/parallel_reduce.h" -#include "module_base/timer.h" -#include "module_psi/psi.h" -#include "module_elecstate/cal_dm.h" -#include "module_elecstate/module_dm/cal_dm_psi.h" - -// force due to the overlap matrix. -// need energy density matrix here. -template<> -void Force_LCAO::cal_fedm( - const bool isforce, - const bool isstress, - ForceStressArrays &fsr, // mohan add 2024-06-16 - const UnitCell& ucell, - const elecstate::DensityMatrix* dm, - const psi::Psi* psi, - const Parallel_Orbitals& pv, - const elecstate::ElecState *pelec, - ModuleBase::matrix &foverlap, - ModuleBase::matrix& soverlap, - const K_Vectors* kv, - Record_adj* ra) -{ - ModuleBase::TITLE("Force_LCAO","cal_fedm"); - ModuleBase::timer::tick("Force_LCAO","cal_fedm"); - - const int nspin = PARAM.inp.nspin; - const int nbands = PARAM.inp.nbands; - const int nlocal = PARAM.globalv.nlocal; - - ModuleBase::matrix wg_ekb; - wg_ekb.create(nspin, nbands); - - for(int is=0; iswg(is,ib) * pelec->ekb(is, ib); - } - } - - // construct a DensityMatrix for Gamma-Only - elecstate::DensityMatrix edm(&pv, nspin); - -#ifdef __PEXSI - if (PARAM.inp.ks_solver == "pexsi") - { - auto pes = dynamic_cast*>(pelec); - for (int ik = 0; ik < nspin; ik++) - { - edm.set_DMK_pointer(ik, pes->get_DM()->pexsi_EDM[ik]); - } - - } - else -#endif - { - elecstate::cal_dm_psi(edm.get_paraV_pointer(), wg_ekb, psi[0], edm); - } - - - for(int i=0; i=0 && nu>=0) - { - const int index = mu * pv.ncol + nu; - double sum = 0.0; - for(int is=0; is -#include - -#ifdef __DEEPKS -#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" -#endif - -#ifdef _OPENMP -#include -#endif - -#include "module_hamilt_lcao/hamilt_lcaodft/record_adj.h" - -template <> -void Force_LCAO>::cal_fedm(const bool isforce, - const bool isstress, - ForceStressArrays& fsr, - const UnitCell& ucell, - const elecstate::DensityMatrix, double>* dm, - const psi::Psi>* psi, - const Parallel_Orbitals& pv, - const elecstate::ElecState* pelec, - ModuleBase::matrix& foverlap, - ModuleBase::matrix& soverlap, - const K_Vectors* kv, - Record_adj* ra) -{ - ModuleBase::TITLE("Force_LCAO", "cal_fedm"); - ModuleBase::timer::tick("Force_LCAO", "cal_fedm"); - - const int nspin = PARAM.inp.nspin; - const int nbands = PARAM.inp.nbands; - - // construct a DensityMatrix object - elecstate::DensityMatrix, double> edm(kv, &pv, nspin); - - //-------------------------------------------- - // calculate the energy density matrix here. - //-------------------------------------------- - - ModuleBase::matrix wg_ekb; - wg_ekb.create(kv->get_nks(), nbands); - ModuleBase::Memory::record("Force::wg_ekb", sizeof(double) * kv->get_nks() * nbands); -#ifdef _OPENMP -#pragma omp parallel for collapse(2) schedule(static, 1024) -#endif - for (int ik = 0; ik < kv->get_nks(); ik++) - { - for (int ib = 0; ib < nbands; ib++) - { - wg_ekb(ik, ib) = pelec->wg(ik, ib) * pelec->ekb(ik, ib); - } - } - std::vector edm_k(kv->get_nks()); - - // use the original formula (Hamiltonian matrix) to calculate energy density matrix - if (dm->EDMK.size()) - { -#ifdef _OPENMP -#pragma omp parallel for schedule(static) -#endif - for (int ik = 0; ik < kv->get_nks(); ++ik) - { - edm.set_DMK_pointer(ik, dm->EDMK[ik].c); - } - } - else - { - // cal_dm_psi - elecstate::cal_dm_psi(edm.get_paraV_pointer(), wg_ekb, psi[0], edm); - } - - // cal_dm_2d - edm.init_DMR(*ra, &ucell); - edm.cal_DMR(); - edm.sum_DMR_spin(); - //-------------------------------------------- - // summation \sum_{i,j} E(i,j)*dS(i,j) - // BEGIN CALCULATION OF FORCE OF EACH ATOM - //-------------------------------------------- - int total_irr = 0; -#ifdef _OPENMP -#pragma omp parallel - { - int num_threads = omp_get_num_threads(); - ModuleBase::matrix local_soverlap(3, 3); - int local_total_irr = 0; -#else - ModuleBase::matrix& local_soverlap = soverlap; - int& local_total_irr = total_irr; -#endif - - ModuleBase::Vector3 tau1, dtau, tau2; - -#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* foverlap_iat; - if (isforce) - { - foverlap_iat = &foverlap(iat, 0); - } - -#ifdef _OPENMP - // using local stack to avoid false sharing in multi-threaded case - double foverlap_temp[3] = {0.0, 0.0, 0.0}; - if (num_threads > 1) - { - foverlap_iat = foverlap_temp; - } -#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 - hamilt::BaseMatrix* tmp_matrix = edm.get_DMR_pointer(1)->find_matrix(iat1, iat2, Rx, Ry, Rz); - if (tmp_matrix == nullptr) - { - continue; - } - int row_ap = pv.atom_begin_row[iat1]; - int col_ap = pv.atom_begin_col[iat2]; - // get DMR - for (int mu = 0; mu < pv.get_row_size(iat1); ++mu) - { - for (int nu = 0; nu < pv.get_col_size(iat2); ++nu) - { - // here do not sum over spin due to edm.sum_DMR_spin(); - double edm2d1 = tmp_matrix->get_value(mu, nu); - double edm2d2 = 2.0 * edm2d1; - - if (isforce) - { - foverlap_iat[0] -= edm2d2 * fsr.DSloc_Rx[irr]; - foverlap_iat[1] -= edm2d2 * fsr.DSloc_Ry[irr]; - foverlap_iat[2] -= edm2d2 * fsr.DSloc_Rz[irr]; - } - if (isstress) - { - for (int ipol = 0; ipol < 3; ipol++) - { - local_soverlap(0, ipol) += edm2d1 * fsr.DSloc_Rx[irr] * fsr.DH_r[irr * 3 + ipol]; - if (ipol < 1) - { - continue; - } - local_soverlap(1, ipol) += edm2d1 * fsr.DSloc_Ry[irr] * fsr.DH_r[irr * 3 + ipol]; - if (ipol < 2) - { - continue; - } - local_soverlap(2, ipol) += edm2d1 * fsr.DSloc_Rz[irr] * fsr.DH_r[irr * 3 + ipol]; - } - } - //} - ++local_total_irr; - ++irr; - } // end kk - } // end jj - } // end cb -#ifdef _OPENMP - if (isforce && num_threads > 1) - { - foverlap(iat, 0) += foverlap_iat[0]; - foverlap(iat, 1) += foverlap_iat[1]; - foverlap(iat, 2) += foverlap_iat[2]; - } -#endif - } // end iat -#ifdef _OPENMP -#pragma omp critical(cal_foverlap_k_reduce) - { - total_irr += local_total_irr; - if (isstress) - { - for (int ipol = 0; ipol < 3; ipol++) - { - soverlap(0, ipol) += local_soverlap(0, ipol); - if (ipol < 1) - { - continue; - } - soverlap(1, ipol) += local_soverlap(1, ipol); - if (ipol < 2) - { - continue; - } - soverlap(2, ipol) += local_soverlap(2, ipol); - } - } - } - } -#endif - - if (isstress) - { - StressTools::stress_fill(ucell.lat0, ucell.omega, soverlap); - } - - 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::fedm_k", "irr!=pv.nnr"); - } - - ModuleBase::timer::tick("Force_LCAO", "cal_fedm"); - return; -} diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/ftvnl_dphi_gamma.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/ftvnl_dphi_gamma.cpp deleted file mode 100644 index f639f4c511..0000000000 --- a/source/module_hamilt_lcao/hamilt_lcaodft/ftvnl_dphi_gamma.cpp +++ /dev/null @@ -1,73 +0,0 @@ -#include "FORCE.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" -#include "module_parameter/parameter.h" -#include "module_cell/module_neighbor/sltk_grid_driver.h" -#include -#include "module_base/timer.h" - -template<> -void Force_LCAO::cal_ftvnl_dphi( - const elecstate::DensityMatrix* dm, - const Parallel_Orbitals& pv, - const UnitCell& ucell, - ForceStressArrays& fsr, - const bool isforce, - const bool isstress, - ModuleBase::matrix& ftvnl_dphi, - ModuleBase::matrix& stvnl_dphi, - Record_adj* ra) -{ - ModuleBase::TITLE("Force_LCAO","cal_ftvnl_dphi"); - ModuleBase::timer::tick("Force_LCAO","cal_ftvnl_dphi"); - - const int nlocal = PARAM.globalv.nlocal; - const int nspin = PARAM.inp.nspin; - - for(int i=0; i= 0 && nu >= 0 ) - { - const int index = mu * pv.ncol + nu; - //contribution from deriv of AO's in T+VNL term - - double sum = 0.0; - for(int is=0; isget_DMK(is+1, 0, nu, mu); - } - sum *= 2.0; - - if(isforce) - { - ftvnl_dphi(iat,0) += sum * fsr.DHloc_fixed_x[index]; - ftvnl_dphi(iat,1) += sum * fsr.DHloc_fixed_y[index]; - ftvnl_dphi(iat,2) += sum * fsr.DHloc_fixed_z[index]; - } - if(isstress) - { - stvnl_dphi(0,0) += sum/2.0 * fsr.DHloc_fixed_11[index]; - stvnl_dphi(0,1) += sum/2.0 * fsr.DHloc_fixed_12[index]; - stvnl_dphi(0,2) += sum/2.0 * fsr.DHloc_fixed_13[index]; - stvnl_dphi(1,1) += sum/2.0 * fsr.DHloc_fixed_22[index]; - stvnl_dphi(1,2) += sum/2.0 * fsr.DHloc_fixed_23[index]; - stvnl_dphi(2,2) += sum/2.0 * fsr.DHloc_fixed_33[index]; - } - } - } - } - - if(isstress) - { - StressTools::stress_fill(ucell.lat0, ucell.omega, stvnl_dphi); - } - - ModuleBase::timer::tick("Force_LCAO","cal_ftvnl_dphi"); - return; -} diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/ftvnl_dphi_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/ftvnl_dphi_k.cpp deleted file mode 100644 index f8f5a463c3..0000000000 --- a/source/module_hamilt_lcao/hamilt_lcaodft/ftvnl_dphi_k.cpp +++ /dev/null @@ -1,170 +0,0 @@ -#include "FORCE.h" - -#include "module_parameter/parameter.h" -#include -#include - -#include "module_base/memory.h" -#include "module_base/parallel_reduce.h" -#include "module_base/timer.h" -#include "module_base/tool_threading.h" -#include "module_cell/module_neighbor/sltk_grid_driver.h" -#include "module_elecstate/cal_dm.h" -#include "module_elecstate/module_dm/cal_dm_psi.h" -#include "module_elecstate/elecstate_lcao.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" -#include "module_io/write_HS.h" - -#ifdef __DEEPKS -#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" -#endif - -#ifdef _OPENMP -#include -#endif - - -template<> -void Force_LCAO>::cal_ftvnl_dphi( - const elecstate::DensityMatrix, double>* dm, - const Parallel_Orbitals& pv, - const UnitCell& ucell, - ForceStressArrays& fsr, - const bool isforce, - const bool isstress, - ModuleBase::matrix& ftvnl_dphi, - ModuleBase::matrix& stvnl_dphi, - Record_adj* ra) -{ - ModuleBase::TITLE("Force_LCAO", "cal_ftvnl_dphi"); - ModuleBase::timer::tick("Force_LCAO", "cal_ftvnl_dphi"); - - const int nspin_DMR = (PARAM.inp.nspin == 2) ? 2 : 1; - - int total_irr = 0; - // get the adjacent atom's information. - -#ifdef _OPENMP -#pragma omp parallel - { - int num_threads = omp_get_num_threads(); - ModuleBase::matrix local_stvnl_dphi(3, 3); - int local_total_irr = 0; -#pragma omp for schedule(dynamic) -#else - ModuleBase::matrix& local_stvnl_dphi = stvnl_dphi; - int& local_total_irr = total_irr; -#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); - // - int irr = pv.nlocstart[iat]; - const int start1 = ucell.itiaiw2iwt(T1, I1, 0); - double* ftvnl_dphi_iat; - if (isforce) - { - ftvnl_dphi_iat = &ftvnl_dphi(iat, 0); - } -#ifdef _OPENMP - // using local stack to avoid false sharing in multi-threaded case - double ftvnl_dphi_temp[3] = {0.0, 0.0, 0.0}; - if (num_threads > 1) - { - ftvnl_dphi_iat = ftvnl_dphi_temp; - } -#endif - 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) - { - // get value from dm - 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) - { - ftvnl_dphi_iat[0] += dm2d2 * fsr.DHloc_fixedR_x[irr]; - ftvnl_dphi_iat[1] += dm2d2 * fsr.DHloc_fixedR_y[irr]; - ftvnl_dphi_iat[2] += dm2d2 * fsr.DHloc_fixedR_z[irr]; - } - if (isstress) - { - local_stvnl_dphi(0, 0) -= dm2d1 * fsr.stvnl11[irr]; - local_stvnl_dphi(0, 1) -= dm2d1 * fsr.stvnl12[irr]; - local_stvnl_dphi(0, 2) -= dm2d1 * fsr.stvnl13[irr]; - local_stvnl_dphi(1, 1) -= dm2d1 * fsr.stvnl22[irr]; - local_stvnl_dphi(1, 2) -= dm2d1 * fsr.stvnl23[irr]; - local_stvnl_dphi(2, 2) -= dm2d1 * fsr.stvnl33[irr]; - } - //} - ++local_total_irr; - ++irr; - } // end kk - } // end jj - } // end cb -#ifdef _OPENMP - if (isforce && num_threads > 1) - { - ftvnl_dphi(iat, 0) += ftvnl_dphi_iat[0]; - ftvnl_dphi(iat, 1) += ftvnl_dphi_iat[1]; - ftvnl_dphi(iat, 2) += ftvnl_dphi_iat[2]; - } -#endif - } // end iat -#ifdef _OPENMP -#pragma omp critical(cal_ftvnl_dphi_k_reduce) - { - total_irr += local_total_irr; - if (isstress) - { - stvnl_dphi(0, 0) += local_stvnl_dphi(0, 0); - stvnl_dphi(0, 1) += local_stvnl_dphi(0, 1); - stvnl_dphi(0, 2) += local_stvnl_dphi(0, 2); - stvnl_dphi(1, 1) += local_stvnl_dphi(1, 1); - stvnl_dphi(1, 2) += local_stvnl_dphi(1, 2); - stvnl_dphi(2, 2) += local_stvnl_dphi(2, 2); - } - } - } -#endif - assert(total_irr == pv.nnr); - - if (isstress) - { - StressTools::stress_fill(ucell.lat0, ucell.omega, stvnl_dphi); - } - - ModuleBase::timer::tick("Force_LCAO", "cal_ftvnl_dphi"); - return; -} diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/fvl_dphi_gamma.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/fvl_dphi_gamma.cpp deleted file mode 100644 index 5f0c97f63b..0000000000 --- a/source/module_hamilt_lcao/hamilt_lcaodft/fvl_dphi_gamma.cpp +++ /dev/null @@ -1,75 +0,0 @@ -#include "FORCE.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" -#include "module_parameter/parameter.h" -#include "module_base/timer.h" - -template<> -void Force_LCAO::cal_fvl_dphi( - const bool isforce, - const bool isstress, - const elecstate::Potential* pot_in, - TGint::type& gint, - ModuleBase::matrix& fvl_dphi, - ModuleBase::matrix& svl_dphi) -{ - ModuleBase::TITLE("Force_LCAO","cal_fvl_dphi"); - ModuleBase::timer::tick("Force_LCAO","cal_fvl_dphi"); - - const int nspin = PARAM.inp.nspin; - - fvl_dphi.zero_out(); - svl_dphi.zero_out(); - - for(int is=0; isget_effective_v(is); - const double* vofk_eff1 = nullptr; - - if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5) - { - vofk_eff1 = pot_in->get_effective_vofk(is); - - Gint_inout inout( - is, - vr_eff1, - vofk_eff1, - isforce, - isstress, - &fvl_dphi, - &svl_dphi, - Gint_Tools::job_type::force_meta); - - gint.cal_gint(&inout); - } - else - { - Gint_inout inout( - is, - vr_eff1, - isforce, - isstress, - &fvl_dphi, - &svl_dphi, - Gint_Tools::job_type::force); - - gint.cal_gint(&inout); - } - } - - if(isstress) - { - for(int i=0;i<3;i++) - { - for(int j=0;j<3;j++) - { - if(i -#include - -#include "module_base/memory.h" -#include "module_base/parallel_reduce.h" -#include "module_base/timer.h" -#include "module_base/tool_threading.h" -#include "module_cell/module_neighbor/sltk_grid_driver.h" -#include "module_elecstate/cal_dm.h" -#include "module_elecstate/module_dm/cal_dm_psi.h" -#include "module_elecstate/elecstate_lcao.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" -#include "module_io/write_HS.h" - -#ifdef __DEEPKS -#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" -#endif - -#ifdef _OPENMP -#include -#endif - - -// calculate the force due to < phi | Vlocal | dphi > -template<> -void Force_LCAO>::cal_fvl_dphi( - const bool isforce, - const bool isstress, - const elecstate::Potential* pot_in, - TGint>::type& gint, - ModuleBase::matrix& fvl_dphi, - ModuleBase::matrix& svl_dphi) -{ - ModuleBase::TITLE("Force_LCAO", "cal_fvl_dphi"); - if (!isforce && !isstress) - { - return; - } - - ModuleBase::timer::tick("Force_LCAO", "cal_fvl_dphi"); - - const int nspin = PARAM.inp.nspin; - - for (int is = 0; is < nspin; ++is) - { - const double* vr_eff1 = pot_in->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_in->get_effective_vofk(is); - } - - //-------------------------------- - // Grid integration here. - //-------------------------------- - // fvl_dphi can not be set to zero here if Vna is used - if (isstress || isforce) - { - // meta GGA functionals - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) - { - Gint_inout inout(is, - vr_eff1, - vofk_eff1, - isforce, - isstress, - &fvl_dphi, - &svl_dphi, - Gint_Tools::job_type::force_meta); - - gint.cal_gint(&inout); - } - else - { - Gint_inout inout(is, - vr_eff1, - isforce, - isstress, - &fvl_dphi, - &svl_dphi, - Gint_Tools::job_type::force); - - gint.cal_gint(&inout); - } - } - } - - if (isstress) - { - StressTools::stress_fill(-1.0, GlobalC::ucell.omega, svl_dphi); - } - - ModuleBase::timer::tick("Force_LCAO", "cal_fvl_dphi"); - return; -} diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_gamma.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_gamma.cpp index f387f0f923..aef53810ad 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_gamma.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_gamma.cpp @@ -3,6 +3,7 @@ #include "module_parameter/parameter.h" #include "module_cell/module_neighbor/sltk_grid_driver.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_hamilt_lcao/hamilt_lcaodft/stress_tools.h" #include diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_k.cpp index fe36856c65..b57ad57d73 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_k.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_k.cpp @@ -10,6 +10,7 @@ #include "module_elecstate/module_dm/cal_dm_psi.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/write_HS.h" +#include "module_hamilt_lcao/hamilt_lcaodft/stress_tools.h" #include #include diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress.h b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress.h new file mode 100644 index 0000000000..e4a7c4fd0b --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress.h @@ -0,0 +1,69 @@ +#pragma once +#include "module_basis/module_nao/two_center_bundle.h" +#include "module_elecstate/module_dm/density_matrix.h" +#include "module_hamilt_lcao/module_gint/gint_gamma.h" +#include "module_hamilt_lcao/module_gint/gint_k.h" +#include "module_elecstate/potentials/potential_new.h" +#include "module_cell/unitcell.h" +#include "module_hamilt_lcao/hamilt_lcaodft/stress_tools.h" +#ifndef TGINT_H +#define TGINT_H +template +struct TGint; +template <> struct TGint { using type = Gint_Gamma; }; +template <> struct TGint> { using type = Gint_k; }; +#endif + +/// 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-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 + ModuleBase::matrix& s, ///< [out] stress + const elecstate::DensityMatrix& dm, ///< [in] density matrix or energy density matrix + const UnitCell& ucell, ///< [in] unit cell + const elecstate::Potential* pot, ///< [in] potential on grid + typename TGint::type& gint, ///< [in] Gint object + const bool& isforce, + const bool& isstress, + const bool& set_dmr_gint = true); +} +#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_center2.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress_center2.cpp new file mode 100644 index 0000000000..ea49096da5 --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/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 = PARAM.inp.nspin; + const int nlocal = PARAM.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_center2_template.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress_center2_template.hpp new file mode 100644 index 0000000000..e0eedffaad --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress_center2_template.hpp @@ -0,0 +1,128 @@ +#pragma once +#include "pulay_force_stress.h" +#include "module_base/timer.h" +#include "module_parameter/parameter.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 = (PARAM.inp.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_gint.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress_gint.hpp new file mode 100644 index 0000000000..feb567eca2 --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress_gint.hpp @@ -0,0 +1,40 @@ +#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" +#include "module_parameter/parameter.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 = PARAM.inp.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/stress_tools.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/stress_tools.cpp new file mode 100644 index 0000000000..9a62eae1a5 --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/stress_tools.cpp @@ -0,0 +1,20 @@ +#include "stress_tools.h" +namespace StressTools +{ +void stress_fill(const double& lat0_, const double& omega_, ModuleBase::matrix& stress_matrix) +{ + assert(omega_ > 0.0); + double weight = lat0_ / omega_; + for (int i = 0; i < 3; ++i) + { + for (int j = 0; j < 3; ++j) + { + if (j > i) + { + stress_matrix(j, i) = stress_matrix(i, j); + } + stress_matrix(i, j) *= weight; + } + } +} +} // namespace StressTools \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/stress_tools.h b/source/module_hamilt_lcao/hamilt_lcaodft/stress_tools.h new file mode 100644 index 0000000000..4ae6ed31ee --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/stress_tools.h @@ -0,0 +1,8 @@ +#pragma once +#include "module_base/matrix.h" +// this namespace used to store global function for some stress operation +namespace StressTools +{ +// set upper matrix to whole matrix +void stress_fill(const double& lat0_, const double& omega_, ModuleBase::matrix& stress_matrix); +} // namespace StressTools \ No newline at end of file