Merge branch 'develop' into test/source-base-coverage

Author: mohanchen

docs/advanced/acceleration/cuda.md

@@ -44,7 +44,7 @@ The ABACUS program will automatically determine whether the current ELPA support
 ## Run with the GPU support by editing the INPUT script:
 
 In `INPUT` file we need to set the input parameter [device](../input_files/input-main.md#device) to `gpu`. If this parameter is not set, ABACUS will try to determine if there are available GPUs.
-- Set `ks_solver`: For the PW basis, CG, BPCG and Davidson methods are supported on GPU; set the input parameter [ks_solver](../input_files/input-main.md#ks_solver) to `cg`, `bpcg` or `dav`. For the LCAO basis, `cusolver`, `cusolvermp` and `elpa` is supported on GPU.
+- Set `ks_solver`: For the PW basis, CG, BPCG, Davidson, and Davidson subspace methods are supported on GPU; set the input parameter [ks_solver](../input_files/input-main.md#ks_solver) to `cg`, `bpcg`, `dav`, or `dav_subspace`. For the LCAO basis, `cusolver`, `cusolvermp`, and `elpa` are supported on GPU.
 - **single-card**: ABACUS allows for single-GPU acceleration. You can run ABACUS without any MPI process by command `abacus`, and `ks_solver cusolver` is recommended for the LCAO basis. *note: avoid using `mpirun -n 1 abacus`*.
 - **multi-cards**: ABACUS allows for multi-GPU acceleration. If you have multiple GPU cards, you can run ABACUS with several MPI processes, and each process will utilize one GPU card. For example, the command `mpirun -n 2 abacus` will by default launch two GPUs for computation. If you only have one card, this command will only start one GPU. *note: the number of MPI processes SHOULD be equal to the number of GPU cards, unless you are using MPS in your computer.*
 
@@ -58,4 +58,4 @@ PW basis:
 
 LCAO basis:
 - Unless there is a specific reason, avoid using multiple GPUs, as it can be slower than using a single GPU. This is because the generalized eigenvalue solution of the LCAO basis set will incur additional communication overhead when calculated on multiple cards. When the memory limit of a GPU card makes it insufficient to complete the task, it is recommended to use multiple cards for calculation.
-- When using elpa on GPUs, some ELPA internal logs will be output.
+- When using elpa on GPUs, some ELPA internal logs will be output.

docs/advanced/input_files/input-main.md

@@ -8,6 +8,7 @@
   - [System variables](#system-variables)
     - [suffix](#suffix)
     - [ntype](#ntype)
+    - [cell\_replica](#cell_replica)
     - [calculation](#calculation)
     - [esolver\_type](#esolver_type)
     - [symmetry](#symmetry)
@@ -356,6 +357,8 @@
     - [md\_restart](#md_restart)
     - [md\_restartfreq](#md_restartfreq)
     - [md\_dumpfreq](#md_dumpfreq)
+    - [md\_neighbor\_skin](#md_neighbor_skin)
+    - [md\_out\_force](#md_out_force)
     - [dump\_force](#dump_force)
     - [dump\_vel](#dump_vel)
     - [dump\_virial](#dump_virial)
@@ -601,6 +604,12 @@
 - **Description**: Number of different atom species in the calculation.
 - **Default**: 0
 
+### cell_replica
+
+- **Type**: Three Integers
+- **Description**: Replicate the input STRU by Na, Nb, and Nc along its lattice vectors for distributed MDCell workflows. This parameter is only used for classical potentials or machine-learned interatomic potentials. The default is 1 1 1, which preserves the input structure.
+- **Default**: 1 1 1
+
 ### calculation
 
 - **Type**: String
@@ -3504,15 +3513,28 @@
 ### md_restartfreq
 
 - **Type**: Integer
-- **Description**: The output frequency of OUT.{suffix}/STRIU/, which are used to restart molecular dynamics calculations, see md_restart in detail.
+- **Description**: The output frequency of OUT.{suffix}/STRU_MD_*, which are used to restart molecular dynamics calculations, see md_restart in detail. Set to 0 to disable MD restart output.
 - **Default**: 5
 
 ### md_dumpfreq
 
 - **Type**: Integer
-- **Description**: The output frequency of OUT.${suffix}/MD_dump in molecular dynamics calculations, which including the information of lattices and atoms.
+- **Description**: The output frequency of OUT.${suffix}/MD_dump in molecular dynamics calculations, which includes lattice and atomic information. Set to 0 to disable MD_dump output.
 - **Default**: 1
 
+### md_neighbor_skin
+
+- **Type**: Real
+- **Description**: The extra neighbor-list radius in Angstrom for MDCell molecular dynamics. This parameter is only used for classical potentials or machine-learned interatomic potentials. A positive value reuses the cutoff-plus-skin candidate list until an atom has moved by half this distance; 0 rebuilds the list every force evaluation.
+- **Default**: 0.0
+- **Unit**: Angstrom
+
+### md_out_force
+
+- **Type**: Boolean
+- **Description**: Whether to output the TOTAL-FORCE table in OUT.${suffix}/running_md.log for MDCell molecular dynamics. This does not affect force calculation or molecular dynamics integration.
+- **Default**: True
+
 ### dump_force
 
 - **Type**: Boolean
@@ -3535,8 +3557,8 @@
 
 - **Type**: Integer
 - **Description**: The random seed to initialize random numbers used in molecular dynamics calculations.
-  - < 0: No srand() function is called.
-  - >= 0: The function srand(md_seed) is called.
+  - < 0: Each MPI rank uses the default seed 1 plus its rank.
+  - >= 0: Each MPI rank uses md_seed plus its rank.
 - **Default**: -1
 
 ### md_tfreq

docs/parameters.yaml

@@ -17,6 +17,14 @@ parameters:
     default_value: "0"
     unit: ""
     availability: ""
+  - name: cell_replica
+    category: System variables
+    type: Three Integers
+    description: |
+      Replicate the input STRU by Na, Nb, and Nc along its lattice vectors for distributed MDCell workflows. This parameter is only used for classical potentials or machine-learned interatomic potentials. The default is 1 1 1, which preserves the input structure.
+    default_value: 1 1 1
+    unit: ""
+    availability: ""
   - name: calculation
     category: System variables
     type: String
@@ -1449,18 +1457,34 @@ parameters:
     category: Molecular dynamics
     type: Integer
     description: |
-      The output frequency of OUT.{suffix}/STRIU/, which are used to restart molecular dynamics calculations, see md_restart in detail.
+      The output frequency of OUT.{suffix}/STRU_MD_*, which are used to restart molecular dynamics calculations, see md_restart in detail. Set to 0 to disable MD restart output.
     default_value: "5"
     unit: ""
     availability: ""
   - name: md_dumpfreq
     category: Molecular dynamics
     type: Integer
     description: |
-      The output frequency of OUT.${suffix}/MD_dump in molecular dynamics calculations, which including the information of lattices and atoms.
+      The output frequency of OUT.${suffix}/MD_dump in molecular dynamics calculations, which includes lattice and atomic information. Set to 0 to disable MD_dump output.
     default_value: "1"
     unit: ""
     availability: ""
+  - name: md_neighbor_skin
+    category: Molecular dynamics
+    type: Real
+    description: |
+      The extra neighbor-list radius in Angstrom for MDCell molecular dynamics. This parameter is only used for classical potentials or machine-learned interatomic potentials. A positive value reuses the cutoff-plus-skin candidate list until an atom has moved by half this distance; 0 rebuilds the list every force evaluation.
+    default_value: "0.0"
+    unit: Angstrom
+    availability: ""
+  - name: md_out_force
+    category: Molecular dynamics
+    type: Boolean
+    description: |
+      Whether to output the TOTAL-FORCE table in OUT.${suffix}/running_md.log for MDCell molecular dynamics. This does not affect force calculation or molecular dynamics integration.
+    default_value: "True"
+    unit: ""
+    availability: ""
   - name: dump_force
     category: Molecular dynamics
     type: Boolean
@@ -1490,8 +1514,8 @@ parameters:
     type: Integer
     description: |
       The random seed to initialize random numbers used in molecular dynamics calculations.
-      * < 0: No srand() function is called.
-      * >= 0: The function srand(md_seed) is called.
+      * < 0: Each MPI rank uses the default seed 1 plus its rank.
+      * >= 0: Each MPI rank uses md_seed plus its rank.
     default_value: "-1"
     unit: ""
     availability: ""

python/pyabacus/CONTRIBUTING.md

@@ -375,7 +375,7 @@ To contribute to the `pyabacus` project, follow these steps:
 
        @staticmethod
        def sphbes_zeros(l: int, n: int, zeros: NDArray[np.float64]) -> None: 
-           super().sphbes_zeros(l, n, zeros)
+           _Sphbes.sphbes_zeros(l, n, zeros)
    ```
 
 ## Conclusion

python/pyabacus/src/pyabacus/ModuleBase/_module_base.py

@@ -48,7 +48,7 @@ def dsphbesj(self, *args, **kwargs):
 
     @staticmethod
     def sphbes_zeros(l: int, n: int, zeros: NDArray[np.float64]) -> None: 
-        super().sphbes_zeros(l, n, zeros)
+        _Sphbes.sphbes_zeros(l, n, zeros)
 
 class Integral(_Integral):
     def __init__(self) -> None: 
@@ -81,7 +81,7 @@ def Simpson_Integral_0toall(
         rab: NDArray[np.float64], 
         asum: NDArray[np.float64]
     ) -> None: 
-        super().Simpson_Integral_0toall(mesh, func, rab, asum)
+        _Integral.Simpson_Integral_0toall(mesh, func, rab, asum)
 
     @staticmethod
     def Simpson_Integral_alltoinf(
@@ -90,7 +90,7 @@ def Simpson_Integral_alltoinf(
         rab: NDArray[np.float64], 
         asum: NDArray[np.float64]
     ) -> None: 
-        super().Simpson_Integral_alltoinf(mesh, func, rab, asum)
+        _Integral.Simpson_Integral_alltoinf(mesh, func, rab, asum)
 
     @overload
     @staticmethod

python/pyabacus/tests/test_base_math.py

@@ -13,6 +13,23 @@ def test_sphbes():
     assert s.sphbesj(0, 0.0) == 1.0
     assert s.sphbesj(1, np.array([0.0]), 1, 1, np.zeros(1)) == None
 
+def test_static_array_wrappers():
+    zeros = np.zeros(2)
+    base.Sphbes.sphbes_zeros(0, 2, zeros)
+    np.testing.assert_allclose(zeros, [np.pi, 2 * np.pi])
+
+    mesh = 3
+    func = np.ones(mesh)
+    rab = np.ones(mesh)
+    integral_from_zero = np.zeros(mesh)
+    integral_to_infinity = np.zeros(mesh)
+
+    base.Integral.Simpson_Integral_0toall(mesh, func, rab, integral_from_zero)
+    base.Integral.Simpson_Integral_alltoinf(mesh, func, rab, integral_to_infinity)
+
+    np.testing.assert_allclose(integral_from_zero, [0.0, 1.0, 2.0])
+    np.testing.assert_allclose(integral_to_infinity, [2.0, 1.0, 0.0])
+
 def test_sbt():
     sbt = base.SphericalBesselTransformer()
 

source/Makefile.Objects

@@ -772,7 +772,7 @@ OBJS_MODULE_RI=conv_coulomb_pot_k.o\
       singular_value.o\
 
 OBJS_PARALLEL=parallel_common.o\
-    communication_domain.o\
+    parallel_cell.o\
     parallel_global.o\
     parallel_comm.o\
     parallel_grid.o\

source/source_base/CMakeLists.txt

@@ -51,7 +51,7 @@ add_library(
     tool_title.cpp
     ylm.cpp
     parallel_common.cpp
-    communication_domain.cpp
+    parallel_cell.cpp
     parallel_global.cpp
     parallel_comm.cpp
     parallel_reduce.cpp

source/source_base/module_container/ATen/core/tensor.h

@@ -390,25 +390,30 @@ class Tensor {
      *
      * @return The pointer to the specified row.
      *
-     * @note This function assumes the tensor is treated as a matrix, where each row
-     * is a contiguous block of memory.
-     * If the row index is out of bounds, the behavior is undefined.
+     * @note This function supports rank-1 and rank-2 tensors. For a rank-1 tensor,
+     * the index selects an element. For a rank-2 tensor, the index selects a row.
+     * An invalid rank or index throws std::invalid_argument.
      */
     template <typename T>
-    T* inner_most_ptr(const int &index) const {
-        if (shape_.ndim() > 2) {
-            throw std::invalid_argument("Invalid call, inner_most_ptr only support tensor rank <= 2!");
+    T* inner_most_ptr(const int& index) const
+    {
+        const unsigned int rank = shape_.ndim();
+        if (rank == 0 || rank > 2)
+        {
+            throw std::invalid_argument("Invalid call, inner_most_ptr only supports tensor ranks 1 and 2!");
         }
-        if (index > shape_.dim_size(static_cast<int>(shape_.ndim() - 2))) {
-            throw std::invalid_argument("Invalid index, index of the inner-most must less than the inner-most shape size!");
+        const int64_t outer_size = shape_.dim_size(0);
+        if (index < 0 || static_cast<int64_t>(index) >= outer_size)
+        {
+            throw std::invalid_argument("Invalid index, inner_most_ptr index is out of bounds!");
         }
-        if (shape_.ndim() == 1) {
+        if (rank == 1)
+        {
             return data<T>() + index;
         }
-        return data<T>() + index * shape_.dim_size(static_cast<int>(shape_.ndim()) - 1);
+        return data<T>() + index * shape_.dim_size(1);
     }
 
-
     /**
      * @brief Equality comparison operator for tensors.
      *

source/source_base/module_container/test/tensor_test.cpp

@@ -230,6 +230,27 @@ TEST(Tensor, GetValueAndInnerMostPtr) {
     EXPECT_EQ(row_ptr[3], 12);
 }
 
+TEST(Tensor, InnerMostPtrBounds)
+{
+    container::Tensor vector(container::DataType::DT_INT, container::DeviceType::CpuDevice, {4});
+    std::vector<int> values = {1, 2, 3, 4};
+    memcpy(vector.data<int>(), values.data(), sizeof(int) * values.size());
+
+    auto element_ptr = vector.inner_most_ptr<int>(2);
+    EXPECT_EQ(*element_ptr, 3);
+    EXPECT_THROW(vector.inner_most_ptr<int>(-1), std::invalid_argument);
+    EXPECT_THROW(vector.inner_most_ptr<int>(4), std::invalid_argument);
+
+    container::Tensor matrix(container::DataType::DT_INT, container::DeviceType::CpuDevice, {2, 2});
+    EXPECT_THROW(matrix.inner_most_ptr<int>(2), std::invalid_argument);
+
+    container::Tensor empty(container::DataType::DT_INT, container::DeviceType::CpuDevice, container::TensorShape());
+    EXPECT_THROW(empty.inner_most_ptr<int>(0), std::invalid_argument);
+
+    container::Tensor rank_three(container::DataType::DT_INT, container::DeviceType::CpuDevice, {1, 1, 1});
+    EXPECT_THROW(rank_three.inner_most_ptr<int>(0), std::invalid_argument);
+}
+
 TEST(Tensor, ReshapeDeathTest) {
     ::testing::FLAGS_gtest_death_test_style = "threadsafe";
     container::Tensor t(container::DataType::DT_FLOAT, container::DeviceType::CpuDevice, {2, 3, 4});
@@ -545,4 +566,4 @@ TEST(Tensor, Accessor) {
     }
 }
 
-} // namespace container
+} // namespace container