Test: add module_io/bessel_basis unit test

source/module_io/bessel_basis.cpp

@@ -1,5 +1,5 @@
-#include "bessel_basis.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
+#include "bessel_basis.h"
 #include "module_base/parallel_common.h"
 #include "module_base/math_integral.h"
 #include "module_base/math_sphbes.h"
@@ -25,9 +25,11 @@ void Bessel_Basis::init(
 	const bool &smooth,
 	const double &sigma,
 	const double &rcut_in,
-	const double &tol_in,	
+	const double &tol_in,
+	const UnitCell& ucell,
 	const double &dk,
-	const double &dr)
+	const double &dr
+	)
 {
 	ModuleBase::TITLE("Bessel_Basis", "init");
 	this->Dk = dk;
@@ -63,14 +65,13 @@ void Bessel_Basis::init(
 	if( start_from_file )
 	{
 		// setup C4
-		this->allocate_C4(ntype, lmax_in, GlobalC::ucell.nmax, Ecut_number);
-
+		this->allocate_C4(ntype, lmax_in, ucell.nmax, Ecut_number, ucell);
 		// check tolerence
-		this->readin_C4("INPUTs", ntype, ecut, rcut, Ecut_number, tolerence);
+		this->readin_C4("INPUTs", ntype, ecut, rcut, Ecut_number, tolerence, ucell);
 #ifdef __MPI
 		Parallel_Common::bcast_double( C4.ptr, C4.getSize() );
 #endif
-		this->init_Faln(ntype, lmax_in, GlobalC::ucell.nmax, Ecut_number);
+		this->init_Faln(ntype, lmax_in, ucell.nmax, Ecut_number, ucell);
 	}
 
 	return;
@@ -124,7 +125,8 @@ void Bessel_Basis::init_Faln(
 	const int &ntype,
 	const int &lmax,
 	const int &nmax,
-	const int &ecut_number)
+	const int &ecut_number,
+	const UnitCell& ucell)
 {
 	ModuleBase::TITLE("Bessel_Basis","init_Faln");
 	ModuleBase::timer::tick("Spillage","init_Faln");
@@ -135,9 +137,9 @@ void Bessel_Basis::init_Faln(
 	this->nwfc = 0;
 	for(int it=0; it<ntype; it++)
 	{
-		for(int il=0; il<GlobalC::ucell.atoms[it].nwl+1; il++)
+		for(int il=0; il<ucell.atoms[it].nwl+1; il++)
 		{
-			for(int in=0; in<GlobalC::ucell.atoms[it].l_nchi[il]; in++)
+			for(int in=0; in<ucell.atoms[it].l_nchi[il]; in++)
 			{
 				for(int ie=0; ie<ecut_number; ie++)
 				{
@@ -340,7 +342,8 @@ void Bessel_Basis::readin_C4(
 	const int &ecut,
 	const int &rcut,
 	const int &ecut_number,
-	const double &tolerence)
+	const double &tolerence,
+	const UnitCell& ucell)
 {
 	ModuleBase::TITLE("Bessel_Basis","readin_C4");
 
@@ -361,9 +364,9 @@ void Bessel_Basis::readin_C4(
 		{
 			std::string filec4;
 			ifs >> filec4;
-			for(int il=0; il< GlobalC::ucell.atoms[it].nwl+1; il++)
+			for(int il=0; il< ucell.atoms[it].nwl+1; il++)
 			{
-				for(int in=0; in< GlobalC::ucell.atoms[it].l_nchi[il]; in++)
+				for(int in=0; in< ucell.atoms[it].l_nchi[il]; in++)
 				{
 					//for tests
 					//std::cout << "\n" << std::setw(5) << it << std::setw(5) << il << std::setw(5) << in;
@@ -448,17 +451,18 @@ void Bessel_Basis::allocate_C4(
 	const int &ntype,
 	const int &lmax,
 	const int &nmax,
-	const int &ecut_number)
+	const int &ecut_number,
+	const UnitCell& ucell)
 {
 	ModuleBase::TITLE("Bessel_Basis","allocate_C4");
 
 	this->C4.create(ntype, lmax+1, nmax, ecut_number);
 
 	for(int it=0; it<ntype; it++)
 	{
-		for(int il=0; il<GlobalC::ucell.atoms[it].nwl+1; il++)
+		for(int il=0; il<ucell.atoms[it].nwl+1; il++)
 		{
-			for(int in=0; in<GlobalC::ucell.atoms[it].l_nchi[il]; in++)
+			for(int in=0; in<ucell.atoms[it].l_nchi[il]; in++)
 			{
 				for(int ie=0; ie<ecut_number; ie++)
 				{

source/module_io/bessel_basis.h

@@ -9,6 +9,8 @@
 #include "../module_base/global_variable.h"
 #include "../module_base/realarray.h"
 
+#include "../module_cell/unitcell.h"
+
 //==========================================================
 // CLASS :
 // NAME :  Bessel_Basis
@@ -19,15 +21,20 @@ class Bessel_Basis
 	Bessel_Basis();
 	~Bessel_Basis();
 
-	//--------------------------------------------------------------------------
-	// used for a specific group of C4 coefficients.
-	// mohan updated 2021-01-04, 
-	// mohan added a new input parameter lmax_in,
-	// if we only generate numerical atomic orbitals based on
-	// spherical Bessel functions, lmax_in = ucell.lmax
-	// However, if we want to generate spherical Bessel functions for descriptor,
-	// then the lmax_in is controlled by user.
-	//--------------------------------------------------------------------------
+	/// @brief Initialization of Bessel function related matrices.
+	/// @details Used for a specific group of C4 coefficients. 2021-01-04, mohan added a new input parameter lmax_in, if we only generate numerical atomic orbitals based on spherical Bessel functions, lmax_in = ucell.lmax. However, if we want to generate Spherical Bessel functions (SBF) for descriptor, then the lmax_in is controlled by user.
+	/// @note This function is called in module_io/numerical_basis.cpp and module_io/numerical_descriptor.cpp
+	/// @param start_from_file whether read C4 coefficients stored in external files
+	/// @param ecutwfc cutoff for numerical atomic orbitals
+	/// @param ntype atom types
+	/// @param lmax_in maximal angular momentum for numerical orbitals
+	/// @param smooth whether smooth SBFs when perform integration to calculate value of matrix element of TableOne. For details, see J. Phys.: Condens. Matter 22 (2010) 445501
+	/// @param sigma stddev of Gaussian function for smoothing SBFs
+	/// @param rcut_in cutoff radius for SBFs
+	/// @param tol_in accurancy control for SBFs
+	/// @param dk kspace grid
+	/// @param dr realspace grid
+	/// @param ucell UnitCell class object, ucell.nmax will be used in this function
 	void init( 
 		const bool start_from_file,
 		const double &ecutwfc,
@@ -37,68 +44,137 @@ class Bessel_Basis
 		const double &sigma,
 		const double &rcut_in,
 		const double &tol_in,
+		const UnitCell& ucell,
 		const double &dk = 0.01,
-		const double &dr = 0.01);
-
+		const double &dr = 0.01
+		);
+	/// @brief return number of SBFs used for one `chi` (see details for more information)
+	/// @details atomic orbital is constructed always with not only one set of SBFs. For different sets, they are marked with different `chi`(s), similar with concept of contracted GTOs. For one `chi`, it is 'q' the summation index, and q is in SBFs like: j_l(q*r), where l is the order of SBF.
+	/// @return number of SBFs
 	const int& get_ecut_number() const { return Ecut_number;}
 
-	// get value from interpolation
+	/// @brief Cubic spline interpolation for matrix Faln
+	/// @param it atom type index
+	/// @param l angular momentum
+	/// @param ic chi index
+	/// @param gnorm norm of G+k vector
+	/// @return interpolated value
 	double Polynomial_Interpolation(const int &it, const int &l, const int &ic, const double &gnorm)const;
+	/// @brief Cubic spline interpolation for matrix TableOne
+	/// @param l angular momentum
+	/// @param ie q index (see explanation in note of function BesselBasis::get_ecut_number())
+	/// @param gnorm norm of G+k vector
+	/// @return interpolated value
 	double Polynomial_Interpolation2(const int &l, const int &ie, const double &gnorm)const;
 
-	// get ecut : get energy cut off, which is used to truncate spherical Bessel function Jlq. 
-	// get rcut : cut off radius(angstrom) of radial spherical Bessel function Jlq.
+
+	/// @brief  get energy cutoff, which is used to truncate SBF Jlq. 
+	/// @param  
+	/// @return energy cutoff in Ry
 	const double &get_ecut(void) const {return ecut;}
+	/// @brief cutoff radius of radial SBF Jlq.
+	/// @param  
+	/// @return cutoff radius in a.u.
 	const double &get_rcut(void) const {return rcut;}
+
 	const double &get_tolerence(void) const {return tolerence;}
 
-	// get_smooth and get_sigma. mohan add 2009-08-28
-	// in this case, the Jlq are not the true Jlq.
+
+	/// @brief check if SBFs are smoothed (mohan add 2009-08-28)
+	/// @attention in this case, the Jlq are not the true Jlq.
+	/// @param  
+	/// @return boolean whether SBFs are smoothed
 	const bool &get_smooth(void) const {return smooth;}
+	/// @brief get sigma the stddev (standard deviation) used in smooth function (Gaussian function)
+	/// @param  
+	/// @return stddev of smooth function
 	const double &get_sigma(void) const {return sigma;}
 
 private:
-	// the most important array to calculate spillage
-	ModuleBase::realArray Faln;// (ntype, max_l+1, max_n, ik)
+	/// @brief the most important array to calculate spillage, has dimension (ntype, lmax+1, max_n, nk)
+	ModuleBase::realArray Faln;
 
-	// Coefficients to be optimized!
+	/// @brief Coefficients to be optimized!
 	ModuleBase::realArray C4;
 
+	/// @brief matrix whose elements are int{dr r^2 j_l(qr)*j_l(kr)}, has dimension (lmax+1, nq, nk)
 	ModuleBase::realArray TableOne;
+
+	/// @brief mesh of k vector, k is in j_l(k*r)
 	int kmesh;
+	/// @brief grid of k
 	double Dk;
+	/// @brief number of q vector, q is in j_l(q*r)
 	int Ecut_number;
+	/// @brief Cutoff radius (in a.u.) of SBFs, for any SBF j_l(qr), r>=rcut, j_l(q*r) = 0 (if not smoothed)
 	double rcut;
+	/// @brief energy cutoff for determining kmesh and number of SBFs
 	double ecut;
 	double tolerence;
-
-	bool smooth; // mohan add 2009-01-18
+	/// @brief whether smooth SBFs around cutoff radius, resulting in non-zero values. For importance of smooth of SBFs, see J. Phys.: Condens. Matter 22 (2010) 445501, eqn 6. (mohan add 2009-01-18)
+	bool smooth;
+	/// @brief stddev of smooth function (Gaussian function, centered at rcut)
 	double sigma;
 
-	//========================================================
-	// MEMBER FUNCTIONS :
-	// NAME : readin ( read in parameters )
-	//========================================================
+	/// @brief Allocate memory for C4 matrix and initialize all elements to one.
+	/// @param ntype number of atom types
+	/// @param lmax maximal angular momentum of localized orbitals
+	/// @param nmax maximal principal quantum number of localized orbitals
+	/// @param ecut_number number of SBFs
 	void allocate_C4(
 		const int &ntype,
 		const int &lmax, 
 		const int &nmax,
-		const int &ecut_number);
-
+		const int &ecut_number,
+		const UnitCell& ucell
+		);
+
+	/// @brief Read C4 from external file. Presently an O(N^2) search algorithm is used. A HTML parser is needed in the future to improve performance.
+	/// @param name name of external file where C4-stored file information is contained
+	/// @param ntype number of atom types
+	/// @param ecut energy cutoff
+	/// @param rcut cutoff radius
+	/// @param ecut_number number of SBFs
+	/// @param tolerence accurancy of SBFs, here only used for consistency check
 	void readin_C4(
 		const std::string &name,
 		const int &ntype,
 		const int &ecut,
 		const int &rcut,
 		const int &ecut_number,
-		const double &tolerence);
+		const double &tolerence,
+		const UnitCell& ucell
+		);
 
 	void init_TableOne(void);
 
-	void init_Faln( const int &ntype,const int &lmax,const int &nmax,const int &ecut_number);
-	int nwfc; // number of localized wave functions
-
-	// init table, used for output overlap Q.
+	/// @brief calculate F_{aln}(it, il, in, ik) = sum_{ie}{C4(it, il, in, ie)*TableOne(il, ie, ik)}, where TableOne is overlap integral between two spherical bessel functions (jle(r) and jlk(r))
+	/// @param ntype number of atomtype
+	/// @param lmax maximal angular momentum
+	/// @param nmax maximal chi
+	/// @param ecut_number number of SBFs
+	void init_Faln(
+		const int &ntype,
+		const int &lmax,
+		const int &nmax,
+		const int &ecut_number,
+		const UnitCell& ucell
+		);
+
+	/// @brief number of localized wave functions
+	int nwfc;
+
+	/// @brief calculate element value of TableOne matrix
+	/// @details (be called in Bessel_Basis::init(), used for outputing overlap Q matrix) initialize the table whose matrix element is the result of integral int{dr r^2 jle(r)*jlk(r)}, TableOne has three subscript (l, ie, ik), the first runs over orbitals' angular momentum and ie, ik run over ecut_number and kmesh SBFs
+	/// @param smooth_in whether jle(r) SBF is smoothed by a Gaussian function
+	/// @param sigma_in stddev for controlling smearing of Gaussian function for smoothing jle(r)
+	/// @param ecutwfc planewave kinetic energy cutoff for controlling kspace sampling
+	/// @param rcut cutoff radius of SBFs
+	/// @param dr realspace grid
+	/// @param dk kspace grid
+	/// @param lmax maximal angular momentum for SBFs
+	/// @param ecut_number number of SBFs
+	/// @param tolerence accurancy of SBFs
 	void init_TableOne(
 		const bool smooth_in,
 		const double &sigma_in,

source/module_io/numerical_basis.cpp

@@ -43,7 +43,9 @@ void Numerical_Basis::start_from_file_k(const int& ik, ModuleBase::ComplexMatrix
 			INPUT.bessel_nao_smooth,
 			INPUT.bessel_nao_sigma,
 			INPUT.bessel_nao_rcut,
-			INPUT.bessel_nao_tolerence );
+			INPUT.bessel_nao_tolerence,
+            GlobalC::ucell
+            );
         this->mu_index = this->init_mu_index();
         this->init_label = true;
     }
@@ -71,7 +73,9 @@ void Numerical_Basis::output_overlap(const psi::Psi<std::complex<double>>& psi,
 			INPUT.bessel_nao_smooth,
 			INPUT.bessel_nao_sigma,
 			INPUT.bessel_nao_rcut,
-			INPUT.bessel_nao_tolerence );
+			INPUT.bessel_nao_tolerence,
+            GlobalC::ucell
+            );
         this->mu_index = this->init_mu_index();
         this->init_label = true;
     }

source/module_io/numerical_descriptor.cpp

@@ -53,7 +53,9 @@ void Numerical_Descriptor::output_descriptor(const psi::Psi<std::complex<double>
 		INPUT.bessel_descriptor_smooth,
 		INPUT.bessel_descriptor_sigma,
 		rcut_in,
-		tol_in );
+		tol_in,
+        GlobalC::ucell
+        );
 	this->nmax = Numerical_Descriptor::bessel_basis.get_ecut_number();
     this->init_mu_index();
     this->init_label = true;

source/module_io/test/CMakeLists.txt

@@ -111,3 +111,9 @@ AddTest(
   TARGET io_parse_args
   SOURCES parse_args_test.cpp
 )
+
+AddTest(
+  TARGET io_bessel_basis_test
+  LIBS ${math_libs} base device
+  SOURCES bessel_basis_test.cpp ../bessel_basis.cpp
+)

source/module_io/test/INPUTs

@@ -0,0 +1,3 @@
+<FILE>
+./support/BesselBasis_UnitTest_C4_AtomType0.html
+</FILE>