Skip to content

Commit 3790620

Browse files
maki49maki49
committed
initial version of writing H and dH terms (except exx)
1 parent 67a0e88 commit 3790620

23 files changed

Lines changed: 1955 additions & 4 deletions

File tree

‎source/source_io/CMakeLists.txt‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -82,10 +82,17 @@ if(ENABLE_LCAO)
8282
module_mulliken/output_mulliken.cpp
8383
module_ml/io_npz.cpp
8484
module_hs/cal_pLpR.cpp
85+
module_dhs/write_dH.cpp
86+
module_dhs/write_dH_t.cpp
87+
module_dhs/write_dH_vnl.cpp
88+
module_dhs/write_dH_vl.cpp
89+
module_dhs/write_dH_vh.cpp
90+
module_dhs/write_dH_vxc.cpp
8591
)
8692
list(APPEND objects_advanced
8793
module_unk/unk_overlap_lcao.cpp
8894
module_hs/write_HS_R.cpp
95+
module_hs/write_H_terms.cpp
8996
module_hs/write_HS_sparse.cpp
9097
module_hs/single_R_io.cpp
9198
module_hs/cal_r_overlap_R.cpp

‎source/source_io/module_ctrl/ctrl_scf_lcao.cpp‎

Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,8 @@
1111
#include "../module_unk/berryphase.h" // use berryphase
1212
#include "../module_hs/cal_pLpR.h" // use AngularMomentumCalculator()
1313
#include "source_io/module_hs/output_mat_sparse.h" // use ModuleIO::output_mat_sparse()
14+
#include "source_io/module_dhs/write_dH.h" // use ModuleIO::write_dH_components()
15+
#include "source_io/module_hs/write_H_terms.h" // use ModuleIO::write_h_*
1416
#include "../module_hs/write_HS_R.h" // use ModuleIO::write_hsr()
1517
#include "../module_mulliken/output_mulliken.h" // use cal_mag()
1618
#include "../module_wannier/to_wannier90_lcao.h" // use toWannier90_LCAO
@@ -248,6 +250,62 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell,
248250
p_ham_tk,
249251
&dftu);
250252

253+
//------------------------------------------------------------------
254+
//! 7c) Output dH components (dT/dR, dV^NL/dR, dV^L/dR, dV^H/dR, dV^XC/dR)
255+
//------------------------------------------------------------------
256+
{
257+
WriteDHParams dh_params;
258+
dh_params.ucell = &ucell;
259+
dh_params.gd = &gd;
260+
dh_params.pv = &pv;
261+
dh_params.two_center_bundle = &two_center_bundle;
262+
dh_params.orb = &orb;
263+
dh_params.kv = &kv;
264+
dh_params.v_eff = &pelec->pot->get_eff_v();
265+
dh_params.iat2iwt = ucell.get_iat2iwt();
266+
dh_params.nat = ucell.nat;
267+
dh_params.nspin = inp.nspin;
268+
dh_params.istep = istep;
269+
dh_params.gamma_only = gamma_only;
270+
dh_params.append = out_app_flag;
271+
if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2)
272+
{
273+
dh_params.dmR = dm->get_DMR_pointer(1);
274+
}
275+
ModuleIO::write_dH_components(dh_params);
276+
}
277+
278+
279+
//------------------------------------------------------------------
280+
//! 7d) Output H components (T, Vnl, Vl, Vh, Vxc)
281+
//------------------------------------------------------------------
282+
{
283+
if (inp.out_mat_h_t[0])
284+
{
285+
ModuleIO::write_h_t(ucell, gd, pv, two_center_bundle, orb, kv, nspin, istep, out_app_flag,
286+
ucell.get_iat2iwt(), ucell.nat);
287+
}
288+
if (inp.out_mat_h_vnl[0])
289+
{
290+
ModuleIO::write_h_vnl(ucell, gd, pv, two_center_bundle, orb, kv, nspin, istep, out_app_flag,
291+
ucell.get_iat2iwt(), ucell.nat);
292+
}
293+
if (inp.out_mat_h_vl[0])
294+
{
295+
ModuleIO::write_h_vl(ucell, gd, pv, orb, pelec->pot, nspin, istep, out_app_flag,
296+
ucell.get_iat2iwt(), ucell.nat, kv);
297+
}
298+
if (inp.out_mat_h_vh[0])
299+
{
300+
ModuleIO::write_h_vh(ucell, gd, pv, orb, pelec->charge, pw_rho, nspin, istep, out_app_flag,
301+
ucell.get_iat2iwt(), ucell.nat, kv);
302+
}
303+
if (inp.out_mat_h_vxc[0])
304+
{
305+
ModuleIO::write_h_vxc(ucell, gd, pv, orb, pelec->charge, pw_rho->nrxx, nspin, istep,
306+
out_app_flag, ucell.get_iat2iwt(), ucell.nat, kv);
307+
}
308+
}
251309
//------------------------------------------------------------------
252310
//! 8) Output kinetic matrix
253311
//------------------------------------------------------------------

‎source/source_io/module_dhs/write_dH.cpp‎

Lines changed: 142 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,142 @@
1+
#include "write_dH.h"
2+
3+
#include "source_base/timer.h"
4+
#include "source_io/module_hs/write_HS_R.h"
5+
#include "source_io/module_output/ucell_io.h"
6+
#include "source_io/module_parameter/parameter.h"
7+
#include "source_lcao/module_hcontainer/hcontainer_funcs.h"
8+
#include "source_lcao/module_hcontainer/output_hcontainer.h"
9+
10+
#include <fstream>
11+
#include <iomanip>
12+
13+
namespace ModuleIO
14+
{
15+
16+
bool any_dh_term_enabled()
17+
{
18+
return PARAM.inp.out_mat_dh_t[0] || PARAM.inp.out_mat_dh_vl[0] || PARAM.inp.out_mat_dh_vnl[0]
19+
|| PARAM.inp.out_mat_dh_vh[0] || PARAM.inp.out_mat_dh_vxc[0];
20+
}
21+
22+
void write_dH_components(WriteDHParams& params)
23+
{
24+
ModuleBase::TITLE("ModuleIO", "write_dH_components");
25+
ModuleBase::timer::start("ModuleIO", "write_dH_components");
26+
27+
if (!any_dh_term_enabled())
28+
{
29+
ModuleBase::timer::end("ModuleIO", "write_dH_components");
30+
return;
31+
}
32+
33+
GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl;
34+
GlobalV::ofs_running << " | |" << std::endl;
35+
GlobalV::ofs_running << " | #Print out dH/dR components# |" << std::endl;
36+
GlobalV::ofs_running << " | |" << std::endl;
37+
GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl;
38+
39+
if (PARAM.inp.out_mat_dh_t[0])
40+
{
41+
write_dH_t(params);
42+
}
43+
44+
if (PARAM.inp.out_mat_dh_vnl[0])
45+
{
46+
write_dH_vnl(params);
47+
}
48+
49+
if (PARAM.inp.out_mat_dh_vl[0])
50+
{
51+
write_dH_vl(params);
52+
}
53+
54+
if (PARAM.inp.out_mat_dh_vh[0])
55+
{
56+
write_dH_vh(params);
57+
}
58+
59+
if (PARAM.inp.out_mat_dh_vxc[0])
60+
{
61+
write_dH_vxc(params);
62+
}
63+
64+
if (PARAM.inp.out_mat_dh[0] && any_dh_term_enabled())
65+
{
66+
write_dH_sum(params);
67+
}
68+
69+
ModuleBase::timer::end("ModuleIO", "write_dH_components");
70+
}
71+
72+
bool write_dH_sum(WriteDHParams& params)
73+
{
74+
ModuleBase::TITLE("ModuleIO", "write_dH_sum");
75+
ModuleBase::timer::start("ModuleIO", "write_dH_sum");
76+
77+
const UnitCell& ucell = *params.ucell;
78+
const Parallel_Orbitals& pv = *params.pv;
79+
const int nat = ucell.nat;
80+
const int nspin = params.nspin;
81+
82+
for (int ispin = 0; ispin < (nspin == 2 ? 2 : 1); ispin++)
83+
{
84+
hamilt::HContainer<double> dH_sum_x(&pv);
85+
hamilt::HContainer<double> dH_sum_y(&pv);
86+
hamilt::HContainer<double> dH_sum_z(&pv);
87+
dH_sum_x.set_zero();
88+
dH_sum_y.set_zero();
89+
dH_sum_z.set_zero();
90+
91+
const int nbasis = dH_sum_x.get_nbasis();
92+
93+
#ifdef __MPI
94+
Parallel_Orbitals serialV;
95+
serialV.init(nbasis, nbasis, nbasis, pv.comm());
96+
serialV.set_serial(nbasis, nbasis);
97+
serialV.set_atomic_trace(params.iat2iwt, nat, nbasis);
98+
99+
hamilt::HContainer<double> dH_x_s(&serialV);
100+
hamilt::HContainer<double> dH_y_s(&serialV);
101+
hamilt::HContainer<double> dH_z_s(&serialV);
102+
hamilt::gatherParallels(dH_sum_x, &dH_x_s, 0);
103+
hamilt::gatherParallels(dH_sum_y, &dH_y_s, 0);
104+
hamilt::gatherParallels(dH_sum_z, &dH_z_s, 0);
105+
106+
if (GlobalV::MY_RANK == 0)
107+
#endif
108+
{
109+
std::string fname_x = ModuleIO::dhr_gen_fname("dhrx", ispin, params.append, params.istep);
110+
std::string fname_y = ModuleIO::dhr_gen_fname("dhry", ispin, params.append, params.istep);
111+
std::string fname_z = ModuleIO::dhr_gen_fname("dhrz", ispin, params.append, params.istep);
112+
113+
if (PARAM.inp.calculation == "md" && !PARAM.inp.out_app_flag)
114+
{
115+
fname_x = PARAM.globalv.global_matrix_dir + fname_x;
116+
fname_y = PARAM.globalv.global_matrix_dir + fname_y;
117+
fname_z = PARAM.globalv.global_matrix_dir + fname_z;
118+
}
119+
else
120+
{
121+
fname_x = PARAM.globalv.global_out_dir + fname_x;
122+
fname_y = PARAM.globalv.global_out_dir + fname_y;
123+
fname_z = PARAM.globalv.global_out_dir + fname_z;
124+
}
125+
126+
#ifdef __MPI
127+
ModuleIO::write_hcontainer_csr(fname_x, &ucell, 8, &dH_x_s, params.istep, ispin, nspin, "dH_sum");
128+
ModuleIO::write_hcontainer_csr(fname_y, &ucell, 8, &dH_y_s, params.istep, ispin, nspin, "dH_sum");
129+
ModuleIO::write_hcontainer_csr(fname_z, &ucell, 8, &dH_z_s, params.istep, ispin, nspin, "dH_sum");
130+
#else
131+
ModuleIO::write_hcontainer_csr(fname_x, &ucell, 8, &dH_sum_x, params.istep, ispin, nspin, "dH_sum");
132+
ModuleIO::write_hcontainer_csr(fname_y, &ucell, 8, &dH_sum_y, params.istep, ispin, nspin, "dH_sum");
133+
ModuleIO::write_hcontainer_csr(fname_z, &ucell, 8, &dH_sum_z, params.istep, ispin, nspin, "dH_sum");
134+
#endif
135+
}
136+
}
137+
138+
ModuleBase::timer::end("ModuleIO", "write_dH_sum");
139+
return true;
140+
}
141+
142+
} // namespace ModuleIO

‎source/source_io/module_dhs/write_dH.h‎

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
#ifndef WRITE_DH_H
2+
#define WRITE_DH_H
3+
4+
#include "source_basis/module_nao/two_center_bundle.h"
5+
#include "source_cell/klist.h"
6+
#include "source_cell/module_neighbor/sltk_grid_driver.h"
7+
#include "source_lcao/LCAO_domain.h"
8+
#include "source_lcao/module_hcontainer/hcontainer.h"
9+
10+
#include <vector>
11+
12+
namespace ModuleIO
13+
{
14+
15+
struct WriteDHParams
16+
{
17+
const UnitCell* ucell = nullptr;
18+
const Grid_Driver* gd = nullptr;
19+
const Parallel_Orbitals* pv = nullptr;
20+
const TwoCenterBundle* two_center_bundle = nullptr;
21+
const LCAO_Orbitals* orb = nullptr;
22+
const K_Vectors* kv = nullptr;
23+
const ModuleBase::matrix* v_eff = nullptr;
24+
const int* iat2iwt = nullptr;
25+
int nat = 0;
26+
int nspin = 1;
27+
int istep = 0;
28+
bool gamma_only = false;
29+
bool append = false;
30+
const hamilt::HContainer<double>* dmR = nullptr;
31+
};
32+
33+
bool any_dh_term_enabled();
34+
35+
void write_dH_components(WriteDHParams& params);
36+
37+
bool write_dH_t(WriteDHParams& params);
38+
39+
bool write_dH_vnl(WriteDHParams& params);
40+
41+
bool write_dH_vl(WriteDHParams& params);
42+
43+
bool write_dH_vh(WriteDHParams& params);
44+
45+
bool write_dH_vxc(WriteDHParams& params);
46+
47+
bool write_dH_sum(WriteDHParams& params);
48+
49+
} // namespace ModuleIO
50+
51+
#endif

‎source/source_io/module_dhs/write_dH_t.cpp‎

Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
1+
#include "source_base/timer.h"
2+
#include "source_io/module_hs/write_HS_R.h"
3+
#include "source_io/module_output/ucell_io.h"
4+
#include "source_io/module_parameter/parameter.h"
5+
#include "source_lcao/module_hcontainer/hcontainer_funcs.h"
6+
#include "source_lcao/module_hcontainer/output_hcontainer.h"
7+
#include "source_lcao/module_operator_lcao/ekinetic.h"
8+
#include "source_lcao/module_operator_lcao/operator_force_stress_utils.h"
9+
#include "write_dH.h"
10+
11+
namespace ModuleIO
12+
{
13+
14+
bool write_dH_t(WriteDHParams& params)
15+
{
16+
ModuleBase::TITLE("ModuleIO", "write_dH_t");
17+
ModuleBase::timer::start("ModuleIO", "write_dH_t");
18+
19+
const UnitCell& ucell = *params.ucell;
20+
const Grid_Driver& gd = *params.gd;
21+
const Parallel_Orbitals& pv = *params.pv;
22+
const TwoCenterBundle& two_center_bundle = *params.two_center_bundle;
23+
const LCAO_Orbitals& orb = *params.orb;
24+
25+
const int nat = ucell.nat;
26+
const int nspin = params.nspin;
27+
const std::vector<double>& orb_cutoff = orb.cutoffs();
28+
29+
for (int ispin = 0; ispin < (nspin == 2 ? 2 : 1); ispin++)
30+
{
31+
hamilt::HContainer<double> dT_x(&pv);
32+
hamilt::HContainer<double> dT_y(&pv);
33+
hamilt::HContainer<double> dT_z(&pv);
34+
35+
hamilt::EKinetic<hamilt::OperatorLCAO<double, double>> tmp_ekinetic(nullptr,
36+
params.kv->kvec_d,
37+
nullptr,
38+
&ucell,
39+
orb_cutoff,
40+
&gd,
41+
two_center_bundle.kinetic_orb.get());
42+
43+
tmp_ekinetic.cal_dH(&dT_x, &dT_y, &dT_z);
44+
45+
const int nbasis = dT_x.get_nbasis();
46+
47+
#ifdef __MPI
48+
Parallel_Orbitals serialV;
49+
serialV.init(nbasis, nbasis, nbasis, pv.comm());
50+
serialV.set_serial(nbasis, nbasis);
51+
serialV.set_atomic_trace(params.iat2iwt, nat, nbasis);
52+
53+
hamilt::HContainer<double> dT_x_serial(&serialV);
54+
hamilt::HContainer<double> dT_y_serial(&serialV);
55+
hamilt::HContainer<double> dT_z_serial(&serialV);
56+
hamilt::gatherParallels(dT_x, &dT_x_serial, 0);
57+
hamilt::gatherParallels(dT_y, &dT_y_serial, 0);
58+
hamilt::gatherParallels(dT_z, &dT_z_serial, 0);
59+
60+
if (GlobalV::MY_RANK == 0)
61+
#endif
62+
{
63+
std::string fname_x = ModuleIO::dhr_gen_fname("dtrx", ispin, params.append, params.istep);
64+
std::string fname_y = ModuleIO::dhr_gen_fname("dtry", ispin, params.append, params.istep);
65+
std::string fname_z = ModuleIO::dhr_gen_fname("dtrz", ispin, params.append, params.istep);
66+
67+
if (PARAM.inp.calculation == "md" && !PARAM.inp.out_app_flag)
68+
{
69+
fname_x = PARAM.globalv.global_matrix_dir + fname_x;
70+
fname_y = PARAM.globalv.global_matrix_dir + fname_y;
71+
fname_z = PARAM.globalv.global_matrix_dir + fname_z;
72+
}
73+
else
74+
{
75+
fname_x = PARAM.globalv.global_out_dir + fname_x;
76+
fname_y = PARAM.globalv.global_out_dir + fname_y;
77+
fname_z = PARAM.globalv.global_out_dir + fname_z;
78+
}
79+
80+
#ifdef __MPI
81+
ModuleIO::write_hcontainer_csr(fname_x, &ucell, 8, &dT_x_serial, params.istep, ispin, nspin, "dT");
82+
ModuleIO::write_hcontainer_csr(fname_y, &ucell, 8, &dT_y_serial, params.istep, ispin, nspin, "dT");
83+
ModuleIO::write_hcontainer_csr(fname_z, &ucell, 8, &dT_z_serial, params.istep, ispin, nspin, "dT");
84+
#else
85+
ModuleIO::write_hcontainer_csr(fname_x, &ucell, 8, &dT_x, params.istep, ispin, nspin, "dT");
86+
ModuleIO::write_hcontainer_csr(fname_y, &ucell, 8, &dT_y, params.istep, ispin, nspin, "dT");
87+
ModuleIO::write_hcontainer_csr(fname_z, &ucell, 8, &dT_z, params.istep, ispin, nspin, "dT");
88+
#endif
89+
}
90+
}
91+
92+
ModuleBase::timer::end("ModuleIO", "write_dH_t");
93+
return true;
94+
}
95+
96+
} // namespace ModuleIO

0 commit comments

Comments
 (0)