diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 3895a65b9f..925298c983 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -655,7 +655,7 @@ OBJS_LCAO=evolve_elec.o\ FORCE_k.o\ stress_tools.o\ edm.o\ - pulay_force_stress_center2.o\ + pulay_fs_center2.o\ grid_init.o\ spar_dh.o\ spar_exx.o\ diff --git a/source/source_basis/module_nao/two_center_bundle.h b/source/source_basis/module_nao/two_center_bundle.h index 7aa6676483..b2708a9dad 100644 --- a/source/source_basis/module_nao/two_center_bundle.h +++ b/source/source_basis/module_nao/two_center_bundle.h @@ -1,5 +1,5 @@ -#ifndef W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_MODULE_BASIS_MODULE_NAO_TWO_CENTER_BUNDLE_H -#define W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_MODULE_BASIS_MODULE_NAO_TWO_CENTER_BUNDLE_H +#ifndef TWO_CENTER_BUNDLE_H +#define TWO_CENTER_BUNDLE_H #include "source_basis/module_ao/ORB_read.h" #include "source_basis/module_nao/two_center_integrator.h" diff --git a/source/source_esolver/esolver_double_xc.cpp b/source/source_esolver/esolver_double_xc.cpp index c6947b8687..db48bc0a8f 100644 --- a/source/source_esolver/esolver_double_xc.cpp +++ b/source/source_esolver/esolver_double_xc.cpp @@ -6,7 +6,6 @@ #include "source_lcao/module_deepks/LCAO_deepks.h" #include "source_lcao/module_deepks/LCAO_deepks_interface.h" #endif -//-----force& stress------------------- #include "source_lcao/FORCE_STRESS.h" //-----HSolver ElecState Hamilt-------- @@ -16,6 +15,7 @@ #include "source_hsolver/hsolver_lcao.h" #include "source_io/module_parameter/parameter.h" #include "source_io/write_HS.h" // use ModuleIO::write_hsk() +#include "source_lcao/setup_deepks.h" // use deepks, mohan add 2025-10-10 namespace ModuleESolver { @@ -388,6 +388,8 @@ void ESolver_DoubleXC::cal_force(UnitCell& ucell, ModuleBase::matrix& fo // set as base functional Temporarily XC_Functional::set_xc_type(PARAM.inp.deepks_out_base); + this->deepks.dpks_out_type = "base"; // for deepks method + fsl.getForceStress(ucell, PARAM.inp.cal_force, PARAM.inp.cal_stress, @@ -408,12 +410,10 @@ void ESolver_DoubleXC::cal_force(UnitCell& ucell, ModuleBase::matrix& fo this->kv, this->pw_rho, this->solvent, -#ifdef __MLALGO - this->deepks.ld, - "base", -#endif + this->deepks, this->exx_nao, &ucell.symm); + // restore to original xc XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index c9753b3e8f..b7601da811 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -7,10 +7,6 @@ #include "source_estate/module_charge/symmetry_rho.h" #include "source_lcao/LCAO_domain.h" // need DeePKS_init #include "source_lcao/module_dftu/dftu.h" -#ifdef __MLALGO -#include "source_lcao/module_deepks/LCAO_deepks.h" -#include "source_lcao/module_deepks/LCAO_deepks_interface.h" -#endif #include "source_lcao/FORCE_STRESS.h" #include "source_estate/elecstate_lcao.h" #include "source_lcao/hamilt_lcao.h" @@ -203,6 +199,8 @@ void ESolver_KS_LCAO::cal_force(UnitCell& ucell, ModuleBase::matrix& for Force_Stress_LCAO fsl(this->RA, ucell.nat); + deepks.dpks_out_type = "tot"; // for deepks method + fsl.getForceStress(ucell, PARAM.inp.cal_force, PARAM.inp.cal_stress, @@ -223,10 +221,7 @@ void ESolver_KS_LCAO::cal_force(UnitCell& ucell, ModuleBase::matrix& for this->kv, this->pw_rho, this->solvent, -#ifdef __MLALGO - this->deepks.ld, - "tot", -#endif + this->deepks, this->exx_nao, &ucell.symm); @@ -337,43 +332,6 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const } #endif elecstate::setup_dm(ucell, estate, this->psi, this->chr, iter, exx_two_level_step); - -/* - if (iter == 1 && exx_two_level_step == 0) - { - std::cout << " WAVEFUN -> CHARGE " << std::endl; - - // calculate the density matrix using read in wave functions - // and then calculate the charge density on grid. - - estate->skip_weights = true; - elecstate::calculate_weights(estate->ekb, - estate->wg, - estate->klist, - estate->eferm, - estate->f_en, - estate->nelec_spin, - estate->skip_weights); - - elecstate::calEBand(estate->ekb, estate->wg, estate->f_en); - elecstate::cal_dm_psi(estate->DM->get_paraV_pointer(), estate->wg, *this->psi, *(estate->DM)); - estate->DM->cal_DMR(); - - estate->psiToRho(*this->psi); - estate->skip_weights = false; - - elecstate::cal_ux(ucell); - - //! update the potentials by using new electron charge density - estate->pot->update_from_charge(&this->chr, &ucell); - - //! compute the correction energy for metals - estate->f_en.descf = estate->cal_delta_escf(); - } - -*/ - - } #ifdef __EXX @@ -410,11 +368,6 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const { this->p_hamilt->refresh(); } - // if (iter == 1 && istep == 0) - // { - // // initialize DMR - // this->deepks.ld.init_DMR(ucell, orb_, this->pv, this->gd); - // } #endif if (PARAM.inp.vl_in_h) @@ -591,11 +544,8 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& // control the output related to the finished iteration ModuleIO::ctrl_iter_lcao(ucell, PARAM.inp, this->kv, estate, this->pv, this->gd, this->psi, this->chr, this->p_chgmix, - hamilt_lcao, this->orb_, -#ifdef __MLALGO - this->deepks.ld, -#endif - exx_nao, iter, istep, conv_esolver, this->scf_ene_thr); + hamilt_lcao, this->orb_, this->deepks, + this->exx_nao, iter, istep, conv_esolver, this->scf_ene_thr); } diff --git a/source/source_esolver/lcao_before_scf.cpp b/source/source_esolver/lcao_before_scf.cpp index e40203a581..cf09a63943 100644 --- a/source/source_esolver/lcao_before_scf.cpp +++ b/source/source_esolver/lcao_before_scf.cpp @@ -6,9 +6,6 @@ #include "source_cell/module_neighbor/sltk_grid_driver.h" #include "source_io/module_parameter/parameter.h" #include "source_estate/elecstate_tools.h" -#ifdef __MLALGO -#include "source_lcao/module_deepks/LCAO_deepks.h" -#endif #include "source_estate/elecstate_lcao.h" #include "source_lcao/LCAO_domain.h" #include "source_lcao/module_operator_lcao/op_exx_lcao.h" @@ -140,25 +137,10 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) ); } -#ifdef __MLALGO + // 9) for each ionic step, the overlap must be rebuilt // since it depends on ionic positions - if (PARAM.globalv.deepks_setorb) - { - const Parallel_Orbitals* pv = &this->pv; - // allocate , phialpha is different every ion step, so it is allocated here - DeePKS_domain::allocate_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, pv, this->deepks.ld.phialpha); - // build and save at beginning - DeePKS_domain::build_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, - pv, *(two_center_bundle_.overlap_orb_alpha), this->deepks.ld.phialpha); - - if (PARAM.inp.deepks_out_unittest) - { - DeePKS_domain::check_phialpha(PARAM.inp.cal_force, ucell, orb_, - this->gd, pv, this->deepks.ld.phialpha, GlobalV::MY_RANK); - } - } -#endif + this->deepks.build_overlap(ucell, orb_, pv, gd, *(two_center_bundle_.overlap_orb_alpha), PARAM.inp); // 10) prepare sc calculation if (PARAM.inp.sc_mag_switch) diff --git a/source/source_esolver/lcao_others.cpp b/source/source_esolver/lcao_others.cpp index 8237c3a700..98cababdd2 100644 --- a/source/source_esolver/lcao_others.cpp +++ b/source/source_esolver/lcao_others.cpp @@ -15,9 +15,6 @@ #include "source_base/timer.h" #include "source_cell/module_neighbor/sltk_atom_arrange.h" #include "source_cell/module_neighbor/sltk_grid_driver.h" -#ifdef __MLALGO -#include "source_lcao/module_deepks/LCAO_deepks.h" -#endif #include "source_lcao/LCAO_domain.h" #include "source_lcao/module_operator_lcao/op_exx_lcao.h" #include "source_lcao/module_operator_lcao/operator_lcao.h" @@ -242,35 +239,11 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) ); } -#ifdef __MLALGO + // for each ionic step, the overlap must be rebuilt // since it depends on ionic positions - if (PARAM.globalv.deepks_setorb) - { - const Parallel_Orbitals* pv = &this->pv; - // allocate , phialpha is different every ion step, so it is allocated here - DeePKS_domain::allocate_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, pv, this->deepks.ld.phialpha); - // build and save at beginning - DeePKS_domain::build_phialpha(PARAM.inp.cal_force, - ucell, - orb_, - this->gd, - pv, - *(two_center_bundle_.overlap_orb_alpha), - this->deepks.ld.phialpha); + this->deepks.build_overlap(ucell, orb_, pv, gd, *(two_center_bundle_.overlap_orb_alpha), PARAM.inp); - if (PARAM.inp.deepks_out_unittest) - { - DeePKS_domain::check_phialpha(PARAM.inp.cal_force, - ucell, - orb_, - this->gd, - pv, - this->deepks.ld.phialpha, - GlobalV::MY_RANK); - } - } -#endif if (PARAM.inp.sc_mag_switch) { spinconstrain::SpinConstrain& sc = spinconstrain::SpinConstrain::getScInstance(); diff --git a/source/source_io/berryphase.h b/source/source_io/berryphase.h index d4befd298d..1685ccdefd 100644 --- a/source/source_io/berryphase.h +++ b/source/source_io/berryphase.h @@ -4,7 +4,6 @@ #ifdef __LCAO #include "unk_overlap_lcao.h" #endif -//#include "source_basis/module_pw/pw_basis.h" #include "source_basis/module_pw/pw_basis_k.h" #include "source_cell/klist.h" #include "source_psi/psi.h" diff --git a/source/source_io/ctrl_iter_lcao.cpp b/source/source_io/ctrl_iter_lcao.cpp index cb0dd3cca0..d0b0c1bd0a 100644 --- a/source/source_io/ctrl_iter_lcao.cpp +++ b/source/source_io/ctrl_iter_lcao.cpp @@ -21,9 +21,7 @@ void ctrl_iter_lcao(UnitCell& ucell, // unit cell * Charge_Mixing* p_chgmix, // charge mixing * hamilt::HamiltLCAO* p_hamilt, // hamiltonian * LCAO_Orbitals &orb, // orbital info * -#ifdef __MLALGO - LCAO_Deepks& ld, -#endif + Setup_DeePKS &deepks, Exx_NAO &exx_nao, int &iter, const int istep, @@ -65,7 +63,7 @@ void ctrl_iter_lcao(UnitCell& ucell, // unit cell * { if (iter % inp.deepks_out_freq_elec == 0 ) { - std::shared_ptr> ld_shared_ptr(&ld, [](LCAO_Deepks*) {}); + std::shared_ptr> ld_shared_ptr(&deepks.ld, [](LCAO_Deepks*) {}); LCAO_Deepks_Interface deepks_interface(ld_shared_ptr); deepks_interface.out_deepks_labels(pelec->f_en.etot, kv.get_nks(), @@ -91,9 +89,7 @@ template void ctrl_iter_lcao(UnitCell& ucell, // unit cell * Charge_Mixing* p_chgmix, // charge mixing * hamilt::HamiltLCAO* p_hamilt, // hamiltonian * LCAO_Orbitals &orb, // orbital info * -#ifdef __MLALGO - LCAO_Deepks& ld, -#endif + Setup_DeePKS &deepks, Exx_NAO &exx_nao, int &iter, const int istep, @@ -112,9 +108,7 @@ template void ctrl_iter_lcao, double>(UnitCell& ucell, // u Charge_Mixing* p_chgmix, // charge mixing * hamilt::HamiltLCAO, double>* p_hamilt, // hamiltonian * LCAO_Orbitals &orb, // orbital info * -#ifdef __MLALGO - LCAO_Deepks>& ld, -#endif + Setup_DeePKS> &deepks, Exx_NAO> &exx_nao, int &iter, const int istep, @@ -133,9 +127,7 @@ template void ctrl_iter_lcao, std::complex>(UnitCel Charge_Mixing* p_chgmix, // charge mixing * hamilt::HamiltLCAO, std::complex>* p_hamilt, // hamiltonian * LCAO_Orbitals &orb, // orbital info * -#ifdef __MLALGO - LCAO_Deepks>& ld, -#endif + Setup_DeePKS> &deepks, Exx_NAO> &exx_nao, int &iter, const int istep, diff --git a/source/source_io/ctrl_iter_lcao.h b/source/source_io/ctrl_iter_lcao.h index 94c3c4364e..c94afac122 100644 --- a/source/source_io/ctrl_iter_lcao.h +++ b/source/source_io/ctrl_iter_lcao.h @@ -9,6 +9,7 @@ #include "source_estate/module_charge/charge_mixing.h" // use charge mixing #include "source_lcao/hamilt_lcao.h" // use hamilt::HamiltLCAO #include "source_lcao/setup_exx.h" // mohan add 20251008 +#include "source_lcao/setup_deepks.h" // mohan add 20251010 namespace ModuleIO { @@ -25,9 +26,7 @@ void ctrl_iter_lcao(UnitCell& ucell, // unit cell * Charge_Mixing* p_chgmix, // charge mixing * hamilt::HamiltLCAO* p_hamilt, // hamiltonian * LCAO_Orbitals &orb, // orbital info * -#ifdef __MLALGO - LCAO_Deepks& ld, -#endif + Setup_DeePKS &deepks, Exx_NAO &exx_nao, int &iter, const int istep, diff --git a/source/source_lcao/CMakeLists.txt b/source/source_lcao/CMakeLists.txt index 30ea5a3118..1831ac6522 100644 --- a/source/source_lcao/CMakeLists.txt +++ b/source/source_lcao/CMakeLists.txt @@ -22,7 +22,7 @@ if(ENABLE_LCAO) module_operator_lcao/td_pot_hybrid.cpp module_operator_lcao/dspin_lcao.cpp module_operator_lcao/dftu_lcao.cpp - pulay_force_stress_center2.cpp + pulay_fs_center2.cpp FORCE_STRESS.cpp FORCE_gamma.cpp FORCE_k.cpp diff --git a/source/source_lcao/FORCE.h b/source/source_lcao/FORCE.h index c5649246d8..e659863b56 100644 --- a/source/source_lcao/FORCE.h +++ b/source/source_lcao/FORCE.h @@ -13,6 +13,7 @@ #include "source_lcao/module_gint/gint_gamma.h" #include "source_lcao/module_gint/gint_k.h" #include "source_psi/psi.h" +#include "source_lcao/setup_deepks.h" #ifndef TGINT_H #define TGINT_H @@ -63,11 +64,9 @@ class Force_LCAO ModuleBase::matrix& stvnl_dphi, ModuleBase::matrix& svnl_dbeta, ModuleBase::matrix& svl_dphi, -#ifdef __MLALGO ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, - LCAO_Deepks& ld, -#endif + Setup_DeePKS& deepks, typename TGint::type& gint, const TwoCenterBundle& two_center_bundle, const LCAO_Orbitals& orb, diff --git a/source/source_lcao/FORCE_STRESS.cpp b/source/source_lcao/FORCE_STRESS.cpp index 3188d26310..ef455260f8 100644 --- a/source/source_lcao/FORCE_STRESS.cpp +++ b/source/source_lcao/FORCE_STRESS.cpp @@ -51,10 +51,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, const K_Vectors& kv, ModulePW::PW_Basis* rhopw, surchem& solvent, -#ifdef __MLALGO - LCAO_Deepks& ld, - const std::string& dpks_out_type, -#endif + Setup_DeePKS& deepks, Exx_NAO &exx_nao, ModuleSymmetry::Symmetry* symm) { @@ -82,39 +79,25 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, ModuleBase::matrix fewalds; ModuleBase::matrix fcc; ModuleBase::matrix fscc; -#ifdef __MLALGO ModuleBase::matrix fvnl_dalpha; // deepks -#endif fvl_dphi.create(nat, 3); // must do it now, update it later, noted by zhengdy if (isforce) { fcs.create(nat, 3); - foverlap.create(nat, 3); - ftvnl_dphi.create(nat, 3); - fvnl_dbeta.create(nat, 3); - fvl_dvl.create(nat, 3); - fewalds.create(nat, 3); - fcc.create(nat, 3); - fscc.create(nat, 3); -#ifdef __MLALGO + foverlap.create(nat, 3); // overlap force + ftvnl_dphi.create(nat, 3); // pulay force of NAO + fvnl_dbeta.create(nat, 3); // pulay force of non-local projectors + fvl_dvl.create(nat, 3); // force from local potentials + fewalds.create(nat, 3); // Ewald force + fcc.create(nat, 3); // force due to core correction + fscc.create(nat, 3); // force due to self-consistent field fvnl_dalpha.create(nat, 3); // deepks -#endif // calculate basic terms in Force, same method with PW base - this->calForcePwPart(ucell, - fvl_dvl, - fewalds, - fcc, - fscc, - pelec->f_en.etxc, - pelec->vnew, - pelec->vnew_exist, - pelec->charge, - rhopw, - locpp, - sf); + this->calForcePwPart(ucell, fvl_dvl, fewalds, fcc, fscc, pelec->f_en.etxc, + pelec->vnew, pelec->vnew_exist, pelec->charge, rhopw, locpp, sf); } // total stress : ModuleBase::matrix scs @@ -127,9 +110,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, ModuleBase::matrix stvnl_dphi; ModuleBase::matrix svnl_dbeta; ModuleBase::matrix svl_dphi; -#ifdef __MLALGO ModuleBase::matrix svnl_dalpha; // deepks -#endif //! stress if (isstress) @@ -145,67 +126,33 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, stvnl_dphi.create(3, 3); svnl_dbeta.create(3, 3); svl_dphi.create(3, 3); -#ifdef __MLALGO svnl_dalpha.create(3, 3); -#endif + // calculate basic terms in Stress, similar method with PW base - this->calStressPwPart(ucell, - sigmadvl, - sigmahar, - sigmaewa, - sigmacc, - sigmaxc, - pelec->f_en.etxc, - pelec->charge, - rhopw, - locpp, - sf); + this->calStressPwPart(ucell, sigmadvl, sigmahar, sigmaewa, sigmacc, + sigmaxc, pelec->f_en.etxc, pelec->charge, rhopw, locpp, sf); } //! atomic forces from integration (4 terms) - this->integral_part(PARAM.globalv.gamma_only_local, - isforce, - isstress, - ucell, - gd, - fsr, - pelec, - psi, - foverlap, - ftvnl_dphi, - fvnl_dbeta, - fvl_dphi, - soverlap, - stvnl_dphi, - svnl_dbeta, - svl_dphi, -#ifdef __MLALGO - fvnl_dalpha, - svnl_dalpha, - ld, -#endif - gint_gamma, - gint_k, - two_center_bundle, - orb, - pv, - kv); + this->integral_part(PARAM.globalv.gamma_only_local, isforce, isstress, + ucell, gd, fsr, pelec, psi, foverlap, ftvnl_dphi, + fvnl_dbeta, fvl_dphi, soverlap, stvnl_dphi, svnl_dbeta, + svl_dphi, fvnl_dalpha, svnl_dalpha, deepks, gint_gamma, + gint_k, two_center_bundle, orb, pv, kv); + // calculate force and stress for Nonlocal part if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) { hamilt::NonlocalNew> tmp_nonlocal(nullptr, - kv.kvec_d, - nullptr, - &ucell, - orb.cutoffs(), - &gd, - two_center_bundle.overlap_orb_beta.get()); + kv.kvec_d, nullptr, &ucell, orb.cutoffs(), &gd, two_center_bundle.overlap_orb_beta.get()); const auto* dm_p = dynamic_cast*>(pelec)->get_DM(); + if (PARAM.inp.nspin == 2) { const_cast*>(dm_p)->switch_dmr(1); } + const hamilt::HContainer* dmr = dm_p->get_DMR_pointer(1); tmp_nonlocal.cal_force_stress(isforce, isstress, dmr, fvnl_dbeta, svnl_dbeta); if (PARAM.inp.nspin == 2) @@ -216,12 +163,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, else if (PARAM.inp.nspin == 4) { hamilt::NonlocalNew, std::complex>> tmp_nonlocal( - nullptr, - kv.kvec_d, - nullptr, - &ucell, - orb.cutoffs(), - &gd, + nullptr, kv.kvec_d, nullptr, &ucell, orb.cutoffs(), &gd, two_center_bundle.overlap_orb_beta.get()); // calculate temporary complex DMR for nonlocal force&stress @@ -491,43 +433,9 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, this->forceSymmetry(ucell, fcs, symm); } -#ifdef __MLALGO - // DeePKS force - if (PARAM.inp.deepks_out_labels) // not parallelized yet - { - if (PARAM.inp.deepks_out_base == "none" - || (PARAM.inp.deepks_out_base != "none" && dpks_out_type == "tot") ) - { - const std::string file_ftot = PARAM.globalv.global_out_dir - + (PARAM.inp.deepks_out_labels == 1 ? "deepks_ftot.npy" : "deepks_force.npy"); - LCAO_deepks_io::save_matrix2npy(file_ftot, fcs, GlobalV::MY_RANK); // Hartree/Bohr, F_tot + // compute forces using the DeePKS model + deepks.write_forces(fcs, fvnl_dalpha, PARAM.inp); - if (PARAM.inp.deepks_out_labels == 1) - { - // this base only considers subtracting the deepks_scf part - const std::string file_fbase = PARAM.globalv.global_out_dir + "deepks_fbase.npy"; - if (PARAM.inp.deepks_scf) - { - LCAO_deepks_io::save_matrix2npy(file_fbase, - fcs - fvnl_dalpha, - GlobalV::MY_RANK); // Hartree/Bohr, F_base - } - else - { - LCAO_deepks_io::save_matrix2npy(file_fbase, fcs, GlobalV::MY_RANK); // no scf, F_base=F_tot - } - } - } - if (PARAM.inp.deepks_out_base != "none") - { - // output fcs as tot or base in another dir - // this base considers changing xc functional to base functional - const std::string file_f = PARAM.globalv.global_deepks_label_elec_dir + - (dpks_out_type == "tot" ? "ftot.npy" : "fbase.npy"); - LCAO_deepks_io::save_matrix2npy(file_f, fcs, GlobalV::MY_RANK); - } - } -#endif // print Rydberg force or not bool ry = false; if (istestf) @@ -693,58 +601,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, symm->symmetrize_mat3(scs, ucell.lat); } // end symmetry -#ifdef __MLALGO - if (PARAM.inp.deepks_out_labels == 1) - { - if (PARAM.inp.deepks_out_base == "none" - || (PARAM.inp.deepks_out_base != "none" && dpks_out_type == "tot") ) - { - const std::string file_stot = PARAM.globalv.global_out_dir + "deepks_stot.npy"; - LCAO_deepks_io::save_matrix2npy(file_stot, - scs, - GlobalV::MY_RANK, - ucell.omega, - 'U'); // change to energy unit Ry when printing, S_tot; - - // this base only considers subtracting the deepks_scf part - const std::string file_sbase = PARAM.globalv.global_out_dir + "deepks_sbase.npy"; - if (PARAM.inp.deepks_scf) - { - LCAO_deepks_io::save_matrix2npy(file_sbase, - scs - svnl_dalpha, - GlobalV::MY_RANK, - ucell.omega, - 'U'); // change to energy unit Ry when printing, S_base; - } - else - { - LCAO_deepks_io::save_matrix2npy(file_sbase, - scs, - GlobalV::MY_RANK, - ucell.omega, - 'U'); // sbase = stot - } - } - if (PARAM.inp.deepks_out_base != "none") - { - // output scs as tot or base in another dir - // this base considers changing xc functional to base functional - const std::string file_s = PARAM.globalv.global_deepks_label_elec_dir + - (dpks_out_type == "tot" ? "stot.npy" : "sbase.npy"); - LCAO_deepks_io::save_matrix2npy(file_s, - scs, - GlobalV::MY_RANK, - ucell.omega, - 'U'); // change to energy unit Ry when printing, S_tot; - } - } - else if (PARAM.inp.deepks_out_labels == 2) - { - const std::string file_stot = PARAM.globalv.global_out_dir + "deepks_stress.npy"; - LCAO_deepks_io::save_matrix2npy(file_stot, scs, GlobalV::MY_RANK, ucell.omega, - 'F'); // flat mode - } -#endif + deepks.write_stress(scs, svnl_dalpha, ucell.omega, PARAM.inp); // print Rydberg stress or not bool ry = false; @@ -872,11 +729,9 @@ void Force_Stress_LCAO::integral_part(const bool isGammaOnly, ModuleBase::matrix& stvnl_dphi, ModuleBase::matrix& svnl_dbeta, ModuleBase::matrix& svl_dphi, -#if __MLALGO ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, - LCAO_Deepks& ld, -#endif + Setup_DeePKS& deepks, Gint_Gamma& gint_gamma, // mohan add 2024-04-01 Gint_k& gint_k, // mohan add 2024-04-01 const TwoCenterBundle& two_center_bundle, @@ -885,30 +740,11 @@ void Force_Stress_LCAO::integral_part(const bool isGammaOnly, const K_Vectors& kv) { - flk.ftable(isforce, - isstress, - fsr, // mohan add 2024-06-15 - ucell, - gd, - psi, - pelec, - foverlap, - ftvnl_dphi, - fvnl_dbeta, - fvl_dphi, - soverlap, - stvnl_dphi, - svnl_dbeta, - svl_dphi, -#if __MLALGO - fvnl_dalpha, - svnl_dalpha, - ld, -#endif - gint_gamma, - two_center_bundle, - orb, - pv); + flk.ftable(isforce, isstress, fsr, ucell, gd, psi, pelec, + foverlap, ftvnl_dphi, fvnl_dbeta, fvl_dphi, + soverlap, stvnl_dphi, svnl_dbeta, svl_dphi, + fvnl_dalpha, svnl_dalpha, deepks, gint_gamma, + two_center_bundle, orb, pv); return; } @@ -929,11 +765,9 @@ void Force_Stress_LCAO>::integral_part(const bool isGammaOn ModuleBase::matrix& stvnl_dphi, ModuleBase::matrix& svnl_dbeta, ModuleBase::matrix& svl_dphi, -#if __MLALGO ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, - LCAO_Deepks>& ld, -#endif + Setup_DeePKS>& deepks, Gint_Gamma& gint_gamma, Gint_k& gint_k, const TwoCenterBundle& two_center_bundle, @@ -941,32 +775,11 @@ void Force_Stress_LCAO>::integral_part(const bool isGammaOn const Parallel_Orbitals& pv, const K_Vectors& kv) { - flk.ftable(isforce, - isstress, - fsr, // mohan add 2024-06-16 - ucell, - gd, - psi, - pelec, - foverlap, - ftvnl_dphi, - fvnl_dbeta, - fvl_dphi, - soverlap, - stvnl_dphi, - svnl_dbeta, - svl_dphi, -#if __MLALGO - fvnl_dalpha, - svnl_dalpha, - ld, -#endif - gint_k, - two_center_bundle, - orb, - pv, - &kv, - this->RA); + flk.ftable(isforce, isstress, fsr, ucell, gd, psi, pelec, + foverlap, ftvnl_dphi, fvnl_dbeta, fvl_dphi, + soverlap, stvnl_dphi, svnl_dbeta, svl_dphi, + fvnl_dalpha, svnl_dalpha, deepks, gint_k, + two_center_bundle, orb, pv, &kv, this->RA); return; } @@ -985,30 +798,20 @@ void Force_Stress_LCAO::calStressPwPart(UnitCell& ucell, const Structure_Factor& sf) { ModuleBase::TITLE("Force_Stress_LCAO", "calStressPwPart"); - //-------------------------------------------------------- + // local pseudopotential stress: - // use charge density; plane wave; local pseudopotential; - //-------------------------------------------------------- sc_pw.stress_loc(ucell, sigmadvl, rhopw, locpp.vloc, &sf, 0, chr); - //-------------------------------------------------------- // hartree term - //-------------------------------------------------------- sc_pw.stress_har(ucell, sigmahar, rhopw, 0, chr); - //-------------------------------------------------------- // ewald stress: use plane wave only. - //-------------------------------------------------------- sc_pw.stress_ewa(ucell, sigmaewa, rhopw, 0); // remain problem - //-------------------------------------------------------- // stress due to core correlation. - //-------------------------------------------------------- sc_pw.stress_cc(sigmacc, rhopw, ucell, &sf, 0, locpp.numeric, chr); - //-------------------------------------------------------- // stress due to self-consistent charge. - //-------------------------------------------------------- for (int i = 0; i < 3; i++) { sigmaxc(i, i) = -etxc / ucell.omega; @@ -1027,21 +830,9 @@ void Force_Stress_LCAO::forceSymmetry(const UnitCell& ucell, ModuleBase::matr double d1, d2, d3; for (int iat = 0; iat < ucell.nat; iat++) { - ModuleBase::Mathzone::Cartesian_to_Direct(fcs(iat, 0), - fcs(iat, 1), - fcs(iat, 2), - ucell.a1.x, - ucell.a1.y, - ucell.a1.z, - ucell.a2.x, - ucell.a2.y, - ucell.a2.z, - ucell.a3.x, - ucell.a3.y, - ucell.a3.z, - d1, - d2, - d3); + ModuleBase::Mathzone::Cartesian_to_Direct(fcs(iat, 0), fcs(iat, 1), fcs(iat, 2), + ucell.a1.x, ucell.a1.y, ucell.a1.z, ucell.a2.x, ucell.a2.y, ucell.a2.z, + ucell.a3.x, ucell.a3.y, ucell.a3.z, d1, d2, d3); fcs(iat, 0) = d1; fcs(iat, 1) = d2; @@ -1050,21 +841,9 @@ void Force_Stress_LCAO::forceSymmetry(const UnitCell& ucell, ModuleBase::matr symm->symmetrize_vec3_nat(fcs.c); for (int iat = 0; iat < ucell.nat; iat++) { - ModuleBase::Mathzone::Direct_to_Cartesian(fcs(iat, 0), - fcs(iat, 1), - fcs(iat, 2), - ucell.a1.x, - ucell.a1.y, - ucell.a1.z, - ucell.a2.x, - ucell.a2.y, - ucell.a2.z, - ucell.a3.x, - ucell.a3.y, - ucell.a3.z, - d1, - d2, - d3); + ModuleBase::Mathzone::Direct_to_Cartesian(fcs(iat, 0), fcs(iat, 1), fcs(iat, 2), + ucell.a1.x, ucell.a1.y, ucell.a1.z, ucell.a2.x, ucell.a2.y, ucell.a2.z, + ucell.a3.x, ucell.a3.y, ucell.a3.z, d1, d2, d3); fcs(iat, 0) = d1; fcs(iat, 1) = d2; diff --git a/source/source_lcao/FORCE_STRESS.h b/source/source_lcao/FORCE_STRESS.h index 5febf0a1aa..5ef6648621 100644 --- a/source/source_lcao/FORCE_STRESS.h +++ b/source/source_lcao/FORCE_STRESS.h @@ -17,6 +17,7 @@ #include "source_lcao/module_gint/gint_gamma.h" #include "source_lcao/module_gint/gint_k.h" #include "source_lcao/setup_exx.h" // for exx, mohan add 20251008 +#include "source_lcao/setup_deepks.h" // for deepks, mohan add 20251010 template @@ -51,10 +52,7 @@ class Force_Stress_LCAO const K_Vectors& kv, ModulePW::PW_Basis* rhopw, surchem& solvent, -#ifdef __MLALGO - LCAO_Deepks& ld, - const std::string& dpks_out_type, -#endif + Setup_DeePKS &deepks, Exx_NAO &exx_nao, ModuleSymmetry::Symmetry* symm); @@ -95,11 +93,9 @@ class Force_Stress_LCAO ModuleBase::matrix& stvnl_dphi, ModuleBase::matrix& svnl_dbeta, ModuleBase::matrix& svl_dphi, -#if __MLALGO ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, - LCAO_Deepks& ld, -#endif + Setup_DeePKS& deepks, Gint_Gamma& gint_gamma, Gint_k& gint_k, const TwoCenterBundle& two_center_bundle, diff --git a/source/source_lcao/FORCE_gamma.cpp b/source/source_lcao/FORCE_gamma.cpp index 96f2537900..ffd17393bb 100644 --- a/source/source_lcao/FORCE_gamma.cpp +++ b/source/source_lcao/FORCE_gamma.cpp @@ -12,7 +12,7 @@ #include "source_cell/module_neighbor/sltk_grid_driver.h" //GridD #include "source_estate/elecstate_lcao.h" #include "source_lcao/LCAO_domain.h" -#include "source_lcao/pulay_force_stress.h" +#include "source_lcao/pulay_fs.h" #include "source_io/write_HS.h" template <> @@ -167,11 +167,9 @@ void Force_LCAO::ftable(const bool isforce, ModuleBase::matrix& stvnl_dphi, ModuleBase::matrix& svnl_dbeta, ModuleBase::matrix& svl_dphi, -#ifdef __MLALGO ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, - LCAO_Deepks& ld, -#endif + Setup_DeePKS& deepks, TGint::type& gint, const TwoCenterBundle& two_center_bundle, const LCAO_Orbitals& orb, @@ -232,16 +230,16 @@ void Force_LCAO::ftable(const bool isforce, { // No need to update E_delta here since it have been done in LCAO_Deepks_Interface in after_scf const int nks = 1; - DeePKS_domain::cal_f_delta(ld.dm_r, + DeePKS_domain::cal_f_delta(deepks.ld.dm_r, ucell, orb, gd, *this->ParaV, nks, kv->kvec_d, - ld.phialpha, - ld.gedm, - ld.inl_index, + deepks.ld.phialpha, + deepks.ld.gedm, + deepks.ld.inl_index, fvnl_dalpha, isstress, svnl_dalpha); diff --git a/source/source_lcao/FORCE_k.cpp b/source/source_lcao/FORCE_k.cpp index 4e65f639e2..ba5b58c89c 100644 --- a/source/source_lcao/FORCE_k.cpp +++ b/source/source_lcao/FORCE_k.cpp @@ -9,7 +9,7 @@ #include "source_estate/elecstate_lcao.h" #include "source_estate/module_dm/cal_dm_psi.h" #include "source_lcao/LCAO_domain.h" -#include "source_lcao/pulay_force_stress.h" +#include "source_lcao/pulay_fs.h" #include "source_pw/module_pwdft/global.h" #include "source_io/write_HS.h" #include "source_io/module_parameter/parameter.h" @@ -204,11 +204,9 @@ void Force_LCAO>::ftable(const bool isforce, ModuleBase::matrix& stvnl_dphi, ModuleBase::matrix& svnl_dbeta, ModuleBase::matrix& svl_dphi, -#ifdef __MLALGO ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, - LCAO_Deepks>& ld, -#endif + Setup_DeePKS>& deepks, TGint>::type& gint, const TwoCenterBundle& two_center_bundle, const LCAO_Orbitals& orb, @@ -232,53 +230,40 @@ void Force_LCAO>::ftable(const bool isforce, 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. PulayForceStress::cal_pulay_fs( - foverlap, - soverlap, + 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); + 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 + // 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. // vl_dphi - PulayForceStress::cal_pulay_fs(fvl_dphi, - svl_dphi, - *dm, - ucell, - pelec->pot, - gint, - isforce, - isstress, + PulayForceStress::cal_pulay_fs(fvl_dphi, svl_dphi, *dm, ucell, + pelec->pot, gint, isforce, isstress, false /*reset dm to gint*/); #ifdef __MLALGO if (PARAM.inp.deepks_scf) { // No need to update E_delta since it have been done in LCAO_Deepks_Interface in after_scf - DeePKS_domain::cal_f_delta>(ld.dm_r, + DeePKS_domain::cal_f_delta>(deepks.ld.dm_r, ucell, orb, gd, pv, kv->get_nks(), kv->kvec_d, - ld.phialpha, - ld.gedm, - ld.inl_index, + deepks.ld.phialpha, + deepks.ld.gedm, + deepks.ld.inl_index, fvnl_dalpha, isstress, svnl_dalpha); diff --git a/source/source_lcao/pulay_force_stress_center2.cpp b/source/source_lcao/pulay_force_stress_center2.cpp deleted file mode 100644 index 9929285665..0000000000 --- a/source/source_lcao/pulay_force_stress_center2.cpp +++ /dev/null @@ -1,132 +0,0 @@ -#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("Forces", "cal_pulay_fs"); - ModuleBase::timer::tick("Forces", "cal_pulay_fs"); - - 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("Forces", "cal_pulay_fs"); - } - - 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); - } - -} diff --git a/source/source_lcao/pulay_force_stress.h b/source/source_lcao/pulay_fs.h similarity index 97% rename from source/source_lcao/pulay_force_stress.h rename to source/source_lcao/pulay_fs.h index f89444ec1e..78094449fc 100644 --- a/source/source_lcao/pulay_force_stress.h +++ b/source/source_lcao/pulay_fs.h @@ -65,5 +65,5 @@ namespace PulayForceStress 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 +#include "pulay_fs_temp.hpp" +#include "pulay_fs_gint.hpp" diff --git a/source/source_lcao/pulay_fs_center2.cpp b/source/source_lcao/pulay_fs_center2.cpp new file mode 100644 index 0000000000..25bc2cdc66 --- /dev/null +++ b/source/source_lcao/pulay_fs_center2.cpp @@ -0,0 +1,144 @@ +#include "pulay_fs.h" + +template<> // gamma-only, provided xy +void PulayForceStress::cal_pulay_fs( + ModuleBase::matrix& force, + ModuleBase::matrix& stress, + 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("Forces", "cal_pulay_fs"); + ModuleBase::timer::tick("Forces", "cal_pulay_fs"); + + 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) + { + force(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) + { + stress(i, j) += sums * dHSxy[ij++][index]; + } + } + } + } + } + } + + if (isstress) + { + StressTools::stress_fill(ucell.lat0, ucell.omega, stress); + } + + ModuleBase::timer::tick("Forces", "cal_pulay_fs"); +} + + +template<> //multi-k, provided xy +void PulayForceStress::cal_pulay_fs( + ModuleBase::matrix& force, + ModuleBase::matrix& stress, + 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(force, stress, dm, ucell, pv, dHSx, dHSxy, + nullptr, isforce, isstress, ra, + factor_force, factor_stress, stress_func); +} + + +template<> // multi-k, provided x +void PulayForceStress::cal_pulay_fs( + ModuleBase::matrix& force, + ModuleBase::matrix& stress, + 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(force, stress, dm, ucell, pv, dHSx, + nullptr, dtau, isforce, isstress, ra, + factor_force, factor_stress, stress_func); +} + diff --git a/source/source_lcao/pulay_force_stress_gint.hpp b/source/source_lcao/pulay_fs_gint.hpp similarity index 98% rename from source/source_lcao/pulay_force_stress_gint.hpp rename to source/source_lcao/pulay_fs_gint.hpp index f310b31fad..f8bdb37d48 100644 --- a/source/source_lcao/pulay_force_stress_gint.hpp +++ b/source/source_lcao/pulay_fs_gint.hpp @@ -1,5 +1,5 @@ #pragma once -#include "pulay_force_stress.h" +#include "pulay_fs.h" #include "source_lcao/stress_tools.h" #include "source_hamilt/module_xc/xc_functional.h" #include "source_io/module_parameter/parameter.h" @@ -62,4 +62,4 @@ namespace PulayForceStress if (isstress) { StressTools::stress_fill(-1.0, ucell.omega, s); } } -} \ No newline at end of file +} diff --git a/source/source_lcao/pulay_force_stress_center2_template.hpp b/source/source_lcao/pulay_fs_temp.hpp similarity index 99% rename from source/source_lcao/pulay_force_stress_center2_template.hpp rename to source/source_lcao/pulay_fs_temp.hpp index a035b5aae1..fb5b728a87 100644 --- a/source/source_lcao/pulay_force_stress_center2_template.hpp +++ b/source/source_lcao/pulay_fs_temp.hpp @@ -1,5 +1,5 @@ #pragma once -#include "pulay_force_stress.h" +#include "pulay_fs.h" #include "source_base/timer.h" #include "source_io/module_parameter/parameter.h" namespace PulayForceStress diff --git a/source/source_lcao/setup_deepks.cpp b/source/source_lcao/setup_deepks.cpp index 8608e65095..460ce7bf76 100644 --- a/source/source_lcao/setup_deepks.cpp +++ b/source/source_lcao/setup_deepks.cpp @@ -1,5 +1,6 @@ #include "source_lcao/setup_deepks.h" #include "source_lcao/LCAO_domain.h" +#include "source_io/module_parameter/parameter.h" // use parameter template Setup_DeePKS::Setup_DeePKS(){} @@ -7,6 +8,38 @@ Setup_DeePKS::Setup_DeePKS(){} template Setup_DeePKS::~Setup_DeePKS(){} + +template +void Setup_DeePKS::build_overlap( + const UnitCell &ucell, + const LCAO_Orbitals &orb, + const Parallel_Orbitals &pv, + const Grid_Driver &gd, + TwoCenterIntegrator &overlap_orb_alpha, + const Input_para &inp) +{ +#ifdef __MLALGO + // 9) for each ionic step, the overlap must be rebuilt + // since it depends on ionic positions + if (PARAM.globalv.deepks_setorb) + { + // allocate , phialpha is different every ion step, so it is allocated here + DeePKS_domain::allocate_phialpha(inp.cal_force, ucell, orb, gd, &pv, this->ld.phialpha); + + // build and save at beginning + DeePKS_domain::build_phialpha(inp.cal_force, ucell, orb, gd, + &pv, overlap_orb_alpha, this->ld.phialpha); + + if (inp.deepks_out_unittest) + { + DeePKS_domain::check_phialpha(inp.cal_force, ucell, orb, + gd, &pv, this->ld.phialpha, GlobalV::MY_RANK); + } + } +#endif +} + + template void Setup_DeePKS::before_runner(const UnitCell& ucell, // unitcell const int nks, // number of k points @@ -28,5 +61,116 @@ void Setup_DeePKS::before_runner(const UnitCell& ucell, // unitcell #endif } + +template +void Setup_DeePKS::write_forces( + const ModuleBase::matrix &fcs, + const ModuleBase::matrix &fvnl_dalpha, + const Input_para &inp) +{ +#ifdef __MLALGO + // DeePKS force + if (inp.deepks_out_labels) // not parallelized yet + { + if (inp.deepks_out_base == "none" + || (inp.deepks_out_base != "none" && this->dpks_out_type == "tot") ) + { + const std::string file_ftot = PARAM.globalv.global_out_dir + + (inp.deepks_out_labels == 1 ? "deepks_ftot.npy" : "deepks_force.npy"); + LCAO_deepks_io::save_matrix2npy(file_ftot, fcs, GlobalV::MY_RANK); // Hartree/Bohr, F_tot + + if (inp.deepks_out_labels == 1) + { + // this base only considers subtracting the deepks_scf part + const std::string file_fbase = PARAM.globalv.global_out_dir + "deepks_fbase.npy"; + if (inp.deepks_scf) + { + LCAO_deepks_io::save_matrix2npy(file_fbase, + fcs - fvnl_dalpha, + GlobalV::MY_RANK); // Hartree/Bohr, F_base + } + else + { + LCAO_deepks_io::save_matrix2npy(file_fbase, fcs, GlobalV::MY_RANK); // no scf, F_base=F_tot + } + } + } + if (inp.deepks_out_base != "none") + { + // output fcs as tot or base in another dir + // this base considers changing xc functional to base functional + const std::string file_f = PARAM.globalv.global_deepks_label_elec_dir + + (dpks_out_type == "tot" ? "ftot.npy" : "fbase.npy"); + LCAO_deepks_io::save_matrix2npy(file_f, fcs, GlobalV::MY_RANK); + } + } +#endif + +} + + +template +void Setup_DeePKS::write_stress( + const ModuleBase::matrix &scs, + const ModuleBase::matrix &svnl_dalpha, + const double &omega, + const Input_para &inp) +{ +#ifdef __MLALGO + if (inp.deepks_out_labels == 1) + { + assert(omega>0.0); + + if (inp.deepks_out_base == "none" + || (inp.deepks_out_base != "none" && this->dpks_out_type == "tot") ) + { + const std::string file_stot = PARAM.globalv.global_out_dir + "deepks_stot.npy"; + LCAO_deepks_io::save_matrix2npy(file_stot, + scs, + GlobalV::MY_RANK, + omega, + 'U'); // change to energy unit Ry when printing, S_tot; + + // this base only considers subtracting the deepks_scf part + const std::string file_sbase = PARAM.globalv.global_out_dir + "deepks_sbase.npy"; + if (inp.deepks_scf) + { + LCAO_deepks_io::save_matrix2npy(file_sbase, + scs - svnl_dalpha, + GlobalV::MY_RANK, + omega, + 'U'); // change to energy unit Ry when printing, S_base; + } + else + { + LCAO_deepks_io::save_matrix2npy(file_sbase, + scs, + GlobalV::MY_RANK, + omega, + 'U'); // sbase = stot + } + } + if (inp.deepks_out_base != "none") + { + // output scs as tot or base in another dir + // this base considers changing xc functional to base functional + const std::string file_s = PARAM.globalv.global_deepks_label_elec_dir + + (this->dpks_out_type == "tot" ? "stot.npy" : "sbase.npy"); + LCAO_deepks_io::save_matrix2npy(file_s, + scs, + GlobalV::MY_RANK, + omega, + 'U'); // change to energy unit Ry when printing, S_tot; + } + } + else if (inp.deepks_out_labels == 2) + { + const std::string file_stot = PARAM.globalv.global_out_dir + "deepks_stress.npy"; + LCAO_deepks_io::save_matrix2npy(file_stot, scs, GlobalV::MY_RANK, omega, + 'F'); // flat mode + } +#endif +} + template class Setup_DeePKS; template class Setup_DeePKS>; diff --git a/source/source_lcao/setup_deepks.h b/source/source_lcao/setup_deepks.h index 6299ad581b..5e9f05ee33 100644 --- a/source/source_lcao/setup_deepks.h +++ b/source/source_lcao/setup_deepks.h @@ -5,6 +5,8 @@ #include "source_io/module_parameter/input_parameter.h" // Input_para #include "source_basis/module_ao/parallel_orbitals.h" // parallel orbitals #include "source_basis/module_ao/ORB_read.h" // orb +#include "source_basis/module_nao/two_center_integrator.h" // overlap_orb_alpha +#include "source_cell/module_neighbor/sltk_grid_driver.h" // grid driver #ifdef __MLALGO #include "source_lcao/module_deepks/LCAO_deepks.h" // deepks @@ -23,13 +25,34 @@ class Setup_DeePKS LCAO_Deepks ld; #endif + std::string dpks_out_type; + void before_runner( - const UnitCell& ucell, // unitcell + const UnitCell &ucell, // unitcell const int nks, // k points const LCAO_Orbitals &orb, // orbital info Parallel_Orbitals &pv, // parallel orbitals const Input_para &inp); + void build_overlap( + const UnitCell &ucell, + const LCAO_Orbitals &orb, + const Parallel_Orbitals &pv, + const Grid_Driver &gd, + TwoCenterIntegrator &overlap_orb_alpha, + const Input_para &inp); + + void write_forces( + const ModuleBase::matrix &fcs, + const ModuleBase::matrix &fvnl_dalpha, + const Input_para &inp); + + void write_stress( + const ModuleBase::matrix &scs, + const ModuleBase::matrix &svnl_dalpha, + const double &omega, + const Input_para &inp); + }; diff --git a/tests/01_PW/209_PW_DFTHALF/INPUT b/tests/01_PW/209_PW_DFTHALF/INPUT index 9ab7dce297..a1f94218d7 100644 --- a/tests/01_PW/209_PW_DFTHALF/INPUT +++ b/tests/01_PW/209_PW_DFTHALF/INPUT @@ -10,20 +10,19 @@ pseudo_dir ../../PP_ORB #Parameters (2.Iteration) ecutwfc 50 -scf_thr 1e-9 -scf_nmax 100 - +scf_thr 1e-9 +scf_nmax 100 #Parameters (3.Basis) -basis_type pw +basis_type pw #Parameters (4.Smearing) -smearing_method gauss -smearing_sigma 0.002 +smearing_method gauss +smearing_sigma 0.002 #Parameters (5.Mixing) -mixing_type broyden -mixing_beta 0.7 +mixing_type broyden +mixing_beta 0.7 dft_functional PBE