Merge branch 'develop' into develop

Author: mohanchen

.github/workflows/test.yml

@@ -17,11 +17,17 @@ jobs:
       volumes:
         - /tmp/ccache:/github/home/.ccache
     steps:
-      - name: Checkout
+      - name: Checkout repository
         uses: actions/checkout@v5
         with:
-          submodules: recursive
           fetch-depth: 0
+          # We will handle submodules manually after fixing ownership
+          submodules: 'false'
+
+      - name: Take ownership of the workspace and update submodules
+        run: |
+          sudo chown -R $(whoami) .
+          git submodule update --init --recursive
 
       - name: Install CI tools
         run: |

docs/advanced/input_files/input-main.md

@@ -1487,15 +1487,20 @@ These variables are used to control the geometry relaxation.
 
 ### relax_method
 
-- **Type**: String
+- **Type**: Vector of string
 - **Description**: The methods to do geometry optimization.
+  the first element:
   - cg: using the conjugate gradient (CG) algorithm. Note that there are two implementations of the conjugate gradient (CG) method, see [relax_new](#relax_new).
-  - bfgs: using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm.
-  - bfgs_trad: using the traditional Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm.
+  - bfgs : using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm.
+  - lbfgs: using the Limited-memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) algorithm.
   - cg_bfgs: using the CG method for the initial steps, and switching to BFGS method when the force convergence is smaller than [relax_cg_thr](#relax_cg_thr).
   - sd: using the steepest descent (SD) algorithm.
   - fire: the Fast Inertial Relaxation Engine method (FIRE), a kind of molecular-dynamics-based relaxation algorithm, is implemented in the molecular dynamics (MD) module. The algorithm can be used by setting [calculation](#calculation) to `md` and [md_type](#md_type) to `fire`. Also ionic velocities should be set in this case. See [fire](../md.md#fire) for more details.
-- **Default**: cg
+
+  the second element:
+  when the first element is bfgs, if the second parameter is 1, it indicates the use of the new BFGS algorithm; if the second parameter is not 1, it indicates the use of the old BFGS algorithm.
+- **Default**: cg 1
+- **Note**:In the 3.10-LTS version, the type of this parameter is std::string. It can be set to "cg","bfgs","cg_bfgs","bfgs_trad","lbfgs","sd","fire".
 
 ### relax_new
 

python/pyabacus/CONTRIBUTING.md

@@ -189,7 +189,7 @@ list(APPEND _diago
     ${HSOLVER_PATH}/diago_david.cpp
     ${HSOLVER_PATH}/diag_const_nums.cpp
     ${HSOLVER_PATH}/diago_iter_assist.cpp
-    ${HSOLVER_PATH}/kernels/dngvd_op.cpp
+    ${HSOLVER_PATH}/kernels/hegvd_op.cpp
     ${HSOLVER_PATH}/kernels/bpcg_kernel_op.cpp
     ${BASE_PATH}/kernels/math_kernel_op.cpp
     ${BASE_PATH}/kernels/math_kernel_op_vec.cpp

python/pyabacus/src/hsolver/CMakeLists.txt

@@ -9,7 +9,7 @@ list(APPEND _diago
     ${HSOLVER_PATH}/diago_pxxxgvx.cpp
 
 
-    ${HSOLVER_PATH}/kernels/dngvd_op.cpp
+    ${HSOLVER_PATH}/kernels/hegvd_op.cpp
     ${HSOLVER_PATH}/kernels/bpcg_kernel_op.cpp
     # dependency
     ${BASE_PATH}/kernels/math_kernel_op.cpp

python/pyabacus/src/hsolver/py_diago_cg.hpp

@@ -145,7 +145,7 @@ class PyDiagoCG
             std::copy(hpsi_ptr, hpsi_ptr + nvec * ld_psi, hpsi_out.data<std::complex<double>>());
         };
 
-        auto subspace_func = [] (const ct::Tensor& psi_in, ct::Tensor& psi_out) { /*do nothing*/ };
+        auto subspace_func = [](const ct::Tensor& psi_in, ct::Tensor& psi_out, const bool S_orth) { /*do nothing*/ };
 
         auto spsi_func = [this] (const ct::Tensor& psi_in, ct::Tensor& spsi_out) {
             const auto ndim = psi_in.shape().ndim();

python/pyabacus/src/hsolver/py_diago_dav_subspace.hpp

@@ -144,7 +144,6 @@ class PyDiagoDavSubspace
             dav_ndim, 
             tol, 
             max_iter, 
-            need_subspace, 
             comm_info,
             diag_subspace,
             nb2d

source/CMakeLists.txt

@@ -35,7 +35,7 @@ list(APPEND device_srcs
   source_pw/module_pwdft/kernels/meta_op.cpp
   source_pw/module_stodft/kernels/hpsi_norm_op.cpp
   source_basis/module_pw/kernels/pw_op.cpp
-  source_hsolver/kernels/dngvd_op.cpp
+  source_hsolver/kernels/hegvd_op.cpp
   source_hsolver/kernels/bpcg_kernel_op.cpp
   source_estate/kernels/elecstate_op.cpp
 
@@ -70,7 +70,7 @@ if(USE_CUDA)
     source_pw/module_stodft/kernels/cuda/hpsi_norm_op.cu
     source_pw/module_pwdft/kernels/cuda/onsite_op.cu
     source_basis/module_pw/kernels/cuda/pw_op.cu
-    source_hsolver/kernels/cuda/dngvd_op.cu
+    source_hsolver/kernels/cuda/hegvd_op.cu
     source_hsolver/kernels/cuda/bpcg_kernel_op.cu
     source_estate/kernels/cuda/elecstate_op.cu
 
@@ -101,7 +101,7 @@ if(USE_ROCM)
     source_pw/module_pwdft/kernels/rocm/onsite_op.hip.cu
     source_pw/module_stodft/kernels/rocm/hpsi_norm_op.hip.cu
     source_basis/module_pw/kernels/rocm/pw_op.hip.cu
-    source_hsolver/kernels/rocm/dngvd_op.hip.cu
+    source_hsolver/kernels/rocm/hegvd_op.hip.cu
     source_hsolver/kernels/rocm/bpcg_kernel_op.hip.cu
     source_estate/kernels/rocm/elecstate_op.hip.cu
 

source/Makefile.Objects

@@ -267,6 +267,7 @@ OBJS_ESOLVER=esolver.o\
     esolver_lj.o\
     esolver_dp.o\
     esolver_of.o\
+    esolver_of_tddft.o\
     esolver_of_tool.o\
     esolver_of_interface.o\
     pw_others.o\
@@ -361,6 +362,7 @@ OBJS_HAMILT_OF=kedf_tf.o\
     kedf_xwm.o\
     kedf_lkt.o\
     kedf_manager.o\
+    evolve_ofdft.o\
 
 OBJS_HAMILT_LCAO=hamilt_lcao.o\
     operator_lcao.o\
@@ -396,7 +398,7 @@ OBJS_HSOLVER=diago_cg.o\
     hsolver_lcaopw.o\
     hsolver_pw_sdft.o\
     diago_iter_assist.o\
-    dngvd_op.o\
+    hegvd_op.o\
     bpcg_kernel_op.o\
     diag_const_nums.o\
     diag_hs_para.o\

source/source_base/kernels/dsp/dsp_connector.cpp

@@ -6,9 +6,9 @@
 extern "C"
 {
 #define complex_double ignore_complex_double
-#include <mt_hthread_blas.h> // MTBLAS_TRANSPOSE etc
+#include <mt_hthread_blas.h> // include faster mtblas kernels
 #undef complex_double
-#include <mtblas_interface.h> // gemm
+#include <mtblas_interface.h> // include normal mtblas kernels that automatically operate memory, but slower.
 }
 namespace mtfunc
 {
@@ -22,45 +22,42 @@ void dspDestoryHandle(int id)
 {
     hthread_dev_close(id);
     std::cout << " ** DSP closed on cluster " << id << " **" << std::endl;
-} // Close dsp cluster at the end
+} // Close dsp cluster at the end of the program
 
-MTBLAS_TRANSPOSE convertBLASTranspose(const char* blasTrans)
+// MTBlas secretly removed its MTBLAS_TRANSPOSE data type and used the original CBLAS_TRANSPOSE. So this function is modified.
+
+CBLAS_TRANSPOSE convertBLASTranspose(const char* blasTrans)
 {
     switch (blasTrans[0])
     {
     case 'N':
     case 'n':
-        return MtblasNoTrans;
+        return CblasNoTrans;
     case 'T':
     case 't':
-        return MtblasTrans;
+        return CblasTrans;
     case 'C':
     case 'c':
-        return MtblasConjTrans;
+        return CblasConjTrans;
     default:
         std::cout << "Invalid BLAS transpose parameter!! Use default instead." << std::endl;
-        return MtblasNoTrans;
+        return CblasNoTrans;
     }
-} // Used to convert normal transpost char to mtblas transpose flag
+} // Used to convert normal transpost char to cblas transpose flag
 
 void* malloc_ht(size_t bytes, int cluster_id)
 {
-    // std::cout << "MALLOC " << cluster_id;
     void* ptr = hthread_malloc((int)cluster_id, bytes, HT_MEM_RW);
-    // std::cout << ptr << " SUCCEED" << std::endl;;
     return ptr;
-}
+} // Malloc on dsp. Used to replace original malloc
+
 
-// Used to replace original malloc
 
 void free_ht(void* ptr)
 {
-    // std::cout << "FREE " << ptr;
     hthread_free(ptr);
-    // std::cout << " FREE SUCCEED" << std::endl;
-}
+} // Free on dsp. Used to replace original free
 
-// Used to replace original free
 
 void sgemm_mt_(const char* transa,
                const char* transb,
@@ -77,7 +74,7 @@ void sgemm_mt_(const char* transa,
                const int* ldc,
                int cluster_id)
 {
-    mtblas_sgemm(MTBLAS_ORDER::MtblasColMajor,
+    mtblas_sgemm(CBLAS_ORDER::CblasColMajor,
                  convertBLASTranspose(transa),
                  convertBLASTranspose(transb),
                  *m,
@@ -109,7 +106,7 @@ void dgemm_mt_(const char* transa,
                const int* ldc,
                int cluster_id)
 {
-    mtblas_dgemm(MTBLAS_ORDER::MtblasColMajor,
+    mtblas_dgemm(CBLAS_ORDER::CblasColMajor,
                  convertBLASTranspose(transa),
                  convertBLASTranspose(transb),
                  *m,
@@ -141,7 +138,7 @@ void zgemm_mt_(const char* transa,
                const int* ldc,
                int cluster_id)
 {
-    mtblas_zgemm(MTBLAS_ORDER::MtblasColMajor,
+    mtblas_zgemm(CBLAS_ORDER::CblasColMajor,
                  convertBLASTranspose(transa),
                  convertBLASTranspose(transb),
                  *m,
@@ -173,7 +170,7 @@ void cgemm_mt_(const char* transa,
                const int* ldc,
                int cluster_id)
 {
-    mtblas_cgemm(MTBLAS_ORDER::MtblasColMajor,
+    mtblas_cgemm(CBLAS_ORDER::CblasColMajor,
                  convertBLASTranspose(transa),
                  convertBLASTranspose(transb),
                  *m,
@@ -207,7 +204,7 @@ void sgemm_mth_(const char* transa,
                 const int* ldc,
                 int cluster_id)
 {
-    mt_hthread_sgemm(MTBLAS_ORDER::MtblasColMajor,
+    mt_hthread_sgemm(CBLAS_ORDER::CblasColMajor,
                      convertBLASTranspose(transa),
                      convertBLASTranspose(transb),
                      *m,
@@ -239,7 +236,7 @@ void dgemm_mth_(const char* transa,
                 const int* ldc,
                 int cluster_id)
 {
-    mt_hthread_dgemm(MTBLAS_ORDER::MtblasColMajor,
+    mt_hthread_dgemm(CBLAS_ORDER::CblasColMajor,
                      convertBLASTranspose(transa),
                      convertBLASTranspose(transb),
                      *m,
@@ -275,7 +272,7 @@ void zgemm_mth_(const char* transa,
     *alp = *alpha;
     std::complex<double>* bet = (std::complex<double>*)malloc_ht(sizeof(std::complex<double>), cluster_id);
     *bet = *beta;
-    mt_hthread_zgemm(MTBLAS_ORDER::MtblasColMajor,
+    mt_hthread_zgemm(CBLAS_ORDER::CblasColMajor,
                      convertBLASTranspose(transa),
                      convertBLASTranspose(transb),
                      *m,
@@ -314,7 +311,7 @@ void cgemm_mth_(const char* transa,
     std::complex<float>* bet = (std::complex<float>*)malloc_ht(sizeof(std::complex<float>), cluster_id);
     *bet = *beta;
 
-    mt_hthread_cgemm(MTBLAS_ORDER::MtblasColMajor,
+    mt_hthread_cgemm(CBLAS_ORDER::CblasColMajor,
                      convertBLASTranspose(transa),
                      convertBLASTranspose(transb),
                      *m,

source/source_base/module_container/ATen/kernels/cuda/lapack.cu

@@ -88,7 +88,7 @@ struct lapack_potrf<T, DEVICE_GPU> {
 };
 
 template <typename T>
-struct lapack_dnevd<T, DEVICE_GPU> {
+struct lapack_heevd<T, DEVICE_GPU> {
     using Real = typename GetTypeReal<T>::type;
     void operator()(
         const char& jobz,
@@ -97,12 +97,12 @@ struct lapack_dnevd<T, DEVICE_GPU> {
         const int& dim,
         Real* eigen_val)
     {
-        cuSolverConnector::dnevd(cusolver_handle, jobz, uplo, dim, Mat, dim, eigen_val);
+        cuSolverConnector::heevd(cusolver_handle, jobz, uplo, dim, Mat, dim, eigen_val);
     }
 };
 
 template <typename T>
-struct lapack_dngvd<T, DEVICE_GPU> {
+struct lapack_hegvd<T, DEVICE_GPU> {
     using Real = typename GetTypeReal<T>::type;
     void operator()(
         const int& itype,
@@ -113,7 +113,7 @@ struct lapack_dngvd<T, DEVICE_GPU> {
         const int& dim,
         Real* eigen_val)
     {
-        cuSolverConnector::dngvd(cusolver_handle, itype, jobz, uplo, dim, Mat_A, dim, Mat_B, dim, eigen_val);
+        cuSolverConnector::hegvd(cusolver_handle, itype, jobz, uplo, dim, Mat_A, dim, Mat_B, dim, eigen_val);
     }
 };
 
@@ -175,15 +175,15 @@ template struct lapack_potrf<double, DEVICE_GPU>;
 template struct lapack_potrf<std::complex<float>,  DEVICE_GPU>;
 template struct lapack_potrf<std::complex<double>, DEVICE_GPU>;
 
-template struct lapack_dnevd<float,  DEVICE_GPU>;
-template struct lapack_dnevd<double, DEVICE_GPU>;
-template struct lapack_dnevd<std::complex<float>,  DEVICE_GPU>;
-template struct lapack_dnevd<std::complex<double>, DEVICE_GPU>;
+template struct lapack_heevd<float,  DEVICE_GPU>;
+template struct lapack_heevd<double, DEVICE_GPU>;
+template struct lapack_heevd<std::complex<float>,  DEVICE_GPU>;
+template struct lapack_heevd<std::complex<double>, DEVICE_GPU>;
 
-template struct lapack_dngvd<float,  DEVICE_GPU>;
-template struct lapack_dngvd<double, DEVICE_GPU>;
-template struct lapack_dngvd<std::complex<float>,  DEVICE_GPU>;
-template struct lapack_dngvd<std::complex<double>, DEVICE_GPU>;
+template struct lapack_hegvd<float,  DEVICE_GPU>;
+template struct lapack_hegvd<double, DEVICE_GPU>;
+template struct lapack_hegvd<std::complex<float>,  DEVICE_GPU>;
+template struct lapack_hegvd<std::complex<double>, DEVICE_GPU>;
 
 template struct lapack_getrf<float,  DEVICE_GPU>;
 template struct lapack_getrf<double, DEVICE_GPU>;