Skip to content

Commit 36bdd4f

Browse files
YuLiu98YuLiu98
authored
Feature: read rho in G space from binary file in qe format (#3303)
* Feature: read in charge density from binary file * Test: add unit tests for read_rhog * Refactor: convert changes about pw_basis * Test: add tests for read_rhog * Test: add charge-density.dat
1 parent d7492b7 commit 36bdd4f

File tree

15 files changed

+523
-73
lines changed

15 files changed

+523
-73
lines changed

‎source/Makefile.Objects‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -413,6 +413,7 @@ OBJS_IO=input.o\
413413
print_info.o\
414414
read_cube.o\
415415
read_rho.o\
416+
read_rhog.o\
416417
restart.o\
417418
binstream.o\
418419
to_wannier90.o\

‎source/module_elecstate/module_charge/charge_init.cpp‎

Lines changed: 90 additions & 73 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@
1414
#include "module_hamilt_pw/hamilt_pwdft/global.h"
1515
#include "module_hamilt_pw/hamilt_pwdft/parallel_grid.h"
1616
#include "module_io/rho_io.h"
17+
#include "module_io/rhog_io.h"
1718
#ifdef USE_PAW
1819
#include "module_cell/module_paw/paw_cell.h"
1920
#endif
@@ -26,99 +27,115 @@ void Charge::init_rho(elecstate::efermi& eferm_iout, const ModuleBase::ComplexMa
2627
bool read_error = false;
2728
if (GlobalV::init_chg == "file")
2829
{
29-
GlobalV::ofs_running << " try to read charge from file : ";
30-
for (int is = 0; is < GlobalV::NSPIN; ++is)
30+
GlobalV::ofs_running << " try to read charge from file : " << std::endl;
31+
32+
// try to read charge from binary file first, which is the same as QE
33+
// liuyu 2023-12-05
34+
std::stringstream binary;
35+
binary << GlobalV::global_readin_dir << "charge-density.dat";
36+
if (ModuleIO::read_rhog(binary.str(), rhopw, rhog))
3137
{
32-
std::stringstream ssc;
33-
ssc << GlobalV::global_readin_dir << "SPIN" << is + 1 << "_CHG.cube";
34-
GlobalV::ofs_running << ssc.str() << std::endl;
35-
double& ef_tmp = eferm_iout.get_ef(is);
36-
if (ModuleIO::read_rho(
37-
#ifdef __MPI
38-
&(GlobalC::Pgrid),
39-
#endif
40-
is,
41-
GlobalV::NSPIN,
42-
ssc.str(),
43-
this->rho[is],
44-
this->rhopw->nx,
45-
this->rhopw->ny,
46-
this->rhopw->nz,
47-
ef_tmp,
48-
&(GlobalC::ucell),
49-
this->prenspin))
50-
{
51-
GlobalV::ofs_running << " Read in the charge density: " << ssc.str() << std::endl;
52-
}
53-
else if (is > 0)
54-
{
55-
if (prenspin == 1)
56-
{
57-
GlobalV::ofs_running << " Didn't read in the charge density but autoset it for spin " << is + 1
58-
<< std::endl;
59-
for (int ir = 0; ir < this->rhopw->nrxx; ir++)
60-
{
61-
this->rho[is][ir] = 0.0;
62-
}
63-
}
64-
//
65-
else if (prenspin == 2)
66-
{ // read up and down , then rearrange them.
67-
if (is == 1)
68-
{
69-
std::cout << "Incomplete charge density file!" << std::endl;
70-
read_error = true;
71-
break;
72-
}
73-
else if (is == 2)
74-
{
75-
GlobalV::ofs_running << " Didn't read in the charge density but would rearrange it later. "
76-
<< std::endl;
77-
}
78-
else if (is == 3)
79-
{
80-
GlobalV::ofs_running << " rearrange charge density " << std::endl;
81-
for (int ir = 0; ir < this->rhopw->nrxx; ir++)
82-
{
83-
this->rho[3][ir] = this->rho[0][ir] - this->rho[1][ir];
84-
this->rho[0][ir] = this->rho[0][ir] + this->rho[1][ir];
85-
this->rho[1][ir] = 0.0;
86-
this->rho[2][ir] = 0.0;
87-
}
88-
}
89-
}
90-
}
91-
else
38+
GlobalV::ofs_running << " Read in the charge density: " << binary.str() << std::endl;
39+
for (int is = 0; is < GlobalV::NSPIN; ++is)
9240
{
93-
read_error = true;
94-
break;
41+
rhopw->recip2real(rhog[is], rho[is]);
9542
}
9643
}
97-
98-
if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
44+
else
9945
{
100-
for (int is = 0; is < GlobalV::NSPIN; is++)
46+
for (int is = 0; is < GlobalV::NSPIN; ++is)
10147
{
10248
std::stringstream ssc;
103-
ssc << GlobalV::global_readin_dir << "SPIN" << is + 1 << "_TAU.cube";
104-
GlobalV::ofs_running << " try to read kinetic energy density from file : " << ssc.str() << std::endl;
105-
// mohan update 2012-02-10, sunliang update 2023-03-09
49+
ssc << GlobalV::global_readin_dir << "SPIN" << is + 1 << "_CHG.cube";
50+
double& ef_tmp = eferm_iout.get_ef(is);
10651
if (ModuleIO::read_rho(
10752
#ifdef __MPI
10853
&(GlobalC::Pgrid),
10954
#endif
11055
is,
11156
GlobalV::NSPIN,
11257
ssc.str(),
113-
this->kin_r[is],
58+
this->rho[is],
11459
this->rhopw->nx,
11560
this->rhopw->ny,
11661
this->rhopw->nz,
117-
eferm_iout.ef,
62+
ef_tmp,
11863
&(GlobalC::ucell),
11964
this->prenspin))
12065
{
121-
GlobalV::ofs_running << " Read in the kinetic energy density: " << ssc.str() << std::endl;
66+
GlobalV::ofs_running << " Read in the charge density: " << ssc.str() << std::endl;
67+
}
68+
else if (is > 0)
69+
{
70+
if (prenspin == 1)
71+
{
72+
GlobalV::ofs_running << " Didn't read in the charge density but autoset it for spin " << is + 1
73+
<< std::endl;
74+
for (int ir = 0; ir < this->rhopw->nrxx; ir++)
75+
{
76+
this->rho[is][ir] = 0.0;
77+
}
78+
}
79+
//
80+
else if (prenspin == 2)
81+
{ // read up and down , then rearrange them.
82+
if (is == 1)
83+
{
84+
std::cout << "Incomplete charge density file!" << std::endl;
85+
read_error = true;
86+
break;
87+
}
88+
else if (is == 2)
89+
{
90+
GlobalV::ofs_running << " Didn't read in the charge density but would rearrange it later. "
91+
<< std::endl;
92+
}
93+
else if (is == 3)
94+
{
95+
GlobalV::ofs_running << " rearrange charge density " << std::endl;
96+
for (int ir = 0; ir < this->rhopw->nrxx; ir++)
97+
{
98+
this->rho[3][ir] = this->rho[0][ir] - this->rho[1][ir];
99+
this->rho[0][ir] = this->rho[0][ir] + this->rho[1][ir];
100+
this->rho[1][ir] = 0.0;
101+
this->rho[2][ir] = 0.0;
102+
}
103+
}
104+
}
105+
}
106+
else
107+
{
108+
read_error = true;
109+
break;
110+
}
111+
}
112+
113+
if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
114+
{
115+
for (int is = 0; is < GlobalV::NSPIN; is++)
116+
{
117+
std::stringstream ssc;
118+
ssc << GlobalV::global_readin_dir << "SPIN" << is + 1 << "_TAU.cube";
119+
GlobalV::ofs_running << " try to read kinetic energy density from file : " << ssc.str()
120+
<< std::endl;
121+
// mohan update 2012-02-10, sunliang update 2023-03-09
122+
if (ModuleIO::read_rho(
123+
#ifdef __MPI
124+
&(GlobalC::Pgrid),
125+
#endif
126+
is,
127+
GlobalV::NSPIN,
128+
ssc.str(),
129+
this->kin_r[is],
130+
this->rhopw->nx,
131+
this->rhopw->ny,
132+
this->rhopw->nz,
133+
eferm_iout.ef,
134+
&(GlobalC::ucell),
135+
this->prenspin))
136+
{
137+
GlobalV::ofs_running << " Read in the kinetic energy density: " << ssc.str() << std::endl;
138+
}
122139
}
123140
}
124141
}

‎source/module_io/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@ list(APPEND objects
1313
print_info.cpp
1414
read_cube.cpp
1515
read_rho.cpp
16+
read_rhog.cpp
1617
restart.cpp
1718
binstream.cpp
1819
write_wfc_pw.cpp

0 commit comments

Comments
 (0)