diff --git a/ABACUS.develop/source/src_lcao/ORB_control.cpp b/ABACUS.develop/source/src_lcao/ORB_control.cpp index 9be615436f..de445a7d32 100644 --- a/ABACUS.develop/source/src_lcao/ORB_control.cpp +++ b/ABACUS.develop/source/src_lcao/ORB_control.cpp @@ -57,5 +57,9 @@ void ORB_control::clear_after_ions() TITLE("ORB_control","clear_after_ions"); UOT.MOT.Destroy_Table(); UOT.tbeta.Destroy_Table_Beta(); + //caoyu add 2021-03-18 + if (INPUT.out_descriptor && BASIS_TYPE == "lcao") { + UOT.talpha.Destroy_Table_Alpha(); + } return; } diff --git a/ABACUS.develop/source/src_lcao/ORB_gen_tables.cpp b/ABACUS.develop/source/src_lcao/ORB_gen_tables.cpp index cda7c3d225..043948466f 100644 --- a/ABACUS.develop/source/src_lcao/ORB_gen_tables.cpp +++ b/ABACUS.develop/source/src_lcao/ORB_gen_tables.cpp @@ -36,6 +36,17 @@ void ORB_gen_tables::gen_tables( const int &job0 ) ORB.get_dR(),// delta R, for making radial table ORB.get_dk() ); // delta k, for integration in k space + //caoyu add 2021-03-18 + if (INPUT.out_descriptor && BASIS_TYPE == "lcao") { + talpha.allocate( + ORB.get_ntype(),// number of atom types + ORB.get_lmax(),// max L used to calculate overlap + ORB.get_kmesh(), // kpoints, for integration in k space + ORB.get_Rmax(),// max value of radial table + ORB.get_dR(),// delta R, for making radial table + ORB.get_dk()); // delta k, for integration in k space + } + // OV: overlap MOT.init_OV_Tpair(); MOT.init_OV_Opair(); @@ -44,6 +55,12 @@ void ORB_gen_tables::gen_tables( const int &job0 ) tbeta.init_NL_Tpair(); tbeta.init_NL_Opair(); // add 2009-5-8 + //caoyu add 2021-03-18 + // DS: Descriptor + if (INPUT.out_descriptor && BASIS_TYPE == "lcao") { + talpha.init_DS_Opair(); + talpha.init_DS_2Lplus1(); + } //========================================= // (2) init Ylm Coef @@ -59,6 +76,12 @@ void ORB_gen_tables::gen_tables( const int &job0 ) MOT.init_Table(job0); tbeta.init_Table_Beta( MOT.pSB );// add 2009-5-8 + //caoyu add 2021-03-18 + if (INPUT.out_descriptor && BASIS_TYPE == "lcao") { + talpha.init_Table_Alpha(MOT.pSB); + talpha.print_Table_DSR(); + } + //========================================= // (3) make Gaunt coefficients table //========================================= @@ -792,4 +815,193 @@ double ORB_gen_tables::get_distance( const Vector3 &R1, const Vector3lat0; } +//caoyu add 2021-03-17 +void ORB_gen_tables::snap_psialpha( + double olm[], + const int& job, + const Vector3& R1, + const int& T1, + const int& L1, + const int& m1, + const int& N1, + const Vector3& R2, + const int& T2, + const int& L2, + const int& m2, + const int& N2, + complex* olm1, + const int is) const +{ + + if (job != 0 && job != 1) + { + WARNING_QUIT("ORB_gen_tables::snap_psialpha", "job must be equal to 0 or 1!"); + } + + Numerical_Orbital::set_position(R1, R2); + assert(this->lat0 > 0.0); + + // (1) get distance between R1 and R2 (a.u.) + // judge if there exist overlap + double distance = Numerical_Orbital::get_distance() * this->lat0; + + const double Rcut1 = ORB.Phi[T1].getRcut(); + const double Rcut2 = ORB.Alpha[0].getRcut(); + + if (job == 0) ZEROS(olm, 1); + else if (job == 1) ZEROS(olm, 3); + + if (distance > (Rcut1 + Rcut2)) return; + + //if distance == 0 + //\int psi(r) psi(r-R) dr independent of R if R == 0 + //distance += tiny1 avoid overflow during calculation + const double tiny1 = 1e-12; + const double tiny2 = 1e-10; + if (distance < tiny1) distance += tiny1; + + // (2) if there exist overlap, calculate the mesh number + // between two atoms + const int rmesh = this->talpha.get_rmesh(Rcut1, Rcut2); + + // (3) Find three dimension of 'Table_DS' + // dim1 : type pairs, equal to T1 here + // dim2 : radial orbital pairs, + // dim3 : find lmax between T1 and T2, and get lmax*2+1 + const int dim1 = T1; + int dim2 = this->talpha.DS_Opair(dim1, L1, L2, N1, N2); + int dim3 = this->talpha.DS_2Lplus1[T1]; + + //Gaunt Index + const int gindex1 = L1 * L1 + m1; + const int gindex2 = L2 * L2 + m2; + + // Peize Lin change rly, grly 2016-08-26 + vector rly; + vector> grly; + + double arr_dR[3]; + arr_dR[0] = Numerical_Orbital::getX() * this->lat0; + arr_dR[1] = Numerical_Orbital::getY() * this->lat0; + arr_dR[2] = Numerical_Orbital::getZ() * this->lat0; + + //double xdr = arr_dR[0] / distance; + //double ydr = arr_dR[1] / distance; + //double zdr = arr_dR[2] / distance; + + //======================= + // *r**l*Ylm_real + // include its derivations + //======================= + if (job == 0) + { + Ylm::rl_sph_harm(dim3 - 1, arr_dR[0], arr_dR[1], arr_dR[2], rly); + } + else + { + Ylm::grad_rl_sph_harm(dim3 - 1, arr_dR[0], arr_dR[1], arr_dR[2], rly, grly); + } + + for (int L = 0; L < dim3; L++) //maxL = dim3-1 + { + //=========================================================== + // triangle rule for L and sum of L, L1, L2 should be even + //=========================================================== + int AL = L1 + L2; + int SL = abs(L1 - L2); + + if ((L > AL) || (L < SL) || ((L - SL) % 2 == 1)) continue; + + double Interp_Slm = 0.0; + double Interp_dSlm = 0.0; + double tmpOlm0 = 0.0; + double tmpOlm1 = 0.0; + // prefactor + double i_exp = pow(-1.0, (L1 - L2 - L) / 2); + double rl = pow(distance, L); + + if (distance > tiny2) + { + Interp_Slm = i_exp * Mathzone::Polynomial_Interpolation( + talpha.Table_DSR[0][dim1][dim2][L], rmesh, MOT.dr, distance); + Interp_Slm /= rl; + } + else // distance = 0.0; + { + Interp_Slm = i_exp * talpha.Table_DSR[0][dim1][dim2][L][0]; + } + + if (job == 1)//calculate the derivative. + { + if (distance > tiny2) + { + Interp_dSlm = i_exp * Mathzone::Polynomial_Interpolation( + talpha.Table_DSR[1][dim1][dim2][L], rmesh, MOT.dr, distance); + Interp_dSlm = Interp_dSlm / pow(distance, L) - Interp_Slm * L / distance; + } + else + { + Interp_dSlm = 0.0; + } + } + + for (int m = 0; m < 2 * L + 1; m++) + { + int gindex = L * L + m; + // double tmpGaunt1 = MGT.Get_Gaunt_SH(L1, m1, L2, m2, L, m); + double tmpGaunt = MGT.Gaunt_Coefficients(gindex1, gindex2, gindex); + + tmpOlm0 = Interp_Slm * tmpGaunt; + + if (job == 1) + { + tmpOlm1 = Interp_dSlm * tmpGaunt; + } + + switch (job) + { + case 0: // calculate overlap. + { + if (NSPIN != 4) olm[0] += tmpOlm0 * rly[MGT.get_lm_index(L, m)]; + else if (olm1 != NULL) + { + olm1[0] += tmpOlm0 * rly[MGT.get_lm_index(L, m)]; + olm1[1] += 0;//tmpOlm0 * (tmp(0,0)+tmp(0,1)); + olm1[2] += 0;//tmpOlm0 * (tmp(1,0)+tmp(1,1)); + olm1[3] += tmpOlm0 * rly[MGT.get_lm_index(L, m)]; + + } + else + { + WARNING_QUIT("ORB_gen_tables::snap_psialpha", "something wrong!"); + + } + + /* + if( abs ( tmpOlm0 * rly[ MGT.get_lm_index(L, m) ] ) > 1.0e-3 ) + { + cout << " L=" << L << " m=" << m << " tmpOlm0=" << tmpOlm0 + << " rly=" << rly[ MGT.get_lm_index(L, m) ] + << " r=" << olm[0] + << endl; + } + */ + break; + } + case 1: // calculate gradient. + { + for (int ir = 0; ir < 3; ir++) + { + olm[ir] += tmpOlm0 * grly[MGT.get_lm_index(L, m)][ir] + + tmpOlm1 * rly[MGT.get_lm_index(L, m)] * arr_dR[ir] / distance; + } + break; + } + default: break; + } + }//!m + } + + return; +} diff --git a/ABACUS.develop/source/src_lcao/ORB_gen_tables.h b/ABACUS.develop/source/src_lcao/ORB_gen_tables.h index 8a257b88e6..2fe9619ee2 100644 --- a/ABACUS.develop/source/src_lcao/ORB_gen_tables.h +++ b/ABACUS.develop/source/src_lcao/ORB_gen_tables.h @@ -10,6 +10,7 @@ #include "src_global/ylm.h" #include "ORB_table_beta.h" #include "ORB_table_phi.h" +#include "ORB_table_alpha.h" //caoyu add 2020-3-18 //------------------------------------ // used to be 'Use_Overlap_Table', @@ -63,6 +64,24 @@ class ORB_gen_tables complex *nlm1=NULL, const int is=0)const; + //caoyu add 2021-03-17 + //job = 0 for vnl matrix elements + //job = 1 for its derivatives + void snap_psialpha( + double olm[], + const int& job, + const Vector3& R1, + const int& I1, + const int& l1, + const int& m1, + const int& n1, + const Vector3& R2, + const int& I2, + const int& l2, + const int& m2, + const int& n2, + complex* olm1 = NULL, + const int is = 0)const; // set as public because in hamilt_linear, // we need to destroy the tables: SR,TR,NR @@ -72,6 +91,7 @@ class ORB_gen_tables // if we want to add table for descriptors, // we should consider here -- mohan 2021-02-09 + ORB_table_alpha talpha; //caoyu add 2021-03-17 private: diff --git a/ABACUS.develop/source/src_lcao/ORB_read.cpp b/ABACUS.develop/source/src_lcao/ORB_read.cpp index e915d6ef2c..a33e946989 100644 --- a/ABACUS.develop/source/src_lcao/ORB_read.cpp +++ b/ABACUS.develop/source/src_lcao/ORB_read.cpp @@ -12,6 +12,7 @@ LCAO_Orbitals::LCAO_Orbitals() this->nchimax = 0;// this initialzied must specified this->Phi = new Numerical_Orbital[1]; this->Beta = new Numerical_Nonlocal[1]; + this->Alpha = new Numerical_Orbital[1]; this->nproj = new int[1]; this->nprojmax = 0; @@ -28,6 +29,7 @@ LCAO_Orbitals::~LCAO_Orbitals() } delete[] Phi; delete[] Beta; + delete[] Alpha; delete[] nproj; } @@ -224,6 +226,10 @@ void LCAO_Orbitals::Read_Orbitals(void) //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> if (INPUT.out_descriptor && BASIS_TYPE == "lcao") //condition: descriptor in lcao line { + + delete[] this->Alpha; + this->Alpha = new Numerical_Orbital[1]; //not related to atom type -- remain to be discussed + this->Read_Descriptor(); } @@ -993,6 +999,7 @@ void LCAO_Orbitals::Read_Descriptor() //read descriptor basis //read lmax and nchi[l] int lmax = 0; + int nchimax = 0; int nchi[10]={0}; char word[80]; if (MY_RANK == 0) @@ -1011,15 +1018,17 @@ void LCAO_Orbitals::Read_Descriptor() //read descriptor basis for (int l = 0; l <= lmax; l++) { in >> word >> word >> word >> nchi[l]; - this->nchimax_d = std::max(this->nchimax_d, nchi[l]); + nchimax = std::max(nchimax, nchi[l]); } } #ifdef __MPI Parallel_Common::bcast_int(lmax); + Parallel_Common::bcast_int(nchimax); Parallel_Common::bcast_int(nchi, lmax + 1); #endif this->lmax_d = lmax; + this->nchimax_d = nchimax; // calculate total number of chi int total_nchi = 0; for (int l = 0; l <= lmax; l++) @@ -1027,7 +1036,7 @@ void LCAO_Orbitals::Read_Descriptor() //read descriptor basis total_nchi += nchi[l]; } //OUT(ofs_running,"Total number of chi(l,n)",total_nchi); - this->Alpha.phiLN = new Numerical_Orbital_Lm[total_nchi]; + this->Alpha[0].phiLN = new Numerical_Orbital_Lm[total_nchi]; int meshr; // number of mesh points int meshr_read; @@ -1178,7 +1187,7 @@ void LCAO_Orbitals::Read_Descriptor() //read descriptor basis delete[] inner; ofs_running << setw(12) << unit << endl; - this->Alpha.phiLN[count].set_orbital_info( + this->Alpha[0].phiLN[count].set_orbital_info( "H", //any type 1, //any type L, //angular momentum L @@ -1203,7 +1212,7 @@ void LCAO_Orbitals::Read_Descriptor() //read descriptor basis } } in.close(); - this->Alpha.set_orbital_info( + this->Alpha[0].set_orbital_info( 1, // any type "H", // any label lmax, diff --git a/ABACUS.develop/source/src_lcao/ORB_read.h b/ABACUS.develop/source/src_lcao/ORB_read.h index f2495936eb..d39e449f68 100644 --- a/ABACUS.develop/source/src_lcao/ORB_read.h +++ b/ABACUS.develop/source/src_lcao/ORB_read.h @@ -63,7 +63,7 @@ class LCAO_Orbitals //caoyu add 2021-3-10 // descriptor bases, saved as one-type atom orbital - Numerical_Orbital Alpha; + Numerical_Orbital* Alpha; // initialized in input.cpp double ecutwfc; diff --git a/ABACUS.develop/source/src_lcao/ORB_table_alpha.cpp b/ABACUS.develop/source/src_lcao/ORB_table_alpha.cpp new file mode 100644 index 0000000000..d0f639d28d --- /dev/null +++ b/ABACUS.develop/source/src_lcao/ORB_table_alpha.cpp @@ -0,0 +1,468 @@ +//caoyu add 2021-03-17 + +#include "ORB_table_alpha.h" +#include "ORB_read.h" +#include "../src_pw/global.h" +#include +#include "../src_ri/exx_abfs.h" +#include "../src_io/winput.h" + +double ORB_table_alpha::dr = -1.0; + +ORB_table_alpha::ORB_table_alpha() +{ + destroy_nr = false; + + ntype = 0; + lmax = 0; + kmesh = 0; + Rmax = 0.0; + dr = 0.0; + dk = 0.0; + + nlm = 0; + Rmesh = 0; + + kpoint = new double[1]; + r = new double[1]; + rab = new double[1]; + kab = new double[1]; + DS_2Lplus1 = new int[1]; +} + +ORB_table_alpha::~ORB_table_alpha() +{ + delete[] kpoint; + delete[] r; + delete[] rab; + delete[] kab; + delete[] DS_2Lplus1; +} + +void ORB_table_alpha::allocate +( + const int& ntype_in, + const int& lmax_in, + const int& kmesh_in, + const double& Rmax_in, + const double& dr_in, + const double& dk_in +) +{ + TITLE("ORB_table_alpha", "allocate"); + + this->ntype = ntype_in;// type of elements. + this->lmax = lmax_in; + this->kmesh = kmesh_in; + this->Rmax = Rmax_in; + this->dr = dr_in; + this->dk = dk_in; + + assert(ntype > 0); + assert(lmax >= 0); + assert(kmesh > 0.0); + assert(Rmax >= 0.0); + assert(dr > 0.0); + assert(dk > 0.0); + + // calculated from input parameters + this->nlm = (2 * lmax + 1) * (2 * lmax + 1); + this->Rmesh = static_cast(Rmax / dr) + 4; + if (Rmesh % 2 == 0) + { + ++Rmesh; + } + + // OUT(ofs_running,"lmax",lmax); + // OUT(ofs_running,"Rmax (Bohr)",Rmax); + // OUT(ofs_running,"dr (Bohr)",dr); + // OUT(ofs_running,"dk",dk); + // OUT(ofs_running,"nlm",nlm); + // OUT(ofs_running,"kmesh",kmesh); + + delete[] kpoint; + delete[] r; + kpoint = new double[kmesh]; + r = new double[Rmesh]; + + delete[] rab; + delete[] kab; + kab = new double[kmesh]; + rab = new double[Rmesh]; + + for (int ik = 0; ik < kmesh; ik++) + { + kpoint[ik] = ik * dk_in; + kab[ik] = dk_in; + } + + for (int ir = 0; ir < Rmesh; ir++) + { + r[ir] = ir * dr; + rab[ir] = dr; + } + + // OUT(ofs_running,"allocate kpoint, r, rab, kab","Done"); + return; +} + + +int ORB_table_alpha::get_rmesh(const double& R1, const double& R2) +{ + int rmesh = static_cast((R1 + R2) / ORB_table_alpha::dr) + 5; + //mohan update 2009-09-08 +1 ==> +5 + //considering interpolation or so on... + if (rmesh % 2 == 0) rmesh++; + + if (rmesh <= 0) + { + ofs_warning << "\n R1 = " << R1 << " R2 = " << R2; + ofs_warning << "\n rmesh = " << rmesh; + WARNING_QUIT("ORB_table_alpha::get_rmesh", "rmesh <= 0"); + } + return rmesh; +} + + + +void ORB_table_alpha::cal_S_PhiAlpha_R( + Sph_Bessel_Recursive::D2* pSB, // mohan add 2021-03-06 + const int& l, + const Numerical_Orbital_Lm& n1, + const Numerical_Orbital_Lm& n2, + const int& rmesh, + double* rs, + double* drs) +{ + timer::tick("ORB_table_alpha", "S_PhiAlpha_R"); + + assert(kmesh > 0); + + //start calc + double* k1_dot_k2 = new double[kmesh]; + + for (int ik = 0; ik < kmesh; ik++) + { + k1_dot_k2[ik] = n1.getPsi_k(ik) * n2.getPsi_k(ik); + } + + //previous version + double* integrated_func = new double[kmesh]; + + const vector>& jlm1 = pSB->get_jlx()[l - 1]; + const vector>& jl = pSB->get_jlx()[l]; + const vector>& jlp1 = pSB->get_jlx()[l + 1]; + for (int ir = 0; ir < rmesh; ir++) + { + ZEROS(integrated_func, kmesh); + double temp = 0.0; + + for (int ik = 0; ik < kmesh; ik++) + { + integrated_func[ik] = jl[ir][ik] * k1_dot_k2[ik]; + } + // Call simpson integration + Mathzone::Simpson_Integral(kmesh, integrated_func, kab, temp); + rs[ir] = temp * FOUR_PI; + + //drs + double temp1, temp2; + + if (l > 0) + { + for (int ik = 0; ik < kmesh; ik++) + { + integrated_func[ik] = jlm1[ir][ik] * k1_dot_k2[ik] * kpoint[ik]; + } + + Mathzone::Simpson_Integral(kmesh, integrated_func, kab, temp1); + } + + + for (int ik = 0; ik < kmesh; ik++) + { + integrated_func[ik] = jlp1[ir][ik] * k1_dot_k2[ik] * kpoint[ik]; + } + + Mathzone::Simpson_Integral(kmesh, integrated_func, kab, temp2); + + if (l == 0) + { + drs[ir] = -FOUR_PI * temp2; + } + else + { + drs[ir] = FOUR_PI * (temp1 * l - (l + 1) * temp2) / (2.0 * l + 1); + } + } + + //liaochen modify on 2010/4/22 + //special case for R=0 + //we store Slm(R) / R**l at the fisrt point, rather than Slm(R) + if (l > 0) + { + ZEROS(integrated_func, kmesh); + double temp = 0.0; + + for (int ik = 0; ik < kmesh; ik++) + { + integrated_func[ik] = k1_dot_k2[ik] * pow(kpoint[ik], l); + } + + // Call simpson integration + Mathzone::Simpson_Integral(kmesh, integrated_func, kab, temp); + rs[0] = FOUR_PI / Mathzone_Add1::dualfac(2 * l + 1) * temp; + } + + delete[] integrated_func; + + + delete[] k1_dot_k2; + timer::tick("ORB_table_alpha", "S_PhiAlpha_R"); + return; +} + + +void ORB_table_alpha::init_Table_Alpha(Sph_Bessel_Recursive::D2* pSB) +{ + TITLE("ORB_table_alpha", "init_Table_Alpha"); + timer::tick("ORB_table_alpha", "init_Table_Alpha", 'D'); + + // (1) allocate 1st dimension ( overlap, derivative) + this->Table_DSR = new double**** [2]; + // (2) allocate 2nd dimension ( overlap, derivative) + this->Table_DSR[0] = new double*** [this->ntype]; + this->Table_DSR[1] = new double*** [this->ntype]; + + // <1Phi|2Alpha> + for (int T1 = 0; T1 < ntype; T1++) // type 1 is orbital + { + const int Lmax1 = ORB.Phi[T1].getLmax(); + const int Lmax2 = ORB.Alpha[0].getLmax(); + const int lmax_now = std::max(Lmax1, Lmax2); + int L2plus1 = 2 * lmax_now + 1; + //------------------------------------------------------------- + // how many + // here we count all possible psi with (L,N) index for type T1. + //------------------------------------------------------------- + const int pairs_chi = ORB.Phi[T1].getTotal_nchi() * ORB.Alpha[0].getTotal_nchi(); + + if (pairs_chi == 0)continue; + + // init 2nd dimension + this->Table_DSR[0][T1] = new double** [pairs_chi]; + this->Table_DSR[1][T1] = new double** [pairs_chi]; + + const double Rcut1 = ORB.Phi[T1].getRcut(); + for (int L1 = 0; L1 < Lmax1 + 1; L1++) + { + for (int N1 = 0; N1 < ORB.Phi[T1].getNchi(L1); N1++) + { + for (int L2 = 0; L2 < Lmax2 + 1; L2++) + { + for (int N2 = 0; N2 < ORB.Alpha[0].getNchi(L2); N2++) + { + // get the second index. + const int Opair = this->DS_Opair(T1, L1, L2, N1, N2); + + // init 3rd dimension + this->Table_DSR[0][T1][Opair] = new double* [L2plus1]; + this->Table_DSR[1][T1][Opair] = new double* [L2plus1]; + + const double Rcut1 = ORB.Phi[T1].getRcut(); + const double Rcut2 = ORB.Alpha[0].getRcut(); + assert(Rcut1 > 0.0 && Rcut1 < 100); + assert(Rcut2 > 0.0 && Rcut2 < 100); + + const int rmesh = this->get_rmesh(Rcut1, Rcut2); + assert(rmesh < this->Rmesh); + + //L=|L1-L2|,|L1-L2|+2,...,L1+L2 + const int SL = abs(L1 - L2); + const int AL = L1 + L2; + + for (int L = 0; L < L2plus1; L++) + { + //Allocation + this->Table_DSR[0][T1][Opair][L] = new double[rmesh]; + this->Table_DSR[1][T1][Opair][L] = new double[rmesh]; + + Memory::record("ORB_table_alpha", "Table_DSR", + 2 * this->ntype * pairs_chi * rmesh, "double"); + + //for those L whose Gaunt Coefficients = 0, we + //assign every element in Table_DSR as zero + if ((L > AL) || (L < SL) || ((L - SL) % 2 == 1)) + { + ZEROS(Table_DSR[0][T1][Opair][L], rmesh); + ZEROS(Table_DSR[1][T1][Opair][L], rmesh); + + continue; + } + + this->cal_S_PhiAlpha_R( + pSB, // mohan add 2021-03-06 + L, + ORB.Phi[T1].PhiLN(L1, N1), + ORB.Alpha[0].PhiLN(L2, N2), // mohan update 2011-03-07 + rmesh, + this->Table_DSR[0][T1][Opair][L], + this->Table_DSR[1][T1][Opair][L]); + }// end L2plus1 + }// end N2 + }// end L2 + }// end N1 + }// end L1 + }// end T1 + destroy_nr = true; + + + // OUT(ofs_running,"allocate non-local potential matrix","Done"); + timer::tick("ORB_table_alpha", "init_Table_Alpha", 'D'); + return; +} + + +void ORB_table_alpha::Destroy_Table_Alpha(void) +{ + if (!destroy_nr) return; + + const int ntype = ORB.get_ntype(); + for (int ir = 0; ir < 2; ir++) + { + for (int T1 = 0; T1 < ntype; T1++) + { + const int Lmax1 = ORB.Phi[T1].getLmax(); + const int Lmax2 = ORB.Alpha[0].getLmax(); + const int lmax_now = std::max(Lmax1, Lmax2); + const int pairs = ORB.Phi[T1].getTotal_nchi() * ORB.Alpha[0].getTotal_nchi(); + + // mohan fix bug 2011-03-30 + if (pairs == 0) continue; + for (int dim2 = 0; dim2 < pairs; dim2++) + { + for (int L = 0; L < 2*lmax_now + 1; L++) + { + delete[] Table_DSR[ir][T1][dim2][L]; + } + delete[] Table_DSR[ir][T1][dim2]; + } + delete[] Table_DSR[ir][T1]; + } + delete[] Table_DSR[ir]; + } + delete[] Table_DSR; + return; +} + +void ORB_table_alpha::init_DS_2Lplus1(void) +{ + TITLE("Make_Overlap_Table", "init_DS_2Lplus1"); + assert(this->ntype > 0); + delete[] DS_2Lplus1; + DS_2Lplus1=new int[ntype]; // 2Lmax+1 for each T1 + + int index = 0; + for (int T1 = 0; T1 < ntype; T1++) + { + this->DS_2Lplus1[T1] = max(ORB.Phi[T1].getLmax(), ORB.Alpha[0].getLmax()) * 2 + 1; + } + return; +} + +void ORB_table_alpha::init_DS_Opair(void) +{ + const int lmax = ORB.get_lmax(); + const int nchimax = ORB.get_nchimax(); + const int lmax_d = ORB.get_lmax_d(); + const int nchimax_d = ORB.get_nchimax_d(); + assert(lmax + 1 > 0); + assert(lmax_d + 1 > 0); + assert(nchimax > 0); + assert(nchimax_d > 0); + + this->DS_Opair.create(this->ntype, lmax+1, lmax_d+1, nchimax, nchimax_d); + + // <1psi|2beta> + // 1. orbital + for (int T1 = 0; T1 < ntype; T1++) //alpha is not related to atom type ! + { + int index = 0; + for (int L1 = 0; L1 < ORB.Phi[T1].getLmax() + 1; L1++) + { + for (int N1 = 0; N1 < ORB.Phi[T1].getNchi(L1); N1++) + { + for (int L2 = 0; L2 < ORB.Alpha[0].getLmax() + 1; L2++) + { + for (int N2 = 0; N2 < ORB.Alpha[0].getNchi(L2); N2++) + { + this->DS_Opair(T1, L1, L2, N1, N2) = index; + ++index; + } + } + } + } + } + return; +} + +//caoyu add 2021-03-20 +void ORB_table_alpha::print_Table_DSR(void) +{ + TITLE("ORB_table_alpha", "print_Table_DSR"); + NEW_PART("Overlap table S between lcao orbital and descriptor basis : S_{I_mu_alpha}"); + + ofstream ofs; + stringstream ss; + // the parameter 'winput::spillage_outdir' is read from INPUTw. + ss << winput::spillage_outdir << "/" << "S_I_mu_alpha.dat"; + if (MY_RANK == 0) + { + ofs.open(ss.str().c_str()); + } + + for (int T1 = 0; T1 < this->ntype; T1++) //T1 + { + const int Lmax1 = ORB.Phi[T1].getLmax(); + const int Lmax2 = ORB.Alpha[0].getLmax(); + for (int L1 = 0; L1 < Lmax1 + 1; L1++) + { + for (int N1 = 0; N1 < ORB.Phi[T1].getNchi(L1); N1++) + { + for (int L2 = 0; L2 < Lmax2 + 1; L2++) + { + for (int N2 = 0; N2 < ORB.Alpha[0].getNchi(L2); N2++) + { + const int Opair = this->DS_Opair(T1, L1, L2, N1, N2); //Opair + ofs <