@@ -193,14 +193,35 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
193193 this ->paraMat_ , this ->nk , this ->kv .kvec_d , this ->ucell , this ->gd , this ->orb_cutoff_ );
194194 LR_Util::print_DMR (dm_trans, this ->ucell .nat , " dm_trans of istate " + std::to_string (istate));
195195
196- elecstate::DensityMatrix<T, T> relaxed_diff_dm = // T+D(Z), (R) can be complex
197- LR_Util::build_dm_from_dmk<T, T>(
198- // LR_Util::operator+(
199- 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_ )
200- + 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_ )
201- ,// ),
202- this ->paraMat_ , this ->nk , this ->kv .kvec_d , this ->ucell , this ->gd , this ->orb_cutoff_ );
203- LR_Util::print_DMR (relaxed_diff_dm, this ->ucell .nat , " relaxed_diff_dm of istate " + std::to_string (istate));
196+ // difference density matrix
197+ std::vector<ct::Tensor> dm_diff_k = 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_ );
198+ // std::cout << "dm_diff_k T(k) before symmetrization, istate " + std::to_string(istate) << std::endl;
199+ // LR_Util::print_value(dm_diff_k[0].data<T>(), this->paraMat_.get_col_size(), this->paraMat_.get_row_size());
200+ for (auto & d : dm_diff_k) { LR_Util::matsym (d.data <T>(), this ->nbasis , this ->paraMat_ ); } // symmetrize
201+ // std::cout << "dm_diff_k T(k) after symmetrization, istate " + std::to_string(istate) << std::endl;
202+ // LR_Util::print_value(dm_diff_k[0].data<T>(), this->paraMat_.get_col_size(), this->paraMat_.get_row_size());
203+
204+ std::vector<ct::Tensor> dm_relaxed_k = 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_ );
205+ // std::cout << "dm_relaxed_k Z(k) before symmetrization, istate " + std::to_string(istate) << std::endl;
206+ // LR_Util::print_value(dm_relaxed_k[0].data<T>(), this->paraMat_.get_col_size(), this->paraMat_.get_row_size());
207+ for (auto & d : dm_relaxed_k) { LR_Util::matsym (d.data <T>(), this ->nbasis , this ->paraMat_ ); } // symmetrize
208+ // std::cout << "dm_relaxed_k Z(k) after symmetrization, istate " + std::to_string(istate) << std::endl;
209+ // LR_Util::print_value(dm_relaxed_k[0].data<T>(), this->paraMat_.get_col_size(), this->paraMat_.get_row_size());
210+ // relaxed difference density matrix
211+ std::vector<ct::Tensor> relaxed_diff_dm_k = dm_diff_k + dm_relaxed_k;
212+ elecstate::DensityMatrix<T, T> relaxed_diff_dm =
213+ LR_Util::build_dm_from_dmk<T, T>(relaxed_diff_dm_k,
214+ this ->paraMat_ , this ->nk , this ->kv .kvec_d , this ->ucell , this ->gd , this ->orb_cutoff_ , /* symmetrize=*/ false );
215+ // LR_Util::print_DMR(relaxed_diff_dm, this->ucell.nat, "relaxed_diff_dm T+Z (Z symmetrized) of istate " + std::to_string(istate));
216+
217+ // elecstate::DensityMatrix<T, T> relaxed_diff_dm = // T+D(Z), (R) can be complex
218+ // LR_Util::build_dm_from_dmk<T, T>(
219+ // // LR_Util::operator+(
220+ // cal_dm_diff_pb las(this->X[ispin].template data<T>() + offset, this->paraX_[ispin], c, this->paraC_, this->nbasis, this->nocc[ispin], this->nvirt[ispin], this->paraMat_)
221+ // + 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_)
222+ // ,// ),
223+ // this->paraMat_, this->nk, this->kv.kvec_d, this->ucell, this->gd, this->orb_cutoff_);
224+ // LR_Util::print_DMR(relaxed_diff_dm, this->ucell.nat, "relaxed_diff_dm of istate " + std::to_string(istate));
204225 elecstate::DensityMatrix<T, double > relaxed_diff_dm_real (&this ->paraMat_ , 1 , this ->kv .kvec_d , this ->nk );
205226 LR_Util::initialize_DMR (relaxed_diff_dm_real, this ->paraMat_ , this ->ucell , this ->gd , this ->orb_cutoff_ );
206227 LR_Util::get_DMR_real_imag_part (relaxed_diff_dm, relaxed_diff_dm_real, this ->ucell .nat , ' R' );
@@ -325,11 +346,11 @@ void LR::ESolver_LR<T, TR>::test_force()
325346 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 );
326347 // 2 for spin in f->v. Spin in v->f is already multiplied in the singlet Hartree factor 2.
327348 // / ======================================= END test 2 =========================================
328- if ( this -> nbasis == 2 && ucell. nat == 2 )
329- {
330- // /========================== test 3: H2 SZ 4-center gradients =========================
331- lr_force.cal_H2_sz_center4_grad_hxc (orb_cutoff_);
332- }
349+ // /========================== test 3: H2 SZ 4-center gradients =========================
350+ // if (this->nbasis == 2 && ucell.nat == 2)
351+ // {
352+ // lr_force.cal_H2_sz_center4_grad_hxc(orb_cutoff_);
353+ // }
333354 // exit(0);
334355}
335356
0 commit comments