Skip to content

Commit 95d79bc

Browse files
maki49maki49
committed
add test: 4-center integral derivatives
1 parent e8c24fb commit 95d79bc

File tree

7 files changed

+112
-51
lines changed

7 files changed

+112
-51
lines changed

‎source/module_elecstate/module_dm/density_matrix.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -166,6 +166,7 @@ class DensityMatrix
166166
* please make sure the size of TK* is correct
167167
*/
168168
void set_DMK_pointer(const int ik, TK* DMK_in);
169+
void set_DMK_vector(const int ik, const std::vector<TK>& v) { this->_DMK[ik] = v; }
169170

170171
/**
171172
* @brief get pointer of paraV

‎source/module_lr/Grad/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@ CVCX/CVCX_parallel.cpp
88
CVCX/CVCX_serial.cpp
99
xc/pot_grad_xc.cpp
1010
force/lr_force.cpp
11+
force/lr_force_test.cpp
1112
multipliers/cal_edm_from_multipliers.cpp
1213
esolver_lr_grad.cpp
1314
)

‎source/module_lr/Grad/esolver_lr_grad.cpp‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -191,6 +191,7 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
191191
const elecstate::DensityMatrix<T, T>& dm_trans = // D(X) complex
192192
LR_Util::build_dm_from_dmk<T, T>(dm_trans_k,
193193
this->paraMat_, this->nk, this->kv.kvec_d, this->ucell, this->gd, this->orb_cutoff_);
194+
LR_Util::print_DMR(dm_trans, this->ucell.nat, "dm_trans of istate " + std::to_string(istate));
194195

195196
elecstate::DensityMatrix<T, T> relaxed_diff_dm = // T+D(Z), (R) can be complex
196197
LR_Util::build_dm_from_dmk<T, T>(
@@ -199,6 +200,7 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
199200
+ 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_)
200201
,// ),
201202
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));
202204
elecstate::DensityMatrix<T, double> relaxed_diff_dm_real(&this->paraMat_, 1, this->kv.kvec_d, this->nk);
203205
LR_Util::initialize_DMR(relaxed_diff_dm_real, this->paraMat_, this->ucell, this->gd, this->orb_cutoff_);
204206
LR_Util::get_DMR_real_imag_part(relaxed_diff_dm, relaxed_diff_dm_real, this->ucell.nat, 'R');
@@ -323,7 +325,11 @@ void LR::ESolver_LR<T, TR>::test_force()
323325
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);
324326
// 2 for spin in f->v. Spin in v->f is already multiplied in the singlet Hartree factor 2.
325327
/// ======================================= END test 2 =========================================
326-
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+
}
327333
// exit(0);
328334
}
329335

‎source/module_lr/Grad/force/cal_hs_grad.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@
55
#include "module_basis/module_nao/two_center_bundle.h"
66
#include "module_hamilt_lcao/module_hcontainer/hcontainer.h"
77

8-
void filter_adjs_by_rcut(const UnitCell& ucell,
8+
inline void filter_adjs_by_rcut(const UnitCell& ucell,
99
const int iat0,
1010
AdjacentAtomInfo& adjs)
1111
{
@@ -63,7 +63,7 @@ hamilt::HContainer<TR> build_hcontainer_local_op(const UnitCell& ucell, const Gr
6363
// }
6464

6565

66-
std::vector<hamilt::HContainer<double>> cal_hs_grad(const char job,
66+
inline std::vector<hamilt::HContainer<double>> cal_hs_grad(const char job,
6767
const UnitCell& ucell,
6868
const Parallel_Orbitals& pv,
6969
const Grid_Driver& gd,

‎source/module_lr/Grad/force/lr_force.cpp‎

Lines changed: 0 additions & 46 deletions
Original file line numberDiff line numberDiff line change
@@ -123,52 +123,6 @@ namespace LR
123123
return PulayForceStress::cal_pulay_fs(dm_trans, this->ucell_, &pot_hxc, *this->gint_);
124124
}
125125

126-
template<typename TK>
127-
ModuleBase::matrix LR_Force<TK>::cal_force_overlap_edm(const elecstate::DensityMatrix<TK, double>& edm)
128-
{
129-
// const double* dS[3] = { dSloc_x, dSloc_y, dSloc_z };
130-
std::vector<hamilt::HContainer<double>> dS = cal_hs_grad('S', this->ucell_, this->pv_, this->gd_, this->two_center_bundle_);
131-
// test: output dS
132-
// std::cout << "dS in 3 directions:\n";
133-
// for (int i = 0;i < 3;++i) { LR_Util::print_HR(dS.at(i), this->ucell_.nat, "dS" + std::to_string(i)); }
134-
ModuleBase::matrix foverlap = PulayForceStress::cal_pulay_fs(edm, this->ucell_, dS, -1.);
135-
if (PARAM.inp.test_force)
136-
{
137-
ModuleIO::print_force(GlobalV::ofs_running, this->ucell_, "OVERLAP FORCE (eV/Angstrom)", foverlap, false);
138-
}
139-
return foverlap;
140-
}
141-
142-
template<typename TK>
143-
ModuleBase::matrix LR_Force<TK>::reproduce_force_gs(
144-
const elecstate::DensityMatrix<TK, double>& dm_gs,
145-
const elecstate::DensityMatrix<TK, double>& edm_gs)
146-
{
147-
this->gint_->reset_DMRGint(PARAM.inp.nspin);
148-
// local + Hartree + xc term, including Hellmann-Feynman and Pulay
149-
ModuleBase::matrix f_gs_hf_pulay = cal_force_hamilt_gs_dm_relaxed_diff(dm_gs, dm_gs); // pw+vnl+t_dphi+vl_dphi
150-
// edm term
151-
ModuleBase::matrix f_nonortho = cal_force_overlap_edm(edm_gs); // overlap
152-
this->gint_->reset_DMRGint(1);
153-
return f_gs_hf_pulay + f_nonortho;
154-
}
155-
156-
template<typename TK>
157-
ModuleBase::matrix LR_Force<TK>::reproduce_force_gs_loc(
158-
const elecstate::DensityMatrix<TK, double>& dm_gs,
159-
const elecstate::Potential& pot_gs)
160-
{
161-
this->gint_->reset_DMRGint(PARAM.inp.nspin);
162-
const Charge chr_gs = dm_to_charge(dm_gs);
163-
// local pp (Pulay) + Hartree + xc (grid integration)
164-
ModuleBase::matrix fvl_dphi(this->ucell_.nat, 3);
165-
ModuleBase::matrix stress_tmp; // no use now, only for passing into interfaces
166-
PulayForceStress::cal_pulay_fs(dm_gs.get_DMR_vector().size()/*nspin*/, fvl_dphi, stress_tmp,
167-
dm_gs, this->ucell_, &pot_gs, *this->gint_, true, false);
168-
this->gint_->reset_DMRGint(1);
169-
return fvl_dphi;
170-
}
171-
172126
#ifdef __EXX
173127
template<typename TK>
174128
ModuleBase::matrix LR_Force<TK>::cal_force_exx_dm_trans(

‎source/module_lr/Grad/force/lr_force.h‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -70,6 +70,9 @@ namespace LR
7070
ModuleBase::matrix reproduce_force_gs_loc(const elecstate::DensityMatrix<TK, double>& dm_gs,
7171
const elecstate::Potential& pot_gs);
7272

73+
/// (i^x j | kl)
74+
void cal_H2_sz_center4_grad_hxc(const std::vector<double>& orb_cutoffs);
75+
7376
protected:
7477
const UnitCell& ucell_;
7578
const std::vector<ModuleBase::Vector3<double>>& kvec_d_;
@@ -89,7 +92,5 @@ namespace LR
8992
Charge dm_to_charge(const elecstate::DensityMatrix<TK, double>& dm);
9093
elecstate::Potential dm_to_hxc_potential(const elecstate::DensityMatrix<TK, double>& dm);
9194
elecstate::Potential local_potential();
92-
// probably move frome the ground state?
93-
// void build_dHS()
9495
};
9596
}

‎source/module_lr/Grad/force/lr_force_test.cpp‎

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
#include "lr_force.h"
2+
#include "cal_hs_grad.h"
3+
#include "pulay_force_hcontainer.h"
4+
#include "module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress.h" // only for gint terms
5+
#include "module_lr/utils/lr_util_hcontainer.h"
6+
namespace LR
7+
{
8+
template<typename TK>
9+
ModuleBase::matrix LR_Force<TK>::cal_force_overlap_edm(const elecstate::DensityMatrix<TK, double>& edm)
10+
{
11+
// const double* dS[3] = { dSloc_x, dSloc_y, dSloc_z };
12+
std::vector<hamilt::HContainer<double>> dS = cal_hs_grad('S', this->ucell_, this->pv_, this->gd_, this->two_center_bundle_);
13+
// test: output dS
14+
// std::cout << "dS in 3 directions:\n";
15+
// for (int i = 0;i < 3;++i) { LR_Util::print_HR(dS.at(i), this->ucell_.nat, "dS" + std::to_string(i)); }
16+
ModuleBase::matrix foverlap = PulayForceStress::cal_pulay_fs(edm, this->ucell_, dS, -1.);
17+
if (PARAM.inp.test_force)
18+
{
19+
ModuleIO::print_force(GlobalV::ofs_running, this->ucell_, "OVERLAP FORCE (eV/Angstrom)", foverlap, false);
20+
}
21+
return foverlap;
22+
}
23+
24+
template<typename TK>
25+
ModuleBase::matrix LR_Force<TK>::reproduce_force_gs(
26+
const elecstate::DensityMatrix<TK, double>& dm_gs,
27+
const elecstate::DensityMatrix<TK, double>& edm_gs)
28+
{
29+
this->gint_->reset_DMRGint(PARAM.inp.nspin);
30+
// local + Hartree + xc term, including Hellmann-Feynman and Pulay
31+
ModuleBase::matrix f_gs_hf_pulay = cal_force_hamilt_gs_dm_relaxed_diff(dm_gs, dm_gs); // pw+vnl+t_dphi+vl_dphi
32+
// edm term
33+
ModuleBase::matrix f_nonortho = cal_force_overlap_edm(edm_gs); // overlap
34+
this->gint_->reset_DMRGint(1);
35+
return f_gs_hf_pulay + f_nonortho;
36+
}
37+
38+
template<typename TK>
39+
ModuleBase::matrix LR_Force<TK>::reproduce_force_gs_loc(
40+
const elecstate::DensityMatrix<TK, double>& dm_gs,
41+
const elecstate::Potential& pot_gs)
42+
{
43+
this->gint_->reset_DMRGint(PARAM.inp.nspin);
44+
const Charge chr_gs = dm_to_charge(dm_gs);
45+
// local pp (Pulay) + Hartree + xc (grid integration)
46+
ModuleBase::matrix fvl_dphi(this->ucell_.nat, 3);
47+
ModuleBase::matrix stress_tmp; // no use now, only for passing into interfaces
48+
PulayForceStress::cal_pulay_fs(dm_gs.get_DMR_vector().size()/*nspin*/, fvl_dphi, stress_tmp,
49+
dm_gs, this->ucell_, &pot_gs, *this->gint_, true, false);
50+
this->gint_->reset_DMRGint(1);
51+
return fvl_dphi;
52+
}
53+
54+
template<typename TK>
55+
void LR_Force<TK>::cal_H2_sz_center4_grad_hxc(const std::vector<double>& orb_cutoffs)
56+
{
57+
GlobalV::ofs_running << " ==== Test H2_SZ_CENTER4_HXC (d_x(i) j | kl) ====" << std::endl;
58+
const std::vector<ModuleBase::Vector3<double>>& kvd_test = { ModuleBase::Vector3<double>(0.0, 0.0, 0.0) };
59+
auto init_dm_eff = [&, this](const int i, const int j) -> elecstate::DensityMatrix<TK, double>
60+
{
61+
std::vector<TK> dm_2d(4, 0.0);
62+
std::cout<<"i<<1 + j =" << ((i<<1) + j) << std::endl;
63+
dm_2d[i*2+j] = 1.0;
64+
elecstate::DensityMatrix<TK, double> dm(&this->pv_, 1, kvd_test, 1);
65+
dm.set_DMK_pointer(0, dm_2d.data());
66+
LR_Util::initialize_DMR(dm, this->pv_, this->ucell_, this->gd_, orb_cutoffs);
67+
dm.cal_DMR();
68+
return dm;
69+
};
70+
ModuleBase::matrix stress_tmp; // dummy
71+
this->gint_->reset_DMRGint(1);
72+
73+
for (auto&& i : { 0, 1 })
74+
for (auto&& j : { 0, 1 })
75+
{
76+
elecstate::DensityMatrix<TK, double> dm_ij = init_dm_eff(i, j);
77+
LR_Util::print_DMR(dm_ij, ucell_.nat, "DMR_" + std::to_string(i) + std::to_string(j));
78+
for (auto&& k : { 0, 1 })
79+
for (auto&& l : { 0, 1 })
80+
{
81+
// 1. build dm(kl)
82+
elecstate::DensityMatrix<TK, double> dm_kl = init_dm_eff(k, l);
83+
// 2. build charge & potential
84+
elecstate::Potential pot_hxc_kl = dm_to_hxc_potential(dm_kl);
85+
// 3. cal force
86+
ModuleBase::matrix fvl_dphi(this->ucell_.nat, 3);
87+
PulayForceStress::cal_pulay_fs(1/*nspin*/, fvl_dphi, stress_tmp,
88+
dm_ij, this->ucell_, &pot_hxc_kl, *this->gint_, true, false);
89+
ModuleIO::print_force(GlobalV::ofs_running, this->ucell_,
90+
"H2_SZ_CENTER4_HXC_" + std::to_string(i) + std::to_string(j) + std::to_string(k) + std::to_string(l) + " FORCE (eV/Angstrom)",
91+
fvl_dphi, false);
92+
}
93+
}
94+
}
95+
}
96+
97+
template class LR::LR_Force<double>;
98+
template class LR::LR_Force<std::complex<double>>;

0 commit comments

Comments
 (0)