Merge pull request #235 from PeizeLin/dev_LibRI_2

1. add Exx_LRI::cal_exx_stress()

source/module_ri/Exx_LRI.h

@@ -26,7 +26,7 @@ class Exx_LRI
 private:
 	using TA = int;
 	using Tcell = int;
-	static const size_t Ndim = 3;
+	static constexpr std::size_t Ndim = 3;
 	using TC = std::array<Tcell,Ndim>;
 	using TAC = std::pair<TA,TC>;
 	using TatomR = std::array<double,Ndim>;		// tmp
@@ -38,10 +38,12 @@ class Exx_LRI
 	void cal_exx_ions();
 	void cal_exx_elec(const Local_Orbital_Charge &loc, const Parallel_Orbitals &pv);
 	void cal_exx_force();
+	void cal_exx_stress();
 
 	std::vector< std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>> Hexxs;
 	Tdata Eexx;
-	ModuleBase::matrix Fexx;
+	ModuleBase::matrix force_exx;
+	ModuleBase::matrix stress_exx;
 
 	void write_Hexxs(const std::string &file_name) const;
 	void read_Hexxs(const std::string &file_name);

source/module_ri/Exx_LRI.hpp

@@ -239,21 +239,46 @@ void Exx_LRI<Tdata>::cal_exx_force()
 
 	this->exx_lri.set_csm_threshold(this->info.cauchy_grad_threshold);
 
-	this->Fexx.create(GlobalC::ucell.nat, Ndim);
+	this->force_exx.create(GlobalC::ucell.nat, Ndim);
 	for(int is=0; is<GlobalV::NSPIN; ++is)
 	{
-		this->exx_lri.cal_Fs({"","",std::to_string(is),"",""});
-		for(std::size_t ix=0; ix<Ndim; ++ix)
-			for(const auto &force_item : this->exx_lri.force[ix])
-				this->Fexx(force_item.first, ix) += std::real(force_item.second);
+		this->exx_lri.cal_force({"","",std::to_string(is),"",""});
+		for(std::size_t idim=0; idim<Ndim; ++idim)
+			for(const auto &force_item : this->exx_lri.force[idim])
+				this->force_exx(force_item.first, idim) += std::real(force_item.second);
 	}
 
 	const double SPIN_multiple = std::map<int,double>{{1,2}, {2,1}, {4,1}}.at(GlobalV::NSPIN);				// why?
 	const double frac = -2 * SPIN_multiple;		// why?
-	this->Fexx *= frac;
+	this->force_exx *= frac;
 	ModuleBase::timer::tick("Exx_LRI", "cal_exx_force");
 }
 
+
+template<typename Tdata>
+void Exx_LRI<Tdata>::cal_exx_stress()
+{
+	ModuleBase::TITLE("Exx_LRI","cal_exx_stress");
+	ModuleBase::timer::tick("Exx_LRI", "cal_exx_stress");
+
+	this->exx_lri.set_csm_threshold(this->info.cauchy_grad_threshold);
+
+	this->stress_exx.create(Ndim, Ndim);
+	for(int is=0; is<GlobalV::NSPIN; ++is)
+	{
+		this->exx_lri.cal_stress({"","",std::to_string(is),"",""});
+		for(std::size_t idim0=0; idim0<Ndim; ++idim0)
+			for(std::size_t idim1=0; idim1<Ndim; ++idim1)
+				this->stress_exx(idim0,idim1) += std::real(this->exx_lri.stress(idim0,idim1));
+	}
+
+	const double SPIN_multiple = std::map<int,double>{{1,2}, {2,1}, {4,1}}.at(GlobalV::NSPIN);				// why?
+	const double frac = 2 * SPIN_multiple / GlobalC::ucell.omega * GlobalC::ucell.lat0;		// why?
+	this->stress_exx *= frac;
+
+	ModuleBase::timer::tick("Exx_LRI", "cal_exx_stress");
+}
+
 template<typename Tdata>
 void Exx_LRI<Tdata>::write_Hexxs(const std::string &file_name) const
 {

source/module_xc/xc_functional_vxc.cpp

@@ -180,8 +180,8 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional::v_xc_libxc(
         const double * const * const rho_in,
         const double * const rho_core_in)
 {
-    ModuleBase::TITLE("XC_Functional","v_xc");
-    ModuleBase::timer::tick("XC_Functional","v_xc");
+    ModuleBase::TITLE("XC_Functional","v_xc_libxc");
+    ModuleBase::timer::tick("XC_Functional","v_xc_libxc");
 
     const int nspin = GlobalV::NSPIN;
 
@@ -385,7 +385,7 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional::v_xc_libxc(
 
     finish_func(funcs);
 
-    ModuleBase::timer::tick("XC_Functional","v_xc");
+    ModuleBase::timer::tick("XC_Functional","v_xc_libxc");
     return std::make_tuple( etxc, vtxc, std::move(v) );
 }
 
@@ -409,8 +409,8 @@ tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional::v_xc_m
 	const double * const rho_core_in,
 	const double * const * const kin_r_in)
 {
-    ModuleBase::TITLE("XC_Functional","v_xc");
-    ModuleBase::timer::tick("XC_Functional","v_xc");
+    ModuleBase::TITLE("XC_Functional","v_xc_meta");
+    ModuleBase::timer::tick("XC_Functional","v_xc_meta");
 
     if(GlobalV::NSPIN==4)
     {

source/src_lcao/FORCE_STRESS.cpp

@@ -233,20 +233,37 @@ void Force_Stress_LCAO::getForceStress(
 			}
 		}
 	}
-	//Force contribution from exx
 #ifdef __EXX
-	ModuleBase::matrix fexx;
+	//Force and Stress contribution from exx
+	ModuleBase::matrix force_exx;
+	ModuleBase::matrix stress_exx;
 	if( GlobalC::exx_info.info_global.cal_exx )
 	{
-		if(GlobalV::GAMMA_ONLY_LOCAL)
+		if(isforce)
 		{
-			GlobalC::exx_lri_double.cal_exx_force();
-			fexx = GlobalC::exx_info.info_global.hybrid_alpha * GlobalC::exx_lri_double.Fexx;
+			if(GlobalV::GAMMA_ONLY_LOCAL)
+			{
+				GlobalC::exx_lri_double.cal_exx_force();
+				force_exx = GlobalC::exx_info.info_global.hybrid_alpha * GlobalC::exx_lri_double.force_exx;
+			}
+			else
+			{
+				GlobalC::exx_lri_complex.cal_exx_force();
+				force_exx = GlobalC::exx_info.info_global.hybrid_alpha * GlobalC::exx_lri_complex.force_exx;
+			}
 		}
-		else
+		if(isstress)
 		{
-			GlobalC::exx_lri_complex.cal_exx_force();
-			fexx = GlobalC::exx_info.info_global.hybrid_alpha * GlobalC::exx_lri_complex.Fexx;
+			if(GlobalV::GAMMA_ONLY_LOCAL)
+			{
+				GlobalC::exx_lri_double.cal_exx_stress();
+				stress_exx = GlobalC::exx_info.info_global.hybrid_alpha * GlobalC::exx_lri_double.stress_exx;
+			}
+			else
+			{
+				GlobalC::exx_lri_complex.cal_exx_stress();
+				stress_exx = GlobalC::exx_info.info_global.hybrid_alpha * GlobalC::exx_lri_complex.stress_exx;
+			}
 		}
 	}
 #endif
@@ -282,7 +299,7 @@ void Force_Stress_LCAO::getForceStress(
 				// Force contribution from exx
 				if( GlobalC::exx_info.info_global.cal_exx )
 				{
-					fcs(iat,i) += fexx(iat,i);
+					fcs(iat,i) += force_exx(iat,i);
 				}
 #endif
 				//VDW force of vdwd2 or vdwd3
@@ -517,6 +534,13 @@ void Force_Stress_LCAO::getForceStress(
 				{
 					scs(i, j) += stress_dftu(i, j);
 				}
+#ifdef __EXX
+				// Stress contribution from exx
+				if( GlobalC::exx_info.info_global.cal_exx )
+				{
+					scs(i,j) += stress_exx(i,j);
+				}
+#endif
 			}
 		}
 
@@ -763,6 +787,7 @@ void Force_Stress_LCAO::calForcePwPart(
 	ModuleBase::matrix &fcc,
 	ModuleBase::matrix &fscc)
 {
+	ModuleBase::TITLE("Force_Stress_LCAO","calForcePwPart");
 	//--------------------------------------------------------
 	// local pseudopotential force:
 	// use charge density; plane wave; local pseudopotential;
@@ -863,6 +888,7 @@ void Force_Stress_LCAO::calStressPwPart(
 	ModuleBase::matrix& sigmaxc
 )
 {
+    ModuleBase::TITLE("Force_Stress_LCAO","calStressPwPart");
 	//--------------------------------------------------------
 	// local pseudopotential stress:
 	// use charge density; plane wave; local pseudopotential;

source/src_pw/forces.cpp

@@ -430,6 +430,7 @@ void Forces::print(const std::string& name, const ModuleBase::matrix& f, bool ry
 
 void Forces::cal_force_loc(ModuleBase::matrix& forcelc, ModulePW::PW_Basis* rho_basis)
 {
+    ModuleBase::TITLE("Forces", "cal_force_loc");
     ModuleBase::timer::tick("Forces", "cal_force_loc");
 
     std::complex<double>* aux = new std::complex<double>[rho_basis->nmaxgr];
@@ -480,6 +481,7 @@ void Forces::cal_force_loc(ModuleBase::matrix& forcelc, ModulePW::PW_Basis* rho_
 #include "H_Ewald_pw.h"
 void Forces::cal_force_ew(ModuleBase::matrix& forceion, ModulePW::PW_Basis* rho_basis)
 {
+    ModuleBase::TITLE("Forces", "cal_force_ew");
     ModuleBase::timer::tick("Forces", "cal_force_ew");
 
     double fact = 2.0;
@@ -643,6 +645,7 @@ void Forces::cal_force_ew(ModuleBase::matrix& forceion, ModulePW::PW_Basis* rho_
 
 void Forces::cal_force_cc(ModuleBase::matrix& forcecc, ModulePW::PW_Basis* rho_basis)
 {
+    ModuleBase::TITLE("Forces", "cal_force_cc");
     // recalculate the exchange-correlation potential.
 
     ModuleBase::matrix v(GlobalV::NSPIN, rho_basis->nrxx);
@@ -924,6 +927,7 @@ void Forces::cal_force_nl(ModuleBase::matrix& forcenl, const psi::Psi<complex<do
 
 void Forces::cal_force_scc(ModuleBase::matrix& forcescc, ModulePW::PW_Basis* rho_basis)
 {
+    ModuleBase::TITLE("Forces", "cal_force_scc");
     std::complex<double>* psic = new std::complex<double>[rho_basis->nmaxgr];
 
     if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4)

source/src_pw/stress_func_cc.cpp

@@ -7,6 +7,7 @@
 //NLCC term, need to be tested
 void Stress_Func::stress_cc(ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, const bool is_pw)
 {
+    ModuleBase::TITLE("Stress_Func","stress_cc");
 	ModuleBase::timer::tick("Stress_Func","stress_cc");
 
 	double fact=1.0;

source/src_pw/stress_func_ewa.cpp

@@ -6,7 +6,8 @@
 //calcualte the Ewald stress term in PW and LCAO
 void Stress_Func::stress_ewa(ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, const bool is_pw)
 {
-    ModuleBase::timer::tick("Stress_Func","stress_ew");
+    ModuleBase::TITLE("Stress_Func","stress_ewa");
+    ModuleBase::timer::tick("Stress_Func","stress_ewa");
 
     double charge=0;
     for(int it=0; it < GlobalC::ucell.ntype; it++)
@@ -152,7 +153,7 @@ void Stress_Func::stress_ewa(ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_
 	delete[] r2;
 	delete[] irr;
 	// this->print(GlobalV::ofs_running, "ewald stress", stression);
-	ModuleBase::timer::tick("Stress_Func","stress_ew");
+	ModuleBase::timer::tick("Stress_Func","stress_ewa");
 
 	return;
 }

source/src_pw/stress_func_gga.cpp

@@ -6,6 +6,7 @@
 //calculate the GGA stress correction in PW and LCAO
 void Stress_Func::stress_gga(ModuleBase::matrix& sigma) 
 {
+    ModuleBase::TITLE("Stress_Func","stress_gga");
 	ModuleBase::timer::tick("Stress_Func","stress_gga");
 
 	int func_type = XC_Functional::get_func_type();

source/src_pw/stress_func_har.cpp

@@ -7,6 +7,7 @@
 //calculate the Hartree part in PW or LCAO base
 void Stress_Func::stress_har(ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, const bool is_pw)
 {
+    ModuleBase::TITLE("Stress_Func","stress_har");
 	ModuleBase::timer::tick("Stress_Func","stress_har");
 	double shart;
 

source/src_pw/stress_func_kin.cpp

@@ -4,6 +4,9 @@
 //calculate the kinetic stress in PW base
 void Stress_Func::stress_kin(ModuleBase::matrix& sigma, const psi::Psi<complex<double>>* psi_in)
 {
+    ModuleBase::TITLE("Stress_Func","stress_kin");
+	ModuleBase::timer::tick("Stress_Func","stress_kin");
+
 	double **gk;
 	gk=new double* [3];
 	int npw;
@@ -131,5 +134,6 @@ void Stress_Func::stress_kin(ModuleBase::matrix& sigma, const psi::Psi<complex<d
 	delete[] gk[2];
 	delete[] gk;
 
+	ModuleBase::timer::tick("Stress_Func","stress_kin");
 	return;
 }

source/src_pw/stress_func_loc.cpp

@@ -6,6 +6,7 @@
 //calculate local pseudopotential stress in PW or VL_dVL stress in LCAO
 void Stress_Func::stress_loc(ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, const bool is_pw)
 {
+    ModuleBase::TITLE("Stress_Func","stress_loc");
     ModuleBase::timer::tick("Stress_Func","stress_loc");
 
     double *dvloc = new double[rho_basis->npw];

source/src_pw/stress_func_print.cpp

@@ -65,7 +65,6 @@ void Stress_Func::printstress_total(const ModuleBase::matrix& scs, bool ry)
 		unit_transform = ModuleBase::RYDBERG_SI / pow(ModuleBase::BOHR_RADIUS_SI,3) * 1.0e-8;
 	}
 //	std::cout.setf(ios::fixed);
-
 
 	//GlobalV::ofs_running << std::setiosflags(ios::right);
  	GlobalV::ofs_running << std::setprecision(6) << std::setiosflags(ios::showpos) << std::setiosflags(ios::fixed) << std::endl;