Refactor: unified two-center integration interface

source/module_basis/module_nao/two_center_bundle.cpp

@@ -1,18 +1,21 @@
 #include "module_basis/module_nao/two_center_bundle.h"
+
+#include "module_base/global_variable.h"
 #include "module_base/memory.h"
+#include "module_base/parallel_common.h"
 #include "module_base/ylm.h"
 #include "module_basis/module_nao/real_gaunt_table.h"
+
 #include <memory>
-#include "module_base/parallel_common.h"
-#include "module_base/global_variable.h"
 
 void TwoCenterBundle::build_orb(int ntype, const std::string* file_orb0)
 {
     std::vector<std::string> file_orb(ntype);
     if (GlobalV::MY_RANK == 0)
     {
-        std::transform(file_orb0, file_orb0 + ntype, file_orb.begin(),
-                [](const std::string& file) { return GlobalV::global_orbital_dir + file; });
+        std::transform(file_orb0, file_orb0 + ntype, file_orb.begin(), [](const std::string& file) {
+            return GlobalV::global_orbital_dir + file;
+        });
     }
 #ifdef __MPI
     Parallel_Common::bcast_string(file_orb.data(), ntype);
@@ -48,7 +51,7 @@ void TwoCenterBundle::build_alpha(int ndesc, std::string* file_desc0)
 
 void TwoCenterBundle::build_orb_onsite(int ntype, double radius)
 {
-    if(GlobalV::onsite_radius > 0)
+    if (GlobalV::onsite_radius > 0)
     {
         orb_onsite_ = std::unique_ptr<RadialCollection>(new RadialCollection);
         orb_onsite_->build(orb_.get(), GlobalV::onsite_radius);
@@ -60,23 +63,28 @@ void TwoCenterBundle::tabulate()
     ModuleBase::SphericalBesselTransformer sbt(true);
     orb_->set_transformer(sbt);
     beta_->set_transformer(sbt);
-    if (alpha_) alpha_->set_transformer(sbt);
-    if (orb_onsite_) orb_onsite_->set_transformer(sbt);
+    if (alpha_)
+        alpha_->set_transformer(sbt);
+    if (orb_onsite_)
+        orb_onsite_->set_transformer(sbt);
 
     //================================================================
     //              build two-center integration tables
     //================================================================
     // set up a universal radial grid
     double rmax = std::max(orb_->rcut_max(), beta_->rcut_max());
-    if (alpha_) rmax = std::max(rmax, alpha_->rcut_max());
+    if (alpha_)
+        rmax = std::max(rmax, alpha_->rcut_max());
     double dr = 0.01;
     double cutoff = 2.0 * rmax;
     int nr = static_cast<int>(rmax / dr) + 1;
 
     orb_->set_uniform_grid(true, nr, cutoff, 'i', true);
     beta_->set_uniform_grid(true, nr, cutoff, 'i', true);
-    if (alpha_) alpha_->set_uniform_grid(true, nr, cutoff, 'i', true);
-    if (orb_onsite_) orb_onsite_->set_uniform_grid(true, nr, cutoff, 'i', true);
+    if (alpha_)
+        alpha_->set_uniform_grid(true, nr, cutoff, 'i', true);
+    if (orb_onsite_)
+        orb_onsite_->set_uniform_grid(true, nr, cutoff, 'i', true);
 
     // build TwoCenterIntegrator objects
     kinetic_orb = std::unique_ptr<TwoCenterIntegrator>(new TwoCenterIntegrator);
@@ -106,13 +114,95 @@ void TwoCenterBundle::tabulate()
 
     ModuleBase::Memory::record("RealGauntTable", RealGauntTable::instance().memory());
 
-    // init Ylm (this shall be done by Ylm automatically! to be done later...)
-    ModuleBase::Ylm::set_coefficients();
+    sbt.clear();
+}
+
+void TwoCenterBundle::tabulate(const double lcao_ecut,
+                               const double lcao_dk,
+                               const double lcao_dr,
+                               const double lcao_rmax)
+{
+    ModuleBase::SphericalBesselTransformer sbt(true);
+    orb_->set_transformer(sbt);
+    beta_->set_transformer(sbt);
+    if (alpha_)
+        alpha_->set_transformer(sbt);
+    if (orb_onsite_)
+        orb_onsite_->set_transformer(sbt);
+
+    //================================================================
+    //              build two-center integration tables
+    //================================================================
+
+    // old formula for the number of k-space grid points
+    int nk = static_cast<int>(sqrt(lcao_ecut) / lcao_dk) + 4;
+    nk += 1 - nk % 2; // make nk odd
+
+    std::vector<double> kgrid(nk);
+    for (int ik = 0; ik < nk; ++ik)
+    {
+        kgrid[ik] = ik * lcao_dk;
+    }
+
+    orb_->set_grid(false, nk, kgrid.data(), 't');
+    beta_->set_grid(false, nk, kgrid.data(), 't');
+    if (alpha_)
+    {
+        alpha_->set_grid(false, nk, kgrid.data(), 't');
+    }
+    if (orb_onsite_)
+    {
+        orb_onsite_->set_grid(false, nk, kgrid.data(), 't');
+    }
+
+    // "st" stands for overlap (s) and kinetic (t)
+    const double cutoff_st = std::min(lcao_rmax, 2.0 * orb_->rcut_max());
+    const int nr_st = static_cast<int>(cutoff_st / lcao_dr) + 5;
+
+    kinetic_orb = std::unique_ptr<TwoCenterIntegrator>(new TwoCenterIntegrator);
+    kinetic_orb->tabulate(*orb_, *orb_, 'T', nr_st, cutoff_st);
+    ModuleBase::Memory::record("TwoCenterTable: Kinetic", kinetic_orb->table_memory());
+
+    overlap_orb = std::unique_ptr<TwoCenterIntegrator>(new TwoCenterIntegrator);
+    overlap_orb->tabulate(*orb_, *orb_, 'S', nr_st, cutoff_st);
+    ModuleBase::Memory::record("TwoCenterTable: Overlap", overlap_orb->table_memory());
+
+    // overlap between orbital and beta (for nonlocal potential)
+    const double cutoff_nl = std::min(lcao_rmax, orb_->rcut_max() + beta_->rcut_max());
+    const int nr_nl = static_cast<int>(cutoff_nl / lcao_dr) + 5;
+    overlap_orb_beta = std::unique_ptr<TwoCenterIntegrator>(new TwoCenterIntegrator);
+    overlap_orb_beta->tabulate(*orb_, *beta_, 'S', nr_nl, cutoff_nl);
+    ModuleBase::Memory::record("TwoCenterTable: Nonlocal", overlap_orb_beta->table_memory());
+
+    // overlap between orbital and deepks projector
+    if (alpha_)
+    {
+        const double cutoff_alpha = std::min(lcao_rmax, orb_->rcut_max() + alpha_->rcut_max());
+        const int nr_alpha = static_cast<int>(cutoff_alpha / lcao_dr) + 5;
+        overlap_orb_alpha = std::unique_ptr<TwoCenterIntegrator>(new TwoCenterIntegrator);
+        overlap_orb_alpha->tabulate(*orb_, *alpha_, 'S', nr_alpha, cutoff_alpha);
+        ModuleBase::Memory::record("TwoCenterTable: Descriptor", overlap_orb_beta->table_memory());
+    }
+
+    // overlap between orbital and "onsite orbital" (for DFT+U)
+    if (orb_onsite_)
+    {
+        const double cutoff_onsite = std::min(lcao_rmax, orb_->rcut_max() + orb_onsite_->rcut_max());
+        const int nr_onsite = static_cast<int>(cutoff_onsite / lcao_dr) + 5;
+        overlap_orb_onsite = std::unique_ptr<TwoCenterIntegrator>(new TwoCenterIntegrator);
+        overlap_orb_onsite->tabulate(*orb_, *orb_onsite_, 'S', nr_onsite, cutoff_onsite);
+    }
+
+    ModuleBase::Memory::record("RealGauntTable", RealGauntTable::instance().memory());
 
     sbt.clear();
 }
 
-void TwoCenterBundle::to_LCAO_Orbitals(LCAO_Orbitals& ORB, const double lcao_ecut, const double lcao_dk, const double lcao_dr, const double lcao_rmax) const
+void TwoCenterBundle::to_LCAO_Orbitals(LCAO_Orbitals& ORB,
+                                       const double lcao_ecut,
+                                       const double lcao_dk,
+                                       const double lcao_dr,
+                                       const double lcao_rmax) const
 {
     ORB.ntype = orb_->ntype();
     ORB.lmax = orb_->lmax();
@@ -121,7 +211,7 @@ void TwoCenterBundle::to_LCAO_Orbitals(LCAO_Orbitals& ORB, const double lcao_ecu
     ORB.dR = lcao_dr;
     ORB.Rmax = lcao_rmax; // lcao_rmax, see ORB_control.cpp
     ORB.dr_uniform = 0.001;
-    
+
     // Due to algorithmic difference in the spherical Bessel transform
     // (FFT vs. Simpson's integration), k grid of FFT is not appropriate
     // for Simpson's integration. The k grid for Simpson's integration is
@@ -130,17 +220,16 @@ void TwoCenterBundle::to_LCAO_Orbitals(LCAO_Orbitals& ORB, const double lcao_ecu
     ORB.ecutwfc = lcao_ecut;
     ORB.dk = lcao_dk;
 
-	if(ORB.ecutwfc < 20)
-	{
-		ORB.kmesh = static_cast<int>( 2 * sqrt(ORB.ecutwfc) / ORB.dk )  + 4;
-	}
-	else
-	{
-		ORB.kmesh = static_cast<int>( sqrt(ORB.ecutwfc) / ORB.dk )  + 4;
-	}
+    if (ORB.ecutwfc < 20)
+    {
+        ORB.kmesh = static_cast<int>(2 * sqrt(ORB.ecutwfc) / ORB.dk) + 4;
+    }
+    else
+    {
+        ORB.kmesh = static_cast<int>(sqrt(ORB.ecutwfc) / ORB.dk) + 4;
+    }
     ORB.kmesh += 1 - ORB.kmesh % 2;
 
-
     delete[] ORB.Phi;
     ORB.Phi = new Numerical_Orbital[orb_->ntype()];
     for (int itype = 0; itype < orb_->ntype(); ++itype)

source/module_basis/module_nao/two_center_bundle.h

@@ -1,8 +1,8 @@
-#ifndef TWO_CENTER_BUNDLE_H
-#define TWO_CENTER_BUNDLE_H
+#ifndef W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_MODULE_BASIS_MODULE_NAO_TWO_CENTER_BUNDLE_H
+#define W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_MODULE_BASIS_MODULE_NAO_TWO_CENTER_BUNDLE_H
 
-#include "module_basis/module_nao/two_center_integrator.h"
 #include "module_basis/module_ao/ORB_read.h"
+#include "module_basis/module_nao/two_center_integrator.h"
 
 #include <memory>
 #include <string>
@@ -21,6 +21,11 @@ class TwoCenterBundle
 
     void tabulate();
 
+    // Unlike the tabulate() above, this overload function computes
+    // two-center integration table by direct integration with Simpson's
+    // rule, which was the algorithm used prior to v3.3.4.
+    void tabulate(const double lcao_ecut, const double lcao_dk, const double lcao_dr, const double lcao_rmax);
+
     /**
      * @brief Overwrites the content of a LCAO_Orbitals object (e.g. GlobalC::ORB)
      * with the current object.
@@ -31,8 +36,7 @@ class TwoCenterBundle
                           const double lcao_ecut,
                           const double lcao_dk,
                           const double lcao_dr,
-                          const double lcao_rmax
-                          ) const;
+                          const double lcao_rmax) const;
 
     std::unique_ptr<TwoCenterIntegrator> kinetic_orb;
     std::unique_ptr<TwoCenterIntegrator> overlap_orb;

source/module_basis/module_nao/two_center_table.cpp

@@ -1,15 +1,15 @@
 #include "module_basis/module_nao/two_center_table.h"
 
+#include "module_base/constants.h"
+#include "module_base/cubic_spline.h"
+#include "module_base/math_integral.h"
+
 #include <algorithm>
 #include <cstring>
 #include <iostream>
 #include <limits>
 #include <numeric>
 
-#include "module_base/constants.h"
-#include "module_base/math_integral.h"
-#include "module_base/cubic_spline.h"
-
 void TwoCenterTable::build(const RadialCollection& bra,
                            const RadialCollection& ket,
                            const char op,
@@ -37,8 +37,13 @@ void TwoCenterTable::build(const RadialCollection& bra,
     std::for_each(rgrid_, rgrid_ + nr_, [this, dr](double& r) { r = (&r - rgrid_) * dr; });
 
     // index the table by generating a map from (itype1, l1, izeta1, itype2, l2, izeta2, l) to a row index
-    index_map_.resize({bra.ntype(), bra.lmax() + 1, bra.nzeta_max(),
-                       ket.ntype(), ket.lmax() + 1, ket.nzeta_max(), bra.lmax() + ket.lmax() + 1});
+    index_map_.resize({bra.ntype(),
+                       bra.lmax() + 1,
+                       bra.nzeta_max(),
+                       ket.ntype(),
+                       ket.lmax() + 1,
+                       ket.nzeta_max(),
+                       bra.lmax() + ket.lmax() + 1});
     std::fill(index_map_.data<int>(), index_map_.data<int>() + index_map_.NumElements(), -1);
 
     ntab_ = 0;
@@ -61,8 +66,8 @@ const double* TwoCenterTable::table(const int itype1,
 #ifdef __DEBUG
     assert(is_present(itype1, l1, izeta1, itype2, l2, izeta2, l));
 #endif
-    return deriv ? dtable_.inner_most_ptr<double>(index_map_.get_value<int>(itype1, l1, izeta1, itype2, l2, izeta2, l)):
-                    table_.inner_most_ptr<double>(index_map_.get_value<int>(itype1, l1, izeta1, itype2, l2, izeta2, l));
+    return deriv ? dtable_.inner_most_ptr<double>(index_map_.get_value<int>(itype1, l1, izeta1, itype2, l2, izeta2, l))
+                 : table_.inner_most_ptr<double>(index_map_.get_value<int>(itype1, l1, izeta1, itype2, l2, izeta2, l));
 }
 
 void TwoCenterTable::lookup(const int itype1,
@@ -82,12 +87,14 @@ void TwoCenterTable::lookup(const int itype1,
 
     if (R > rmax())
     {
-        if (val) *val = 0.0;
-        if (dval) *dval = 0.0;
+        if (val)
+            *val = 0.0;
+        if (dval)
+            *dval = 0.0;
         return;
     }
 
-    const double*  tab = table(itype1, l1, izeta1, itype2, l2, izeta2, l, false);
+    const double* tab = table(itype1, l1, izeta1, itype2, l2, izeta2, l, false);
     const double* dtab = table(itype1, l1, izeta1, itype2, l2, izeta2, l, true);
     ModuleBase::CubicSpline::eval(nr_, rgrid_, tab, dtab, 1, &R, val, dval);
 }
@@ -100,7 +107,9 @@ int& TwoCenterTable::table_index(const NumericalRadial* it1, const NumericalRadi
 void TwoCenterTable::cleanup()
 {
     op_ = '\0';
+    ntab_ = 0;
     nr_ = 0;
+    rmax_ = 0.0;
     delete[] rgrid_;
     rgrid_ = nullptr;
 
@@ -137,9 +146,7 @@ double TwoCenterTable::dfact(int l) const
     return result;
 }
 
-void TwoCenterTable::two_center_loop(const RadialCollection& bra, 
-                                     const RadialCollection& ket, 
-                                     looped_func f)
+void TwoCenterTable::two_center_loop(const RadialCollection& bra, const RadialCollection& ket, looped_func f)
 {
     for (int l = 0; l <= bra.lmax() + ket.lmax(); ++l)
         for (int l1 = 0; l1 <= bra.lmax(); ++l1)
@@ -179,7 +186,7 @@ void TwoCenterTable::_tabulate(const NumericalRadial* it1, const NumericalRadial
     //
     // See the developer's document for more details.
     double dr = rmax_ / (nr_ - 1);
-    if ( l > 0 )
+    if (l > 0)
     {
         // divide S(R) by R^l (except the R=0 point)
         std::for_each(&tab[1], tab + nr_, [&](double& val) { val /= std::pow(dr * (&val - tab), l); });
@@ -195,37 +202,37 @@ void TwoCenterTable::_tabulate(const NumericalRadial* it1, const NumericalRadial
         int op_exp = l;
         switch (op_)
         {
-        case 'S': op_exp += 2;
-                  break;
-        case 'T': op_exp += 4;
-                  break;
-        default: ; // currently not supposed to happen
+        case 'S':
+            op_exp += 2;
+            break;
+        case 'T':
+            op_exp += 4;
+            break;
+        default:; // currently not supposed to happen
         }
 
         for (int ik = 0; ik != nk; ++ik)
         {
-            fk[ik] = it1->kvalue(ik) * it2->kvalue(ik) 
-                    * std::pow(kgrid[ik], op_exp);
+            fk[ik] = it1->kvalue(ik) * it2->kvalue(ik) * std::pow(kgrid[ik], op_exp);
         }
 
-        tab[0] = ModuleBase::Integral::simpson(nk, fk, &h[1]) 
-                * ModuleBase::FOUR_PI / dfact(2 * l + 1);
+        tab[0] = ModuleBase::Integral::simpson(nk, fk, &h[1]) * ModuleBase::FOUR_PI / dfact(2 * l + 1);
 
         delete[] fk;
         delete[] h;
     }
 
-    // The derivative table stores the derivative of S(R)/R^l or T(R)/R^l 
+    // The derivative table stores the derivative of S(R)/R^l or T(R)/R^l
     // instead of bare dS(R)/dR or dT(R)/dR, which simplifies further calculation.
     //
     // The derivatives are computed from a cubic spline interpolation rather
     // than two spherical Bessel transforms. By doing so, we achieve a good
     // consistency between the table and its derivative during interpolation.
     using ModuleBase::CubicSpline;
-    CubicSpline::build(
-            nr_, rgrid_, table_.inner_most_ptr<double>(itab),
-            {CubicSpline::BoundaryType::first_deriv, 0.0},
-            {CubicSpline::BoundaryType::first_deriv, 0.0},
-            dtable_.inner_most_ptr<double>(itab));
+    CubicSpline::build(nr_,
+                       rgrid_,
+                       table_.inner_most_ptr<double>(itab),
+                       {CubicSpline::BoundaryType::first_deriv, 0.0},
+                       {CubicSpline::BoundaryType::first_deriv, 0.0},
+                       dtable_.inner_most_ptr<double>(itab));
 }
-