fix nspin=1

Author: maki49

source/source_esolver/esolver_lr_lcao_tddft.cpp

@@ -845,7 +845,7 @@ void ModuleESolver::ESolver_LR<T, TR>::init_pot(const Charge& chg_gs)
         // while the S1 integrand indexes them as if there were 1 -- it does not even read a
         // consistent spin combination. Use `ST::S2_gs` there, which is exactly half of S2_singlet,
         // matching the `K_Hxc(singlet) = 2 * pot_hxc_gs` convention of the gradient operators.
-        const ST st_gs = (nspin == 1) ? ST::S1 : (oshell ? ST::S2_updown : ST::S2_gs);
+        const ST st_gs = (nspin == 1) ? ST::S1_gs : (oshell ? ST::S2_updown : ST::S2_gs);
         // `pot_hxc_gs` supplies the $g^{xc}$ of both $W^c$ and the force term, for either spin.
         // Those two call sites are guarded by `has_local_xc(xc_kernel)` -- the *LR* kernel name --
         // so an `xc_kernel rpa` run never touches them however local `dft_functional` is.
@@ -918,7 +918,8 @@ void ModuleESolver::ESolver_LR<T, TR>::read_ks_chg(Charge& chg_gs)
     for (int is = 0; is < this->nspin; ++is)
     {
         std::stringstream ssc;
-        ssc << this->in_dir << "chgs" << is + 1 << ".cube";
+        if (this->nspin == 1) { ssc << this->in_dir << "chg.cube"; }
+        else { ssc << this->in_dir << "chgs" << is + 1 << ".cube"; }
         GlobalV::ofs_running << ssc.str() << std::endl;
         if (ModuleIO::read_vdata_palgrid(Pgrid,
             GlobalV::MY_RANK,

source/source_lcao/module_lr/Grad/force/lr_force.cpp

@@ -7,6 +7,9 @@
 // #include "source_lcao/module_lr/utils/lr_util_hcontainer.h"
 namespace LR
 {
+    /// `dm_gs` carries the ground-state occupations, at nspin=1 they are 2 for fully occupied bands.
+    inline double gs_dm_channel_factor() { return (PARAM.inp.nspin == 1) ? 0.5 : 1.0; }
+
     template<typename TK>
     Charge LR_Force<TK>::dm_to_charge(const elecstate::DensityMatrix<TK, double>& dm)
     {
@@ -128,7 +131,8 @@ namespace LR
             std::vector<const double*> vr_eff = { v_lin.c };
             ModuleGint::cal_gint_fvl(1, vr_eff, dm_gs.get_DMR_vector(), true, false, &fhxc_dvhxc, &stress_tmp);
         }
-        if(!reproduce_gs) {fhxc_dvhxc *= 2;} // for the two channels of the ground-state dm. 
+        if(!reproduce_gs) {fhxc_dvhxc *= 2;} // for the two channels of the ground-state dm.
+        fhxc_dvhxc *= gs_dm_channel_factor();
 
         // 4. kinetic (Pulay)
         std::vector<hamilt::HContainer<double>> dT = cal_hs_grad('T', this->ucell_, this->pv_, this->gd_, this->two_center_bundle_);
@@ -193,6 +197,7 @@ namespace LR
         ModuleBase::matrix stress_tmp;
         std::vector<const double*> vr_eff = { v2.c };
         ModuleGint::cal_gint_fvl(1, vr_eff, dm_gs.get_DMR_vector(), true, false, &f, &stress_tmp);
+        f *= gs_dm_channel_factor();
         return f;
     }
 

source/source_lcao/module_lr/potentials/pot_hxc_lrtd.cpp

@@ -59,10 +59,10 @@ namespace LR
         // Hartree
         switch (this->spin_type_)
         {
-        case SpinType::S1: case SpinType::S2_updown: case SpinType::S2_gs:
+        case SpinType::S1_gs: case SpinType::S2_updown: case SpinType::S2_gs:
             v_eff += elecstate::H_Hartree_pw::v_hartree(ucell, const_cast<ModulePW::PW_Basis*>(&this->rho_basis_), 1, rho);
             break;
-        case SpinType::S2_singlet:
+        case SpinType::S1: case SpinType::S2_singlet:
             v_eff += 2 * elecstate::H_Hartree_pw::v_hartree(ucell, const_cast<ModulePW::PW_Basis*>(&this->rho_basis_), 1, rho);
             break;
         default:
@@ -90,11 +90,16 @@ namespace LR
             switch (s)
             {
             case SpinType::S1:
-                funcs[s] = [this, &fxc](FXC_PARA_TYPE)->void
+            case SpinType::S1_gs:
+            {
+                // S1_gs is exactly half of S1 (see the SpinType doc in the header).
+                const double prefac = (s == SpinType::S1_gs) ? 1.0 : 2.0;
+                funcs[s] = [this, &fxc, prefac](FXC_PARA_TYPE)->void
                     {
-                        for (int ir = 0;ir < nrxx;++ir) { v_eff(0, ir) += ModuleBase::e2 * fxc.v2rho2.at(ir) * rho[ir]; }
+                        for (int ir = 0;ir < nrxx;++ir) { v_eff(0, ir) += ModuleBase::e2 * prefac * fxc.v2rho2.at(ir) * rho[ir]; }
                     };
                 break;
+            }
             case SpinType::S2_singlet:
             case SpinType::S2_gs:
             {
@@ -140,7 +145,10 @@ namespace LR
             switch (s)
             {
             case SpinType::S1:
-                funcs[s] = [this, &fxc](FXC_PARA_TYPE)->void
+            case SpinType::S1_gs:
+            {
+                const double prefac = (s == SpinType::S1_gs) ? 1.0 : 2.0;   // S1_gs is exactly half of S1.
+                funcs[s] = [this, &fxc, prefac](FXC_PARA_TYPE)->void
                     {
                         // test: output drho
                         // double thr = 1e-1;
@@ -176,9 +184,10 @@ namespace LR
                             vxc_tmp[ir] += (fxc.v2rho2.at(ir) * rho[ir]
                                 + fxc.v2rhosigma_2drho.at(ir) * drho.at(ir));
                         }
-                        BlasConnector::axpy(nrxx, ModuleBase::e2, vxc_tmp.data(), 1, v_eff.c, 1);
+                        BlasConnector::axpy(nrxx, ModuleBase::e2 * prefac, vxc_tmp.data(), 1, v_eff.c, 1);
                     };
                 break;
+            }
             case SpinType::S2_singlet:
             case SpinType::S2_gs:
             {

source/source_lcao/module_lr/potentials/pot_hxc_lrtd.h

@@ -10,22 +10,22 @@ namespace LR
     class PotHxcLR : public PotLRBase
     {
     public:
-        /// S1: K^Hartree + K^xc
+        /// S1: the nspin=1 singlet LR kernel, 2*K^Hartree + 2*K^xc.
+        ///     The factor 2 on *both* terms is what makes it equal to S2_singlet. 
+        /// S1_gs: the nspin=1 counterpart of S2_gs, i.e. the *ground-state* Hxc kernel
+        ///     K^Hartree + K^xc = S1 / 2. This is what `pot_hxc_gs` needs at nspin=1.
         /// S2_singlet: 2*K^Hartree + K^xc_{upup} + K^xc_{updown}
         /// S2_triplet: K^xc_{upup} - K^xc_{updown}
         /// S2_updown: K^Hartree + (K^xc_{upup}, K^xc_{updown},  K^xc_{downup} or K^xc_{downdown}), according to `ispin_op` (for spin-polarized systems)
-        /// S2_gs: the nspin=2 counterpart of S1, i.e. the *ground-state* Hxc kernel
-        ///     K^Hartree + (K^xc_{upup} + K^xc_{updown})/2 = S2_singlet / 2.
-        ///     Used for `pot_hxc_gs` in LR gradients, where the convention is
-        ///     `K_Hxc(singlet) = 2 * pot_hxc_gs` (see `cal_multiplier_w_from_z.h`).
+        /// S2_gs: the nspin=2 *ground-state* Hxc kernel K^Hartree + (K^xc_{upup} + K^xc_{updown})/2 = S2_singlet / 2.
+        ///     Used for `pot_hxc_gs` in LR gradients, where the convention is `K_Hxc(singlet) = 2 * pot_hxc_gs`.
         ///     The 1/2 on the xc part is not a convention but the chain rule: the derivative is
         ///     taken w.r.t. the *total* density matrix, and $\partial v_u/\partial\rho =
-        ///     (f_{uu}+f_{ud})/2$ because $\rho_u=\rho_d=\rho/2$. The Hartree part needs no
-        ///     halving, which is exactly why S1 and S2_gs share the same Hartree weight.
+        ///     (f_{uu}+f_{ud})/2$ because $\rho_u=\rho_d=\rho/2$. The Hartree part needs no halving.
         ///     Do NOT use S1 here when nspin=2: `KernelXC` is built with `PARAM.inp.nspin`, so the
         ///     kernel arrays carry 3 spin components per grid point while the S1 integrand indexes
         ///     them as if there were 1.
-        enum SpinType { S1 = 0, S2_singlet = 1, S2_triplet = 2, S2_updown = 3, S2_gs = 4 };
+        enum SpinType { S1 = 0, S2_singlet = 1, S2_triplet = 2, S2_updown = 3, S2_gs = 4, S1_gs = 5 };
         /// XCType here is to determin the method of integration from kernel to potential, not the way calculating the kernel
         enum XCType { None = 0, LDA = 1, GGA = 2, HYB_GGA = 4 };
         /// constructor building exchange-correlation kernel

source/source_lcao/module_lr/potentials/xc_kernel.cpp

@@ -223,12 +223,17 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl
             xc_gga_fxc(&func, nrxx, rho.data(), sigma.data(), v2rho2_tmp.data(), v2rhosigma_tmp.data(), v2sigma2_tmp.data());
             // std::cout << "max element of v2sigma2_tmp: " << *std::max_element(v2sigma2_tmp.begin(), v2sigma2_tmp.end()) << std::endl;
             // std::cout << "rho corresponding to max element of v2sigma2_tmp: " << rho[(std::max_element(v2sigma2_tmp.begin(), v2sigma2_tmp.end()) - v2sigma2_tmp.begin()) / 6] << std::endl;
-            // cut off by sgn
-            cutoff_grid_data_spin2(vrho_tmp, sgn);
-            cutoff_grid_data_spin2(vsigma_tmp, sgn);
-            cutoff_grid_data_spin2(v2rho2_tmp, sgn);
-            cutoff_grid_data_spin2(v2rhosigma_tmp, sgn);
-            cutoff_grid_data_spin2(v2sigma2_tmp, sgn);
+            // cut off by sgn. nspin=2 only: `cutoff_grid_data_spin2` assumes >1 component per
+            // grid point (it asserts on it), and at nspin=1 there is exactly one, for which both
+            // of its `for_each` ranges are empty -- the cutoff is a no-op anyway. 
+            if (nspin == 2)
+            {
+                cutoff_grid_data_spin2(vrho_tmp, sgn);
+                cutoff_grid_data_spin2(vsigma_tmp, sgn);
+                cutoff_grid_data_spin2(v2rho2_tmp, sgn);
+                cutoff_grid_data_spin2(v2rhosigma_tmp, sgn);
+                cutoff_grid_data_spin2(v2sigma2_tmp, sgn);
+            }
             if (need_kxc)
             {
                 xc_gga_kxc(&func, nrxx, rho.data(), sigma.data(),
@@ -503,23 +508,36 @@ void LR::KernelXC::build_gxc_coef(GxcCoef& dst, const bool triplet, const int& n
 
     if (nspin == 1)
     {
-        // Single component everywhere; all the weight sums collapse to 1 (section 3).
+        // Single component everywhere; all the weight sums collapse to 1.
+        //
+        // ... and then the whole set is scaled by 4 to reach the *singlet* normalization, the same
+        // one `nspin=2` produces below and the only one the rest of the gradient code knows about.
+        // libxc's unpolarized derivatives are taken w.r.t. the TOTAL density, so for a closed shell
+        // (rho_u = rho_d = rho/2) an n-th derivative is 2^(n-1) smaller than the singlet spin
+        // combination: d^2E/drho^2 = (f_uu+f_ud)/2 -- which is why `SpinType::S1` carries a 2 --
+        // and d^3E/drho^3 = (g_uuu + 3g_uud)/4, while the nspin=2 branch below builds
+        // a_s2 = g_uuu + 2g_uud + g_udd = g_uuu + 3g_uud. Hence 4 here, 2 there.
+        //
+        // Measured on H2/SZ/LDA: without it the GXC DMTRANS force was 0.08549 eV/Ang against the
+        // nspin=2 singlet's 0.17097 -- a factor 2, being 1/4 from this and 2 from the `dm_gs`
+        // channel convention (`gs_dm_channel_factor` in `lr_force.cpp`).
+        constexpr double to_singlet = 4.;
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 4096)
 #endif
         for (int i = 0;i < nrxx;++i)
         {
-            dst.a_s2[i] = v3r3[i];
+            dst.a_s2[i] = to_singlet * v3r3[i];
             if (!is_gga) { continue; }
-            dst.a_st[i] = v3r2s[i] * 4.;
-            dst.a_t2[i] = v3rs2[i] * 4.;
-            dst.a_q[i] = v2rs[i] * 2.;
-            dst.c_s[i] = v2rs[i] * 4.;
-            dst.c_t[i] = v2s2[i] * 8.;
-            dst.e_s2[i] = v3r2s[i] * 2.;
-            dst.e_st[i] = v3rs2[i] * 8.;
-            dst.e_t2[i] = v3s3[i] * 8.;
-            dst.e_q[i] = v2s2[i] * 4.;
+            dst.a_st[i] = to_singlet * v3r2s[i] * 4.;
+            dst.a_t2[i] = to_singlet * v3rs2[i] * 4.;
+            dst.a_q[i] = to_singlet * v2rs[i] * 2.;
+            dst.c_s[i] = to_singlet * v2rs[i] * 4.;
+            dst.c_t[i] = to_singlet * v2s2[i] * 8.;
+            dst.e_s2[i] = to_singlet * v3r2s[i] * 2.;
+            dst.e_st[i] = to_singlet * v3rs2[i] * 8.;
+            dst.e_t2[i] = to_singlet * v3s3[i] * 8.;
+            dst.e_q[i] = to_singlet * v2s2[i] * 4.;
         }
         return;
     }