Skip to content

Commit 51657ad

Browse files
JGHan7JGHan7
committed
last save befor move cal_V_ TV/hartree/XC from rdmft.cpp to cal_V_rdmft.cpp
1 parent 73a4eae commit 51657ad

File tree

5 files changed

+622
-277
lines changed

5 files changed

+622
-277
lines changed

‎abaInstall_HZWpara.sh‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
rm -rf build
22

3-
cmake -B build -DCMAKE_INSTALL_PREFIX=/public1/home/t6s000394/jghan/software/abacus-develop/rdmft-abacus/ -DCMAKE_CXX_COMPILER=icpx -DELPA_DIR=/public1/home/t6s000394/jghan/software/elpa-2021.11/ -DCEREAL_INCLUDE_DIR=/public1/home/t6s000394/jghan/software/cereal/cereal-1.3.2/include
4-
#cmake -B build -DCMAKE_INSTALL_PREFIX=/public1/home/t6s000394/jghan/software/abacus-develop/rdmft-abacus/ -DCMAKE_CXX_COMPILER=icpx -DMPI_CXX_COMPILER=mpiicpc -DELPA_DIR=/public1/home/t6s000394/jghan/software/elpa-2021.11/ -DLibxc_DIR=/public1/home/t6s000394/jghan/software/libxc/ -DLIBRI_DIR=/public1/home/t6s000394/jghan/software/LibRI -DLIBCOMM_DIR=/public1/home/t6s000394/jghan/software/LibComm -DCEREAL_INCLUDE_DIR=/public1/home/t6s000394/jghan/software/cereal/cereal-1.3.2/include
3+
#cmake -B build -DCMAKE_INSTALL_PREFIX=/public1/home/t6s000394/jghan/software/abacus-develop/rdmft-abacus/ -DCMAKE_CXX_COMPILER=icpx -DELPA_DIR=/public1/home/t6s000394/jghan/software/elpa-2021.11/ -DCEREAL_INCLUDE_DIR=/public1/home/t6s000394/jghan/software/cereal/cereal-1.3.2/include
4+
cmake -B build -DCMAKE_INSTALL_PREFIX=/public1/home/t6s000394/jghan/software/abacus-develop/rdmft-abacus/ -DCMAKE_CXX_COMPILER=icpx -DMPI_CXX_COMPILER=mpiicpc -DELPA_DIR=/public1/home/t6s000394/jghan/software/elpa-2021.11/ -DLibxc_DIR=/public1/home/t6s000394/jghan/software/libxc/ -DLIBRI_DIR=/public1/home/t6s000394/jghan/software/LibRI -DLIBCOMM_DIR=/public1/home/t6s000394/jghan/software/LibComm -DCEREAL_INCLUDE_DIR=/public1/home/t6s000394/jghan/software/cereal/cereal-1.3.2/include
55

66
#cmake --build build -j 52 2>job.err
77
cmake --build build -j 8

‎source/Makefile.Objects‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -736,3 +736,4 @@ OBJS_TENSOR=tensor.o\
736736

737737
OBJS_RDMFT=rdmft.o\
738738
rdmft_tools.o\
739+
cal_V_rdmft.o\

‎source/module_rdmft/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@ if(ENABLE_LCAO)
44
OBJECT
55
rdmft.cpp
66
rdmft_tools.cpp
7+
cal_V_rdmft.cpp
78
)
89
endif()
910

‎source/module_rdmft/cal_V_rdmft.cpp‎

Lines changed: 343 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,343 @@
1+
//==========================================================
2+
// Author: Jingang Han
3+
// DATE : 2024-10-30
4+
//==========================================================
5+
6+
7+
8+
#include "module_hamilt_general/operator.h"
9+
#include "module_hamilt_lcao/module_hcontainer/hcontainer.h"
10+
#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h"
11+
#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h"
12+
#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.h"
13+
#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.h"
14+
#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h"
15+
16+
17+
// used by Exx&LRI
18+
#ifdef __EXX
19+
#include "module_ri/RI_2D_Comm.h"
20+
#include "module_ri/Exx_LRI.h"
21+
#endif
22+
23+
24+
25+
26+
namespace rdmft
27+
{
28+
29+
30+
31+
template <typename TK, typename TR>
32+
void RDMFT<TK, TR>::cal_V_TV()
33+
{
34+
HR_TV->set_zero();
35+
36+
V_ekinetic_potential = new hamilt::EkineticNew<hamilt::OperatorLCAO<TK, TR>>(
37+
hsk_TV,
38+
kv->kvec_d,
39+
HR_TV,
40+
&GlobalC::ucell,
41+
orb->cutoffs(),
42+
&GlobalC::GridD,
43+
two_center_bundle->kinetic_orb.get()
44+
);
45+
46+
V_nonlocal = new hamilt::NonlocalNew<hamilt::OperatorLCAO<TK, TR>>(
47+
hsk_TV,
48+
kv->kvec_d,
49+
HR_TV,
50+
&GlobalC::ucell,
51+
orb->cutoffs(),
52+
&GlobalC::GridD,
53+
two_center_bundle->overlap_orb_beta.get()
54+
);
55+
56+
if( PARAM.inp.gamma_only )
57+
{
58+
V_local = new rdmft::Veff_rdmft<TK,TR>(
59+
GG,
60+
hsk_TV,
61+
kv->kvec_d,
62+
this->pelec->pot,
63+
HR_TV,
64+
&GlobalC::ucell,
65+
orb->cutoffs(),
66+
&GlobalC::GridD,
67+
nspin,
68+
69+
charge,
70+
rho_basis,
71+
vloc,
72+
sf,
73+
"local"
74+
);
75+
}
76+
else
77+
{
78+
V_local = new rdmft::Veff_rdmft<TK,TR>(
79+
GK,
80+
hsk_TV,
81+
kv->kvec_d,
82+
this->pelec->pot,
83+
HR_TV,
84+
&GlobalC::ucell,
85+
orb->cutoffs(),
86+
&GlobalC::GridD,
87+
nspin,
88+
89+
charge,
90+
rho_basis,
91+
vloc,
92+
sf,
93+
"local"
94+
);
95+
}
96+
97+
// update HR_TV in ion-step, now HR_TV has the HR of V_ekinetic + V_nonlcao + V_local
98+
V_ekinetic_potential->contributeHR();
99+
V_nonlocal->contributeHR();
100+
V_local->contributeHR();
101+
102+
}
103+
104+
105+
template <typename TK, typename TR>
106+
void RDMFT<TK, TR>::cal_V_hartree()
107+
{
108+
HR_hartree->set_zero();
109+
110+
if( PARAM.inp.gamma_only )
111+
{
112+
V_hartree = new rdmft::Veff_rdmft<TK,TR>(
113+
GG,
114+
hsk_hartree,
115+
kv->kvec_d,
116+
this->pelec->pot,
117+
HR_hartree,
118+
&GlobalC::ucell,
119+
orb->cutoffs(),
120+
&GlobalC::GridD,
121+
nspin,
122+
123+
charge,
124+
rho_basis,
125+
vloc,
126+
sf,
127+
"hartree"
128+
);
129+
}
130+
else
131+
{
132+
// this can be optimized, use potHartree.update_from_charge()
133+
V_hartree = new rdmft::Veff_rdmft<TK,TR>(
134+
GK,
135+
hsk_hartree,
136+
kv->kvec_d,
137+
this->pelec->pot,
138+
HR_hartree,
139+
&GlobalC::ucell,
140+
orb->cutoffs(),
141+
&GlobalC::GridD,
142+
nspin,
143+
144+
charge,
145+
rho_basis,
146+
vloc,
147+
sf,
148+
"hartree"
149+
);
150+
}
151+
152+
// in gamma only, must calculate HR_hartree before HR_local
153+
// HR_exx_XC get from another way, so don't need to do this
154+
V_hartree->contributeHR();
155+
156+
// // update HR_local in e-step, now HR_TV has the HR of V_ekinetic + V_nonlcao + V_local,
157+
// V_local->contributeHR();
158+
// HR_local->add(*HR_TV); // now HR_local has the HR of V_ekinetic + V_nonlcao + V_local
159+
160+
}
161+
162+
163+
template <typename TK, typename TR>
164+
void RDMFT<TK, TR>::cal_V_XC()
165+
{
166+
// // //test
167+
// DM_XC_pass = DM_XC;
168+
169+
// elecstate::DensityMatrix<TK, double> DM_test(ParaV, nspin, kv->kvec_d, nk_total);
170+
// elecstate::cal_dm_psi(ParaV, wg, wfc, DM_test);
171+
// DM_test.init_DMR(&GlobalC::GridD, &GlobalC::ucell);
172+
// DM_test.cal_DMR();
173+
174+
// // compare DM_XC and DM get in update_charge(or ABACUS)
175+
// std::cout << "\n\ntest DM_XC - DM in ABACUS: \n" << std::endl;
176+
// double DM_XC_minus_DMtest = 0.0;
177+
// for(int ik=0; ik<nk_total; ++ik)
178+
// {
179+
// TK* dmk_pointer = DM_test.get_DMK_pointer(ik);
180+
// for(int iloc=0; iloc<ParaV->nloc; ++iloc)
181+
// {
182+
// double test = std::abs(DM_XC[ik][iloc] - dmk_pointer[iloc]);
183+
// DM_XC_minus_DMtest += test;
184+
// if( test > 1e-16 )
185+
// {
186+
// std::cout << "\nik, iloc, minus[ik][iloc]: " << ik << " " << iloc << " " << test << std::endl;
187+
// }
188+
// }
189+
// }
190+
// std::cout << "\nsum of DM_XC - DM in ABACUS: " << DM_XC_minus_DMtest << std::endl;
191+
192+
if( !only_exx_type )
193+
{
194+
HR_dft_XC->set_zero();
195+
if( PARAM.inp.gamma_only )
196+
{
197+
// this can be optimized, use potXC.update_from_charge()
198+
V_dft_XC = new rdmft::Veff_rdmft<TK,TR>(
199+
GG,
200+
hsk_dft_XC,
201+
kv->kvec_d,
202+
this->pelec->pot,
203+
HR_dft_XC,
204+
&GlobalC::ucell,
205+
orb->cutoffs(),
206+
&GlobalC::GridD,
207+
nspin,
208+
209+
charge,
210+
rho_basis,
211+
vloc,
212+
sf,
213+
"xc",
214+
&etxc,
215+
&vtxc
216+
);
217+
}
218+
else
219+
{
220+
// this can be optimized, use potXC.update_from_charge()
221+
V_dft_XC = new rdmft::Veff_rdmft<TK,TR>(
222+
GK,
223+
hsk_dft_XC,
224+
kv->kvec_d,
225+
this->pelec->pot,
226+
HR_dft_XC,
227+
&GlobalC::ucell,
228+
orb->cutoffs(),
229+
&GlobalC::GridD,
230+
nspin,
231+
232+
charge,
233+
rho_basis,
234+
vloc,
235+
sf,
236+
"xc",
237+
&etxc,
238+
&vtxc
239+
);
240+
}
241+
V_dft_XC->contributeHR();
242+
}
243+
244+
#ifdef __EXX
245+
if(GlobalC::exx_info.info_global.cal_exx)
246+
{
247+
HR_exx_XC->set_zero();
248+
249+
std::vector< std::vector<TK> > DM_XC(nk_total, std::vector<TK>(ParaV->nloc));
250+
get_DM_XC(DM_XC);
251+
// get DM_XC of all k points
252+
if( exx_spacegroup_symmetry )
253+
{
254+
DM_XC = symrot_exx.restore_dm(*this->kv, DM_XC, *ParaV); // class vector could be auto resize()
255+
}
256+
std::vector< const std::vector<TK>* > DM_XC_pointer(DM_XC.size());
257+
for(int ik=0; ik<DM_XC.size(); ++ik) DM_XC_pointer[ik] = &DM_XC[ik];
258+
259+
if (GlobalC::exx_info.info_ri.real_number)
260+
{
261+
// transfer the DM_XC to appropriate format
262+
std::vector<std::map<int,std::map<std::pair<int,std::array<int,3>>,RI::Tensor<double>>>>
263+
Ds_XC_d = std::is_same<TK, double>::value //gamma_only_local
264+
? RI_2D_Comm::split_m2D_ktoR<double>(*kv, DM_XC_pointer, *ParaV, nspin)
265+
: RI_2D_Comm::split_m2D_ktoR<double>(*kv, DM_XC_pointer, *ParaV, nspin, this->exx_spacegroup_symmetry);
266+
267+
// provide the Ds_XC to Vxc_fromRI(V_exx_XC)
268+
if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace)
269+
{
270+
Vxc_fromRI_d->cal_exx_elec(Ds_XC_d, *ParaV, &this->symrot_exx);
271+
}
272+
else
273+
{
274+
Vxc_fromRI_d->cal_exx_elec(Ds_XC_d, *ParaV);
275+
}
276+
277+
// when we doing V_exx_XC.contributeHk(ik), we get HK_XC constructed by the special DM_XC
278+
V_exx_XC = new hamilt::OperatorEXX<hamilt::OperatorLCAO<TK, TR>>(
279+
hsk_exx_XC,
280+
HR_exx_XC,
281+
*kv,
282+
&Vxc_fromRI_d->Hexxs,
283+
nullptr,
284+
hamilt::Add_Hexx_Type::k
285+
);
286+
}
287+
else
288+
{
289+
// transfer the DM_XC to appropriate format
290+
std::vector<std::map<int,std::map<std::pair<int,std::array<int,3>>,RI::Tensor<std::complex<double>>>>>
291+
Ds_XC_c = std::is_same<TK, double>::value //gamma_only_local
292+
? RI_2D_Comm::split_m2D_ktoR<std::complex<double>>(*kv, DM_XC_pointer, *ParaV, nspin)
293+
: RI_2D_Comm::split_m2D_ktoR<std::complex<double>>(*kv, DM_XC_pointer, *ParaV, nspin, this->exx_spacegroup_symmetry);
294+
295+
// // provide the Ds_XC to Vxc_fromRI(V_exx_XC)
296+
if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace)
297+
{
298+
Vxc_fromRI_c->cal_exx_elec(Ds_XC_c, *ParaV, &this->symrot_exx);
299+
}
300+
else
301+
{
302+
Vxc_fromRI_c->cal_exx_elec(Ds_XC_c, *ParaV);
303+
}
304+
305+
// when we doing V_exx_XC.contributeHk(ik), we get HK_XC constructed by the special DM_XC
306+
V_exx_XC = new hamilt::OperatorEXX<hamilt::OperatorLCAO<TK, TR>>(
307+
hsk_exx_XC,
308+
HR_exx_XC,
309+
*kv,
310+
nullptr,
311+
&Vxc_fromRI_c->Hexxs,
312+
hamilt::Add_Hexx_Type::k
313+
);
314+
}
315+
// use hamilt::Add_Hexx_Type::k, not R, contributeHR() should be skipped
316+
// V_exx_XC->contributeHR();
317+
}
318+
#endif
319+
320+
}
321+
322+
323+
324+
325+
326+
327+
328+
329+
330+
331+
332+
333+
334+
335+
336+
337+
338+
339+
340+
341+
342+
343+
}

0 commit comments

Comments
 (0)