diff --git a/python/pyabacus/CMakeLists.txt b/python/pyabacus/CMakeLists.txt index e21f5e1446..5baa420ce8 100644 --- a/python/pyabacus/CMakeLists.txt +++ b/python/pyabacus/CMakeLists.txt @@ -1,18 +1,20 @@ cmake_minimum_required(VERSION 3.15...3.26) +# project settings project( ${SKBUILD_PROJECT_NAME} VERSION ${SKBUILD_PROJECT_VERSION} LANGUAGES CXX) +# find python and pybind11 find_package(Python REQUIRED COMPONENTS Interpreter Development.Module) find_package(pybind11 CONFIG REQUIRED) +# set source path set(ABACUS_SOURCE_DIR "${PROJECT_SOURCE_DIR}/../../source") set(BASE_PATH "${ABACUS_SOURCE_DIR}/module_base") set(NAO_PATH "${ABACUS_SOURCE_DIR}/module_basis/module_nao") set(HSOLVER_PATH "${ABACUS_SOURCE_DIR}/module_hsolver") -set(HAMILT_PATH "${ABACUS_SOURCE_DIR}/module_hamilt_general") set(PSI_PATH "${ABACUS_SOURCE_DIR}/module_psi") set(ENABLE_LCAO ON) list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/../../cmake") @@ -63,11 +65,14 @@ else() endif() endif() -include_directories(${BASE_PATH} +# add include directories +include_directories( + ${BASE_PATH} ${ABACUS_SOURCE_DIR} - ${ABACUS_SOURCE_DIR}/module_base/module_container) + ${ABACUS_SOURCE_DIR}/module_base/module_container + ) -# add library +# add basic libraries set(CMAKE_POSITION_INDEPENDENT_CODE ON) # add base set(BASE_BINARY_DIR "${PROJECT_SOURCE_DIR}/build/base") @@ -80,83 +85,6 @@ set(ORB_BINARY_DIR "${PROJECT_SOURCE_DIR}/build/orb") add_subdirectory(${ABACUS_SOURCE_DIR}/module_basis/module_ao ${ORB_BINARY_DIR}) # set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE) -# add nao shared library -list(APPEND _naos - ${NAO_PATH}/atomic_radials.cpp - ${NAO_PATH}/beta_radials.cpp - ${NAO_PATH}/hydrogen_radials.cpp - ${NAO_PATH}/numerical_radial.cpp - ${NAO_PATH}/pswfc_radials.cpp - ${NAO_PATH}/radial_collection.cpp - ${NAO_PATH}/radial_set.cpp - ${NAO_PATH}/real_gaunt_table.cpp - ${NAO_PATH}/sphbes_radials.cpp - ${NAO_PATH}/two_center_bundle.cpp - ${NAO_PATH}/two_center_integrator.cpp - ${NAO_PATH}/two_center_table.cpp - ${NAO_PATH}/projgen.cpp - # dependency - ${ABACUS_SOURCE_DIR}/module_base/kernels/math_op.cpp - # ${ABACUS_SOURCE_DIR}/module_psi/kernels/psi_memory_op.cpp - ${ABACUS_SOURCE_DIR}/module_base/module_device/memory_op.cpp - ) -add_library(naopack SHARED - ${_naos} - ) -# add diago shared library -list(APPEND _diago - ${HSOLVER_PATH}/diago_dav_subspace.cpp - ${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/math_kernel_op.cpp - # dependency - ${BASE_PATH}/module_device/device.cpp - ${BASE_PATH}/module_device/memory_op.cpp - - ${HAMILT_PATH}/operator.cpp - ${PSI_PATH}/psi.cpp - ) -add_library(diagopack SHARED - ${_diago} - ) -target_link_libraries(diagopack - PRIVATE - ${OpenBLAS_LIBRARIES} - ${LAPACK_LIBRARIES} - ) -# link math_libs -if(MKLROOT) - target_link_libraries(naopack - base - parameter - container - orb - ${math_libs} - MPI::MPI_CXX - OpenMP::OpenMP_CXX - ) -else() - target_link_libraries(naopack - base - parameter - container - orb - ${math_libs} - ) -endif() -# list(APPEND _sources ${_naos} ${_bases}) -list(APPEND _sources - ${PROJECT_SOURCE_DIR}/src/py_abacus.cpp - ${PROJECT_SOURCE_DIR}/src/py_base_math.cpp - ${PROJECT_SOURCE_DIR}/src/py_m_nao.cpp - ${PROJECT_SOURCE_DIR}/src/py_hsolver.cpp - ) -pybind11_add_module(_core MODULE ${_sources}) -target_link_libraries(_core PRIVATE pybind11::headers naopack diagopack) -target_compile_definitions(_core PRIVATE VERSION_INFO=${PROJECT_VERSION}) # set RPATH execute_process( COMMAND "${PYTHON_EXECUTABLE}" -c "from distutils.sysconfig import get_python_lib; print(get_python_lib())" @@ -167,9 +95,9 @@ execute_process( # set package name to pyabacus set(TARGET_PACK pyabacus) set(CMAKE_INSTALL_RPATH "${PYTHON_SITE_PACKAGES}/${TARGET_PACK}") -set_target_properties(_core PROPERTIES INSTALL_RPATH "$ORIGIN") -set_target_properties(naopack PROPERTIES INSTALL_RPATH "$ORIGIN") -set_target_properties(diagopack PROPERTIES INSTALL_RPATH "$ORIGIN") -install(TARGETS _core naopack DESTINATION ${TARGET_PACK}) -install(TARGETS _core diagopack DESTINATION ${TARGET_PACK}) + +# add subdirectories for submodules +add_subdirectory(${PROJECT_SOURCE_DIR}/src/hsolver) +add_subdirectory(${PROJECT_SOURCE_DIR}/src/ModuleBase) +add_subdirectory(${PROJECT_SOURCE_DIR}/src/ModuleNAO) diff --git a/python/pyabacus/CONTRIBUTING.md b/python/pyabacus/CONTRIBUTING.md new file mode 100644 index 0000000000..fbd23ad9ff --- /dev/null +++ b/python/pyabacus/CONTRIBUTING.md @@ -0,0 +1,376 @@ +# Developer Guide + +## Introduction + +Welcome to the `pyabacus` project! This document provides guidelines and instructions for developers who want to contribute to this project. + +`pyabacus` is a Python interface for the ABACUS package. It provides a high-level Python API for interacting with the ABACUS library, allowing users to perform electronic structure calculations and analyze the results using Python. + + + +- [Project structure](#project-structure) + - [Root CMake Configuration](#root-cmake-configuration) + - [Module CMake Configuration](#module-cmake-configuration) +- [Development Process](#development-process) + + + +**If you are new to the project**, please refer to the [README.md](./README.md) file for an overview of the project and its goals. + +**If you are already familiar with the project and want to contribute**, this guide will help you understand the project structure, development process, and best practices for contributing code. + +**If you have any questions or need help**, feel free to reach out to the maintainers or create an issue in the repository. + +**Please feel free to contribute to this guide** by submitting a pull request with any improvements or additional information. + +Let's get started! + +## Project Structure + +The project is organized as follows: + +``` +pyabacus/ +├── CMakeLists.txt +└── src + ├── pyabacus + │ └── {your_module} + │ ├── {interface}.py + │ └── __init__.py + └── {your_module} + ├── {your_code}.cpp + └── CMakeLists.txt +``` + +Our project is built using [pybind11](http://github.com/pybind/pybind11) and [scikit-build-core](https://scikit-build-core.readthedocs.io/) for facilitating the `CMake` build toolchain. So the `CMakeLists.txt` configuration is the key to thoroughly understanding the project structure. + +### Root CMake Configuration + +The `CMakeLists.txt` in root directory is the main configuration file for the pyabacus project. It sets up the project, finds necessary dependencies, configures build options, and includes subdirectories for different modules. Below is a detailed explanation of each section of the file: + +```cmake +cmake_minimum_required(VERSION 3.15...3.26) + +# Project settings +project( + ${SKBUILD_PROJECT_NAME} + VERSION ${SKBUILD_PROJECT_VERSION} + LANGUAGES CXX) +``` +- This section sets the project name, version, and the programming languages used (C++ in this case). The project name and version are obtained from the `SKBUILD_PROJECT_NAME` and `SKBUILD_PROJECT_VERSION` variables, respectively. + +```cmake +# Find Python and pybind11 +find_package(Python REQUIRED COMPONENTS Interpreter Development.Module) +find_package(pybind11 CONFIG REQUIRED) +``` +- This section finds the required Python and pybind11 packages. + +```cmake +# Set source path +set(ABACUS_SOURCE_DIR "${PROJECT_SOURCE_DIR}/../../source") +set(BASE_PATH "${ABACUS_SOURCE_DIR}/module_base") +set(NAO_PATH "${ABACUS_SOURCE_DIR}/module_basis/module_nao") +set(HSOLVER_PATH "${ABACUS_SOURCE_DIR}/module_hsolver") +set(PSI_PATH "${ABACUS_SOURCE_DIR}/module_psi") +set(ENABLE_LCAO ON) +list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/../../cmake") +``` +- This section sets various source paths and configuration options. It defines the paths to different modules and appends the custom CMake module path. + +```cmake +# Add math_libs +if(DEFINED ENV{MKLROOT} AND NOT DEFINED MKLROOT) + set(MKLROOT "$ENV{MKLROOT}") +endif() +if(MKLROOT) + set(MKL_INTERFACE lp64) + set(ENABLE_MPI ON) + if (ENABLE_MPI) + find_package(MPI REQUIRED) + include_directories(${MPI_CXX_INCLUDE_PATH}) + endif() + + set(USE_OPENMP ON) + if(USE_OPENMP) + find_package(OpenMP REQUIRED) + set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") + add_link_options(${OpenMP_CXX_LIBRARIES}) + endif() + find_package(MKL REQUIRED) + add_definitions(-D__MKL) + include_directories(${MKL_INCLUDE} ${MKL_INCLUDE}/fftw) + + if(NOT ENABLE_DEEPKS) + list(APPEND math_libs IntelMKL::MKL) + endif() + + if(CMAKE_CXX_COMPILER_ID MATCHES Intel) + list(APPEND math_libs -lifcore) + endif() +else() + find_package(FFTW3 REQUIRED) + add_compile_definitions(__FFTW3) + find_package(LAPACK REQUIRED) + include_directories(${FFTW3_INCLUDE_DIRS}) + list(APPEND math_libs FFTW3::FFTW3 LAPACK::LAPACK) + + if(ENABLE_LCAO) + find_package(ScaLAPACK REQUIRED) + list(APPEND math_libs ScaLAPACK::ScaLAPACK) + endif() +endif() +``` +- This section configures the math libraries. It checks for the presence of the Intel Math Kernel Library (MKL) and configures it if available. If MKL is not available, it falls back to using FFTW3 and LAPACK. It also configures MPI and OpenMP if enabled. + +```cmake +# Add include directories +include_directories( + ${BASE_PATH} + ${ABACUS_SOURCE_DIR} + ${ABACUS_SOURCE_DIR}/module_base/module_container + ) +``` +- This section adds the necessary include directories for the project. + +```cmake +# Add basic libraries +set(CMAKE_POSITION_INDEPENDENT_CODE ON) +# Add base +set(BASE_BINARY_DIR "${PROJECT_SOURCE_DIR}/build/base") +add_subdirectory(${ABACUS_SOURCE_DIR}/module_base ${BASE_BINARY_DIR}) +# Add parameter +set(PARAMETER_BINARY_DIR "${PROJECT_SOURCE_DIR}/build/parameter") +add_subdirectory(${ABACUS_SOURCE_DIR}/module_parameter ${PARAMETER_BINARY_DIR}) +# Add orb +set(ORB_BINARY_DIR "${PROJECT_SOURCE_DIR}/build/orb") +add_subdirectory(${ABACUS_SOURCE_DIR}/module_basis/module_ao ${ORB_BINARY_DIR}) +``` +- This section sets the position-independent code flag and adds subdirectories for the base, parameter, and orb modules. It specifies the build directories for these modules. + +```cmake +# Set RPATH +execute_process( + COMMAND "${PYTHON_EXECUTABLE}" -c "from distutils.sysconfig import get_python_lib; print(get_python_lib())" + OUTPUT_VARIABLE PYTHON_SITE_PACKAGES + OUTPUT_STRIP_TRAILING_WHITESPACE +) +``` +- This section sets the runtime search path (RPATH) for the Python site-packages directory. It uses a Python command to get the site-packages path and stores it in the `PYTHON_SITE_PACKAGES` variable. + +```cmake +# Set package name to pyabacus +set(TARGET_PACK pyabacus) +set(CMAKE_INSTALL_RPATH "${PYTHON_SITE_PACKAGES}/${TARGET_PACK}") +``` +- This section sets the package name to `pyabacus` and configures the install RPATH to include the Python site-packages directory. + +```cmake +# Add subdirectories for submodules +add_subdirectory(${PROJECT_SOURCE_DIR}/src/hsolver) +add_subdirectory(${PROJECT_SOURCE_DIR}/src/ModuleBase) +add_subdirectory(${PROJECT_SOURCE_DIR}/src/ModuleNAO) +``` +- This section adds subdirectories for modules. Each subdirectory contains its own `CMakeLists.txt` for further configuration. + +By following this structure, the `CMakeLists.txt` file ensures that all necessary dependencies are found, configured, and included in the build process. It also sets up the project environment and includes submodules for different components of the `pyabacus` project. + +### Module CMake Configuration + +I'll show you a `CMakeLists.txt` for example (pyabacus.hsolver) + +```cmake +# Add diago shared library +list(APPEND _diago + ${HSOLVER_PATH}/diago_dav_subspace.cpp + ${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/math_kernel_op.cpp + ${BASE_PATH}/kernels/math_op.cpp + ${BASE_PATH}/module_device/device.cpp + ${BASE_PATH}/module_device/memory_op.cpp + ${PSI_PATH}/psi.cpp +) +add_library(diagopack SHARED ${_diago}) +target_link_libraries(diagopack + base + parameter + container + orb + ${math_libs} + ${OpenBLAS_LIBRARIES} + ${LAPACK_LIBRARIES} +) + +list(APPEND pymodule_hsolver + ${PROJECT_SOURCE_DIR}/src/hsolver/py_hsolver.cpp +) + +# Use pybind11 to add python module +pybind11_add_module(_hsolver_pack MODULE ${pymodule_hsolver}) +# Link your dependencies and pybind11 libraries to your module +target_link_libraries(_hsolver_pack PRIVATE pybind11::headers diagopack) +target_compile_definitions(_hsolver_pack PRIVATE VERSION_INFO=${PROJECT_VERSION}) + +set_target_properties(diagopack PROPERTIES INSTALL_RPATH "$ORIGIN") +set_target_properties(_hsolver_pack PROPERTIES INSTALL_RPATH "$ORIGIN") + +# Install your module package to destination path +install(TARGETS _hsolver_pack diagopack DESTINATION ${TARGET_PACK}/hsolver) +``` + +You can refer to the `CMakeLists.txt` files in other modules for guidance on how to configure your module. + +## Development Process + +To contribute to the `pyabacus` project, follow these steps: + +1. **Check the issues**: + - Look for issues to ensure that you are not working on something that is already in progress. + - If you want to work on a new feature or bug fix, create an issue first to discuss it with the maintainers. + +2. **Create a new folder for your module**: + - If you want to add a new module with pure Python code, create a new folder in the `src/pyabacus` directory. + - If you want to add a new module with C++ code, create a new folder in the `src` directory and a corresponding directory in the `src/pyabacus` directory. + +3. **Write source code using pybind11**: + - Follow the structure of other modules. + - Manage dependencies and installation paths in the `CMakeLists.txt` file. + +3. **Modify `src/pyabacus/__init__.py`**: + - Add the name of your module to the `__submodules__` list and import the module in the `__getattr__` function. + + ```python + from __future__ import annotations + + __submodules__ = ["ModuleBase", "ModuleNAO", "hsolver", "{module_name}"] + + __all__ = list(__submodules__) + + def __getattr__(attr): + if attr == "ModuleBase": + import pyabacus.ModuleBase as ModuleBase + return ModuleBase + elif attr == "ModuleNAO": + import pyabacus.ModuleNAO as ModuleNAO + return ModuleNAO + elif attr == "hsolver": + import pyabacus.hsolver as hsolver + return hsolver + elif attr == '{module_name}': + import pyabacus.{module_name} as {module_name} + return {module_name} + else: + raise AttributeError(f"module {__name__} has no attribute {attr}") + ``` + +4. **Create two files in `src/pyabacus/{module_name}`**: + - `__init__.py`: This file allows Python to recognize the folder as a module. + - `_{module_name}.py`: This file is responsible for designing the Python interface (frontend). + + **Example `__init__.py`**: + + ```python + from __future__ import annotations + from ._{module_name} import * + + __all__ = ["{class_name}", "{func_name}", ...] + ``` + + **Example `_{module_name}.py`**: + + ```python + from .{module_library_name} import {your_class} as _your_class, ... + + """ + Your class should inherit from the corresponding class in the C++ library. + All methods should be overridden to provide type hints and auto-completion. + You can use the `super()` method to call the base class(C++ class) methods. + """ + class {your_class}(_your_class): + def __init__(self) -> None: + super().__init__() + + def foo(self, arg1, arg2, ...) -> RetType: + return super().foo(arg1, arg2, ...) + + def bar(self, arg1, arg2, ...): + super().bar(arg1, arg2, ...) + ``` + + For a class, if you do not declare the interface in the frontend, the IDE will not provide type hints and auto-completion. However, if the interface name matches the name binding in pybind11, it will be overridden. To address this, you can use the method as shown above. + +5. **Handle overloaded functions in C++**: + - Since Python does not support function overloading with different parameters, use the following method: + + ```python + @overload + def foo(self, x: float) -> float: ... + @overload + def foo(self, n: int, x: float, y: float) -> float: ... + + def foo(self, *args, **kwargs): + return super().foo(*args, **kwargs) + ``` + +**Example Python Interface**: + + ```python + class diag_comm_info(_diag_comm_info): + def __init__(self, rank: int, nproc: int): + super().__init__(rank, nproc) + + @property + def rank(self) -> int: + return super().rank + + @property + def nproc(self) -> int: + return super().nproc + + class Sphbes(_Sphbes): + def __init__(self) -> None: + super().__init__() + + @overload + @staticmethod + def sphbesj(l: int, x: float) -> float: ... + @overload + @staticmethod + def sphbesj( + n: int, + r: NDArray[np.float64], + q: int, + l: int, + jl: NDArray[np.float64] + ) -> None: ... + + def sphbesj(self, *args, **kwargs): + return super().sphbesj(*args, **kwargs) + + @overload + @staticmethod + def dsphbesj(l: int, x: float) -> float: ... + @overload + @staticmethod + def dsphbesj( + n: int, + r: NDArray[np.float64], + q: int, + l: int, + djl: NDArray[np.float64] + ) -> None: ... + + def dsphbesj(self, *args, **kwargs): + return super().dsphbesj(*args, **kwargs) + + @staticmethod + def sphbes_zeros(l: int, n: int, zeros: NDArray[np.float64]) -> None: + super().sphbes_zeros(l, n, zeros) + ``` + +## Conclusion + +By following this guide, you can effectively contribute to the `pyabacus` project. Ensure that you follow the structure and conventions outlined here to maintain consistency and readability in the codebase. Happy coding! diff --git a/python/pyabacus/README.md b/python/pyabacus/README.md index f751e52c1b..5634bf8cf2 100644 --- a/python/pyabacus/README.md +++ b/python/pyabacus/README.md @@ -1,19 +1,32 @@ -Build Example: TwoCenterIntegral Section in ABACUS -================================================== +# pyabacus: a Python interface for the ABACUS package -An example project built with [pybind11](https://github.com/pybind/pybind11) -and scikit-build-core. Python 3.7+ (see older commits for older versions of -Python). +`pyabacus` is a Python interface for the ABACUS package, which provides a high-level Python API for interacting with the `ABACUS` library. -Installation ------------- +This project is built using [pybind11](http://github.com/pybind/pybind11) and [scikit-build-core](https://scikit-build-core.readthedocs.io/), so you can easily build the project and use it in your Python environment. + +Now, `pyabacus` provides the following modules: +- `io`: a module for input/output in pyabacus. +- `Cell`: a module for the cell structure to bridge `ModuleNAO` in python module with users' input. +- `ModuleBase`: a module for basic math functions. +- `ModuleNAO`: a module for numerical atomic orbitals (NAO). +- `hsolver`: a module for solving the Hamiltonian. + + + +- [Installation](#installation) +- [CI Examples](#ci-examples) +- [License](#license) +- [Test call](#test-call) + + + +## Installation - Create and activate a new conda env, e.g. `conda create -n myenv python=3.8 & conda activate myenv`. - Clone ABACUS main repository and `cd abacus-develop/python/pyabacus`. - Build pyabacus by `pip install -v .` or install test dependencies & build pyabacus by `pip install .[test]`. (Use `pip install -v .[test] -i https://pypi.tuna.tsinghua.edu.cn/simple` to accelerate installation process.) -CI Examples ------------ +## CI Examples There are examples for CI in `.github/workflows`. A simple way to produces binary "wheels" for all platforms is illustrated in the "wheels.yml" file, @@ -71,15 +84,13 @@ eigenvalues difference: -8.94295749e-12 4.71351846e-11 5.39378986e-10 1.97244101e-08] ``` -License -------- +## License pybind11 is provided under a BSD-style license that can be found in the LICENSE file. By using, distributing, or contributing to this project, you agree to the terms and conditions of this license. -Test call ---------- +## Test call ```python import pyabacus as m diff --git a/python/pyabacus/examples/diago_matrix.py b/python/pyabacus/examples/diago_matrix.py index b7dba3865c..6a3d1eb75a 100644 --- a/python/pyabacus/examples/diago_matrix.py +++ b/python/pyabacus/examples/diago_matrix.py @@ -50,7 +50,7 @@ def calc_eig_scipy(mat_file): if __name__ == '__main__': mat_file = './Si2.mat' - method = ['davidson', 'dav_subspace'] + method = ['dav_subspace', 'davidson'] for m in method: print(f'\n====== Calculating eigenvalues using {m} method... ======') diff --git a/python/pyabacus/src/ModuleBase/CMakeLists.txt b/python/pyabacus/src/ModuleBase/CMakeLists.txt new file mode 100644 index 0000000000..f150cf1b5e --- /dev/null +++ b/python/pyabacus/src/ModuleBase/CMakeLists.txt @@ -0,0 +1,25 @@ +list(APPEND pymodule_base + ${PROJECT_SOURCE_DIR}/src/ModuleBase/py_base_math.cpp + ${BASE_PATH}/kernels/math_op.cpp + ${BASE_PATH}/module_device/memory_op.cpp + ) + +pybind11_add_module(_base_pack MODULE ${pymodule_base}) + +target_link_libraries(_base_pack + PRIVATE + base + parameter + container + orb + ${math_libs} + ${OpenBLAS_LIBRARIES} + ${LAPACK_LIBRARIES} + ) + +target_link_libraries(_base_pack PRIVATE pybind11::headers) +target_compile_definitions(_base_pack PRIVATE VERSION_INFO=${PROJECT_VERSION}) + +set_target_properties(_base_pack PROPERTIES INSTALL_RPATH "$ORIGIN") + +install(TARGETS _base_pack DESTINATION ${TARGET_PACK}/ModuleBase) \ No newline at end of file diff --git a/python/pyabacus/src/py_base_math.cpp b/python/pyabacus/src/ModuleBase/py_base_math.cpp similarity index 97% rename from python/pyabacus/src/py_base_math.cpp rename to python/pyabacus/src/ModuleBase/py_base_math.cpp index 5f87831477..b9803c81cd 100644 --- a/python/pyabacus/src/py_base_math.cpp +++ b/python/pyabacus/src/ModuleBase/py_base_math.cpp @@ -12,10 +12,8 @@ using overload_cast_ = pybind11::detail::overload_cast_impl; void bind_base_math(py::module& m) { - py::module module_base = m.def_submodule("ModuleBase"); - // python binding for class Sphbes - py::class_(module_base, "Sphbes") + py::class_(m, "Sphbes") .def(py::init<>()) .def_static("sphbesj", overload_cast_()(&ModuleBase::Sphbes::sphbesj), "l"_a, "x"_a) .def_static("dsphbesj", overload_cast_()(&ModuleBase::Sphbes::dsphbesj), "l"_a, "x"_a) @@ -65,7 +63,7 @@ void bind_base_math(py::module& m) }); // python binding for class Integral - py::class_(module_base, "Integral") + py::class_(m, "Integral") .def(py::init<>()) .def_static("Simpson_Integral", [](const int mesh, py::array_t func, py::array_t rab, double asum) { py::buffer_info func_info = func.request(); @@ -220,6 +218,13 @@ void bind_base_math(py::module& m) static_cast(x_info.ptr), static_cast(w_info.ptr)); }); - py::class_(module_base, "SphericalBesselTransformer") - .def(py::init<>()); + py::class_(m, "SphericalBesselTransformer") + .def(py::init<>()); +} + +PYBIND11_MODULE(_base_pack, m) +{ + m.doc() = "Submodule for pyabacus: ModuleBase"; + + bind_base_math(m); } \ No newline at end of file diff --git a/python/pyabacus/src/ModuleNAO/CMakeLists.txt b/python/pyabacus/src/ModuleNAO/CMakeLists.txt new file mode 100644 index 0000000000..65e209ca40 --- /dev/null +++ b/python/pyabacus/src/ModuleNAO/CMakeLists.txt @@ -0,0 +1,52 @@ +# add nao shared library +list(APPEND _naos + ${NAO_PATH}/atomic_radials.cpp + ${NAO_PATH}/beta_radials.cpp + ${NAO_PATH}/hydrogen_radials.cpp + ${NAO_PATH}/numerical_radial.cpp + ${NAO_PATH}/pswfc_radials.cpp + ${NAO_PATH}/radial_collection.cpp + ${NAO_PATH}/radial_set.cpp + ${NAO_PATH}/real_gaunt_table.cpp + ${NAO_PATH}/sphbes_radials.cpp + ${NAO_PATH}/two_center_bundle.cpp + ${NAO_PATH}/two_center_integrator.cpp + ${NAO_PATH}/two_center_table.cpp + ${NAO_PATH}/projgen.cpp + # dependency + ${ABACUS_SOURCE_DIR}/module_base/kernels/math_op.cpp + # ${ABACUS_SOURCE_DIR}/module_psi/kernels/psi_memory_op.cpp + ${ABACUS_SOURCE_DIR}/module_base/module_device/memory_op.cpp +) +add_library(naopack SHARED + ${_naos} +) + +# link math_libs +if(MKLROOT) + target_link_libraries(naopack + MPI::MPI_CXX + OpenMP::OpenMP_CXX + ) +endif() + +target_link_libraries(naopack + base + parameter + container + orb + ${math_libs} +) + +list(APPEND pymodule_nao + ${PROJECT_SOURCE_DIR}/src/ModuleNAO/py_m_nao.cpp +) + +pybind11_add_module(_nao_pack MODULE ${pymodule_nao}) +target_link_libraries(_nao_pack PRIVATE pybind11::headers naopack) +target_compile_definitions(_nao_pack PRIVATE VERSION_INFO=${PROJECT_VERSION}) + +set_target_properties(naopack PROPERTIES INSTALL_RPATH "$ORIGIN") +set_target_properties(_nao_pack PROPERTIES INSTALL_RPATH "$ORIGIN") + +install(TARGETS _nao_pack naopack DESTINATION ${TARGET_PACK}/ModuleNAO) \ No newline at end of file diff --git a/python/pyabacus/src/py_m_nao.cpp b/python/pyabacus/src/ModuleNAO/py_m_nao.cpp similarity index 98% rename from python/pyabacus/src/py_m_nao.cpp rename to python/pyabacus/src/ModuleNAO/py_m_nao.cpp index d7c7283614..92c5084d44 100644 --- a/python/pyabacus/src/py_m_nao.cpp +++ b/python/pyabacus/src/ModuleNAO/py_m_nao.cpp @@ -13,12 +13,8 @@ using overload_cast_ = pybind11::detail::overload_cast_impl; void bind_m_nao(py::module& m) { - // Create the submodule for Module NAO - py::module m_nao = m.def_submodule("ModuleNAO"); - m_nao.doc() = "Module for Numerical Atomic Orbitals (NAO) in ABACUS"; - // Bind the RadialCollection class - py::class_(m_nao, "RadialCollection") + py::class_(m, "RadialCollection") .def(py::init<>(), R"pbdoc( A class that holds all numerical radial functions of the same kind. @@ -95,7 +91,7 @@ void bind_m_nao(py::module& m) .def("nchi", overload_cast_()(&RadialCollection::nchi, py::const_), "itype"_a) .def("nchi", overload_cast_<>()(&RadialCollection::nchi, py::const_)); // Bind the TwoCenterIntegrator class - py::class_(m_nao, "TwoCenterIntegrator") + py::class_(m, "TwoCenterIntegrator") .def(py::init<>(), R"pbdoc( A class to compute two-center integrals. @@ -261,7 +257,7 @@ void bind_m_nao(py::module& m) "pvR"_a, "deriv"_a = false); // Bind the NumericalRadial class - py::class_(m_nao, "NumericalRadial") + py::class_(m, "NumericalRadial") .def(py::init<>(), R"pbdoc( A class that represents a numerical radial function. @@ -493,4 +489,11 @@ void bind_m_nao(py::module& m) return py::array_t(self.nk(), kvalue); }) .def_property_readonly("is_fft_compliant", overload_cast_<>()(&NumericalRadial::is_fft_compliant, py::const_)); +} + +PYBIND11_MODULE(_nao_pack, m) +{ + m.doc() = "Module for Numerical Atomic Orbitals (NAO) in ABACUS"; + + bind_m_nao(m); } \ No newline at end of file diff --git a/python/pyabacus/src/hsolver/CMakeLists.txt b/python/pyabacus/src/hsolver/CMakeLists.txt new file mode 100644 index 0000000000..6bf690ecec --- /dev/null +++ b/python/pyabacus/src/hsolver/CMakeLists.txt @@ -0,0 +1,41 @@ +# add diago shared library +list(APPEND _diago + ${HSOLVER_PATH}/diago_dav_subspace.cpp + ${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/math_kernel_op.cpp + # dependency + ${BASE_PATH}/kernels/math_op.cpp + ${BASE_PATH}/module_device/device.cpp + ${BASE_PATH}/module_device/memory_op.cpp + + ${PSI_PATH}/psi.cpp + ) +add_library(diagopack SHARED + ${_diago} + ) +target_link_libraries(diagopack + base + parameter + container + orb + ${math_libs} + ${OpenBLAS_LIBRARIES} + ${LAPACK_LIBRARIES} + ) + +list(APPEND pymodule_hsolver + ${PROJECT_SOURCE_DIR}/src/hsolver/py_hsolver.cpp + ) + +pybind11_add_module(_hsolver_pack MODULE ${pymodule_hsolver}) +target_link_libraries(_hsolver_pack PRIVATE pybind11::headers diagopack) +target_compile_definitions(_hsolver_pack PRIVATE VERSION_INFO=${PROJECT_VERSION}) + +set_target_properties(diagopack PROPERTIES INSTALL_RPATH "$ORIGIN") +set_target_properties(_hsolver_pack PROPERTIES INSTALL_RPATH "$ORIGIN") + +install(TARGETS _hsolver_pack diagopack DESTINATION ${TARGET_PACK}/hsolver) \ No newline at end of file diff --git a/python/pyabacus/src/py_diago_dav_subspace.hpp b/python/pyabacus/src/hsolver/py_diago_dav_subspace.hpp similarity index 100% rename from python/pyabacus/src/py_diago_dav_subspace.hpp rename to python/pyabacus/src/hsolver/py_diago_dav_subspace.hpp diff --git a/python/pyabacus/src/py_diago_david.hpp b/python/pyabacus/src/hsolver/py_diago_david.hpp similarity index 100% rename from python/pyabacus/src/py_diago_david.hpp rename to python/pyabacus/src/hsolver/py_diago_david.hpp diff --git a/python/pyabacus/src/py_hsolver.cpp b/python/pyabacus/src/hsolver/py_hsolver.cpp similarity index 92% rename from python/pyabacus/src/py_hsolver.cpp rename to python/pyabacus/src/hsolver/py_hsolver.cpp index e1c65d876e..2401832f72 100644 --- a/python/pyabacus/src/py_hsolver.cpp +++ b/python/pyabacus/src/hsolver/py_hsolver.cpp @@ -17,14 +17,12 @@ using namespace pybind11::literals; void bind_hsolver(py::module& m) { - py::module hsolver = m.def_submodule("hsolver"); - - py::class_(hsolver, "diag_comm_info") + py::class_(m, "diag_comm_info") .def(py::init(), "rank"_a, "nproc"_a) .def_readonly("rank", &hsolver::diag_comm_info::rank) .def_readonly("nproc", &hsolver::diag_comm_info::nproc); - py::class_(hsolver, "diago_dav_subspace") + py::class_(m, "diago_dav_subspace") .def(py::init(), R"pbdoc( Constructor of diago_dav_subspace, a class for diagonalizing a linear operator using the Davidson-Subspace Method. @@ -59,9 +57,8 @@ void bind_hsolver(py::module& m) The maximum number of iterations. need_subspace : bool Whether to use the subspace function. - diag_ethr : list[float] - A list of float values indicating the thresholds of each band for the diagonalization, - meaning that the corresponding eigenvalue is to be calculated. + diag_ethr : List[float] | None, optional + The list of thresholds of bands, by default None. scf_type : bool Whether to use the SCF type, which is used to determine the convergence criterion. @@ -92,7 +89,7 @@ void bind_hsolver(py::module& m) Get the eigenvalues. )pbdoc"); - py::class_(hsolver, "diago_david") + py::class_(m, "diago_david") .def(py::init(), R"pbdoc( Constructor of diago_david, a class for diagonalizing a linear operator using the Davidson Method. @@ -148,3 +145,10 @@ void bind_hsolver(py::module& m) Get the eigenvalues. )pbdoc"); } + +PYBIND11_MODULE(_hsolver_pack, m) +{ + m.doc() = "Submodule for pyabacus: hsolver"; + + bind_hsolver(m); +} \ No newline at end of file diff --git a/python/pyabacus/src/py_abacus.cpp b/python/pyabacus/src/py_abacus.cpp deleted file mode 100644 index 5874667761..0000000000 --- a/python/pyabacus/src/py_abacus.cpp +++ /dev/null @@ -1,16 +0,0 @@ -#include -#include - -namespace py = pybind11; - -void bind_base_math(py::module& m); -void bind_m_nao(py::module& m); -void bind_hsolver(py::module& m); - -PYBIND11_MODULE(_core, m) -{ - m.doc() = "Python extension for ABACUS built with pybind11 and scikit-build."; - bind_base_math(m); - bind_m_nao(m); - bind_hsolver(m); -} \ No newline at end of file diff --git a/python/pyabacus/src/pyabacus/ModuleBase/_module_base.py b/python/pyabacus/src/pyabacus/ModuleBase/_module_base.py index fe96247e43..e02b9fea6a 100644 --- a/python/pyabacus/src/pyabacus/ModuleBase/_module_base.py +++ b/python/pyabacus/src/pyabacus/ModuleBase/_module_base.py @@ -5,12 +5,15 @@ """ import numpy as np - from typing import overload from numpy.typing import NDArray -class Sphbes: - def __init__(self) -> None: ... +from ._base_pack import Sphbes as _Sphbes, Integral as _Integral, SphericalBesselTransformer as _SphericalBesselTransformer + +class Sphbes(_Sphbes): + def __init__(self) -> None: + super().__init__() + @overload @staticmethod def sphbesj(l: int, x: float) -> float: ... @@ -23,6 +26,10 @@ def sphbesj( l: int, jl: NDArray[np.float64] ) -> None: ... + + def sphbesj(self, *args, **kwargs): + return super().sphbesj(*args, **kwargs) + @overload @staticmethod def dsphbesj(l: int, x: float) -> float: ... @@ -35,11 +42,18 @@ def dsphbesj( l: int, djl: NDArray[np.float64] ) -> None: ... - @staticmethod - def sphbes_zeros(l: int, n: int, zeros: NDArray[np.float64]) -> None: ... -class Integral: - def __init__(self) -> None: ... + def dsphbesj(self, *args, **kwargs): + return super().dsphbesj(*args, **kwargs) + + @staticmethod + def sphbes_zeros(l: int, n: int, zeros: NDArray[np.float64]) -> None: + super().sphbes_zeros(l, n, zeros) + +class Integral(_Integral): + def __init__(self) -> None: + super().__init__() + @overload @staticmethod def Simpson_Integral( @@ -56,20 +70,28 @@ def Simpson_Integral( dr: float, asum: float ) -> float: ... + + def Simpson_Integral(self, *args, **kwargs): + return super().Simpson_Integral(*args, **kwargs) + @staticmethod def Simpson_Integral_0toall( mesh: int, func: NDArray[np.float64], rab: NDArray[np.float64], asum: NDArray[np.float64] - ) -> None: ... + ) -> None: + super().Simpson_Integral_0toall(mesh, func, rab, asum) + @staticmethod def Simpson_Integral_alltoinf( mesh: int, func: NDArray[np.float64], rab: NDArray[np.float64], asum: NDArray[np.float64] - ) -> None: ... + ) -> None: + super().Simpson_Integral_alltoinf(mesh, func, rab, asum) + @overload @staticmethod def simpson( @@ -84,6 +106,10 @@ def simpson( f: NDArray[np.float64], h: NDArray[np.float64], ) -> float: ... + + def simpson(self, *args, **kwargs): + return super().simpson(*args, **kwargs) + @overload @staticmethod def Gauss_Legendre_grid_and_weight( @@ -100,6 +126,10 @@ def Gauss_Legendre_grid_and_weight( x: NDArray[np.float64], w: NDArray[np.float64], ) -> None: ... + + def Gauss_Legendre_grid_and_weight(self, *args, **kwargs): + return super().Gauss_Legendre_grid_and_weight(*args, **kwargs) -class SphericalBesselTransformer: - def __init__(self) -> None: ... \ No newline at end of file +class SphericalBesselTransformer(_SphericalBesselTransformer): + def __init__(self) -> None: + super().__init__() diff --git a/python/pyabacus/src/pyabacus/ModuleNAO/__init__.py b/python/pyabacus/src/pyabacus/ModuleNAO/__init__.py index 6c43ea250f..b12dddb305 100644 --- a/python/pyabacus/src/pyabacus/ModuleNAO/__init__.py +++ b/python/pyabacus/src/pyabacus/ModuleNAO/__init__.py @@ -1 +1,4 @@ -from __future__ import annotations \ No newline at end of file +from __future__ import annotations +from ._module_nao import * + +__all__ = ["RadialCollection", "TwoCenterIntegrator", "NumericalRadial"] \ No newline at end of file diff --git a/python/pyabacus/src/pyabacus/ModuleNAO/_module_nao.py b/python/pyabacus/src/pyabacus/ModuleNAO/_module_nao.py index c27f2ffe5e..b8ed03a1a4 100644 --- a/python/pyabacus/src/pyabacus/ModuleNAO/_module_nao.py +++ b/python/pyabacus/src/pyabacus/ModuleNAO/_module_nao.py @@ -5,13 +5,13 @@ """ import numpy as np - from numpy.typing import NDArray from pyabacus.ModuleBase import SphericalBesselTransformer -from typing import overload, List +from typing import overload, List, Tuple +from ._nao_pack import RadialCollection as _RadialCollection, TwoCenterIntegrator as _TwoCenterIntegrator, NumericalRadial as _NumericalRadial -class RadialCollection: +class RadialCollection(_RadialCollection): def __init__(self) -> None: """ A class that holds all numerical radial functions of the same kind. @@ -20,7 +20,7 @@ def __init__(self) -> None: of numerical atomic orbitals, or all Kleinman-Bylander beta functions from all elements involved in a calculation. """ - pass + super().__init__() def build( self, @@ -31,7 +31,7 @@ def build( """ Builds the collection from (orbital) files. """ - pass + super().build(nfile, file_list, ftype) def set_transformer( self, @@ -41,7 +41,7 @@ def set_transformer( """ Sets a spherical Bessel transformers for all RadialSet objects. """ - pass + super().set_transformer(sbt, update) def set_uniform_grid( self, @@ -54,7 +54,7 @@ def set_uniform_grid( """ Sets a common uniform grid for all RadialSet objects. """ - pass + super().set_uniform_grid(for_r_space, ngrid, cutoff, mode, enable_fft) def set_grid( self, @@ -66,33 +66,408 @@ def set_grid( """ Sets a common grid for all RadialSet objects """ - pass + super().set_grid(for_r_space, ngrid, grid, mode) + + def __call__( + self, + itype: int, + l: int, + izeta: int + ) -> 'NumericalRadial': + return super().__call__(itype, l, izeta) + + def symbol(self, itype: int) -> str: + return super().symbol(itype) - def __call__(self, itype: int, l: int, izeta: int) -> 'NumericalRadial': ... - def symbol(self, itype: int) -> str: ... @property - def ntype(self) -> int: ... - @overload - def lmax(self, itype: int) -> int: ... + def ntype(self) -> int: + return super().ntype + + def lmax(self, itype: int) -> int: + return super().lmax(itype) + @property - def lmax(self) -> int: ... - @overload - def rcut_max(self, itype: int) -> float: ... + def lmax(self) -> int: + return super().lmax + + def rcut_max(self, itype: int) -> float: + return super().rcut_max(itype) + @property - def rcut_max(self) -> float: ... - def nzeta(self, itype: int, l: int) -> int: ... + def rcut_max(self) -> float: + return super().rcut_max + + def nzeta(self, itype: int, l: int) -> int: + return super().nzeta(itype, l) + @overload def nzeta_max(self, itype: int) -> int: ... @overload def nzeta_max(self) -> int: ... + + def nzeta_max(self, *args, **kwargs): + return super().nzeta_max(*args, **kwargs) + @overload def nchi(self, itype: int) -> int: ... @overload def nchi(self) -> int: ... -class TwoCenterIntegrator: ... + def nchi(self, *args, **kwargs): + return super().nchi(*args, **kwargs) + +class TwoCenterIntegrator(_TwoCenterIntegrator): + def __init__(self) -> None: + """ + A class to compute two-center integrals. + + This class computes two-center integrals of the form: + + / + I(R) = | dr phi1(r) (op) phi2(r - R) + / + + as well as their gradients, where op is 1 (overlap) or minus Laplacian (kinetic), and phi1, + phi2 are "atomic-orbital-like" functions of the form: + + phi(r) = chi(|r|) * Ylm(r/|r|) + + where chi is some numerical radial function and Ylm is some real spherical harmonics. + + This class is designed to efficiently compute the two-center integrals + between two "collections" of the above functions with various R, e.g., the + overlap integrals between all numerical atomic orbitals and all + Kleinman-Bylander nonlocal projectors, the overlap & kinetic integrals between all numerical atomic orbitals, etc. + This is done by tabulating the radial part of the integrals on an r-space grid and the real Gaunt coefficients in advance. + """ + super().__init__() + + def tabulate( + self, + bra: 'RadialCollection', + ket: 'RadialCollection', + op: str, + nr: int, + cutoff: float + ) -> None: + """ + Tabulates the radial part of a two-center integral. + + Parameters: + bra (RadialFunctions): The radial functions of the first collection. + ket (RadialFunctions): The radial functions of the second collection. + op (char): Operator, could be 'S' (overlap) or 'T' (kinetic). + nr (int): Number of r-space grid points. + cutoff (float): r-space cutoff radius. + """ + super().tabulate(bra, ket, op, nr, cutoff) + + def calculate( + self, + itype1: int, + l1: int, + izeta1: int, + m1: int, + itype2: int, + l2: int, + izeta2: int, + m2: int, + pvR: NDArray[np.float64], + cal_grad: bool = False + ) -> Tuple[NDArray[np.float64], NDArray[np.float64]]: + """ + Compute the two-center integrals. + + This function calculates the two-center integral + + / + I(R) = | dr phi1(r) (op_) phi2(r - R) + / + + or its gradient by using the tabulated radial part and real Gaunt coefficients. + + Parameters + ---------- + itype1 : int + Element index of orbital 1. + l1 : int + Angular momentum of orbital 1. + izeta1 : int + Zeta number of orbital 1. + m1 : int + Magnetic quantum number of orbital 1. + itype2 : int + Element index of orbital 2. + l2 : int + Angular momentum of orbital 2. + izeta2 : int + Zeta number of orbital 2. + m2 : int + Magnetic quantum number of orbital 2. + pvR : array_like + R2 - R1, the displacement vector between the two centers. + cal_grad : bool, optional + The gradient will not be computed if cal_grad is false. + + Returns + ------- + out_array : array_like + The two-center integral. + grad_out_array : array_like + Gradient of the two-center integral. + """ + return super().calculate(itype1, l1, izeta1, m1, itype2, l2, izeta2, m2, pvR, cal_grad) + + def snap( + self, + itype1: int, + l1: int, + izeta1: int, + m1: int, + itype2: int, + pvR: NDArray[np.float64], + deriv: bool + ) -> List[List[float]]: + """ + Compute a batch of two-center integrals. + + This function calculates the two-center integrals (and optionally their gradients) + between one orbital and all orbitals of a certain type from the other collection. + """ + return super().snap(itype1, l1, izeta1, m1, itype2, pvR, deriv) + +class NumericalRadial(_NumericalRadial): + def __init__(self) -> None: + """ + A class that represents a numerical radial function. + + This class is designed to be the container for the radial part of numerical atomic orbitals, Kleinman-Bylander beta functions, and all other similar numerical radial functions in three-dimensional space, each of which is associated with some angular momentum l and whose r and k space values are related by an l-th order spherical Bessel transform. + + A NumericalRadial object can be initialized by "build", which requires the angular momentum, the number of grid points, the grid and the corresponding values. Grid does not have to be uniform. One can initialize the object in either r or k space. After initialization, one can set the + grid in the other space via set_grid or set_uniform_grid. Values in the other space are automatically computed by a spherical Bessel transform. + """ + super().__init__() + + def build( + self, + l: int, + for_r_space: bool, + ngrid: int, + grid: NDArray[np.float64], + value: NDArray[np.float64], + p: int = 0, + izeta: int = 0, + symbol: str = "", + itype: int = 0, + init_sbt: bool = True + ) -> None: + """ + Initializes the object by providing the grid & values in one space. + + Parameters + ---------- + l : int + Angular momentum. + for_r_space : bool + Specifies whether the input corresponds to r or k space. + ngrid : int + Number of grid points. + grid : array_like + Grid points, must be positive & strictly increasing. + value : array_like + Values on the grid. + p : float + Implicit exponent in input values, see pr_ & pk_. + izeta : int + Multiplicity index of radial functions of the same itype and l. + symbol : str + Chemical symbol. + itype : int + Index for the element in calculation. + init_sbt : bool + If true, internal SphericalBesselTransformer will be initialized. + + Notes + ----- + init_sbt is only useful when the internal SphericalBesselTransformer (sbt_) is null-initialized; The function will NOT reset sbt_ if it is already usable. + """ + super().build(l, for_r_space, ngrid, grid, value, p, izeta, symbol, itype, init_sbt) + + def set_transformer( + self, + sbt: SphericalBesselTransformer, + update: int = 0 + ) -> None: + """ + Sets a SphericalBesselTransformer. + + By default, the class uses an internal SphericalBesselTransformer, but one can optionally use a shared one. This could be beneficial when there are a lot of NumericalRadial objects whose grids have the same size. + + Parameters + ---------- + sbt : SphericalBesselTransformer + An external transformer. + update : int + Specifies whether and how values are recomputed with the new transformer. + Accepted values are: + * 0: does not recompute values; + * 1: calls a forward transform; + * -1: calls a backward transform. + """ + super().set_transformer(sbt, update) + + def set_grid( + self, + for_r_space: bool, + ngrid: int, + grid: NDArray[np.float64], + mode: str = 'i' + ) -> None: + """ + Sets up a grid. + + This function can be used to set up the grid which is absent in "build" (in which case values on the new grid are automatically computed by a spherical Bessel transform) or interpolate on an existing grid to a new grid. + + Parameters + ---------- + for_r_space : bool + Specifies whether to set grid for the r or k space. + ngrid : int + Number of grid points. + grid : array_like + Grid points, must be positive & strictly increasing. + mode : char + Specifies how values are updated, could be 'i' or 't': + - 'i': New values are obtained by interpolating and zero-padding + the existing values from current space. With this option, + it is an error if the designated space does not have a grid; + - 't': New values are obtained via transform from the other space. + With this option, it is an error if the other space does not + have a grid. + """ + super().set_grid(for_r_space, ngrid, grid, mode) + + def set_uniform_grid( + self, + for_r_space: bool, + ngrid: int, + cutoff: float, + mode: str = 'i', + enable_fft: bool = False + ) -> None: + """ + Sets up a uniform grid. + + The functionality of this function is similar to set_grid, except that the new grid is a uniform grid specified by the cutoff and the number of grid points, which are calculated as: + + grid[i] = i * (cutoff / (ngrid - 1)) + + Parameters + ---------- + for_r_space : bool + Specifies whether to set grid for the r or k space. + ngrid : int + Number of grid points. + cutoff : float + The maximum value of the grid, which determines the grid spacing. + enable_fft : bool + If true, this function will not only set up the grid & values in the designated space, but also sets the grid in the other space such that the r & k grids are FFT-compliant (and updates values via a FFT-based spherical Bessel transform). + mode : char + Specifies how values are updated, could be 'i' or 't'. + """ + super().set_uniform_grid(for_r_space, ngrid, cutoff, mode, enable_fft) + + def set_value( + self, + for_r_space: bool, + value: NDArray[np.float64], + p: int + ) -> None: + """ + Updates values on an existing grid. + + This function does not alter the grid; it merely updates values on the existing grid. The number of values to read from "value" is determined by the current number of points in the r or k space (nr_ or nk_). Values of the other space will be automatically updated if they exist. + + Warning + ------- + This function does not check the index bound; use with care! + """ + super().set_value(for_r_space, value, p) + + def wipe( + self, + r_space: bool = True, + k_space: bool = True + ) -> None: + super().wipe(r_space, k_space) + + def normalize(self, for_r_space: bool = True) -> None: + """ + Normalizes the radial function. + + The radial function is normalized such that the integral of the square of the function multiplied by the square of the radial coordinate over the entire space is equal to one: + + ∫ from 0 to +∞ of (x^2 * f(x)^2) dx = 1 + + where x is r or k. The integral is evaluated with Simpson's rule. Values in the other space are updated automatically via a spherical Bessel transform. + """ + super().normalize(for_r_space) -class NumericalRadial: ... + @property + def symbol(self) -> str: + return super().symbol + @property + def itype(self) -> int: + return super().itype + @property + def izeta(self) -> int: + return super().izeta + @property + def l(self) -> int: + return super().l + @property + def nr(self) -> int: + return super().nr + @property + def nk(self) -> int: + return super().nk + @property + def rcut(self) -> float: + return super().rcut + @property + def kcut(self) -> float: + return super().kcut + @property + def rmax(self) -> float: + return super().rmax + @property + def kmax(self) -> float: + return super().kmax + @property + def pr(self) -> float: + return super().pr + @property + def pk(self) -> float: + return super().pk + @property + def sbt(self) -> SphericalBesselTransformer: + return super().sbt + @property + def rgrid(self) -> NDArray[np.float64]: + return super().rgrid + @property + def kgrid(self) -> NDArray[np.float64]: + return super().kgrid + @property + def rvalue(self) -> NDArray[np.float64]: + return super().rvalue + @property + def kvalue(self) -> NDArray[np.float64]: + return super().kvalue + @property + def is_fft_compliant(self) -> bool: + return super().is_fft_compliant + diff --git a/python/pyabacus/src/pyabacus/__init__.py b/python/pyabacus/src/pyabacus/__init__.py index b81562e214..d8a6e22b72 100644 --- a/python/pyabacus/src/pyabacus/__init__.py +++ b/python/pyabacus/src/pyabacus/__init__.py @@ -5,10 +5,10 @@ def __getattr__(attr): if attr == "ModuleBase": - from ._core import ModuleBase + import pyabacus.ModuleBase as ModuleBase return ModuleBase elif attr == "ModuleNAO": - from ._core import ModuleNAO + import pyabacus.ModuleNAO as ModuleNAO return ModuleNAO elif attr == "hsolver": import pyabacus.hsolver as hsolver diff --git a/python/pyabacus/src/pyabacus/hsolver/__init__.py b/python/pyabacus/src/pyabacus/hsolver/__init__.py index e6a2da7ea6..d9d954d9ee 100644 --- a/python/pyabacus/src/pyabacus/hsolver/__init__.py +++ b/python/pyabacus/src/pyabacus/hsolver/__init__.py @@ -1,4 +1,4 @@ from __future__ import annotations from ._hsolver import * -__all__ = ["dav_subspace"] \ No newline at end of file +__all__ = ["diag_comm_info", "dav_subspace", "davidson"] \ No newline at end of file diff --git a/python/pyabacus/src/pyabacus/hsolver/_hsolver.py b/python/pyabacus/src/pyabacus/hsolver/_hsolver.py index 87ad97664e..9dece524a7 100644 --- a/python/pyabacus/src/pyabacus/hsolver/_hsolver.py +++ b/python/pyabacus/src/pyabacus/hsolver/_hsolver.py @@ -1,20 +1,28 @@ -# pyabacus.hsolver +""" +pyabacus.hsolver +================== +Module for Solver for Hamiltonian in ABACUS +""" import numpy as np from numpy.typing import NDArray from typing import Tuple, List, Union, Callable -from .._core import hsolver +from ._hsolver_pack import diag_comm_info as _diag_comm_info +from ._hsolver_pack import diago_dav_subspace, diago_david -class diag_comm_info: - def __init__(self, rank: int, nproc: int) -> None: ... +class diag_comm_info(_diag_comm_info): + def __init__(self, rank: int, nproc: int): + super().__init__(rank, nproc) @property - def rank(self) -> int: ... + def rank(self) -> int: + return super().rank @property - def nproc(self) -> int: ... - + def nproc(self) -> int: + return super().nproc + def dav_subspace( mvv_op: Callable[[NDArray[np.complex128]], NDArray[np.complex128]], init_v: NDArray[np.complex128], @@ -76,11 +84,11 @@ def dav_subspace( if init_v.ndim != 1 or init_v.dtype != np.complex128: init_v = init_v.flatten().astype(np.complex128, order='C') - _diago_obj_dav_subspace = hsolver.diago_dav_subspace(dim, num_eigs) + _diago_obj_dav_subspace = diago_dav_subspace(dim, num_eigs) _diago_obj_dav_subspace.set_psi(init_v) _diago_obj_dav_subspace.init_eigenvalue() - comm_info = hsolver.diag_comm_info(0, 1) + comm_info = diag_comm_info(0, 1) assert dav_ndim > 1, "dav_ndim must be greater than 1." assert dav_ndim * num_eigs < dim * comm_info.nproc, "dav_ndim * num_eigs must be less than dim * comm_info.nproc." @@ -151,13 +159,13 @@ def davidson( if init_v.ndim != 1 or init_v.dtype != np.complex128: init_v = init_v.flatten().astype(np.complex128, order='C') - _diago_obj_dav_subspace = hsolver.diago_david(dim, num_eigs) - _diago_obj_dav_subspace.set_psi(init_v) - _diago_obj_dav_subspace.init_eigenvalue() + _diago_obj_david = diago_david(dim, num_eigs) + _diago_obj_david.set_psi(init_v) + _diago_obj_david.init_eigenvalue() - comm_info = hsolver.diag_comm_info(0, 1) + comm_info = diag_comm_info(0, 1) - _ = _diago_obj_dav_subspace.diag( + _ = _diago_obj_david.diag( mvv_op, pre_condition, dav_ndim, @@ -167,8 +175,8 @@ def davidson( comm_info ) - e = _diago_obj_dav_subspace.get_eigenvalue() - v = _diago_obj_dav_subspace.get_psi() + e = _diago_obj_david.get_eigenvalue() + v = _diago_obj_david.get_psi() return e, v \ No newline at end of file diff --git a/python/pyabacus/tests/test_base_math.py b/python/pyabacus/tests/test_base_math.py index 1d5190521f..1ff5205981 100644 --- a/python/pyabacus/tests/test_base_math.py +++ b/python/pyabacus/tests/test_base_math.py @@ -11,6 +11,7 @@ def test_sphbes(): # test for sphbesj assert s.sphbesj(1, 0.0) == 0.0 assert s.sphbesj(0, 0.0) == 1.0 + assert s.sphbesj(1, np.array([0.0]), 1, 1, np.zeros(1)) == None def test_sbt(): sbt = base.SphericalBesselTransformer()