Skip to content

Commit 27fcd05

Browse files
author
wenfei-li
committed
gint_k : removed obsolete code evaluate_pDMp
1 parent fcd593c commit 27fcd05

File tree

2 files changed

+4
-236
lines changed

2 files changed

+4
-236
lines changed

‎source/src_lcao/gint_k.h‎

Lines changed: 0 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -124,16 +124,6 @@ class Gint_k : public Gint_k_init
124124
double** distance,
125125
const double &delta_r);
126126

127-
//------------------------------------------------------
128-
// in gint_k_rho.cpp
129-
//------------------------------------------------------
130-
// evaluate the <phi0 | Density Matrix | phiR> to get the charge density.
131-
void evaluate_pDMp(
132-
const int &grid_index,
133-
const int &size,
134-
bool** cal_flag,
135-
double*** psir_ylm,
136-
int* vindex);
137127

138128
//------------------------------------------------------
139129
// in gint_k_fvl.cpp

‎source/src_lcao/gint_k_rho.cpp‎

Lines changed: 4 additions & 226 deletions
Original file line numberDiff line numberDiff line change
@@ -217,7 +217,7 @@ inline void cal_band_rho(
217217
std::cout << " size1=" << ia1 << " size2=" << ia1 << std::endl;
218218
std::cout << " iat1=" << iat1 << " iat2=" << iat1 << std::endl;
219219
std::cout << " dR=" << dRx << " " << dRy << " " << dRz << std::endl;
220-
ModuleBase::WARNING_QUIT("gint_k","evaluate_pDMp wrong");
220+
ModuleBase::WARNING_QUIT("gint_k","cal_band_rho wrong");
221221
}
222222
//const int offset=AllOffset[ia1][ia2];
223223
assert(offset < GlobalC::GridT.nad[iat1]);
@@ -256,7 +256,7 @@ inline void cal_band_rho(
256256
std::cout << " size1=" << ia1 << " size2=" << ia1 << std::endl;
257257
std::cout << " iat1=" << iat1 << " iat2=" << iat1 << std::endl;
258258
std::cout << " dR=" << dRx << " " << dRy << " " << dRz << std::endl;
259-
ModuleBase::WARNING_QUIT("gint_k","evaluate_pDMp wrong");
259+
ModuleBase::WARNING_QUIT("gint_k","cal_band_rho wrong");
260260
}
261261
//const int offset=AllOffset[ia1][ia2];
262262
assert(offset < GlobalC::GridT.nad[iat1]);
@@ -320,7 +320,7 @@ inline void cal_band_rho(
320320
std::cout << " dR=" << dRx << " " << dRy << " " << dRz << std::endl;
321321
std::cout << " R1=" << R1x << " " << R1y << " " << R1z << std::endl;
322322
std::cout << " R2=" << R2x << " " << R2y << " " << R2z << std::endl;
323-
ModuleBase::WARNING_QUIT("gint_k","evaluate_pDMp wrong");
323+
ModuleBase::WARNING_QUIT("gint_k","cal_band_rho wrong");
324324
}
325325
assert(offset < GlobalC::GridT.nad[iat1]);
326326

@@ -369,7 +369,7 @@ inline void cal_band_rho(
369369
std::cout << " dR=" << dRx << " " << dRy << " " << dRz << std::endl;
370370
std::cout << " R1=" << R1x << " " << R1y << " " << R1z << std::endl;
371371
std::cout << " R2=" << R2x << " " << R2y << " " << R2z << std::endl;
372-
ModuleBase::WARNING_QUIT("gint_k","evaluate_pDMp wrong");
372+
ModuleBase::WARNING_QUIT("gint_k","cal_band_rho wrong");
373373
}
374374
assert(offset < GlobalC::GridT.nad[iat1]);
375375

@@ -549,226 +549,4 @@ void Gint_k::cal_rho_k(double** DM_R_in)
549549
return;
550550
}
551551

552-
void Gint_k::evaluate_pDMp(
553-
const int &grid_index,
554-
const int &size,
555-
bool** cal_flag,
556-
double*** psir_ylm,
557-
int* vindex)
558-
{
559-
//-----------------------------------------------------
560-
// in order to calculate <i,alpha,R1 | DM_R | j,beta,R2>
561-
//-----------------------------------------------------
562-
double **tchg = new double*[GlobalV::NSPIN];
563-
for(int is=0; is<GlobalV::NSPIN; is++)
564-
{
565-
tchg[is] = new double[GlobalC::pw.bxyz];
566-
ModuleBase::GlobalFunc::ZEROS(tchg[is], GlobalC::pw.bxyz);
567-
}
568-
569-
bool *all_out_of_range = new bool[size];
570-
for(int ia=0; ia<size; ia++)
571-
{
572-
all_out_of_range[ia] = true;
573-
for(int ib=0; ib<GlobalC::pw.bxyz; ib++)
574-
{
575-
if(cal_flag[ib][ia])
576-
{
577-
all_out_of_range[ia] = false;
578-
}
579-
}
580-
}
581-
582-
double *psi1, *psi2;
583-
double *iw1p, *iw2p;
584-
double *end1, *end2;
585-
double *dm, *tchgs;
586-
int ixxx, ixxx2;
587-
int iw1_lo, iw2_lo;
588-
double psi1_2;
589-
590-
// get (i,alpha,R1)
591-
// size: how many atoms in this big cell.
592-
for (int ia1=0; ia1<size; ia1++)
593-
{
594-
if(all_out_of_range[ia1]) continue;
595-
const int mcell_index1 = GlobalC::GridT.bcell_start[grid_index] + ia1;
596-
const int iat = GlobalC::GridT.which_atom[mcell_index1];
597-
const int T1 = GlobalC::ucell.iat2it[iat];
598-
const int I1 = GlobalC::ucell.iat2ia[iat];
599-
const int start1 = GlobalC::ucell.itiaiw2iwt(T1, I1, 0);
600-
Atom *atom1 = &GlobalC::ucell.atoms[T1];
601-
const int nw1 = atom1->nw;
602-
603-
//~~~~~~~~~~~~~~~~
604-
// get cell R1.
605-
//~~~~~~~~~~~~~~~~
606-
const int id1 = GlobalC::GridT.which_unitcell[mcell_index1];
607-
const int R1x = GlobalC::GridT.ucell_index2x[id1];
608-
const int R1y = GlobalC::GridT.ucell_index2y[id1];
609-
const int R1z = GlobalC::GridT.ucell_index2z[id1];
610-
const int DM_start = GlobalC::GridT.nlocstartg[iat];
611-
612-
// get (j,beta,R2)
613-
for (int ia2=0; ia2<size; ia2++)
614-
{
615-
if(all_out_of_range[ia2]) continue;
616-
617-
618-
bool same_flag = false;
619-
for(int ib=0; ib<GlobalC::pw.bxyz; ib++)
620-
{
621-
if(cal_flag[ib][ia1] && cal_flag[ib][ia2])
622-
{
623-
same_flag = true;
624-
break;
625-
}
626-
}
627-
628-
// only go on if the two atoms are in the same meshcell!
629-
if(!same_flag)continue;
630-
631-
const int mcell_index2 = GlobalC::GridT.bcell_start[grid_index] + ia2;
632-
const int iat2 = GlobalC::GridT.which_atom[mcell_index2];
633-
const int T2 = GlobalC::ucell.iat2it[iat2];
634-
635-
if (T2 >= T1)
636-
{
637-
Atom *atom2 = &GlobalC::ucell.atoms[T2];
638-
const int nw2 = atom2->nw;
639-
const int I2 = GlobalC::ucell.iat2ia[iat2];
640-
const int start2 = GlobalC::ucell.itiaiw2iwt(T2, I2, 0);
641-
642-
//~~~~~~~~~~~~~~~~
643-
// get cell R2.
644-
//~~~~~~~~~~~~~~~~
645-
const int id2 = GlobalC::GridT.which_unitcell[mcell_index2];
646-
const int R2x = GlobalC::GridT.ucell_index2x[id2];
647-
const int R2y = GlobalC::GridT.ucell_index2y[id2];
648-
const int R2z = GlobalC::GridT.ucell_index2z[id2];
649-
650-
const int dRx = R1x - R2x;
651-
const int dRy = R1y - R2y;
652-
const int dRz = R1z - R2z;
653-
654-
// get the index from dRx, dRy, dRz.
655-
// in fact I can calculate this once, and update only when meshcell is changed.
656-
const int index = GlobalC::GridT.cal_RindexAtom(dRx, dRy, dRz, iat2);
657-
658-
int offset = -1;
659-
int* find_start = GlobalC::GridT.find_R2[iat];
660-
int* find_end = GlobalC::GridT.find_R2[iat] + GlobalC::GridT.nad[iat];
661-
662-
for(int* find=find_start; find < find_end; find++)
663-
{
664-
//--------------------------------------------------------------
665-
// start positions of adjacent atom of 'iat'
666-
//--------------------------------------------------------------
667-
if( find[0] == index )
668-
{
669-
offset = find - find_start;
670-
break;
671-
}
672-
}
673-
674-
if(offset == -1)
675-
{
676-
std::cout << "== charge ==========================" << std::endl;
677-
std::cout << " grid_index = " << grid_index << std::endl;
678-
std::cout << " index = " << index << std::endl;
679-
std::cout << " size1=" << ia1 << " size2=" << ia2 << std::endl;
680-
std::cout << " iat=" << iat << " iat2=" << iat2 << std::endl;
681-
std::cout << " dR=" << dRx << " " << dRy << " " << dRz << std::endl;
682-
std::cout << " R1=" << R1x << " " << R1y << " " << R1z << std::endl;
683-
std::cout << " R2=" << R2x << " " << R2y << " " << R2z << std::endl;
684-
ModuleBase::WARNING_QUIT("gint_k","evaluate_pDMp wrong");
685-
}
686-
assert(offset < GlobalC::GridT.nad[iat]);
687-
688-
689-
// key variable:
690-
ixxx = DM_start + GlobalC::GridT.find_R2st[iat][offset];
691-
692-
for(int is=0; is<GlobalV::NSPIN; is++)
693-
{
694-
dm = this->DM_R[is];
695-
tchgs = tchg[is];
696-
for(int ib=0; ib<GlobalC::pw.bxyz; ib++)
697-
{
698-
if(cal_flag[ib][ia1] && cal_flag[ib][ia2])
699-
{
700-
psi1 = psir_ylm[ib][ia1];
701-
psi2 = psir_ylm[ib][ia2];
702-
end1 = psi1 + nw1;
703-
end2 = psi2 + nw2;
704-
705-
iw1_lo = GlobalC::GridT.trace_lo[start1]/GlobalV::NPOL;
706-
ixxx2 = ixxx;
707-
//------------------------------------
708-
// circle for wave functions of atom 1.
709-
//------------------------------------
710-
for (iw1p=psi1; iw1p<end1; ++iw1p)
711-
{
712-
iw2_lo = GlobalC::GridT.trace_lo[start2]/GlobalV::NPOL;
713-
// 2.0 counts for the undiagonalized part
714-
psi1_2 = 2.0 * iw1p[0];
715-
//------------------------------------
716-
// circle for wave functions of atom 2.
717-
// dmt: temperary density matrix
718-
//------------------------------------
719-
double* dmt = &dm[ixxx2];
720-
for (iw2p=psi2; iw2p<end2; ++iw2p)
721-
{
722-
if ( iw1_lo > iw2_lo)
723-
{
724-
++iw2_lo;
725-
++dmt;
726-
continue;
727-
}
728-
else if ( iw1_lo < iw2_lo)
729-
{
730-
tchgs[ib] +=
731-
psi1_2
732-
* iw2p[0]
733-
* dmt[0];
734-
}
735-
else// means iw1_lo == iw2_lo
736-
{
737-
tchgs[ib] +=
738-
iw1p[0]
739-
* iw2p[0]
740-
* dmt[0];
741-
}
742-
++iw2_lo;
743-
++dmt;
744-
}// iw2
745-
++iw1_lo;
746-
ixxx2 += nw2;
747-
}// iw
748-
}//end cal_flag
749-
}// end ib
750-
}// end is
751-
}// T
752-
}// ia2
753-
}// ia1
754-
755-
for(int is=0; is<GlobalV::NSPIN; is++)
756-
{
757-
for(int ib=0; ib<GlobalC::pw.bxyz; ib++)
758-
{
759-
GlobalC::CHR.rho[is][vindex[ib]] += tchg[is][ib];
760-
}
761-
}
762-
763-
delete[] all_out_of_range;
764-
765-
for(int is=0; is<GlobalV::NSPIN; is++)
766-
{
767-
delete[] tchg[is];
768-
}
769-
delete[] tchg;
770-
771-
return;
772-
}
773-
774552

0 commit comments

Comments
 (0)