diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 51f4181913..b4933de79a 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -56,6 +56,7 @@ VPATH=./src_global:\ ./module_hamilt_lcao/module_dftu:\ ./module_hamilt_lcao/module_deltaspin:\ ./module_hamilt_lcao/hamilt_lcaodft/operator_lcao:\ +./module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress:\ ./module_hamilt_lcao/module_gint:\ ./module_relax:\ ./module_relax/relax_old:\ @@ -571,10 +572,6 @@ OBJS_LCAO=evolve_elec.o\ stress_tools.o\ edm.o\ pulay_force_stress_center2.o\ - fedm_gamma.o\ - fedm_k.o\ - ftvnl_dphi_gamma.o\ - ftvnl_dphi_k.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 0bf170672f..abb97265e6 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt +++ b/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt @@ -21,10 +21,6 @@ if(ENABLE_LCAO) FORCE_k.cpp stress_tools.cpp edm.cpp - fedm_gamma.cpp - fedm_k.cpp - ftvnl_dphi_gamma.cpp - ftvnl_dphi_k.cpp fvnl_dbeta_gamma.cpp fvnl_dbeta_k.cpp grid_init.cpp diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp index 2c4b45c312..42ea5fedf8 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp @@ -219,16 +219,13 @@ void Force_LCAO::ftable(const bool isforce, 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_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), ucell, pv, dSx, dSxy, isforce, isstress); 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); + //tvnl_dphi PulayForceStress::cal_pulay_fs(ftvnl_dphi, stvnl_dphi, *dm, ucell, pv, dHx, dHxy, isforce, isstress); this->cal_fvnl_dbeta(dm, @@ -242,6 +239,7 @@ void Force_LCAO::ftable(const bool isforce, fvnl_dbeta, svnl_dbeta); + // vl_dphi PulayForceStress::cal_pulay_fs(fvl_dphi, svl_dphi, *dm, ucell, pelec->pot, gint, isforce, isstress, false/*reset dm to gint*/); #ifdef __DEEPKS diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp index c96cbf708d..7692f5cb30 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp @@ -316,20 +316,17 @@ void Force_LCAO>::ftable(const bool isforce, 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_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), ucell, pv, dSx, fsr.DH_r, isforce, isstress, ra, -1.0, 1.0); 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); + // 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, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/fedm_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/fedm_k.cpp deleted file mode 100644 index 62ae74f870..0000000000 --- a/source/module_hamilt_lcao/hamilt_lcaodft/fedm_k.cpp +++ /dev/null @@ -1,245 +0,0 @@ -#include "FORCE.h" -#include "module_base/memory.h" -#include "module_parameter/parameter.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/elecstate_lcao.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 - -#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; -<<<<<<< HEAD - 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(); -======= ->>>>>>> 27869bac3 (extract edm) - //-------------------------------------------- - // 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 = dm.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 72f741818f..0000000000 --- a/source/module_hamilt_lcao/hamilt_lcaodft/ftvnl_dphi_gamma.cpp +++ /dev/null @@ -1,74 +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" -#include "module_hamilt_lcao/hamilt_lcaodft/stress_tools.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 4f7b8398cd..0000000000 --- a/source/module_hamilt_lcao/hamilt_lcaodft/ftvnl_dphi_k.cpp +++ /dev/null @@ -1,171 +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" -#include "module_hamilt_lcao/hamilt_lcaodft/stress_tools.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; -}