Skip to content

Commit 437cc26

Browse files
wenfei-liwenfei-li
authored
Fix: out_dm1 (#1577)
* Fix: out_dm1
1 parent ef7b09c commit 437cc26

File tree

3 files changed

+73
-82
lines changed

3 files changed

+73
-82
lines changed

‎source/module_esolver/esolver_ks_lcao.cpp‎

Lines changed: 22 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -561,6 +561,18 @@ void ESolver_KS_LCAO::afterscf(const int istep)
561561
// 0 means don't need to consider iter,
562562
//--------------------------------------
563563

564+
double** dm2d;
565+
if(this->LOC.out_dm1 == 1)
566+
{
567+
dm2d = new double*[GlobalV::NSPIN];
568+
for (int is = 0; is < GlobalV::NSPIN; is++)
569+
{
570+
dm2d[is] = new double[this->LOC.ParaV->nnr];
571+
ModuleBase::GlobalFunc::ZEROS(dm2d[is], this->LOC.ParaV->nnr);
572+
}
573+
this->LOC.cal_dm_R(this->LOC.dm_k,this->RA,dm2d);
574+
}
575+
564576
for (int is = 0; is < GlobalV::NSPIN; is++)
565577
{
566578
const int precision = 3;
@@ -584,7 +596,7 @@ void ESolver_KS_LCAO::afterscf(const int istep)
584596
this->LOC.write_dm(is, 0, ssd.str(), precision);
585597
if(this->LOC.out_dm1 == 1)
586598
{
587-
this->LOC.write_dm1(is, istep);
599+
this->LOC.write_dm1(is, istep, dm2d);
588600
}
589601
/* Broken, please fix it
590602
if (GlobalV::out_pot == 1) // LiuXh add 20200701
@@ -596,6 +608,15 @@ void ESolver_KS_LCAO::afterscf(const int istep)
596608
*/
597609
}
598610

611+
if(this->LOC.out_dm1 == 1)
612+
{
613+
for (int is = 0; is < GlobalV::NSPIN; is++)
614+
{
615+
delete[] dm2d[is];
616+
}
617+
delete[] dm2d;
618+
}
619+
599620
#ifdef __EXX
600621
if (GlobalC::exx_info.info_global.cal_exx) // Peize Lin add if 2022.11.14
601622
{

‎source/src_io/write_dm.cpp‎

Lines changed: 49 additions & 79 deletions
Original file line numberDiff line numberDiff line change
@@ -124,103 +124,73 @@ void Local_Orbital_Charge::write_dm(
124124

125125
//ofs << "\n " << GlobalV::GAMMA_ONLY_LOCAL << " (GAMMA ONLY LOCAL)" << std::endl;
126126
#ifndef __MPI
127-
if(GlobalV::GAMMA_ONLY_LOCAL)
128-
{
129-
for(int i=0; i<GlobalV::NLOCAL; ++i)
130-
{
131-
for(int j=0; j<GlobalV::NLOCAL; ++j)
132-
{
133-
if(j%8==0) ofs << "\n";
134-
ofs << " " << this->DM[is][i][j];
135-
}
136-
}
137-
}
138-
else
127+
128+
for(int i=0; i<GlobalV::NLOCAL; ++i)
139129
{
140-
ModuleBase::WARNING_QUIT("write_dm","not ready yet");
141-
ofs << " " << GlobalC::GridT.nnrg << " (nnrg)" << std::endl;
142-
for(int i=0; i<GlobalC::GridT.nnrg; ++i)
130+
for(int j=0; j<GlobalV::NLOCAL; ++j)
143131
{
144-
if(i%8==0) ofs << "\n";
145-
ofs << " " << this->DM_R[is][i];
132+
if(j%8==0) ofs << "\n";
133+
ofs << " " << this->DM[is][i][j];
146134
}
147135
}
136+
148137
#else
149-
if(GlobalV::GAMMA_ONLY_LOCAL)
138+
//xiaohui modify 2014-06-18
139+
140+
double* tmp = new double[GlobalV::NLOCAL];
141+
int* count = new int[GlobalV::NLOCAL];
142+
for (int i=0; i<GlobalV::NLOCAL; ++i)
150143
{
151-
//xiaohui modify 2014-06-18
152-
153-
double* tmp = new double[GlobalV::NLOCAL];
154-
int* count = new int[GlobalV::NLOCAL];
155-
for (int i=0; i<GlobalV::NLOCAL; ++i)
144+
// when reduce, there may be 'redundance', we need to count them.
145+
ModuleBase::GlobalFunc::ZEROS(count, GlobalV::NLOCAL);
146+
const int mu = GlobalC::GridT.trace_lo[i];
147+
if (mu >= 0)
156148
{
157-
// when reduce, there may be 'redundance', we need to count them.
158-
ModuleBase::GlobalFunc::ZEROS(count, GlobalV::NLOCAL);
159-
const int mu = GlobalC::GridT.trace_lo[i];
160-
if (mu >= 0)
149+
for (int j=0; j<GlobalV::NLOCAL; ++j)
161150
{
162-
for (int j=0; j<GlobalV::NLOCAL; ++j)
151+
const int nu = GlobalC::GridT.trace_lo[j];
152+
if (nu >= 0)
163153
{
164-
const int nu = GlobalC::GridT.trace_lo[j];
165-
if (nu >= 0)
166-
{
167-
count[j]=1;
168-
}
154+
count[j]=1;
169155
}
170156
}
171-
Parallel_Reduce::reduce_int_all( count, GlobalV::NLOCAL );
157+
}
158+
Parallel_Reduce::reduce_int_all( count, GlobalV::NLOCAL );
172159

173-
// reduce the density matrix for 'i' line.
174-
ModuleBase::GlobalFunc::ZEROS(tmp, GlobalV::NLOCAL);
175-
if (mu >= 0)
160+
// reduce the density matrix for 'i' line.
161+
ModuleBase::GlobalFunc::ZEROS(tmp, GlobalV::NLOCAL);
162+
if (mu >= 0)
163+
{
164+
for (int j=0; j<GlobalV::NLOCAL; j++)
176165
{
177-
for (int j=0; j<GlobalV::NLOCAL; j++)
166+
const int nu = GlobalC::GridT.trace_lo[j];
167+
if (nu >=0)
178168
{
179-
const int nu = GlobalC::GridT.trace_lo[j];
180-
if (nu >=0)
181-
{
182-
tmp[j] = DM[is][mu][nu];
183-
//GlobalV::ofs_running << " dmi=" << i << " j=" << j << " " << DM[is][mu][nu] << std::endl;
184-
}
169+
tmp[j] = DM[is][mu][nu];
170+
//GlobalV::ofs_running << " dmi=" << i << " j=" << j << " " << DM[is][mu][nu] << std::endl;
185171
}
186172
}
187-
Parallel_Reduce::reduce_double_all( tmp, GlobalV::NLOCAL );
173+
}
174+
Parallel_Reduce::reduce_double_all( tmp, GlobalV::NLOCAL );
188175

189-
if(GlobalV::MY_RANK==0)
176+
if(GlobalV::MY_RANK==0)
177+
{
178+
for (int j=0; j<GlobalV::NLOCAL; j++)
190179
{
191-
for (int j=0; j<GlobalV::NLOCAL; j++)
180+
if(j%8==0) ofs << "\n";
181+
if(count[j]>0)
192182
{
193-
if(j%8==0) ofs << "\n";
194-
if(count[j]>0)
195-
{
196-
ofs << " " << tmp[j]/(double)count[j];
197-
}
198-
else
199-
{
200-
ofs << " 0";
201-
}
183+
ofs << " " << tmp[j]/(double)count[j];
184+
}
185+
else
186+
{
187+
ofs << " 0";
202188
}
203189
}
204190
}
205-
delete[] tmp;
206-
delete[] count;
207-
208-
//xiaohui add 2014-06-18
209-
//for(int i=0; i<GlobalV::NLOCAL; ++i)
210-
//{
211-
// for(int j=0; j<GlobalV::NLOCAL; ++j)
212-
// {
213-
// if(j%8==0) ofs << "\n";
214-
// ofs << " " << this->DM[is][i][j];
215-
// }
216-
//}
217-
218-
}
219-
else
220-
{
221-
ofs << " " << GlobalC::GridT.nnrg << " (nnrg)" << std::endl;
222-
ModuleBase::WARNING_QUIT("local_orbital_charge","not ready to output DM_R");
223191
}
192+
delete[] tmp;
193+
delete[] count;
224194
#endif
225195
if(GlobalV::MY_RANK==0)
226196
{
@@ -233,12 +203,12 @@ void Local_Orbital_Charge::write_dm(
233203
return;
234204
}
235205

236-
void Local_Orbital_Charge::write_dm1(const int &is, const int &istep)
206+
void Local_Orbital_Charge::write_dm1(const int &is, const int &istep, double** dm2d)
237207
{
238208
ModuleBase::timer::tick("Local_Orbital_Charge","write_dm");
239209
ModuleBase::TITLE("Local_Orbital_Charge","write_dm");
240210

241-
get_dm_sparse(is);
211+
get_dm_sparse(is, dm2d);
242212
write_dm_sparse(is, istep);
243213
destroy_dm_sparse();
244214

@@ -328,7 +298,7 @@ inline void output_single_R(std::ofstream &ofs, const std::map<size_t, std::map<
328298

329299
}
330300

331-
void Local_Orbital_Charge::get_dm_sparse(const int &is)
301+
void Local_Orbital_Charge::get_dm_sparse(const int &is, double** dm2d)
332302
{
333303
ModuleBase::timer::tick("Local_Orbital_Charge","get_dm_sparse");
334304
ModuleBase::TITLE("Local_Orbital_Charge","get_dm_sparse");
@@ -409,7 +379,7 @@ void Local_Orbital_Charge::get_dm_sparse(const int &is)
409379

410380
if(nu<0)continue;
411381

412-
temp_value_double = this->DM_R[is][index];
382+
temp_value_double = dm2d[is][index];
413383
if (std::abs(temp_value_double) > sparse_threshold)
414384
{
415385
this->DMR_sparse[dR][iw1_all][iw2_all] = temp_value_double;
@@ -500,7 +470,7 @@ void Local_Orbital_Charge::write_dm_sparse(const int &is, const int &istep)
500470

501471
if(GlobalV::DRANK==0)
502472
{
503-
g1.open(ssdm.str().c_str(), ios::app);
473+
g1.open(ssdm.str().c_str());
504474
g1 << "STEP: " << istep << std::endl;
505475
g1 << "Matrix Dimension of DM(R): " << GlobalV::NLOCAL <<std::endl;
506476
g1 << "Matrix number of DM(R): " << output_R_number << std::endl;

‎source/src_lcao/local_orbital_charge.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,7 @@ class Local_Orbital_Charge
5050
static int out_dm1;
5151

5252
void write_dm(const int &is, const int &iter, const std::string &fn, const int &precision);
53-
void write_dm1(const int &is, const int &istep);
53+
void write_dm1(const int &is, const int &istep, double** dm2d);
5454

5555
void read_dm(const int &is, const std::string &fn);
5656

@@ -120,7 +120,7 @@ class Local_Orbital_Charge
120120

121121
std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, double>>> DMR_sparse;
122122

123-
void get_dm_sparse(const int &is);
123+
void get_dm_sparse(const int &is, double** dm2d);
124124
void write_dm_sparse(const int &is, const int &istep);
125125
void destroy_dm_sparse();
126126

0 commit comments

Comments
 (0)