Feature: LR-TDDFT analytical forces (gradients) by Z-vector method and excited-state geometry relaxation

.github/workflows/test.yml

@@ -83,6 +83,41 @@ jobs:
           ./scripts/stage4/install_stage4.sh
           cd ..
 
+      # Keep the image's LibXC version: only enable the KXC needed by LR gradients.
+      - name: Identify LibXC version
+        id: libxc
+        run: echo "version=$(pkg-config --modversion libxc)" >> "${GITHUB_OUTPUT}"
+
+      - name: Cache LibXC with KXC
+        id: libxc-cache
+        uses: actions/cache@v6
+        with:
+          path: toolchain/install/libxc-kxc
+          key: libxc-kxc-v1-${{ runner.os }}-${{ runner.arch }}-${{ steps.compiler.outputs.fingerprint }}-${{ steps.libxc.outputs.version }}
+
+      - name: Build LibXC with KXC
+        if: steps.libxc-cache.outputs.cache-hit != 'true'
+        env:
+          LIBXC_VERSION: ${{ steps.libxc.outputs.version }}
+        run: |
+          wget -q "https://gitlab.com/libxc/libxc/-/archive/${LIBXC_VERSION}/libxc-${LIBXC_VERSION}.tar.gz" -O /tmp/libxc-kxc.tar.gz
+          mkdir -p /tmp/libxc-kxc-source
+          tar -xzf /tmp/libxc-kxc.tar.gz --strip-components=1 -C /tmp/libxc-kxc-source
+          cmake -S /tmp/libxc-kxc-source -B /tmp/libxc-kxc-build -G Ninja \
+            -DCMAKE_BUILD_TYPE=Release \
+            -DCMAKE_INSTALL_PREFIX="${GITHUB_WORKSPACE}/toolchain/install/libxc-kxc" \
+            -DCMAKE_INSTALL_LIBDIR=lib \
+            -DBUILD_SHARED_LIBS=ON -DENABLE_FORTRAN=OFF -DBUILD_TESTING=OFF \
+            -DDISABLE_KXC=OFF -DDISABLE_LXC=ON
+          cmake --build /tmp/libxc-kxc-build -j8
+          cmake --install /tmp/libxc-kxc-build
+
+      - name: Use KXC-enabled LibXC at runtime
+        run: echo "LD_LIBRARY_PATH=${GITHUB_WORKSPACE}/toolchain/install/libxc-kxc/lib:${LD_LIBRARY_PATH:-}" >> "${GITHUB_ENV}"
+
+      - name: Test LR excited-state tracking checker
+        run: python3 tests/integrate/validation_tools/test_es_check.py
+
       - name: Configure
         run: |
           source toolchain/install/setup
@@ -91,6 +126,7 @@ jobs:
             -DCMAKE_EXPORT_COMPILE_COMMANDS=ON \
             -DENABLE_MLALGO=ON \
             -DENABLE_LIBXC=ON \
+            -DLibxc_DIR="${GITHUB_WORKSPACE}/toolchain/install/libxc-kxc/share/cmake/Libxc" \
             -DENABLE_LIBRI=ON \
             -DENABLE_GOOGLEBENCH=OFF \
             -DENABLE_JSON=ON \

source/Makefile.Objects

@@ -131,6 +131,7 @@ ${OBJS_DELTASPIN}\
 ${OBJS_TENSOR}\
 ${OBJS_HSOLVER_PEXSI}\
 ${OBJS_LR}\
+${OBJS_LR_GRAD}\
 ${OBJS_RDMFT}
 
 OBJS_MAIN=main.o\
@@ -447,6 +448,7 @@ OBJS_HAMILT_LCAO=hamilt_lcao.o\
     nonlocal_dh.o\
     nonlocal_fs.o\
     overlap.o\
+    ovlp_block.o\
     overlap_fs.o\
     td_ekinetic_lcao.o\
     td_nonlocal_lcao.o\
@@ -859,6 +861,7 @@ OBJS_MODULE_RI=conv_coulomb_pot_k.o\
       symm_rot_out.o\
       gaussian_abfs.o\
       exx_lri_detail.o\
+      exx_workspace.o\
       ri_util.o\
       singular_value.o\
 
@@ -1072,8 +1075,8 @@ OBJS_TENSOR=tensor.o\
     refcount.o
 
     OBJS_LR=lr_util.o\
-    lr_util_hcontainer.o\
     utils/lr_io.o\
+    utils/lr_sparse.o\
     utils/exciton_plotter.o\
     ao_to_mo_parallel.o\
     ao_to_mo_serial.o\
@@ -1083,11 +1086,45 @@ OBJS_TENSOR=tensor.o\
     operator_lr_exx.o\
     xc_kernel.o\
     pot_hxc_lrtd.o\
+    kernel_pots.o\
+    init_helpers.o\
+    hsolver_lrtd.o\
+    zeq_solver.o\
+    zeq_cg.o\
+    hamilt_ulr.o\
+    lr_density.o\
+    spectrum_mo.o\
+    spectrum_out.o\
+    dm_diff_ser.o\
+    dm_diff_par.o\
     lr_spectrum.o\
     lr_spectrum_velocity.o\
     hamilt_casida.o\
+    esolver_lr_orb.o\
+    esolver_lr_init.o\
+    esolver_lr_cfg.o\
     esolver_lr_lcao_tddft.o\
 
+OBJS_LR_GRAD=lr_force.o\
+    xc_kxc.o\
+    lr_force_aux.o\
+    grad_degen.o\
+    state_group.o\
+    state_track.o\
+    state_ovlp.o\
+    grad_jt.o\
+    gradient_output.o\
+    lr_grad_cs.o\
+    lr_grad_os.o\
+    exx_proj.o\
+    cvcx_serial.o\
+    cvcx_par.o\
+    cal_edm.o\
+    zeqlin_solv.o\
+    pot_grad_xc.o\
+    esolver_lr_grad.o\
+    esolver_lr_rlx.o\
+
 OBJS_RDMFT=rdmft.o\
         rdmft_tools.o\
         rdmft_pot.o\

source/source_base/matrix.cpp

@@ -75,7 +75,12 @@ matrix::matrix( matrix && m_in )
 matrix& matrix::operator=( const matrix & m_in )
 {
     this->create( m_in.nr, m_in.nc, false );
-    memcpy( c, m_in.c, nr*nc*sizeof(double) );
+    // `create` leaves `c` null for an empty matrix, and memcpy's arguments are declared
+    // non-null even for a zero count -- so assigning an empty matrix is undefined behaviour.
+    if( nr && nc )
+    {
+        memcpy( c, m_in.c, nr*nc*sizeof(double) );
+    }
     return *this;
 }
 
@@ -157,6 +162,16 @@ void matrix::create( const int nrow, const int ncol, const bool flag_zero )
     }
 }
 
+/* Unary minus*/
+matrix operator-(const matrix& m1)
+{
+    matrix tm(m1);
+    const int size = m1.nr * m1.nc;
+    for (int i = 0; i < size; i++)
+        tm.c[i] = -tm.c[i];
+    return tm;
+}
+
 /* Adding matrices, as a friend */
 matrix operator+(const matrix &m1, const matrix &m2)
 {

source/source_base/matrix.h

@@ -79,7 +79,7 @@ class matrix
     using type=double;					// Peize Lin add 2022.08.08 for template
 };
 
-
+matrix operator-(const matrix& m1); // unary minus
 matrix operator+(const matrix &m1, const matrix &m2);
 matrix operator-(const matrix &m1, const matrix &m2);
 matrix operator*(const matrix &m1, const matrix &m2);

source/source_base/module_container/base/third_party/lapack.h

@@ -45,6 +45,22 @@ extern "C"
 int ilaenv_(const int* ispec, const char* name, const char* opts,
             const int* n1, const int* n2, const int* n3, const int* n4);
 
+// Symmetric/Hermitian eigenproblems (QR algorithm)
+void dsyev_(const char* jobz, const char* uplo, const int* n, double* a,
+            const int* lda, double* w, double* work, const int* lwork, int* info);
+void zheev_(const char* jobz, const char* uplo, const int* n, std::complex<double>* a,
+            const int* lda, double* w, std::complex<double>* work, const int* lwork,
+            double* rwork, int* info);
+
+// General nonsymmetric eigenproblems
+void dgeev_(const char* jobvl, const char* jobvr, const int* n, double* a, const int* lda,
+            double* wr, double* wi, double* vl, const int* ldvl, double* vr,
+            const int* ldvr, double* work, const int* lwork, int* info);
+void zgeev_(const char* jobvl, const char* jobvr, const int* n, std::complex<double>* a,
+            const int* lda, std::complex<double>* w, std::complex<double>* vl, const int* ldvl,
+            std::complex<double>* vr, const int* ldvr, std::complex<double>* work,
+            const int* lwork, double* rwork, int* info);
+
 // Generalized symmetric-definite eigenproblems (divide-and-conquer)
 void ssygvd_(const int* itype, const char* jobz, const char* uplo, const int* n,
              float* a, const int* lda,
@@ -340,6 +356,37 @@ void hegvx(const int itype, const char jobz, const char range, const char uplo,
 }
 
 
+// Column-major general eigenproblems; preserve the supplied workspace/query policy.
+static inline
+void geev(const char jobvl, const char jobvr, const int n, double* a, const int lda,
+          double* wr, double* wi, double* vl, const int ldvl, double* vr,
+          const int ldvr, double* work, const int lwork, int& info)
+{
+    dgeev_(&jobvl, &jobvr, &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, work, &lwork, &info);
+}
+static inline
+void geev(const char jobvl, const char jobvr, const int n, std::complex<double>* a, const int lda,
+          std::complex<double>* w, std::complex<double>* vl, const int ldvl,
+          std::complex<double>* vr, const int ldvr, std::complex<double>* work,
+          const int lwork, double* rwork, int& info)
+{
+    zgeev_(&jobvl, &jobvr, &n, a, &lda, w, vl, &ldvl, vr, &ldvr, work, &lwork, rwork, &info);
+}
+
+// Column-major symmetric/Hermitian eigenproblems; lwork=-1 queries the workspace.
+static inline
+void heev(const char jobz, const char uplo, const int n, double* a, const int lda,
+          double* w, double* work, const int lwork, int& info)
+{
+    dsyev_(&jobz, &uplo, &n, a, &lda, w, work, &lwork, &info);
+}
+static inline
+void heev(const char jobz, const char uplo, const int n, std::complex<double>* a, const int lda,
+          double* w, std::complex<double>* work, const int lwork, double* rwork, int& info)
+{
+    zheev_(&jobz, &uplo, &n, a, &lda, w, work, &lwork, rwork, &info);
+}
+
 // wrap function of fortran lapack routine zheevx.
 static inline
 void heevx(const char jobz, const char range, const char uplo, const int n,

source/source_base/module_external/blacs_connector.h

@@ -59,6 +59,9 @@ extern "C"
 
     void Czgebs2d(int ConTxt, char *scope, char *top, int m, int n, std::complex<double> *A, int lda);
     void Czgebr2d(int ConTxt, char *scope, char *top, int m, int n, std::complex<double> *A, int lda, int rsrc, int csrc);
+
+    // element-wise sum over `scope`; rdest = -1 leaves the result on every process
+    void Cigsum2d(int ConTxt, char *scope, char *top, int m, int n, int *A, int lda, int rdest, int cdest);
 }
 
 // unified interface for broadcast

source/source_base/module_external/scalapack_connector.h

@@ -118,6 +118,20 @@ extern "C"
             int *ipiv, std::complex<double>* B, const int* ib, const int* jb, const int*descb, const int *info
         );
 
+    void pdgesv_(
+            const int *n, const int *nrhs,
+            double *A, const int *ia, const int *ja, const int *desca,
+            int *ipiv, double* B, const int* ib, const int* jb, const int*descb, int *info
+        );
+
+    void pdpotrs_(const char* uplo, const int* n, const int* nrhs,
+        const double* A, const int* ia, const int* ja, const int* desca,
+        double* B, const int* ib, const int* jb, const int* descb, int* info);
+
+    void pzpotrs_(const char* uplo, const int* n, const int* nrhs,
+        const std::complex<double>* A, const int* ia, const int* ja, const int* desca,
+        std::complex<double>* B, const int* ib, const int* jb, const int* descb, int* info);
+
     void pdsygvx_(const int* itype, const char* jobz, const char* range, const char* uplo,
         const int* n, double* A, const int* ia, const int* ja, const int*desca, double* B, const int* ib, const int* jb, const int*descb,
         const double* vl, const double* vu, const int* il, const int* iu,
@@ -396,34 +410,61 @@ class ScalapackConnector
         pzgesv_(&n, &nrhs, A, &ia, &ja, desca, ipiv, B, &ib, &jb, descb, info);
     }
 
+    static inline
+    void gesv(
+        const int n, const int nrhs,
+        double *A, const int ia, const int ja, const int *desca,
+        int *ipiv, double* B, const int ib, const int jb, const int*descb, int *info)
+    {
+        pdgesv_(&n, &nrhs, A, &ia, &ja, desca, ipiv, B, &ib, &jb, descb, info);
+    }
+
+    static inline
+    void potrs(
+        const char uplo, const int n, const int nrhs,
+        const double* A, const int ia, const int ja, const int* desca,
+        double* B, const int ib, const int jb, const int* descb, int* info)
+    {
+        pdpotrs_(&uplo, &n, &nrhs, A, &ia, &ja, desca, B, &ib, &jb, descb, info);
+    }
+
+    static inline
+    void potrs(
+        const char uplo, const int n, const int nrhs,
+        const std::complex<double>* A, const int ia, const int ja, const int* desca,
+        std::complex<double>* B, const int ib, const int jb, const int* descb, int* info)
+    {
+        pzpotrs_(&uplo, &n, &nrhs, A, &ia, &ja, desca, B, &ib, &jb, descb, info);
+    }
+
     static inline
     void tranu(
         const int m, const int n,
-        const double alpha, double *a, const int ia, const int ja, const int *desca,
+        const double alpha, const double *a, const int ia, const int ja, const int *desca,
         const double beta, double *c, const int ic, const int jc, const int *descc)
     {
         pdtran_(&m, &n, &alpha, a, &ia, &ja, desca, &beta, c, &ic, &jc, descc);
     }
     static inline
     void tranu(
         const int m, const int n,
-        const std::complex<double> alpha, std::complex<double> *a, const int ia, const int ja, const int *desca,
+        const std::complex<double> alpha, const std::complex<double> *a, const int ia, const int ja, const int *desca,
         const std::complex<double> beta, std::complex<double> *c, const int ic, const int jc, const int *descc)
     {
         pztranu_(&m, &n, &alpha, a, &ia, &ja, desca, &beta, c, &ic, &jc, descc);
     }
     static inline
     void tranc(
         const int m, const int n,
-        const double alpha, double *a, const int ia, const int ja, const int *desca,
+        const double alpha, const double *a, const int ia, const int ja, const int *desca,
         const double beta, double *c, const int ic, const int jc, const int *descc)
     {
         pdtran_(&m, &n, &alpha, a, &ia, &ja, desca, &beta, c, &ic, &jc, descc);
     }
     static inline
     void tranc(
         const int m, const int n,
-        const std::complex<double> alpha, std::complex<double> *a, const int ia, const int ja, const int *desca,
+        const std::complex<double> alpha, const std::complex<double> *a, const int ia, const int ja, const int *desca,
         const std::complex<double> beta, std::complex<double> *c, const int ic, const int jc, const int *descc)
     {
         pztranc_(&m, &n, &alpha, a, &ia, &ja, desca, &beta, c, &ic, &jc, descc);

source/source_base/parallel_cell.cpp

@@ -47,6 +47,17 @@ int CommunicationDomain::max(int value) const
     return value;
 }
 
+int CommunicationDomain::sum(int value) const
+{
+#ifdef __MPI
+    if (communicator_ != MPI_COMM_NULL)
+    {
+        MPI_Allreduce(MPI_IN_PLACE, &value, 1, MPI_INT, MPI_SUM, communicator_);
+    }
+#endif
+    return value;
+}
+
 double CommunicationDomain::max(double value) const
 {
 #ifdef __MPI

source/source_base/parallel_cell.h

@@ -18,6 +18,7 @@ class CommunicationDomain
     int rank() const;
     int size() const;
     int max(int value) const;
+    int sum(int value) const;
     double max(double value) const;
 
 private:

source/source_base/parallel_grid.h

@@ -43,7 +43,7 @@ class Parallel_Grid
     int get_nz() const { return ncz; }
     int get_nrxx() const { return nrxx; }
 
-    private:
+	private:
 
     void z_distribution(const int kpar);
 

source/source_esolver/CMakeLists.txt

@@ -20,7 +20,12 @@ if(ENABLE_LCAO)
   list(APPEND objects
       esolver_ks_lcao.cpp
       esolver_ks_lcao_tddft.cpp
+      esolver_lr_orb.cpp
+      esolver_lr_init.cpp
+      esolver_lr_cfg.cpp
       esolver_lr_lcao_tddft.cpp
+      esolver_lr_grad.cpp
+      esolver_lr_rlx.cpp
       esolver_gets.cpp
       lcao_others.cpp
       esolver_dm2rho.cpp