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

@@ -129,6 +129,7 @@ ${OBJS_DELTASPIN}\
 ${OBJS_TENSOR}\
 ${OBJS_HSOLVER_PEXSI}\
 ${OBJS_LR}\
+${OBJS_LR_GRAD}\
 ${OBJS_RDMFT}
 
 OBJS_MAIN=main.o\
@@ -445,6 +446,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\
@@ -857,6 +859,7 @@ OBJS_MODULE_RI=conv_coulomb_pot_k.o\
       symm_rot_out.o\
       gaussian_abfs.o\
       exx_lri_detail.o\
+      exx_lr_ws.o\
       ri_util.o\
       rpa_lri.o\
       rpa_lri_detail.o\
@@ -1063,7 +1066,6 @@ OBJS_TENSOR=tensor.o\
     refcount.o
 
     OBJS_LR=lr_util.o\
-    lr_util_hcontainer.o\
     utils/lr_io.o\
     utils/exciton_plotter.o\
     ao_to_mo_parallel.o\
@@ -1077,8 +1079,31 @@ OBJS_TENSOR=tensor.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

@@ -82,7 +82,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_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,6 +410,33 @@ 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,

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