Fix LRI force with a spectrally truncated Coulomb metric

source/Makefile.Objects

@@ -866,6 +866,9 @@ OBJS_LCAO=evolve_elec.o\
       center2orb_orb22.o\
 
 OBJS_MODULE_RI=conv_coulomb_pot_k.o\
+      invmat_ri.o\
+      lri_cv.o\
+      lri_cv_tools.o\
       exx_abfs_jle.o\
       module_restart/restart_exx_csr.o\
       exx_abfs_io.o\

source/source_io/module_parameter/read_inp_exx_dftu.cpp

@@ -264,7 +264,7 @@ void ReadInput::item_exx()
         item.annotation = "threshold to inverse Vq in abfs for generating Cs";
         item.category = "Exact Exchange (LCAO)";
         item.type = "Real";
-        item.description = "By default, the Coulomb matrix inversion required for obtaining LRI coefficients is performed using LU decomposition. However, this approach may suffer from numerical instabilities when a large set of auxiliary basis functions (ABFs) is employed. When exx_cs_inv_thr > 0, the inversion is instead carried out via matrix diagonalization. Eigenvalues smaller than exx_cs_inv_thr are discarded to improve numerical stability. A relatively safe and commonly recommended value is 1e-5.";
+        item.description = "By default, the Coulomb matrix inversion required for obtaining LRI coefficients is performed using Cholesky factorization. However, this approach may suffer from numerical instabilities when a large set of auxiliary basis functions (ABFs) is employed. When exx_cs_inv_thr > 0, the inversion is instead carried out via matrix diagonalization. Eigenvalues at or below exx_cs_inv_thr times the largest eigenvalue are discarded to improve numerical stability. A relatively safe and commonly recommended value is 1e-5. Force derivatives assume that modes remain on the same side of the cutoff. Eigenvalue crossings can make the truncated model nonsmooth during relaxation or molecular dynamics.";
         item.default_value = "-1";
         item.unit = "";
         read_sync_double(input.exx_cs_inv_thr);

source/source_lcao/module_ri/CMakeLists.txt

@@ -5,6 +5,8 @@ if(ENABLE_LIBRI)
   list(
     APPEND
     objects
+    invmat_ri.cpp
+    lri_cv_tools.cpp
     matrix_orbs11.cpp
     matrix_orbs21.cpp
     matrix_orbs22.cpp
@@ -15,6 +17,7 @@ if(ENABLE_LIBRI)
     list(
       APPEND
       objects
+      lri_cv.cpp
       conv_coulomb_pot_k.cpp
       exx_abfs_ctor_orbs.cpp
       exx_abfs_io.cpp

source/source_lcao/module_ri/ewald_vq.hpp

@@ -6,6 +6,7 @@
 #ifndef EWALD_VQ_HPP
 #define EWALD_VQ_HPP
 
+#include "lri_cv_tools.h"
 #include <RI/comm/mix/Communicate_Tensors_Map_Judge.h>
 #include <RI/distribute/Distribute_Equally.h>
 #include <RI/global/Global_Func-1.h>

source/source_lcao/module_ri/exx_lri.hpp

@@ -6,6 +6,7 @@
 #ifndef EXX_LRI_HPP
 #define EXX_LRI_HPP
 
+#include "lri_cv_tools.h"
 #include "exx_lri.h"
 #include "source_cell/unitcell.h"
 #include "ri_2d_comm.h"

source/source_lcao/module_ri/invmat_ri.hpp → source/source_lcao/module_ri/invmat_ri.cpp

@@ -3,10 +3,9 @@
 // DATE :   2022-08-17
 //=======================
 
-#ifndef INVERSE_MATRIX_HPP
-#define INVERSE_MATRIX_HPP
 
 #include "invmat_ri.h"
+#include "abfs_construct_pca.h"
 #include "source_base/module_external/lapack_connector.h"
 
 #include <algorithm>
@@ -103,7 +102,7 @@ struct InverseMatrixTraits<std::complex<float>>
 };
 
 template <typename T>
-inline void Inverse_Matrix<T>::using_syev(const double& threshold_condition_number)
+void Inverse_Matrix<T>::using_syev(const double& threshold_condition_number)
 {
     using traits = InverseMatrixTraits<T>;
     using val_t = typename traits::value_type;
@@ -143,6 +142,71 @@ inline void Inverse_Matrix<T>::using_syev(const double& threshold_condition_numb
     BlasConnector::gemm(traits::gemm_trans, 'N', eA.nc, eA.nc, ie, 1, eA.c, eA.nc, eA.c, eA.nc, 0, A.ptr(), A.shape[1]);
 }
 
+template <typename T>
+std::array<RI::Tensor<T>, 3> Inverse_Matrix<T>::cal_inverse_with_derivatives(
+    const std::array<RI::Tensor<T>, 3>& derivatives,
+    const double threshold_condition_number)
+{
+    using traits = InverseMatrixTraits<T>;
+    using real_type = typename traits::value_type;
+    const std::size_t n = A.shape[0];
+    std::vector<real_type> eigenvalues(n);
+    int info = 0;
+    traits::syev(A, eigenvalues.data(), info);
+    if (info != 0)
+    {
+        throw std::range_error("Spectral inverse eigensolve failed, info=" + std::to_string(info));
+    }
+    const real_type maximum = *std::max_element(eigenvalues.begin(), eigenvalues.end());
+    const double cutoff = maximum * threshold_condition_number;
+    std::vector<bool> retained(n);
+    std::vector<real_type> inverse_values(n, 0);
+    for (std::size_t i = 0; i < n; ++i)
+    {
+        retained[i] = eigenvalues[i] > cutoff;
+        if (retained[i])
+        {
+            inverse_values[i] = 1 / eigenvalues[i];
+        }
+    }
+    // LAPACK eigenvectors are stored as rows U^H in the row-major tensor.
+    const RI::Tensor<T> vectors = A.copy();
+    const RI::Tensor<T> adjoint = vectors.conjugate().transpose();
+    RI::Tensor<T> weighted = vectors.copy();
+    for (std::size_t i = 0; i < n; ++i)
+    {
+        for (std::size_t j = 0; j < n; ++j)
+        {
+            weighted(i, j) *= inverse_values[i];
+        }
+    }
+    A = adjoint * weighted;
+    std::array<RI::Tensor<T>, 3> result;
+    for (std::size_t axis = 0; axis < 3; ++axis)
+    {
+        RI::Tensor<T> response = vectors * derivatives[axis] * adjoint;
+        for (std::size_t i = 0; i < n; ++i)
+        {
+            for (std::size_t j = 0; j < n; ++j)
+            {
+                real_type divided_difference = 0;
+                if (retained[i] && retained[j])
+                {
+                    divided_difference = -inverse_values[i] * inverse_values[j];
+                }
+                else if (retained[i] != retained[j])
+                {
+                    divided_difference = (inverse_values[i] - inverse_values[j])
+                                         / (eigenvalues[i] - eigenvalues[j]);
+                }
+                response(i, j) *= divided_difference;
+            }
+        }
+        result[axis] = adjoint * response * vectors;
+    }
+    return result;
+}
+
 /*
 void Inverse_Matrix::using_syev( const double &threshold_condition_number )
 {
@@ -259,4 +323,6 @@ void Inverse_Matrix<Tdata>::copy_down_triangle()
             A(i0, i1) = A(i1, i0);
 }
 
-#endif
+// Supported real and complex double-precision RI backends.
+template class Inverse_Matrix<double>;
+template class Inverse_Matrix<std::complex<double>>;

source/source_lcao/module_ri/invmat_ri.h

@@ -7,10 +7,10 @@
 #define ABACUS_SOURCE_LCAO_MODULE_RI_INVMAT_RI_H
 
 
-#include "abfs_construct_pca.h"
 
 #include <RI/global/Tensor.h>
 #include <vector>
+#include <array>
 
 template <typename Tdata>
 class Inverse_Matrix
@@ -23,6 +23,11 @@ class Inverse_Matrix
     };
     void cal_inverse(const Method& method, const double& threshold_condition_number = 0.);
 
+    /// Invert with a fixed-rank spectral cutoff and differentiate the retained subspace.
+    std::array<RI::Tensor<Tdata>, 3> cal_inverse_with_derivatives(
+        const std::array<RI::Tensor<Tdata>, 3>& derivatives,
+        const double threshold_condition_number);
+
     void input(const RI::Tensor<Tdata>& m);
     void input(const std::vector<std::vector<RI::Tensor<Tdata>>>& ms);
     RI::Tensor<Tdata> output() const;
@@ -36,6 +41,5 @@ class Inverse_Matrix
     RI::Tensor<Tdata> A;
 };
 
-#include "invmat_ri.hpp"
 
 #endif // ABACUS_SOURCE_LCAO_MODULE_RI_INVMAT_RI_H

source/source_lcao/module_ri/lri_cv.hpp → source/source_lcao/module_ri/lri_cv.cpp

@@ -3,8 +3,6 @@
 // DATE :   2022-08-17
 //=======================
 
-#ifndef LRI_CV_HPP
-#define LRI_CV_HPP
 
 #include "lri_cv.h"
 #include "lri_cv_tools.h"
@@ -406,11 +404,25 @@ LRI_CV<Tdata>::DPcal_C_dC(
                       DPcal_V(it1, it1, {0,0,0}, {{"writable_Vws",true}})}};
 
             std::vector<std::vector<RI::Tensor<Tdata>>> L;
+            std::array<LRI_CV_Tools::RIBlockMatrix<Tdata>, 3> dL;
             const double cs_inv_thr = this->p_info_ri != nullptr ? this->p_info_ri->Cs_inv_thr : 0.0;
-            if (cs_inv_thr > 0)
+            if (cs_inv_thr > 0 && !flag_finish_dC)
+            {
+                const auto derivative01 = DPcal_dV(it0, it1, R, flags);
+                const auto reverse_derivative = DPcal_dV(it1, it0, Rm, flags);
+                const auto derivative10 = LRI_CV_Tools::negative(reverse_derivative);
+                auto inverse = LRI_CV_Tools::cal_I_dI(V, derivative01, derivative10, cs_inv_thr);
+                L = std::move(inverse.first);
+                dL = std::move(inverse.second);
+            }
+            else if (cs_inv_thr > 0)
+            {
                 L = LRI_CV_Tools::cal_I(V, Inverse_Matrix<Tdata>::Method::syev, cs_inv_thr);
+            }
             else
+            {
                 L = LRI_CV_Tools::cal_I(V);
+            }
 
             const std::vector<RI::Tensor<Tdata>> C = LRI_CV_Tools::mul2(L,A);
             if(flags.at("writable_Cws"))
@@ -438,17 +450,34 @@ LRI_CV<Tdata>::DPcal_C_dC(
                                 this->index_abfs, this->index_lcaos, this->index_lcaos,
                                 Matrix_Orbs21::Matrix_Order::A1BA2))};
 
-                const std::array<RI::Tensor<Tdata>,3> dV_01 = DPcal_dV(it0, it1, R, flags);
-                const std::array<RI::Tensor<Tdata>,3> dV_10 = LRI_CV_Tools::negative(DPcal_dV(it1, it0, Rm, flags));
-
-                std::array<std::vector<RI::Tensor<Tdata>>,3>		// dC = L*(dA-dV*C)
+                std::array<std::vector<RI::Tensor<Tdata>>, 3> dC_tmp;
+                if (cs_inv_thr > 0)
+                {
+                    // C = F(V) A: include rotations between retained and discarded modes.
+                    for (std::size_t axis = 0; axis < 3; ++axis)
+                    {
+                        const std::vector<RI::Tensor<Tdata>> derivative_A = {dA[0][axis], dA[1][axis]};
+                        dC_tmp[axis] = LRI_CV_Tools::mul2(L, derivative_A);
+                        const auto inverse_response = LRI_CV_Tools::mul2(dL[axis], A);
+                        for (std::size_t block = 0; block < 2; ++block)
+                        {
+                            dC_tmp[axis][block] += inverse_response[block];
+                        }
+                    }
+                }
+                else
+                {
+                    const auto dV_01 = DPcal_dV(it0, it1, R, flags);
+                    const auto reverse_derivative = DPcal_dV(it1, it0, Rm, flags);
+                    const auto dV_10 = LRI_CV_Tools::negative(reverse_derivative);
                     dC_tmp = LRI_CV_Tools::mul2(
                             L,
                             LRI_CV_Tools::change_order( LRI_CV_Tools::minus(
                                 dA,
                                 std::vector<std::array<RI::Tensor<Tdata>,3>>{
                                     LRI_CV_Tools::mul1(dV_01, C[1]),
                                     LRI_CV_Tools::mul1(dV_10, C[0])})));
+                }
                 const std::vector<std::array<RI::Tensor<Tdata>,3>>
                     dC = LRI_CV_Tools::change_order(std::move(dC_tmp));
                 if(flags.at("writable_dCws"))
@@ -464,5 +493,6 @@ LRI_CV<Tdata>::DPcal_C_dC(
     } // end else (!(C_read && flag_finish_dC))
 }
 
-
-#endif
+// Supported real and complex double-precision RI backends.
+template class LRI_CV<double>;
+template class LRI_CV<std::complex<double>>;

source/source_lcao/module_ri/lri_cv.h

@@ -45,13 +45,13 @@ class LRI_CV
         const double &kmesh_times,
         std::shared_ptr<ORB_gaunt_table> MGT,
         const bool& init_C);
-    inline std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>
+    std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>
     cal_Vs(
         const UnitCell &ucell,
         const std::vector<TA> &list_A0,
         const std::vector<TAC> &list_A1,
         const std::map<std::string,bool> &flags);						// "writable_Vws"
-    inline std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>>
+    std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>>
     cal_dVs(
         const UnitCell &ucell,
         const std::vector<TA> &list_A0,
@@ -108,16 +108,16 @@ class LRI_CV
         const T_func_cal_Rcut& func_cal_Rcut,
         const T_func_DPcal_data<Tresult>& func_DPcal_data);
 
-    inline double cal_V_Rcut(const int it0, const int it1);
-    inline double cal_C_Rcut(const int it0, const int it1);
+    double cal_V_Rcut(const int it0, const int it1);
+    double cal_C_Rcut(const int it0, const int it1);
 
-    inline RI::Tensor<Tdata>
+    RI::Tensor<Tdata>
     DPcal_V(
         const int it0,
         const int it1,
         const Abfs::Vector3_Order<double> &R,
         const std::map<std::string,bool> &flags);						// "writable_Vws"
-    inline std::array<RI::Tensor<Tdata>,3>
+    std::array<RI::Tensor<Tdata>,3>
     DPcal_dV(
         const int it0,
         const int it1,
@@ -141,6 +141,5 @@ class LRI_CV
         const Tfunc &func_cal_o11);
 };
 
-#include "lri_cv.hpp"
 
 #endif

source/source_lcao/module_ri/lri_cv_tools.cpp

@@ -0,0 +1,51 @@
+#include "lri_cv_tools.h"
+
+#include <complex>
+
+template <typename T>
+std::pair<LRI_CV_Tools::RIBlockMatrix<T>, std::array<LRI_CV_Tools::RIBlockMatrix<T>, 3>>
+LRI_CV_Tools::cal_I_dI(const RIBlockMatrix<T>& metric,
+                      const std::array<RI::Tensor<T>, 3>& derivative01,
+                      const std::array<RI::Tensor<T>, 3>& derivative10,
+                      const double threshold)
+{
+    const std::vector<std::size_t> sizes = {metric[0][0].shape[0], metric[1][1].shape[0]};
+    std::array<RI::Tensor<T>, 3> derivatives;
+    for (std::size_t axis = 0; axis < 3; ++axis)
+    {
+        const RIBlockMatrix<T> blocks = {
+            {RI::Tensor<T>({sizes[0], sizes[0]}), derivative01[axis]},
+            {derivative10[axis], RI::Tensor<T>({sizes[1], sizes[1]})}};
+        Inverse_Matrix<T> pack;
+        pack.input(blocks);
+        derivatives[axis] = pack.output();
+    }
+    Inverse_Matrix<T> inverse;
+    inverse.input(metric);
+    const auto responses = inverse.cal_inverse_with_derivatives(derivatives, threshold);
+    std::array<RIBlockMatrix<T>, 3> block_responses;
+    for (std::size_t axis = 0; axis < 3; ++axis)
+    {
+        Inverse_Matrix<T> pack;
+        pack.input(responses[axis]);
+        block_responses[axis] = pack.output(sizes, sizes);
+    }
+    auto inverse_blocks = inverse.output(sizes, sizes);
+    return std::make_pair(inverse_blocks, block_responses);
+}
+
+template std::pair<LRI_CV_Tools::RIBlockMatrix<double>,
+                   std::array<LRI_CV_Tools::RIBlockMatrix<double>, 3>>
+LRI_CV_Tools::cal_I_dI(
+    const RIBlockMatrix<double>& metric,
+    const std::array<RI::Tensor<double>, 3>& derivative01,
+    const std::array<RI::Tensor<double>, 3>& derivative10,
+    const double threshold);
+
+template std::pair<LRI_CV_Tools::RIBlockMatrix<std::complex<double>>,
+                   std::array<LRI_CV_Tools::RIBlockMatrix<std::complex<double>>, 3>>
+LRI_CV_Tools::cal_I_dI(
+    const RIBlockMatrix<std::complex<double>>& metric,
+    const std::array<RI::Tensor<std::complex<double>>, 3>& derivative01,
+    const std::array<RI::Tensor<std::complex<double>>, 3>& derivative10,
+    const double threshold);

source/source_lcao/module_ri/lri_cv_tools.h

@@ -7,6 +7,7 @@
 #define LRI_CV_TOOLS_H
 
 #include "invmat_ri.h"
+#include "source_cell/unitcell.h"
 #include "abfs_vector3_order.h"
 #include "source_lcao/module_ri/abfs.h"
 
@@ -16,6 +17,7 @@
 #include <cstddef>
 #include <map>
 #include <vector>
+#include <utility>
 
 namespace LRI_CV_Tools
 {
@@ -30,6 +32,16 @@ extern std::vector<std::vector<RI::Tensor<Tdata>>> cal_I(const std::vector<std::
                                                          = Inverse_Matrix<Tdata>::Method::potrf,
                                                          const double& threshold_condition_number = 0.);
 
+template <typename T>
+using RIBlockMatrix = std::vector<std::vector<RI::Tensor<T>>>;
+
+template <typename T>
+std::pair<RIBlockMatrix<T>, std::array<RIBlockMatrix<T>, 3>> cal_I_dI(
+    const RIBlockMatrix<T>& metric,
+    const std::array<RI::Tensor<T>, 3>& derivative01,
+    const std::array<RI::Tensor<T>, 3>& derivative10,
+    const double threshold);
+
 template <typename Tdata>
 inline RI::Tensor<Tdata> transform_Rm(const RI::Tensor<Tdata>& V);
 template <typename Tdata>