Refactor Propagator class implementation into multiple files for better code organization

Author: AsTonyshment

source/Makefile.Objects

@@ -561,6 +561,9 @@ OBJS_LCAO=evolve_elec.o\
       middle_hamilt.o\
       norm_psi.o\
       propagator.o\
+      propagator_cn2.o\
+      propagator_taylor.o\
+      propagator_etrs.o\
       td_velocity.o\
       td_current.o\
       snap_psibeta_half_tddft.o\

source/module_hamilt_lcao/module_tddft/CMakeLists.txt

@@ -6,6 +6,9 @@ if(ENABLE_LCAO)
         middle_hamilt.cpp
         norm_psi.cpp
         propagator.cpp
+        propagator_cn2.cpp
+        propagator_taylor.cpp
+        propagator_etrs.cpp
         upsi.cpp
         td_velocity.cpp
         td_current.cpp

source/module_hamilt_lcao/module_tddft/propagator.cpp

@@ -14,6 +14,91 @@
 
 namespace module_tddft
 {
+//------------------------ Utility function ------------------------//
+#ifdef __MPI
+inline int globalIndex(int localindex, int nblk, int nprocs, int myproc)
+{
+    int iblock, gIndex;
+    iblock = localindex / nblk;
+    gIndex = (iblock * nprocs + myproc) * nblk + localindex % nblk;
+    return gIndex;
+}
+#endif // __MPI
+
+// Auxiliary function: process non-complex types, return value 1.0
+template <typename T>
+inline T init_value(typename std::enable_if<!std::is_same<T, std::complex<float>>::value
+                                            && !std::is_same<T, std::complex<double>>::value>::type* = nullptr)
+{
+    return T(1.0);
+}
+
+// Auxiliary function: process complex types, return value 1.0 + 0.0i
+template <typename T>
+inline T init_value(typename std::enable_if<std::is_same<T, std::complex<float>>::value
+                                            || std::is_same<T, std::complex<double>>::value>::type* = nullptr)
+{
+    return T(1.0, 0.0);
+}
+
+// Create an identity matrix of size n×n
+template <typename T>
+ct::Tensor create_identity_matrix(const int n, ct::DeviceType device = ct::DeviceType::CpuDevice)
+{
+    // Choose the data type of the Tensor
+    ct::DataType data_type;
+    if (std::is_same<T, float>::value)
+    {
+        data_type = ct::DataType::DT_FLOAT;
+    }
+    else if (std::is_same<T, double>::value)
+    {
+        data_type = ct::DataType::DT_DOUBLE;
+    }
+    else if (std::is_same<T, std::complex<float>>::value)
+    {
+        data_type = ct::DataType::DT_COMPLEX;
+    }
+    else if (std::is_same<T, std::complex<double>>::value)
+    {
+        data_type = ct::DataType::DT_COMPLEX_DOUBLE;
+    }
+    else
+    {
+        static_assert(std::is_same<T, float>::value || std::is_same<T, double>::value
+                          || std::is_same<T, std::complex<float>>::value
+                          || std::is_same<T, std::complex<double>>::value,
+                      "Unsupported data type!");
+    }
+
+    ct::Tensor tensor(data_type, device, ct::TensorShape({n, n}));
+    tensor.zero();
+
+    // Set the diagonal elements to 1
+    if (device == ct::DeviceType::CpuDevice)
+    {
+        // For CPU, we can directly access the data
+        T* data_ptr = tensor.data<T>();
+        for (int i = 0; i < n; ++i)
+        {
+            data_ptr[i * n + i] = init_value<T>();
+        }
+    }
+    else if (device == ct::DeviceType::GpuDevice)
+    {
+        // For GPU, we need to use a kernel to set the diagonal elements
+        T* data_ptr = tensor.data<T>();
+        for (int i = 0; i < n; ++i)
+        {
+            T value = init_value<T>();
+            ct::kernels::set_memory<T, ct::DEVICE_GPU>()(data_ptr + i * n + i, value, 1);
+        }
+    }
+
+    return tensor;
+}
+//------------------------ Utility function ------------------------//
+
 class Propagator
 {
   public: