|
| 1 | +#include "../esolver_lrtd_lcao.h" |
| 2 | +#include "module_lr/Grad/multipliers/zeq_solver.h" |
| 3 | +#include "module_lr/Grad/multipliers/cal_edm_from_multipliers.h" |
| 4 | +#include "module_lr/Grad/force/lr_force.h" |
| 5 | +#include "module_elecstate/module_dm/cal_dm_psi.h" |
| 6 | +#include "module_io/output_log.h" |
| 7 | + |
| 8 | +template <typename Tstream> |
| 9 | +inline void print_force(const std::vector<ModuleBase::matrix>& force, Tstream& ofs) |
| 10 | +{ |
| 11 | + const int nstate = force.size(); |
| 12 | + ofs << "Gradients of each excited state: (eV/Angstrom)" << std::endl; |
| 13 | + ofs << std::setw(6) << "state" << std::setw(6) << "atom" |
| 14 | + << std::setw(15) << "x" << std::setw(15) << "y" << std::setw(15) << "z" << std::endl; |
| 15 | + const double fac = ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A; |
| 16 | + for (int i = 0;i < nstate;++i) |
| 17 | + { |
| 18 | + for (int iat = 0;iat < force[i].nr;++iat) |
| 19 | + { |
| 20 | + std::string istate = iat == 0 ? std::to_string(i) : " "; |
| 21 | + ofs << std::setw(6) << istate << std::setw(6) << iat << std::setw(6) << "force"; |
| 22 | + for (int ixyz = 0;ixyz < 3;++ixyz) { ofs << std::setw(15) << force[i](iat, ixyz) * fac; } |
| 23 | + ofs << std::endl; |
| 24 | + } |
| 25 | + } |
| 26 | +} |
| 27 | +template<typename T, typename TR> |
| 28 | +void LR::ESolver_LR<T, TR>::init_pot_groundstate(const Charge& chg_gs) |
| 29 | +{ |
| 30 | + ModuleBase::TITLE("ESolver_LR", "init_pot_gs"); |
| 31 | + std::vector<std::string> pot_register; |
| 32 | + if (PARAM.inp.vl_in_h) |
| 33 | + { |
| 34 | + //! 11) calculate the structure factor |
| 35 | + this->sf.setup_structure_factor(&ucell, Pgrid, this->pw_rhod); |
| 36 | + this->locpp.init_vloc(this->ucell, this->pw_rho); |
| 37 | + pot_register.push_back("local"); |
| 38 | + } |
| 39 | + if(PARAM.inp.vh_in_h) |
| 40 | + { |
| 41 | + pot_register.push_back("hartree"); |
| 42 | + } |
| 43 | + pot_register.push_back("xc"); |
| 44 | + |
| 45 | + // initialize the ground state potential |
| 46 | + this->pot_gs = LR_Util::make_unique<elecstate::Potential>(this->pw_rhod, this->pw_rho, |
| 47 | + &this->ucell, &this->locpp.vloc, &this->sf, &this->solvent, |
| 48 | + &this->etxc_gs, &this->vtxc_gs); |
| 49 | + this->pot_gs.get()->pot_register(pot_register); |
| 50 | + XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); // set XC type of the ground state |
| 51 | + this->pot_gs->init_pot(0, &chg_gs); // call update_from_charge inside |
| 52 | + if (has_local_xc(this->xc_kernel)) |
| 53 | + { |
| 54 | + XC_Functional::set_xc_type(this->xc_kernel); // recover the excited state xc kernel type |
| 55 | + } |
| 56 | + if (PARAM.inp.test_force) |
| 57 | + { |
| 58 | + this->pot_gs_hartree = LR_Util::make_unique<elecstate::Potential>(this->pw_rhod, this->pw_rho, |
| 59 | + &this->ucell, &this->locpp.vloc, &this->sf, &this->solvent, |
| 60 | + &this->etxc_gs, &this->vtxc_gs); |
| 61 | + this->pot_gs_hartree->pot_register({ "hartree" }); |
| 62 | + this->pot_gs_hartree->init_pot(0, &chg_gs); // call update_from_charge inside |
| 63 | + } |
| 64 | +} |
| 65 | + |
| 66 | +template<typename T, typename TR> |
| 67 | +ct::Tensor LR::ESolver_LR<T, TR>::solve_zvector_eqation(const int ispin) |
| 68 | +{ |
| 69 | + ModuleBase::TITLE("ESolver_LR", "cal_force"); |
| 70 | + ModuleBase::timer::tick("ESolver_LR", "solve_zvector_eqation"); |
| 71 | + ct::Tensor Z = LR_Util::newTensor<T>({ this->nstates, this->nloc_per_band }); |
| 72 | + // construct and solve the Z-vector equation |
| 73 | + Z_vector_equation(this->X[ispin].template data<T>(), Z.template data<T>(), |
| 74 | + this->xc_kernel, this->nstates, this->nspin, this->nbasis, this->nocc, this->nvirt, |
| 75 | + this->ucell, orb_cutoff_, this->gd, *this->psi_ks, this->eig_ks, |
| 76 | +#ifdef __EXX |
| 77 | + std::weak_ptr<Exx_LRI<T>>(this->exx_lri), this->exx_info.info_global.hybrid_alpha, |
| 78 | +#endif |
| 79 | + this-> gint_, std::weak_ptr<PotHxcLR>( this-> pot[ispin]), this-> kv, |
| 80 | + this->paraX_, this->paraC_, this->paraMat_, this->spin_types[ispin]); |
| 81 | + ModuleBase::timer::tick("ESolver_LR", "solve_zvector_eqation"); |
| 82 | + return Z; |
| 83 | +} |
| 84 | + |
| 85 | +template<typename T, typename TR> |
| 86 | +std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin) |
| 87 | +{ |
| 88 | + if (PARAM.inp.test_force) { this->test_force(ispin); } |
| 89 | + |
| 90 | + const ct::Tensor& Z = this->solve_zvector_eqation(ispin); |
| 91 | + |
| 92 | + ModuleBase::TITLE("ESolver_LR", "cal_force"); |
| 93 | + ModuleBase::timer::tick("ESolver_LR", "cal_force"); |
| 94 | + const auto& c = LR_Util::get_psi_spin(*this->psi_ks, ispin, this->nk); // wavefunction coefficients of ground state |
| 95 | + |
| 96 | + // calculate the force (the partial gradient of Lagrangian) |
| 97 | + LR_Force<T> lr_force(this->ucell, this->kv.kvec_d, this->paraMat_, *this->pw_rho, this->locpp, this->sf, this->gd, this->gint_, this->two_center_bundle_); |
| 98 | + GlobalV::ofs_running << "Start to calculate excited-state force of " << this->spin_types[ispin] << std::endl; |
| 99 | + // ground state dm for currrent spin (only for test the correctness of the force) |
| 100 | + // elecstate::DensityMatrix<T, T> dm_gs(this->paraMat_, 1, this->kv.kvec_d, this->nk); |
| 101 | + |
| 102 | + // for each state, calculate dm_trans, dm_relaxed_diff, edm and force |
| 103 | + std::vector<ModuleBase::matrix> forces(this->nstates); |
| 104 | + for (int istate = 0;istate < this->nstates;++istate) |
| 105 | + { |
| 106 | + const int offset = istate * this->nloc_per_band; |
| 107 | + // The imag part will be cancelled in the force calculation, so we use double DM(R) to calculate force. |
| 108 | + // But complex transition DM(R) is still used in energy density matrix calculation. |
| 109 | + elecstate::DensityMatrix<T, double> dm_trans_real = // D(X) |
| 110 | + LR_Util::build_dm_from_dmk<T, double>( |
| 111 | + cal_dm_trans_pblas(this->X[ispin].template data<T>() + offset, this->paraX_[ispin], c, this->paraC_, this->nbasis, this->nocc[ispin], this->nvirt[ispin], this->paraMat_), |
| 112 | + this->paraMat_, this->nk, this->kv.kvec_d, this->ucell, this->gd, this->orb_cutoff_); |
| 113 | + |
| 114 | + elecstate::DensityMatrix<T, T> relaxed_diff_dm = // T+D(Z), (R) can be complex |
| 115 | + LR_Util::build_dm_from_dmk<T, T>( |
| 116 | + // LR_Util::operator+( |
| 117 | + cal_dm_diff_pblas(this->X[ispin].template data<T>() + offset, this->paraX_[ispin], c, this->paraC_, this->nbasis, this->nocc[ispin], this->nvirt[ispin], this->paraMat_) |
| 118 | + + cal_dm_trans_pblas(Z.template data<T>() + offset, this->paraX_[ispin], c, this->paraC_, this->nbasis, this->nocc[ispin], this->nvirt[ispin], this->paraMat_) |
| 119 | + ,// ), |
| 120 | + this->paraMat_, this->nk, this->kv.kvec_d, this->ucell, this->gd, this->orb_cutoff_); |
| 121 | + elecstate::DensityMatrix<T, double> relaxed_diff_dm_real(&this->paraMat_, 1, this->kv.kvec_d, this->nk); |
| 122 | + LR_Util::initialize_DMR(relaxed_diff_dm_real, this->paraMat_, this->ucell, this->gd, this->orb_cutoff_); |
| 123 | + LR_Util::get_DMR_real_imag_part(relaxed_diff_dm, relaxed_diff_dm_real, this->ucell.nat, 'R'); |
| 124 | + |
| 125 | + // get edm of type DensityMatrix |
| 126 | + // weak_ptr here is to avoid "could not match 'weak_ptr' against 'shared_ptr'" |
| 127 | + // but why there're no bug in the previous code (HamiltLR and HamiltULF)? |
| 128 | + // seems because those two are classes having constructors |
| 129 | + // but `cal_edm_from_XZ_istate` here is a functions |
| 130 | + std::weak_ptr<PotHxcLR> pot_weak = this->pot[ispin]; |
| 131 | +#ifdef __EXX |
| 132 | + std::weak_ptr<Exx_LRI<T>> exx_lri_weak = this->exx_lri; |
| 133 | +#endif |
| 134 | + const std::vector<ct::Tensor>& edm_k = |
| 135 | + cal_edm_from_XZ_istate(this->X[ispin].template data<T>() + offset, |
| 136 | + Z.template data<T>() + offset, |
| 137 | + this->pelec->ekb.c[ ispin * nstates + istate], |
| 138 | + // pack the following as a struct or use parameter package |
| 139 | + this->eig_ks.c, relaxed_diff_dm, |
| 140 | + c, this->nspin, this->nbasis, this->nocc, this->nvirt, this->ucell, this->orb_cutoff_, |
| 141 | +#ifdef __EXX |
| 142 | + exx_lri_weak, this->exx_info.info_global.hybrid_alpha, |
| 143 | +#endif |
| 144 | + this->gint_, pot_weak, this->kv, this->gd, this->paraX_, this->paraC_, this->paraMat_, |
| 145 | + has_local_xc(this->xc_kernel)); |
| 146 | + elecstate::DensityMatrix<T, double> edm_real = LR_Util::build_dm_from_dmk<T, double>(edm_k, |
| 147 | + this->paraMat_, this->nk, this->kv.kvec_d, this->ucell, this->gd, this->orb_cutoff_); |
| 148 | + |
| 149 | + ModuleBase::matrix force_hxc_dmtrans = lr_force.cal_force_hxc_dmtrans(dm_trans_real, *this->pot[ispin]); |
| 150 | + std::cout << "Force (Hxc-DMTrans term) of state " << istate << ": " << std::endl; |
| 151 | + LR_Util::print_value(force_hxc_dmtrans.c, ucell.nat, 3); |
| 152 | + ModuleBase::matrix force_hamiltgs_relaxed_diff = lr_force.cal_force_hamilt_gs_dm_relaxed_diff(relaxed_diff_dm_real, *pot_gs); |
| 153 | + std::cout << "Force (GS-(T+Z) term) of state " << istate << ": " << std::endl; |
| 154 | + LR_Util::print_value(force_hamiltgs_relaxed_diff.c, ucell.nat, 3); |
| 155 | + ModuleBase::matrix force_overlap_edm = lr_force.cal_force_overlap_edm(edm_real); // "-" sign has been included in the force factor |
| 156 | + std::cout << "Force (Overlap-EDM term) of state " << istate << ": " << std::endl; |
| 157 | + LR_Util::print_value(force_overlap_edm.c, ucell.nat, 3); |
| 158 | + // remaining for EXX |
| 159 | + forces[istate] = force_hxc_dmtrans + force_hamiltgs_relaxed_diff + force_overlap_edm; |
| 160 | + } |
| 161 | + ModuleBase::timer::tick("ESolver_LR", "cal_force"); |
| 162 | + print_force(forces, std::cout); |
| 163 | + print_force(forces, GlobalV::ofs_running); |
| 164 | + return forces; |
| 165 | +} |
| 166 | + |
| 167 | +template<typename T, typename TR> |
| 168 | +void LR::ESolver_LR<T, TR>::test_force(const int ispin) |
| 169 | +{ |
| 170 | + LR_Force<T> lr_force(this->ucell, this->kv.kvec_d, this->paraMat_, *this->pw_rho, this->locpp, this->sf, this->gd, this->gint_, this->two_center_bundle_); |
| 171 | + |
| 172 | + elecstate::DensityMatrix<T, double> dm_gs(&this->paraMat_, this->nspin, this->kv.kvec_d, this->nk); //DX |
| 173 | + elecstate::cal_dm_psi(&this->paraMat_, this->wg_ks, *this->psi_ks, dm_gs); |
| 174 | + LR_Util::initialize_DMR(dm_gs, this->paraMat_, this->ucell, this->gd, this->orb_cutoff_); |
| 175 | + dm_gs.cal_DMR(); |
| 176 | + |
| 177 | + ///========================== test 1: reproduce the force of ground state ========================= |
| 178 | + // energy density matrix of the ground state |
| 179 | + elecstate::DensityMatrix<T, double> edm_gs(&this->paraMat_, this->nspin, this->kv.kvec_d, this->nk); //DX |
| 180 | + ModuleBase::matrix wg_ekb_ks(nspin, nbands); |
| 181 | + std::transform(this->wg_ks.c, this->wg_ks.c + nspin * nbands, this->eig_ks.c, wg_ekb_ks.c, std::multiplies<double>()); |
| 182 | + elecstate::cal_dm_psi(&this->paraMat_, wg_ekb_ks, *this->psi_ks, edm_gs); |
| 183 | + LR_Util::initialize_DMR(edm_gs, this->paraMat_, this->ucell, this->gd, this->orb_cutoff_); |
| 184 | + edm_gs.cal_DMR(); |
| 185 | + // ground-state force |
| 186 | + ModuleBase::matrix force_gs = lr_force.reproduce_force_gs(dm_gs, edm_gs, *this->pot_gs); |
| 187 | + ModuleIO::print_force(GlobalV::ofs_running, this->ucell, "Ground State FORCE (eV/Angstrom)", force_gs, false); |
| 188 | + /// ======================================= END test 1 ========================================= |
| 189 | + ///========================== test 2: reproduce the DX Hartree term ========================= |
| 190 | + ModuleBase::matrix f_hxc_potgs = lr_force.reproduce_force_gs_loc(dm_gs, *this->pot_gs_hartree); |
| 191 | + ModuleIO::print_force(GlobalV::ofs_running, this->ucell, "GS Hartree force calculated by 'cal_pulay_fs' from potential (eV/Angstrom)", f_hxc_potgs, false); |
| 192 | + ModuleBase::matrix f_hxc_potlr = lr_force.cal_force_hxc_dmtrans(dm_gs, *this->pot[ispin]); |
| 193 | + ModuleIO::print_force(GlobalV::ofs_running, this->ucell, "2* GS Hxc force calculated by 'LR_Force' from kernel (eV/Angstrom)", f_hxc_potlr * 2, false); |
| 194 | + // 2 for spin in f->v. Spin in v->f is already multiplied in the singlet Hartree factor 2. |
| 195 | + /// ======================================= END test 2 ========================================= |
| 196 | + |
| 197 | + // exit(0); |
| 198 | +} |
| 199 | + |
| 200 | +template class LR::ESolver_LR<double, double>; |
| 201 | +template class LR::ESolver_LR<std::complex<double>, double>; |
0 commit comments