output 5D overlap table between lcao and descriptor orbitals

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;
 }

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<double> &R1, const Vector3<do
 	return dR.norm() * this->lat0;	
 }
 
+//caoyu add 2021-03-17
+void ORB_gen_tables::snap_psialpha(
+	double olm[],
+	const int& job,
+	const Vector3<double>& R1,
+	const int& T1,
+	const int& L1,
+	const int& m1,
+	const int& N1,
+	const Vector3<double>& R2,
+	const int& T2,
+	const int& L2,
+	const int& m2,
+	const int& N2,
+	complex<double>* 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<double> rly;
+	vector<vector<double>> 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;
+}

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<double> *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<double>& R1,
+		const int& I1,
+		const int& l1,
+		const int& m1,
+		const int& n1,
+		const Vector3<double>& R2,
+		const int& I2,
+		const int& l2,
+		const int& m2,
+		const int& n2,
+		complex<double>* 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:
 

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,23 +1018,25 @@ 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++)
 	{
 		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,

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;