Feature: add ABFS_ORBITAL and ABFS_JLES_ORBITAL in class Exx_Opt_Orb

source/source_cell/bcast_cell.cpp

@@ -112,8 +112,14 @@ namespace unitcell
 
         #ifdef __EXX
         ModuleBase::bcast_data_cereal(GlobalC::exx_info.info_ri.files_abfs,
-                                    MPI_COMM_WORLD,
-                                    0);
+                                      MPI_COMM_WORLD,
+                                      0);
+        ModuleBase::bcast_data_cereal(GlobalC::exx_info.info_opt_abfs.files_abfs,
+                                      MPI_COMM_WORLD,
+                                      0);
+        ModuleBase::bcast_data_cereal(GlobalC::exx_info.info_opt_abfs.files_jles,
+                                      MPI_COMM_WORLD,
+                                      0);
         #endif
         return;
     #endif

source/source_cell/read_atom_species.cpp

@@ -112,6 +112,16 @@ bool read_atom_species(std::ifstream& ifa,
                 std::string ofile;
                 ifa >> ofile;
                 GlobalC::exx_info.info_ri.files_abfs.push_back(ofile);
+                GlobalC::exx_info.info_opt_abfs.files_abfs.push_back(ofile);
+            }
+        }
+        if( ModuleBase::GlobalFunc::SCAN_LINE_BEGIN(ifa, "ABFS_JLES_ORBITAL") )
+        {
+            for(int i=0; i<ntype; i++)
+            {
+                std::string ofile;
+                ifa >> ofile;
+                GlobalC::exx_info.info_opt_abfs.files_jles.push_back(ofile);
             }
         }
     }

source/source_hamilt/module_xc/exx_info.h

@@ -82,6 +82,11 @@ struct Exx_Info
         int abfs_Lmax = 0; // tmp
         double ecut_exx = 60;
         double tolerence = 1E-12;
+        std::vector<std::string> files_jles;
+
+        double pca_threshold = 0;
+        std::vector<std::string> files_abfs;
+
         double kmesh_times = 4;
     };
     Exx_Info_Opt_ABFs info_opt_abfs;

source/source_io/input_conv.cpp

@@ -473,6 +473,7 @@ void Input_Conv::Convert()
         GlobalC::exx_info.info_ri.ccp_rmesh_times = std::stod(PARAM.inp.exx_ccp_rmesh_times);
         GlobalC::exx_info.info_ri.exx_symmetry_realspace = PARAM.inp.exx_symmetry_realspace;
 
+        GlobalC::exx_info.info_opt_abfs.pca_threshold = PARAM.inp.exx_pca_threshold;
         GlobalC::exx_info.info_opt_abfs.abfs_Lmax  = PARAM.inp.exx_opt_orb_lmax;
         GlobalC::exx_info.info_opt_abfs.ecut_exx = PARAM.inp.exx_opt_orb_ecut;
         GlobalC::exx_info.info_opt_abfs.tolerence = PARAM.inp.exx_opt_orb_tolerence;

source/source_io/read_input_item_exx_dftu.cpp

@@ -337,12 +337,6 @@ void ReadInput::item_exx()
         Input_Item item("exx_opt_orb_lmax");
         item.annotation = "the maximum l of the spherical Bessel functions for opt ABFs";
         read_sync_int(input.exx_opt_orb_lmax);
-        item.check_value = [](const Input_Item& item, const Parameter& para) {
-            if (para.input.exx_opt_orb_lmax < 0)
-            {
-                ModuleBase::WARNING_QUIT("ReadInput", "exx_opt_orb_lmax must >= 0");
-            }
-        };
         this->add_item(item);
     }
     {

source/source_io/test_serial/read_input_item_test.cpp

@@ -1403,14 +1403,6 @@ TEST_F(InputTest, Item_test2)
         output = testing::internal::GetCapturedStdout();
         EXPECT_THAT(output, testing::HasSubstr("NOTICE"));
     }
-    { // exx_opt_orb_lmax
-        auto it = find_label("exx_opt_orb_lmax", readinput.input_lists);
-        param.input.exx_opt_orb_lmax = -1;
-        testing::internal::CaptureStdout();
-        EXPECT_EXIT(it->second.check_value(it->second, param), ::testing::ExitedWithCode(1), "");
-        output = testing::internal::GetCapturedStdout();
-        EXPECT_THAT(output, testing::HasSubstr("NOTICE"));
-    }
     { // exx_opt_orb_ecut
         auto it = find_label("exx_opt_orb_ecut", readinput.input_lists);
         param.input.exx_opt_orb_ecut = -1;

source/source_lcao/module_ri/exx_abfs-io.cpp

@@ -36,7 +36,7 @@ std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> Exx_Abfs::IO::constr
 	const double kmesh_times )
 {
 	std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> 
-		&&abfs = construct_abfs( orbs, files_abfs, kmesh_times );
+		abfs = construct_abfs( orbs, files_abfs, kmesh_times );
 
 	assert( abfs.size() == abfs_pre.size() );
 	for( size_t T=0; T!=abfs.size(); ++T )
@@ -73,7 +73,7 @@ std::vector<std::vector<Numerical_Orbital_Lm>> Exx_Abfs::IO::construct_abfs_T(
 
 	std::ifstream ifs( file_name.c_str() );
 	if(!ifs)
-		throw std::runtime_error(" Can't find the abfs ORBITAL file.");	
+		throw std::runtime_error(" Can't find the abfs ORBITAL file " + file_name);
 
 	while( ifs.good() )
 	{
@@ -130,10 +130,22 @@ std::vector<std::vector<Numerical_Orbital_Lm>> Exx_Abfs::IO::construct_abfs_T(
 		{
 			ModuleBase::GlobalFunc::READ_VALUE( ifs, N_size[8] );
 		}
+		else if ( "Lorbital-->"==word )
+		{
+			ModuleBase::GlobalFunc::READ_VALUE( ifs, N_size[9] );
+		}
+		else if ( "Morbital-->"==word )
+		{
+			ModuleBase::GlobalFunc::READ_VALUE( ifs, N_size[10] );
+		}
+		else if ( "Norbital-->"==word )
+		{
+			ModuleBase::GlobalFunc::READ_VALUE( ifs, N_size[11] );
+		}
 		else if ( "END"==word )
 		{
 			break;
-		}		
+		}
 	}
 
 	ModuleBase::CHECK_NAME(ifs, "Mesh");
@@ -169,19 +181,11 @@ std::vector<std::vector<Numerical_Orbital_Lm>> Exx_Abfs::IO::construct_abfs_T(
 	----------------------*/
 	for( size_t L=0; L<=L_size; ++L )
 		if( N_size.find(L) == N_size.end() )
-		{
-			std::stringstream ss;
-			ss<<"Can't find N of L="<<L<<" in "<<file_name;
-			throw std::domain_error(ss.str());
-		}
+			{ throw std::domain_error("Can't find N of L="+std::to_string(L)+" in "+file_name); }
 	for( size_t L=0; L<=L_size; ++L )
 		for( size_t N=0; N!=N_size[L]; ++N )
 			if( psis.find(L)==psis.end() || psis[L].find(N)==psis[L].end() )
-			{
-				std::stringstream ss;
-				ss<<"Can't find abf of L="<<L<<" T="<<T<<" in "<<file_name;
-				throw std::domain_error(ss.str());
-			}
+				{ throw std::domain_error("Can't find abf of L="+std::to_string(L)+" T="+std::to_string(T)+" in "+file_name); }
 
 
 	/*----------------------

source/source_lcao/module_ri/exx_abfs-jle.cpp

@@ -5,6 +5,7 @@
 #include "../../source_cell/unitcell.h"
 #include "../../source_base/mathzone.h"
 #include "../../source_base/math_sphbes.h" // mohan add 2021-05-06
+#include "source_base/tool_title.h"
 
 std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>
 Exx_Abfs::Jle::init_jle(
@@ -13,15 +14,19 @@ Exx_Abfs::Jle::init_jle(
 	const UnitCell& ucell,
 	const LCAO_Orbitals& orb)
 {
+	ModuleBase::TITLE("Exx_Abfs::Jle","init_jle");
 	std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> jle( ucell.ntype );
 
 	for(int T=0; T<ucell.ntype; ++T)
 	{
-		jle[T].resize( info.abfs_Lmax+1 );
+		if(info.abfs_Lmax+1>0)
+			{ jle[T].resize( info.abfs_Lmax+1 ); }
 		for(int L=0; L<=info.abfs_Lmax; ++L)
 		{
 			const size_t ecut_number 
 				= static_cast<size_t>( std::sqrt( info.ecut_exx ) * orb.Phi[T].getRcut() / ModuleBase::PI ); // Rydberg Unit.
+			if(ecut_number<=0)
+				{ continue; }
 
 			jle[T][L].resize( ecut_number );
 

source/source_lcao/module_ri/exx_opt_orb-print.cpp

@@ -1,18 +1,21 @@
 #include "exx_opt_orb.h"
 #include "../../source_pw/module_pwdft/global.h"
 #include "exx_abfs-jle.h"
+#include "source_base/tool_title.h"
 #include <iomanip>
 
 void Exx_Opt_Orb::print_matrix(
 	const Exx_Info::Exx_Info_Opt_ABFs &info,
 	const UnitCell& ucell,
-    const K_Vectors &kv,
-    const std::string& file_name,
+	const K_Vectors &kv,
+	const int Lmax,
+	const std::vector<std::size_t> &ecut_number,
+	const std::string &file_name,
 	const std::vector<RI::Tensor<double>> &matrix_Q, 
 	const std::vector<std::vector<RI::Tensor<double>>> &matrix_S,
 	const RI::Tensor<double> &matrix_V,
-	const size_t TA, const size_t IA, const size_t TB, const size_t IB,
-    const std::vector<double>& orb_cutoff,
+	const std::size_t TA, const std::size_t IA, const std::size_t TB, const std::size_t IB,
+	const std::vector<double>& orb_cutoff,
 	const ModuleBase::Element_Basis_Index::Range &range_jles, 
 	const ModuleBase::Element_Basis_Index::IndexLNM &index_jles) const
 {
@@ -61,45 +64,40 @@ void Exx_Opt_Orb::print_matrix(
 				<< ucell.atoms[TB].tau[IB].z << std::endl;
 		}
 
-		// ecutwfc_jlq determine the jlq corresponding to plane wave calculation.
-		ofs << info.ecut_exx << " ecutwfc" << std::endl; // mohan add 2009-09-08
+		ofs << info.ecut_exx << " ecutwfc" << std::endl;
 
 		// this parameter determine the total number of jlq.
-		ofs << info.ecut_exx << " ecutwfc_jlq" << std::endl;//mohan modify 2009-09-08
+		ofs << info.ecut_exx << " ecutwfc_jlq" << std::endl;
 
 		if(TA==TB)
 			{ ofs << orb_cutoff[TA] << " rcut_Jlq" << std::endl; }
 		else
 			{ ofs << orb_cutoff[TA] << " " << orb_cutoff[TB] << " rcut_Jlq" << std::endl; }
 
-		// mohan add 'smooth' and 'smearing_sigma' 2009-08-28
 		ofs << 0 << " smooth" << std::endl;
 		ofs << 0 << " smearing_sigma" << std::endl;
 
 		ofs << info.tolerence << " tolerence" << std::endl;
 
-		ofs << info.abfs_Lmax << " lmax" << std::endl;
+		ofs << Lmax << " lmax" << std::endl;
 
 		ofs << kv.get_nkstot() << " nks" << std::endl;
 		assert( matrix_V.shape[0]*matrix_V.shape[1] == matrix_V.shape[2]*matrix_V.shape[3] );
 		ofs	<< matrix_V.shape[0]*matrix_V.shape[1] << " nbands" << std::endl;
 
-		auto cal_sum_M = [&range_jles](size_t T) -> size_t
+		auto cal_sum_M = [&range_jles](std::size_t T) -> std::size_t
 		{
-			size_t sum_M = 0;
-			for( size_t L = 0; L!=range_jles[T].size(); ++L )
+			std::size_t sum_M = 0;
+			for( std::size_t L = 0; L!=range_jles[T].size(); ++L )
 				{ sum_M += range_jles[T][L].M; }
 			return sum_M;
 		};
-		const size_t nwfc = (TA==TB && IA==IB) ? cal_sum_M(TA) : cal_sum_M(TA)+cal_sum_M(TB);
+		const std::size_t nwfc = (TA==TB && IA==IB) ? cal_sum_M(TA) : cal_sum_M(TA)+cal_sum_M(TB);
 		ofs	<< nwfc << " nwfc" << std::endl;
 
-		const size_t ecut_numberA = static_cast<size_t>( std::sqrt( info.ecut_exx ) * orb_cutoff[TA] / ModuleBase::PI ); // Rydberg Unit
-		const size_t ecut_numberB = static_cast<size_t>( std::sqrt( info.ecut_exx ) * orb_cutoff[TB] / ModuleBase::PI ); // Rydberg Unit
-		if(TA==TB)
-			{ ofs << ecut_numberA << " ne" << std::endl; }
-		else
-			{ ofs << ecut_numberA << " " << ecut_numberB << " ne" << std::endl; }
+		for(const std::size_t ne : ecut_number)
+			{ ofs << ne << " "; }
+		ofs << "ne" << std::endl;
 
 		ofs << "<WEIGHT_OF_KPOINTS>" << std::endl;
 		for( int ik=0; ik!=kv.get_nkstot(); ++ik )		
@@ -119,18 +117,18 @@ void Exx_Opt_Orb::print_matrix(
 		//  < Psi | jY >
 		//---------------------
 		ofs<< "<OVERLAP_Q>" << std::endl;		
-		for( size_t iw0=0; iw0!=matrix_V.shape[0]; ++iw0 )
+		for( std::size_t iw0=0; iw0!=matrix_V.shape[0]; ++iw0 )
 		{
-			for( size_t iw1=0; iw1!=matrix_V.shape[1]; ++iw1 )
+			for( std::size_t iw1=0; iw1!=matrix_V.shape[1]; ++iw1 )
 			{
-				for( size_t iat=0; iat!=matrix_Q.size(); ++iat )
+				for( std::size_t iat=0; iat!=matrix_Q.size(); ++iat )
 				{
-					const size_t it = (iat==0) ? TA : TB;
-					for( size_t il=0; il!=range_jles[it].size(); ++il )
+					const std::size_t it = (iat==0) ? TA : TB;
+					for( std::size_t il=0; il!=range_jles[it].size(); ++il )
 					{
-						for( size_t im=0; im!=range_jles[it][il].M; ++im )
+						for( std::size_t im=0; im!=range_jles[it][il].M; ++im )
 						{
-							for( size_t iq=0; iq!=range_jles[it][il].N; ++iq )
+							for( std::size_t iq=0; iq!=range_jles[it][il].N; ++iq )
 							{
 								ofs<<matrix_Q[iat]( iw0, iw1, index_jles[it][il][iq][im] )<<"\t"<<0<<std::endl;
 							}
@@ -149,23 +147,23 @@ void Exx_Opt_Orb::print_matrix(
 		//  < jY | jY >
 		//---------------------
 		ofs<< "<OVERLAP_Sq>" <<std::endl;
-		for( size_t iat1=0; iat1!=matrix_S.size(); ++iat1 )
+		for( std::size_t iat1=0; iat1!=matrix_S.size(); ++iat1 )
 		{
-			const size_t it1 = (iat1==0) ? TA : TB;
-			for( size_t il1=0; il1!=range_jles[it1].size(); ++il1 )
+			const std::size_t it1 = (iat1==0) ? TA : TB;
+			for( std::size_t il1=0; il1!=range_jles[it1].size(); ++il1 )
 			{
-				for( size_t im1=0; im1!=range_jles[it1][il1].M; ++im1 )
+				for( std::size_t im1=0; im1!=range_jles[it1][il1].M; ++im1 )
 				{
-					for( size_t iat2=0; iat2!=matrix_S[iat1].size(); ++iat2 )
+					for( std::size_t iat2=0; iat2!=matrix_S[iat1].size(); ++iat2 )
 					{
-						const size_t it2 = (iat2==0) ? TA : TB;
-						for( size_t il2=0; il2!=range_jles[it2].size(); ++il2 )
+						const std::size_t it2 = (iat2==0) ? TA : TB;
+						for( std::size_t il2=0; il2!=range_jles[it2].size(); ++il2 )
 						{
-							for( size_t im2=0; im2!=range_jles[it2][il2].M; ++im2 )
+							for( std::size_t im2=0; im2!=range_jles[it2][il2].M; ++im2 )
 							{
-								for( size_t iq1=0; iq1!=range_jles[it1][il1].N; ++iq1 )
+								for( std::size_t iq1=0; iq1!=range_jles[it1][il1].N; ++iq1 )
 								{
-									for( size_t iq2=0; iq2!=range_jles[it2][il2].N; ++iq2 )
+									for( std::size_t iq2=0; iq2!=range_jles[it2][il2].N; ++iq2 )
 									{
 										ofs<<matrix_S[iat1][iat2]( index_jles[it1][il1][iq1][im1], index_jles[it2][il2][iq2][im2] )*scale<<"\t"<<0<<std::endl;
 									}
@@ -186,13 +184,13 @@ void Exx_Opt_Orb::print_matrix(
 		//  < Psi | Psi >
 		//---------------------	
 		ofs << "<OVERLAP_V>" << std::endl;
-		for( size_t iw0=0; iw0!=matrix_V.shape[0]; ++iw0 )
+		for( std::size_t iw0=0; iw0!=matrix_V.shape[0]; ++iw0 )
 		{
-			for( size_t iw1=0; iw1!=matrix_V.shape[1]; ++iw1 )
+			for( std::size_t iw1=0; iw1!=matrix_V.shape[1]; ++iw1 )
 			{
-				for( size_t iw2=0; iw2!=matrix_V.shape[2]; ++iw2 )
+				for( std::size_t iw2=0; iw2!=matrix_V.shape[2]; ++iw2 )
 				{
-					for( size_t iw3=0; iw3!=matrix_V.shape[3]; ++iw3 )
+					for( std::size_t iw3=0; iw3!=matrix_V.shape[3]; ++iw3 )
 					{
 						ofs<<matrix_V(iw0,iw1,iw2,iw3)*scale<<"\t";
 					}
@@ -203,6 +201,7 @@ void Exx_Opt_Orb::print_matrix(
 		ofs << "</OVERLAP_V>" << std::endl << std::endl;
 	};
 
+	ModuleBase::TITLE("Exx_Opt_Orb","print_matrix");
 	std::ofstream ofs(file_name+"_"+std::to_string(TA)+"_"+std::to_string(IA)+"_"+std::to_string(TB)+"_"+std::to_string(IB));
 	print_header(ofs);
 	ofs<<std::setprecision(15);