XC functional refactor

CMakeLists.txt

@@ -260,7 +260,6 @@ target_link_libraries(${ABACUS_BIN_NAME}
     cell
     symmetry
     md
-    symmetry
     neighbor
     orb
     io
@@ -272,6 +271,7 @@ target_link_libraries(${ABACUS_BIN_NAME}
     pw
     ri
     driver
+    xc
     -lm
 )
 

source/CMakeLists.txt

@@ -5,6 +5,7 @@ add_subdirectory(module_neighbor)
 add_subdirectory(module_orbital)
 add_subdirectory(module_md)
 add_subdirectory(module_deepks)
+add_subdirectory(module_xc)
 add_subdirectory(src_io)
 add_subdirectory(src_ions)
 add_subdirectory(src_lcao)

source/Makefile

@@ -15,6 +15,7 @@ VPATH=./src_global\
 :./module_base\
 :./module_md\
 :./module_deepks\
+:./module_xc\
 :./src_pw\
 :./src_lcao\
 :./src_ions\

source/Makefile.Objects

@@ -19,8 +19,7 @@ input_conv.o\
 run_pw.o\
 run_lcao.o\
 
-OBJS_PW=xc_type.o \
-xc_functional.o\
+OBJS_PW=xc_functional.o\
 xc_3.o \
 vdwd2.o\
 vdwd2_parameters.o\
@@ -265,7 +264,6 @@ write_HS.o\
 write_HS_R.o\
 write_dm.o\
 write_wfc_realspace.o\
-potential_libxc.o \
 potential_libxc_meta.o \
 efield.o \
 magnetism.o\

source/input.cpp

@@ -144,7 +144,7 @@ void Input::Default(void)
 //----------------------------------------------------------
 // electrons / spin
 //----------------------------------------------------------
-	dft_functional = "none";
+	dft_functional = "default";
     nspin = 1;
     nelec = 0.0;
     lmaxmax = 2;
@@ -214,7 +214,6 @@ void Input::Default(void)
 	vl_in_h = 1;
 	vnl_in_h = 1;
 	vh_in_h = 1;
-	vxc_in_h = 1;
 	vion_in_h = 1;
 	test_force = 0;
 	test_stress = 0;
@@ -355,7 +354,6 @@ void Input::Default(void)
 //----------------------------------------------------------
 // exx										//Peize Lin add 2018-06-20
 //----------------------------------------------------------
-	exx_hybrid_type = "no";
 
 	exx_hybrid_alpha = 0.25;
 	exx_hse_omega = 0.11;
@@ -892,10 +890,6 @@ bool Input::Read(const std::string &fn)
         {
             read_value(ifs, vh_in_h);
         }
-        else if (strcmp("vxc_in_h", word) == 0)
-        {
-            read_value(ifs, vxc_in_h);
-        }
         else if (strcmp("vion_in_h", word) == 0)
         {
             read_value(ifs, vion_in_h);
@@ -1490,9 +1484,9 @@ bool Input::Read(const std::string &fn)
 // exx
 // Peize Lin add 2018-06-20
 //----------------------------------------------------------
-        else if (strcmp("exx_hybrid_type", word) == 0)
+        else if (strcmp("dft_functional", word) == 0)
         {
-            read_value(ifs, exx_hybrid_type);
+            read_value(ifs, dft_functional);
         }
         else if (strcmp("exx_hybrid_alpha", word) == 0)
         {
@@ -2059,7 +2053,6 @@ void Input::Bcast()
 	Parallel_Common::bcast_int( vl_in_h );
 	Parallel_Common::bcast_int( vnl_in_h );
 	Parallel_Common::bcast_int( vh_in_h );
-	Parallel_Common::bcast_int( vxc_in_h );
 	Parallel_Common::bcast_int( vion_in_h );
 
 	Parallel_Common::bcast_int( test_force );
@@ -2258,7 +2251,7 @@ void Input::Bcast()
 	Parallel_Common::bcast_double( eps_degauss );
 
 	// Peize Lin add 2018-06-20
-	Parallel_Common::bcast_string( exx_hybrid_type );
+	Parallel_Common::bcast_string( dft_functional );
 	Parallel_Common::bcast_double( exx_hybrid_alpha );
 	Parallel_Common::bcast_double( exx_hse_omega );
 	Parallel_Common::bcast_bool( exx_separate_loop );
@@ -2862,16 +2855,7 @@ void Input::Check(void)
 		}
 	}
 
-	if(exx_hybrid_type!="no" &&
-		exx_hybrid_type!="hf" &&
-		exx_hybrid_type!="pbe0" &&
-		exx_hybrid_type!="hse" &&
-		exx_hybrid_type!="opt_orb")
-	{
-		ModuleBase::WARNING_QUIT("INPUT","exx_hybrid_type must be no or hf or pbe0 or hse or opt_orb");
-	}
-
-	if(exx_hybrid_type=="hf" || exx_hybrid_type=="pbe0" || exx_hybrid_type=="hse")
+	if(dft_functional=="hf" || dft_functional=="pbe0" || dft_functional=="hse")
 	{
 		if(exx_hybrid_alpha<0 || exx_hybrid_alpha>1)
 		{
@@ -2893,7 +2877,7 @@ void Input::Check(void)
 			ModuleBase::WARNING_QUIT("INPUT","exx_distribute_type must be htime or kmeans2 or kmeans1");
 		}
 	}
-	if(exx_hybrid_type=="opt_orb")
+	if(dft_functional=="opt_orb")
 	{
 		if(exx_opt_orb_lmax<0)
 		{

source/input.h

@@ -92,6 +92,7 @@ class Input
 // electrons / spin
 //==========================================================
     std::string dft_functional;	// input DFT functional.
+	bool use_libxc;         // whether to use LIBXC
 	int nspin;				// LDA ; LSDA ; non-linear spin
     double nelec;			// total number of electrons
     int lmaxmax;
@@ -166,7 +167,6 @@ class Input
 	int vnl_in_h;			// calculate the vnl or not.
 
 	int vh_in_h;			// calculate the hartree potential or not
-	int vxc_in_h;			// calculate the xc potential or not
 	int vion_in_h;			// calculate the local ionic potential or not
 	//only relevant when vl_in_h = 1
 
@@ -360,8 +360,6 @@ class Input
 // exx
 // Peize Lin add 2018-06-20
 //==========================================================
-	std::string exx_hybrid_type;		// "no", "hf", "pbe0", "hse"
-
 	double exx_hybrid_alpha;
 	double exx_hse_omega;
 

source/input_conv.cpp

@@ -153,7 +153,6 @@ void Input_Conv::Convert(void)
 	GlobalV::VL_IN_H = INPUT.vl_in_h;
 	GlobalV::VNL_IN_H = INPUT.vnl_in_h;
 	GlobalV::VH_IN_H = INPUT.vh_in_h;
-	GlobalV::VXC_IN_H = INPUT.vxc_in_h;
 	GlobalV::VION_IN_H = INPUT.vion_in_h;
 	GlobalV::TEST_FORCE = INPUT.test_force;
 	GlobalV::TEST_STRESS = INPUT.test_stress;
@@ -419,20 +418,20 @@ void Input_Conv::Convert(void)
 		const std::string command0 = "test -d " + GlobalC::restart.folder + " || mkdir " + GlobalC::restart.folder;
 		if (GlobalV::MY_RANK == 0)
 			system(command0.c_str());
-		if (INPUT.exx_hybrid_type == "no")
+		if (INPUT.dft_functional == "hf" || INPUT.dft_functional == "pbe0" || INPUT.dft_functional == "hse" || INPUT.dft_functional == "opt_orb")
 		{
 			GlobalC::restart.info_save.save_charge = true;
+			GlobalC::restart.info_save.save_H = true;
 		}
 		else
 		{
 			GlobalC::restart.info_save.save_charge = true;
-			GlobalC::restart.info_save.save_H = true;
 		}
 	}
 	if (INPUT.restart_load)
 	{
 		GlobalC::restart.folder = GlobalV::global_out_dir + "restart/";
-		if (INPUT.exx_hybrid_type == "no")
+		if (INPUT.dft_functional == "hf" || INPUT.dft_functional == "pbe0" || INPUT.dft_functional == "hse" || INPUT.dft_functional == "opt_orb")
 		{
 			GlobalC::restart.info_load.load_charge = true;
 		}
@@ -447,28 +446,30 @@ void Input_Conv::Convert(void)
 // about exx, Peize Lin add 2018-06-20
 //----------------------------------------------------------
 #ifdef __LCAO
-	if (INPUT.exx_hybrid_type == "no")
+
+	if (INPUT.dft_functional == "hf")
 	{
-		GlobalC::exx_global.info.hybrid_type = Exx_Global::Hybrid_Type::No;
+		GlobalC::exx_global.info.hybrid_type = Exx_Global::Hybrid_Type::HF;
+	}
+	else if (INPUT.dft_functional == "pbe0")
+	{
+		GlobalC::exx_global.info.hybrid_type = Exx_Global::Hybrid_Type::PBE0;
+	}
+	else if (INPUT.dft_functional == "hse")
+	{
+		GlobalC::exx_global.info.hybrid_type = Exx_Global::Hybrid_Type::HSE;
+	}
+	else if (INPUT.dft_functional == "opt_orb")
+	{
+		GlobalC::exx_global.info.hybrid_type = Exx_Global::Hybrid_Type::Generate_Matrix;
 	}
 	else
 	{
-		if (INPUT.exx_hybrid_type == "hf")
-		{
-			GlobalC::exx_global.info.hybrid_type = Exx_Global::Hybrid_Type::HF;
-		}
-		else if (INPUT.exx_hybrid_type == "pbe0")
-		{
-			GlobalC::exx_global.info.hybrid_type = Exx_Global::Hybrid_Type::PBE0;
-		}
-		else if (INPUT.exx_hybrid_type == "hse")
-		{
-			GlobalC::exx_global.info.hybrid_type = Exx_Global::Hybrid_Type::HSE;
-		}
-		else if (INPUT.exx_hybrid_type == "opt_orb")
-		{
-			GlobalC::exx_global.info.hybrid_type = Exx_Global::Hybrid_Type::Generate_Matrix;
-		}
+		GlobalC::exx_global.info.hybrid_type = Exx_Global::Hybrid_Type::No;
+	}
+
+	if(GlobalC::exx_global.info.hybrid_type != Exx_Global::Hybrid_Type::No)
+	{
 		GlobalC::exx_global.info.hybrid_alpha = INPUT.exx_hybrid_alpha;
 		GlobalC::exx_global.info.hse_omega = INPUT.exx_hse_omega;
 		GlobalC::exx_global.info.separate_loop = INPUT.exx_separate_loop;
@@ -569,16 +570,6 @@ void Input_Conv::Convert(void)
 	// mohan add 2021-02-16
 	berryphase::berry_phase_flag = INPUT.berry_phase;
 
-	// wenfei 2021-7-28
-	if (GlobalV::DFT_FUNCTIONAL == "scan")
-	{
-		if (GlobalV::BASIS_TYPE != "pw")
-		{
-			ModuleBase::WARNING_QUIT("Input_conv", "add metaGGA for pw first");
-		}
-		GlobalV::DFT_META = 1;
-	}
-
 	ModuleBase::timer::tick("Input_Conv", "Convert");
 //-----------------------------------------------
 // caoyu add for DeePKS

source/module_base/global_variable.cpp

@@ -29,8 +29,7 @@ std::string	CALCULATION = "scf";
 int		EFIELD = 0; // 5: add electric field
 int		DIPOLE = 0; // 7: add dipole field
 
-std::string  DFT_FUNCTIONAL = "none";
-bool	DFT_META = 0;
+std::string  DFT_FUNCTIONAL = "default";
 int 	NSPIN = 1; // LDA
 bool	TWO_EFERMI = 0; // two fermi energy, exist only magnetization is fixed.
 int 	CURRENT_SPIN = 0;
@@ -76,7 +75,6 @@ int T_IN_H = 1; // mohan add 2010-11-28
 int VL_IN_H = 1;
 int VNL_IN_H = 1;
 int VH_IN_H = 1;
-int VXC_IN_H = 1;
 int VION_IN_H = 1;
 int ZEEMAN_IN_H = 1;
 double  STRESS_THR = 1.0e-2; //LiuXh add 20180515

source/module_base/global_variable.h

@@ -31,7 +31,6 @@ extern int 		DIPOLE;				// 7 add dipole correction
 
 
 extern std::string   DFT_FUNCTIONAL;		// 6.5 change the DFT functional from input file.
-extern bool		DFT_META;		// whether is meta-GGA
 extern int 		NSPIN;				// 7
 extern bool		TWO_EFERMI; 		// 7.5 mohan add 2011-04-03, two fermi energy, exist if magnetization is fixed.
 extern int 		CURRENT_SPIN;		// 8
@@ -89,7 +88,6 @@ extern int T_IN_H; // 23, calculate T in H or not.
 extern int VL_IN_H; // 24, calculate Vl in H or not.
 extern int VNL_IN_H; // 25, calculate Vnl in H or not.
 extern int VH_IN_H; // 26, calculate Vh in H or not.
-extern int VXC_IN_H; // 27, calculate Vxc in H or not.
 extern int VION_IN_H; // 28, calculate Vion_loc in H or not.
 extern double   STRESS_THR; //LiuXh add 20180515
 

source/module_cell/atom_pseudo.cpp

@@ -197,7 +197,7 @@ void Atom_pseudo::bcast_atom_pseudo2(void)
 //std::string
 	Parallel_Common::bcast_string( psd );
 	Parallel_Common::bcast_string( pp_type );
-	Parallel_Common::bcast_string( dft, 4 );
+	Parallel_Common::bcast_string( xc_func );
 
 	if(GlobalV::MY_RANK!=0)
 	{

source/module_cell/pseudo_nc.cpp

@@ -125,10 +125,7 @@ void pseudo_nc::set_pseudo_h(const Pseudopot_upf &upf)
 	this->tvanp = upf.tvanp;// if USPP
 	this->nlcc = upf.nlcc;// Non linear core corrections( bool ?)
 
-	for(int i=0; i<4; i++)
-	{
-		this->dft[i] = upf.dft[i];
-	}
+	this->xc_func = upf.xc_func;
 
 	this->zv = upf.zp;
 	this->etotps = upf.etotps;
@@ -377,7 +374,7 @@ void pseudo_nc::print_pseudo_h(std::ofstream &ofs)
     ofs << "\n pp_type  " << pp_type;
     ofs << "\n tvanp    " << tvanp;
     ofs << "\n nlcc " << nlcc;
-    ofs << "\n dft  " << dft;
+    ofs << "\n dft  " << xc_func;
     ofs << "\n zv       " << zv;
     ofs << "\n etotps   " << etotps;
     ofs << "\n ecutwfc " << ecutwfc;

source/module_cell/pseudo_nc.h

@@ -22,7 +22,7 @@ class pseudo_nc
 	std::string pp_type;		// Pseudo type ( NC or US )
 	bool tvanp;			// .true. if Ultrasoft
 	bool nlcc;			// Non linear core corrections(bool)
-	std::string dft[4];		// Exch-Corr type
+	std::string xc_func;		// Exch-Corr type
 	int  zv;			// z valence
 	double etotps;		// total energy
 	double ecutwfc;		// suggested cut-off for wfc

source/module_cell/read_cell_pseudopots.cpp

@@ -82,10 +82,7 @@ void UnitCell_pseudo::read_cell_pseudopots(const std::string &pp_dir, std::ofstr
 
 			log << "\n Read in pseudopotential file is " << pseudo_fn[i] << std::endl;
 			ModuleBase::GlobalFunc::OUT(log,"pseudopotential type",atoms[i].pp_type);
-			ModuleBase::GlobalFunc::OUT(log,"functional Ex", atoms[i].dft[0]);
-			ModuleBase::GlobalFunc::OUT(log,"functional Ec", atoms[i].dft[1]);
-			ModuleBase::GlobalFunc::OUT(log,"functional GCEx", atoms[i].dft[2]);
-			ModuleBase::GlobalFunc::OUT(log,"functional GCEc", atoms[i].dft[3]);
+			ModuleBase::GlobalFunc::OUT(log,"exchange-correlation functional", atoms[i].xc_func);
 			ModuleBase::GlobalFunc::OUT(log,"nonlocal core correction", atoms[i].nlcc);
 //			ModuleBase::GlobalFunc::OUT(log,"spin orbital",atoms[i].has_so);
 			ModuleBase::GlobalFunc::OUT(log,"valence electrons", atoms[i].zv);
@@ -100,20 +97,16 @@ void UnitCell_pseudo::read_cell_pseudopots(const std::string &pp_dir, std::ofstr
 		}
 		if(upf.functional_error == 1)
 		{
-			std::cout << "In Pseudopot_upf::read_pseudo_header : input xc functional does not match that in pseudopot file" << std::endl;
+			std::cout << "In Pseudopot_upf::read_pseudo_header : dft_functional from INPUT does not match that in pseudopot file" << std::endl;
 			std::cout << "Please make sure this is what you need" << std::endl;
-			atoms[i].dft[0] = GlobalV::DFT_FUNCTIONAL;
-			transform(atoms[i].dft[0].begin(), atoms[i].dft[0].end(), atoms[i].dft[0].begin(), (::toupper));
-			atoms[i].dft[1].clear();
-			atoms[i].dft[2].clear();
-			atoms[i].dft[3].clear();
+			atoms[i].xc_func = GlobalV::DFT_FUNCTIONAL;
+			transform(atoms[i].xc_func.begin(), atoms[i].xc_func.end(), atoms[i].xc_func.begin(), (::toupper));
 			if(GlobalV::MY_RANK==0)
 			{
-				log << "\n XC functional updated to : " << std::endl;
-				ModuleBase::GlobalFunc::OUT(log,"functional Ex", atoms[i].dft[0]);
-				ModuleBase::GlobalFunc::OUT(log,"functional Ec", atoms[i].dft[1]);
-				ModuleBase::GlobalFunc::OUT(log,"functional GCEx", atoms[i].dft[2]);
-				ModuleBase::GlobalFunc::OUT(log,"functional GCEc", atoms[i].dft[3]);
+				log << "\n In Pseudopot_upf::read_pseudo_header : dft_functional from INPUT does not match that in pseudopot file" << std::endl;
+				log << " Please make sure this is what you need" << std::endl;
+				log << " XC functional updated to : " << GlobalV::DFT_FUNCTIONAL << std::endl;
+				ModuleBase::GlobalFunc::OUT(log,"exchange-correlation functional", atoms[i].xc_func);
 			}
 		}
 

source/module_cell/read_pp.h

@@ -32,7 +32,7 @@ class Pseudopot_upf
 	std::string pp_type;		// header_3 // Pseudo type ( NC or US )
 	bool tvanp;         // header_4 // .true. if Ultrasoft 
 	bool nlcc;		    // header_5 // Non linear core corrections
-	std::string dft[4];		// header_6 // Exch-Corr type
+	std::string xc_func;		// header_6 // Exch-Corr type
 	int  zp;            // header_7 // z valence 
 	double etotps;		// header_8 // total energy
 	double ecutwfc;		// header_9 // suggested cut-off for wfc