Feature(xc): implement gga_grad=3 via SF builtin for noncollinear GGA

Call SF builtin (v_xc_ncgga_sf_builtin) from v_xc for gga_grad=3
noncollinear calculations, which computes LDA and GGA in a single
self-consistent function with per-mu independent div correction.

Use gga_grad=2 per-component gradient decomposition and projected div
in gradcorr for gga_grad=3 as fallback when not using SF builtin.

Scope libxc SF path to gga_grad=2 only.

Author: dyzheng

source/source_hamilt/module_xc/xc_functional_gradcorr.cpp

@@ -12,6 +12,9 @@
 #include "source_base/timer.h"
 #include "source_basis/module_pw/pw_basis_k.h"
 #include "source_io/module_parameter/parameter.h"
+
+#include "xc_functional_gga_noncol_sf_builtin.h"
+
 #include <ATen/core/tensor.h>
 #include <ATen/core/tensor_map.h>
 #include <ATen/core/tensor_types.h>
@@ -21,7 +24,7 @@
 #include "xc_functional_libxc.h"
 #endif
 
-#include "xc_functional_gga_noncol_sf_builtin.h"
+
 
 void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 	const Charge* const chr, ModulePW::PW_Basis* rhopw, const UnitCell *ucell,
@@ -187,7 +190,7 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 		gdr2 = new ModuleBase::Vector3<double>[rhopw->nrxx];
 		if(!is_stress) h2 = new ModuleBase::Vector3<double>[rhopw->nrxx];
 
-		if (gga_grad == 2)
+		if (gga_grad == 2 || gga_grad == 3)
 		{
 			tmp_recip = new std::complex<double>[rhopw->npw];
 			gdr_mag = new ModuleBase::Vector3<double>[rhopw->nrxx];
@@ -491,7 +494,7 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 }
 #endif
 
-	if(!is_stress && PARAM.inp.nspin==4 && gga_grad==2)
+	if(!is_stress && PARAM.inp.nspin==4 && (PARAM.globalv.domag||PARAM.globalv.domag_z) && (gga_grad==2 || gga_grad==3))
 	{
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 1024)

source/source_hamilt/module_xc/xc_functional_libxc_tools.cpp

@@ -235,7 +235,7 @@ std::pair<double,ModuleBase::matrix> XC_Functional_Libxc::convert_vtxc_v(
 
 	if(func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA)
 	{
-		if(PARAM.inp.nspin==4 && PARAM.inp.gga_grad >= 2 && (PARAM.globalv.domag || PARAM.globalv.domag_z))
+		if(PARAM.inp.nspin==4 && PARAM.inp.gga_grad == 2 && (PARAM.globalv.domag || PARAM.globalv.domag_z))
 		{
 			std::vector<double> mag_part_tmp = XC_Functional_Libxc::compute_mag_part_nspin4(nrxx, chr);
 			const std::vector<std::vector<double>> dh = XC_Functional_Libxc::cal_dh_sf(nspin, nrxx, sgn, gdr, vsigma, mag_part_tmp, tpiba, chr);

source/source_hamilt/module_xc/xc_functional_libxc_vxc.cpp

@@ -66,7 +66,7 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(		/
         std::tuple<std::vector<double>,std::vector<double>> rho_amag = XC_Functional_Libxc::convert_rho_amag_nspin4(nspin, nrxx, chr);
         rho = std::get<0>(std::move(rho_amag));
         amag = std::get<1>(std::move(rho_amag));
-        if(PARAM.inp.gga_grad >= 2 && (PARAM.globalv.domag || PARAM.globalv.domag_z))
+        if(PARAM.inp.gga_grad == 2 && (PARAM.globalv.domag || PARAM.globalv.domag_z))
         {
             mag_part = XC_Functional_Libxc::compute_mag_part_nspin4(nrxx, chr);
         }
@@ -76,7 +76,7 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(		/
     std::vector<double> sigma;
     if(is_gga)
     {
-        if(PARAM.inp.nspin==4 && PARAM.inp.gga_grad >= 2 && (PARAM.globalv.domag || PARAM.globalv.domag_z))
+        if(PARAM.inp.nspin==4 && PARAM.inp.gga_grad == 2 && (PARAM.globalv.domag || PARAM.globalv.domag_z))
         {
             gdr = XC_Functional_Libxc::cal_gdr_sf(nspin, nrxx, rho, mag_part, tpiba, chr);
         }
@@ -181,7 +181,7 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(		/
 
     if(4==PARAM.inp.nspin)
     {
-        if(PARAM.inp.gga_grad >= 2 && (PARAM.globalv.domag || PARAM.globalv.domag_z))
+        if(PARAM.inp.gga_grad == 2 && (PARAM.globalv.domag || PARAM.globalv.domag_z))
         {
             v = XC_Functional_Libxc::convert_v_nspin4_sf(nrxx, chr, mag_part, v);
         }

source/source_hamilt/module_xc/xc_functional_vxc.cpp

@@ -9,11 +9,13 @@
 #include "source_io/module_parameter/parameter.h"
 #include "xc_functional.h"
 
+#include "xc_functional_gga_noncol_sf_builtin.h"
+
 #ifdef USE_LIBXC
 #include "xc_functional_libxc.h"
 #endif
 
-#include "xc_functional_gga_noncol_sf_builtin.h"
+
 
 // [etxc, vtxc, v] = XC_Functional::v_xc(...)
 std::tuple<double, double, ModuleBase::matrix> XC_Functional::v_xc(const int& nrxx, // number of real-space grid
@@ -36,6 +38,11 @@ std::tuple<double, double, ModuleBase::matrix> XC_Functional::v_xc(const int& nr
 #endif
     }
 
+    if (PARAM.inp.nspin == 4 && (PARAM.globalv.domag || PARAM.globalv.domag_z) && PARAM.inp.gga_grad == 3)
+    {
+        return ModuleXC::NCGGA_SF_Builtin::v_xc_ncgga_sf_builtin(nrxx, ucell->omega, ucell->tpiba, chr);
+    }
+
     ModuleBase::timer::start("XC_Functional", "v_xc");
 
     //Exchange-Correlation potential Vxc(r) from n(r)