Skip to content

Commit a2bd6bf

Browse files
dyzhengdyzheng
authored
Merge pull request #95 from jingan-181/develop
Output SR sparse matrix separately
2 parents 7f4224d + d9c3be5 commit a2bd6bf

File tree

6 files changed

+291
-13
lines changed

6 files changed

+291
-13
lines changed

‎source/src_io/write_HS.cpp‎

Lines changed: 132 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1145,6 +1145,138 @@ void HS_Matrix::save_HSR_sparse(
11451145
return;
11461146
}
11471147

1148+
void HS_Matrix::save_SR_sparse(
1149+
const double &sparse_threshold,
1150+
const bool &binary,
1151+
const std::string &SR_filename
1152+
)
1153+
{
1154+
ModuleBase::TITLE("HS_Matrix","save_SR_sparse");
1155+
ModuleBase::timer::tick("HS_Matrix","save_SR_sparse");
1156+
1157+
auto &all_R_coor_ptr = GlobalC::LM.all_R_coor;
1158+
auto &SR_sparse_ptr = GlobalC::LM.SR_sparse;
1159+
auto &SR_soc_sparse_ptr = GlobalC::LM.SR_soc_sparse;
1160+
1161+
int total_R_num = all_R_coor_ptr.size();
1162+
int output_R_number = 0;
1163+
int *S_nonzero_num = nullptr;
1164+
1165+
S_nonzero_num = new int[total_R_num];
1166+
ModuleBase::GlobalFunc::ZEROS(S_nonzero_num, total_R_num);
1167+
1168+
int count = 0;
1169+
for (auto &R_coor : all_R_coor_ptr)
1170+
{
1171+
if (GlobalV::NSPIN != 4)
1172+
{
1173+
auto iter = SR_sparse_ptr.find(R_coor);
1174+
if (iter != SR_sparse_ptr.end())
1175+
{
1176+
for (auto &row_loop : iter->second)
1177+
{
1178+
S_nonzero_num[count] += row_loop.second.size();
1179+
}
1180+
}
1181+
}
1182+
else
1183+
{
1184+
auto iter = SR_soc_sparse_ptr.find(R_coor);
1185+
if (iter != SR_soc_sparse_ptr.end())
1186+
{
1187+
for (auto &row_loop : iter->second)
1188+
{
1189+
S_nonzero_num[count] += row_loop.second.size();
1190+
}
1191+
}
1192+
}
1193+
1194+
count++;
1195+
}
1196+
1197+
Parallel_Reduce::reduce_int_all(S_nonzero_num, total_R_num);
1198+
1199+
for (int index = 0; index < total_R_num; ++index)
1200+
{
1201+
if (S_nonzero_num[index] != 0)
1202+
{
1203+
output_R_number++;
1204+
}
1205+
}
1206+
1207+
std::stringstream sss;
1208+
sss << GlobalV::global_out_dir << SR_filename;
1209+
std::ofstream g2;
1210+
1211+
if(GlobalV::DRANK==0)
1212+
{
1213+
if (binary)
1214+
{
1215+
g2.open(sss.str().c_str(), ios::binary);
1216+
g2.write(reinterpret_cast<char *>(&GlobalV::NLOCAL), sizeof(int));
1217+
g2.write(reinterpret_cast<char *>(&output_R_number), sizeof(int));
1218+
}
1219+
else
1220+
{
1221+
g2.open(sss.str().c_str());
1222+
g2 << "Matrix Dimension of S(R): " << GlobalV::NLOCAL <<std::endl;
1223+
g2 << "Matrix number of S(R): " << output_R_number << std::endl;
1224+
}
1225+
}
1226+
1227+
count = 0;
1228+
for (auto &R_coor : all_R_coor_ptr)
1229+
{
1230+
int dRx = R_coor.x;
1231+
int dRy = R_coor.y;
1232+
int dRz = R_coor.z;
1233+
1234+
if (S_nonzero_num[count] == 0)
1235+
{
1236+
count++;
1237+
continue;
1238+
}
1239+
1240+
if (GlobalV::DRANK == 0)
1241+
{
1242+
if (binary)
1243+
{
1244+
g2.write(reinterpret_cast<char *>(&dRx), sizeof(int));
1245+
g2.write(reinterpret_cast<char *>(&dRy), sizeof(int));
1246+
g2.write(reinterpret_cast<char *>(&dRz), sizeof(int));
1247+
g2.write(reinterpret_cast<char *>(&S_nonzero_num[count]), sizeof(int));
1248+
}
1249+
else
1250+
{
1251+
g2 << dRx << " " << dRy << " " << dRz << " " << S_nonzero_num[count] << std::endl;
1252+
}
1253+
}
1254+
1255+
if (GlobalV::NSPIN != 4)
1256+
{
1257+
output_single_R(g2, SR_sparse_ptr[R_coor], sparse_threshold, binary);
1258+
}
1259+
else
1260+
{
1261+
output_soc_single_R(g2, SR_soc_sparse_ptr[R_coor], sparse_threshold, binary);
1262+
}
1263+
1264+
count++;
1265+
1266+
}
1267+
1268+
if(GlobalV::DRANK==0)
1269+
{
1270+
g2.close();
1271+
}
1272+
1273+
delete[] S_nonzero_num;
1274+
S_nonzero_num = nullptr;
1275+
1276+
ModuleBase::timer::tick("HS_Matrix","save_SR_sparse");
1277+
return;
1278+
}
1279+
11481280
void HS_Matrix::output_single_R(std::ofstream &ofs, const std::map<size_t, std::map<size_t, double>> &XR, const double &sparse_threshold, const bool &binary)
11491281
{
11501282
double *line = nullptr;

‎source/src_io/write_HS.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,11 @@ namespace HS_Matrix
2222
const std::string &HR_filename_up,
2323
const std::string &HR_filename_down
2424
);
25+
void save_SR_sparse(
26+
const double &sparse_threshold,
27+
const bool &binary,
28+
const std::string &SR_filename
29+
);
2530
void output_single_R(std::ofstream &ofs, const std::map<size_t, std::map<size_t, double>> &XR, const double &sparse_threshold, const bool &binary);
2631
void output_soc_single_R(std::ofstream &ofs, const std::map<size_t, std::map<size_t, std::complex<double>>> &XR, const double &sparse_threshold, const bool &binary);
2732

‎source/src_io/write_HS_R.cpp‎

Lines changed: 24 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -3,8 +3,16 @@
33
#include "../src_pw/global.h"
44
#include "write_HS.h"
55

6-
7-
void LOOP_ions::output_HS_R(void)
6+
// if 'binary=true', output binary file.
7+
// The 'sparse_threshold' is the accuracy of the sparse matrix.
8+
// If the absolute value of the matrix element is less than or equal to the 'sparse_threshold', it will be ignored.
9+
void LOOP_ions::output_HS_R(
10+
const std::string &SR_filename,
11+
const std::string &HR_filename_up,
12+
const std::string HR_filename_down,
13+
const bool &binary,
14+
const double &sparse_threshold
15+
)
816
{
917
ModuleBase::TITLE("LOOP_ions","output_HS_R");
1018
ModuleBase::timer::tick("LOOP_ions","output_HS_R");
@@ -17,10 +25,6 @@ void LOOP_ions::output_HS_R(void)
1725
r_matrix.out_r_overlap_R(GlobalV::NSPIN);
1826
}
1927

20-
// Parameters for HR and SR output
21-
double sparse_threshold = 1e-10;
22-
bool binary = false; // output binary file
23-
2428
if(GlobalV::NSPIN==1||GlobalV::NSPIN==4)
2529
{
2630
// GlobalC::UHM.calculate_STN_R();
@@ -88,9 +92,6 @@ void LOOP_ions::output_HS_R(void)
8892
}
8993
}
9094

91-
std::string SR_filename = "data-SR-sparse_SPIN0.csr";
92-
std::string HR_filename_up = "data-HR-sparse_SPIN0.csr";
93-
std::string HR_filename_down = "data-HR-sparse_SPIN1.csr";
9495
HS_Matrix::save_HSR_sparse(sparse_threshold, binary, SR_filename, HR_filename_up, HR_filename_down);
9596
GlobalC::UHM.destroy_all_HSR_sparse();
9697

@@ -102,3 +103,17 @@ void LOOP_ions::output_HS_R(void)
102103
ModuleBase::timer::tick("LOOP_ions","output_HS_R");
103104
return;
104105
}
106+
107+
108+
void LOOP_ions::output_SR(const std::string &SR_filename, const bool &binary, const double &sparse_threshold)
109+
{
110+
ModuleBase::TITLE("LOOP_ions","output_SR");
111+
ModuleBase::timer::tick("LOOP_ions","output_SR");
112+
113+
GlobalC::UHM.calculate_SR_sparse(sparse_threshold);
114+
HS_Matrix::save_SR_sparse(sparse_threshold, binary, SR_filename);
115+
GlobalC::UHM.destroy_all_HSR_sparse();
116+
117+
ModuleBase::timer::tick("LOOP_ions","output_SR");
118+
return;
119+
}

‎source/src_lcao/LCAO_hamilt.cpp‎

Lines changed: 120 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -620,10 +620,13 @@ void LCAO_Hamilt::calculate_STN_R_sparse(const int &current_spin, const double &
620620

621621
if(GlobalV::NSPIN!=4)
622622
{
623-
temp_value_double = GlobalC::LM.SlocR[index];
624-
if (std::abs(temp_value_double) > sparse_threshold)
623+
if (current_spin == 0)
625624
{
626-
GlobalC::LM.SR_sparse[dR][iw1_all][iw2_all] = temp_value_double;
625+
temp_value_double = GlobalC::LM.SlocR[index];
626+
if (std::abs(temp_value_double) > sparse_threshold)
627+
{
628+
GlobalC::LM.SR_sparse[dR][iw1_all][iw2_all] = temp_value_double;
629+
}
627630
}
628631

629632
temp_value_double = GlobalC::LM.Hloc_fixedR[index];
@@ -659,6 +662,113 @@ void LCAO_Hamilt::calculate_STN_R_sparse(const int &current_spin, const double &
659662
}
660663

661664

665+
void LCAO_Hamilt::calculate_STN_R_sparse_for_S(const double &sparse_threshold)
666+
{
667+
ModuleBase::TITLE("LCAO_Hamilt","calculate_STN_R_sparse_for_S");
668+
669+
int index = 0;
670+
ModuleBase::Vector3<double> dtau, tau1, tau2;
671+
ModuleBase::Vector3<double> dtau1, dtau2, tau0;
672+
673+
double temp_value_double;
674+
std::complex<double> temp_value_complex;
675+
676+
for(int T1 = 0; T1 < GlobalC::ucell.ntype; ++T1)
677+
{
678+
Atom* atom1 = &GlobalC::ucell.atoms[T1];
679+
for(int I1 = 0; I1 < atom1->na; ++I1)
680+
{
681+
tau1 = atom1->tau[I1];
682+
GlobalC::GridD.Find_atom(GlobalC::ucell, tau1, T1, I1);
683+
Atom* atom1 = &GlobalC::ucell.atoms[T1];
684+
const int start = GlobalC::ucell.itiaiw2iwt(T1,I1,0);
685+
686+
for(int ad = 0; ad < GlobalC::GridD.getAdjacentNum()+1; ++ad)
687+
{
688+
const int T2 = GlobalC::GridD.getType(ad);
689+
const int I2 = GlobalC::GridD.getNatom(ad);
690+
Atom* atom2 = &GlobalC::ucell.atoms[T2];
691+
692+
tau2 = GlobalC::GridD.getAdjacentTau(ad);
693+
dtau = tau2 - tau1;
694+
double distance = dtau.norm() * GlobalC::ucell.lat0;
695+
double rcut = GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ORB.Phi[T2].getRcut();
696+
697+
bool adj = false;
698+
699+
if(distance < rcut) adj = true;
700+
else if(distance >= rcut)
701+
{
702+
for(int ad0 = 0; ad0 < GlobalC::GridD.getAdjacentNum()+1; ++ad0)
703+
{
704+
const int T0 = GlobalC::GridD.getType(ad0);
705+
706+
tau0 = GlobalC::GridD.getAdjacentTau(ad0);
707+
dtau1 = tau0 - tau1;
708+
dtau2 = tau0 - tau2;
709+
710+
double distance1 = dtau1.norm() * GlobalC::ucell.lat0;
711+
double distance2 = dtau2.norm() * GlobalC::ucell.lat0;
712+
713+
double rcut1 = GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max();
714+
double rcut2 = GlobalC::ORB.Phi[T2].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max();
715+
716+
if( distance1 < rcut1 && distance2 < rcut2 )
717+
{
718+
adj = true;
719+
break;
720+
}
721+
}
722+
}
723+
724+
if(adj)
725+
{
726+
const int start2 = GlobalC::ucell.itiaiw2iwt(T2,I2,0);
727+
728+
Abfs::Vector3_Order<int> dR(GlobalC::GridD.getBox(ad).x, GlobalC::GridD.getBox(ad).y, GlobalC::GridD.getBox(ad).z);
729+
730+
for(int ii=0; ii<atom1->nw*GlobalV::NPOL; ii++)
731+
{
732+
const int iw1_all = start + ii;
733+
const int mu = GlobalC::ParaO.trace_loc_row[iw1_all];
734+
735+
if(mu<0)continue;
736+
737+
for(int jj=0; jj<atom2->nw*GlobalV::NPOL; jj++)
738+
{
739+
int iw2_all = start2 + jj;
740+
const int nu = GlobalC::ParaO.trace_loc_col[iw2_all];
741+
742+
if(nu<0)continue;
743+
744+
if(GlobalV::NSPIN!=4)
745+
{
746+
temp_value_double = GlobalC::LM.SlocR[index];
747+
if (std::abs(temp_value_double) > sparse_threshold)
748+
{
749+
GlobalC::LM.SR_sparse[dR][iw1_all][iw2_all] = temp_value_double;
750+
}
751+
}
752+
else
753+
{
754+
temp_value_complex = GlobalC::LM.SlocR_soc[index];
755+
if(std::abs(temp_value_complex) > sparse_threshold)
756+
{
757+
GlobalC::LM.SR_soc_sparse[dR][iw1_all][iw2_all] = temp_value_complex;
758+
}
759+
}
760+
761+
++index;
762+
}
763+
}
764+
}
765+
}
766+
}
767+
}
768+
769+
return;
770+
}
771+
662772
void LCAO_Hamilt::calculate_HSR_sparse(const int &current_spin, const double &sparse_threshold)
663773
{
664774
ModuleBase::TITLE("LCAO_Hamilt","calculate_HSR_sparse");
@@ -685,6 +795,13 @@ void LCAO_Hamilt::calculate_HSR_sparse(const int &current_spin, const double &sp
685795

686796
}
687797

798+
void LCAO_Hamilt::calculate_SR_sparse(const double &sparse_threshold)
799+
{
800+
ModuleBase::TITLE("LCAO_Hamilt","calculate_SR_sparse");
801+
set_R_range_sparse();
802+
calculate_STN_R_sparse_for_S(sparse_threshold);
803+
}
804+
688805
void LCAO_Hamilt::calculat_HR_dftu_sparse(const int &current_spin, const double &sparse_threshold)
689806
{
690807
ModuleBase::TITLE("LCAO_Hamilt","calculat_HR_dftu_sparse");

‎source/src_lcao/LCAO_hamilt.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,9 +26,11 @@ class LCAO_Hamilt
2626
// jingan add 2021-6-4
2727
void set_R_range_sparse();
2828
void calculate_STN_R_sparse(const int &current_spin, const double &sparse_threshold);
29+
void calculate_STN_R_sparse_for_S(const double &sparse_threshold);
2930
void calculat_HR_dftu_sparse(const int &current_spin, const double &sparse_threshold);
3031
void calculat_HR_dftu_soc_sparse(const int &current_spin, const double &sparse_threshold);
3132
void calculate_HSR_sparse(const int &current_spin, const double &sparse_threshold);
33+
void calculate_SR_sparse(const double &sparse_threshold);
3234
void clear_zero_elements(const int &current_spin, const double &sparse_threshold);
3335
void destroy_all_HSR_sparse(void);
3436

‎source/src_lcao/LOOP_ions.h‎

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,14 @@ class LOOP_ions
1717
LOOP_elec LOE;
1818

1919
void opt_ions(void);
20-
void output_HS_R(void); //LiuXh add 2019-07-15
20+
void output_HS_R(
21+
const std::string &SR_filename="data-SR-sparse_SPIN0.csr",
22+
const std::string &HR_filename_up="data-HR-sparse_SPIN0.csr",
23+
const std::string HR_filename_down="data-HR-sparse_SPIN1.csr",
24+
const bool &binary=false,
25+
const double &sparse_threshold=1e-10
26+
); //LiuXh add 2019-07-15, modify in 2021-12-3
27+
void output_SR(const std::string &SR_filename, const bool &binary=false, const double &sparse_threshold=1e-10);
2128

2229
private:
2330

0 commit comments

Comments
 (0)