Merge branch 'develop' of https://github.com/deepmodeling/abacus-develop into hotfix

Author: Qianruipku

README.md

@@ -16,3 +16,5 @@ ABACUS (Atomic-orbital Based Ab-initio Computation at UStc) is an open-source pa
 
 # Online Documentation
 For detailed documentation, please refer to [our documentation website](https://abacus.deepmodeling.com/).
+
+See our [Github Pages](https://mcresearch.github.io/abacus-user-guide/) for more tutorials and developer guides.

docs/advanced/input_files/input-main.md

@@ -243,8 +243,8 @@
     - [exx\_opt\_orb\_ecut](#exx_opt_orb_ecut)
     - [exx\_opt\_orb\_tolerence](#exx_opt_orb_tolerence)
     - [exx\_real\_number](#exx_real_number)
-    - [exx\_symmetry\_realspace](#exx_symmetry_realspace)
     - [rpa\_ccp\_rmesh\_times](#rpa_ccp_rmesh_times)
+    - [exx\_symmetry\_realspace](#exx_symmetry_realspace)
     - [out\_ri\_cv](#out_ri_cv)
   - [Molecular dynamics](#molecular-dynamics)
     - [md\_type](#md_type)
@@ -273,6 +273,9 @@
     - [lj\_epsilon](#lj_epsilon)
     - [lj\_sigma](#lj_sigma)
     - [pot\_file](#pot_file)
+    - [dp\_rescaling](#dp_rescaling)
+    - [dp\_fparam](#dp_fparam)
+    - [dp\_aparam](#dp_aparam)
     - [msst\_direction](#msst_direction)
     - [msst\_vel](#msst_vel)
     - [msst\_vis](#msst_vis)
@@ -422,11 +425,12 @@
     - [nocc](#nocc)
     - [nvirt](#nvirt)
     - [lr\_nstates](#lr_nstates)
+    - [lr\_unrestricted](#lr_unrestricted)
     - [abs\_wavelen\_range](#abs_wavelen_range)
     - [out\_wfc\_lr](#out_wfc_lr)
     - [abs\_broadening](#abs_broadening)
     - [ri\_hartree\_benchmark](#ri_hartree_benchmark)
-    - [aims_nbasis](#aims_nbasis)
+    - [aims\_nbasis](#aims_nbasis)
 
 [back to top](#full-list-of-input-keywords)
 ## System variables
@@ -2908,46 +2912,38 @@ These variables are used to control vdW-corrected related parameters.
 - **Type**: String
 - **Description**: Specifies the method used for Van der Waals (VdW) correction. Available options are:
   - `d2`: [Grimme's D2](https://onlinelibrary.wiley.com/doi/abs/10.1002/jcc.20495) dispersion correction method
-  - `d3_0`: [Grimme's DFT-D3(0)](https://aip.scitation.org/doi/10.1063/1.3382344) dispersion correction method
-  - `d3_bj`: [Grimme's DFTD3(BJ)](https://onlinelibrary.wiley.com/doi/abs/10.1002/jcc.21759) dispersion correction method
+  - `d3_0`: [Grimme's DFT-D3(0)](https://aip.scitation.org/doi/10.1063/1.3382344) dispersion correction method (zero-damping)
+  - `d3_bj`: [Grimme's DFTD3(BJ)](https://onlinelibrary.wiley.com/doi/abs/10.1002/jcc.21759) dispersion correction method (BJ-damping)
   - `none`: no vdW correction
 - **Default**: none
+- **Note**: ABACUS supports automatic setting on DFT-D3 parameters for common functionals after version 3.8.3 (and several develop versions earlier). To benefit from this feature, please specify the parameter `dft_functional` explicitly (for more details on this parameter, please see [dft_functional](#dft_functional)), otherwise the autoset procedure will crash with error message like `cannot find DFT-D3 parameter for XC(***)`. If not satisfied with those in-built parameters, any manually setting on `vdw_s6`, `vdw_s8`, `vdw_a1` and `vdw_a2` will overwrite. 
+- **Special**: There are special cases for functional family wB97 (Omega-B97): if want to use the functional wB97X-D3BJ, one needs to specify the `dft_functional` as `HYB_GGA_WB97X_V` and `vdw_method` as `d3_bj`. If want to use the functional wB97X-D3, specify `dft_functional` as `HYB_GGA_WB97X_D3` and `vdw_method` as `d3_0`.
 
 ### vdw_s6
 
 - **Type**: Real
 - **Availability**: `vdw_method` is set to `d2`, `d3_0`, or `d3_bj`
-- **Description**: This scale factor is used to optimize the interaction energy deviations in van der Waals (vdW) corrected calculations. The recommended values of this parameter are dependent on the chosen vdW correction method and the DFT functional being used. For DFT-D2, the recommended values are 0.75 (PBE), 1.2 (BLYP), 1.05 (B-P86), 1.0 (TPSS), and 1.05 (B3LYP). For DFT-D3, recommended values with different DFT functionals can be found on the [here](https://www.chemiebn.uni-bonn.de/pctc/mulliken-center/software/dft-d3/dft-d3). The default value of this parameter in ABACUS is set to be the recommended value for PBE.
+- **Description**: This scale factor is used to optimize the interaction energy deviations in van der Waals (vdW) corrected calculations. The recommended values of this parameter are dependent on the chosen vdW correction method and the DFT functional being used. For DFT-D2, the recommended values are 0.75 (PBE), 1.2 (BLYP), 1.05 (B-P86), 1.0 (TPSS), and 1.05 (B3LYP). If not set, will use values of PBE functional. For DFT-D3, recommended values with different DFT functionals can be found on the [here](https://github.com/dftd3/simple-dftd3/blob/main/assets/parameters.toml). If not set, will search in ABACUS built-in dataset based on the `dft_functional` keywords. User set value will overwrite the searched value.
 - **Default**:
   - 0.75: if `vdw_method` is set to `d2`
-  - 1.0: if `vdw_method` is set to `d3_0` or `d3_bj`
 
 ### vdw_s8
 
 - **Type**: Real
 - **Availability**: `vdw_method` is set to `d3_0` or `d3_bj`
-- **Description**: This scale factor is relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://www.chemiebn.uni-bonn.de/pctc/mulliken-center/software/dft-d3/dft-d3). The default value of this parameter in ABACUS is set to be the recommended value for PBE.
-- **Default**:
-  - 0.722: if `vdw_method` is set to `d3_0`
-  - 0.7875: if `vdw_method` is set to `d3_bj`
+- **Description**: This scale factor is relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://github.com/dftd3/simple-dftd3/blob/main/assets/parameters.toml). If not set, will search in ABACUS built-in dataset based on the `dft_functional` keywords. User set value will overwrite the searched value.
 
 ### vdw_a1
 
 - **Type**: Real
 - **Availability**: `vdw_method` is set to `d3_0` or `d3_bj`
-- **Description**: This damping function parameter is relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://www.chemiebn.uni-bonn.de/pctc/mulliken-center/software/dft-d3/dft-d3). The default value of this parameter in ABACUS is set to be the recommended value for PBE.
-- **Default**:
-  - 1.217: if `vdw_method` is set to `d3_0`
-  - 0.4289: if `vdw_method` is set to `d3_bj`
+- **Description**: This damping function parameter is relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://github.com/dftd3/simple-dftd3/blob/main/assets/parameters.toml). If not set, will search in ABACUS built-in dataset based on the `dft_functional` keywords. User set value will overwrite the searched value.
 
 ### vdw_a2
 
 - **Type**: Real
 - **Availability**: `vdw_method` is set to `d3_0` or `d3_bj`
-- **Description**: This damping function parameter is only relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://www.chemiebn.uni-bonn.de/pctc/mulliken-center/software/dft-d3/dft-d3). The default value of this parameter in ABACUS is set to be the recommended value for PBE.
-- **Default**:
-  - 1.0: if `vdw_method` is set to `d3_0`
-  - 4.4407: if `vdw_method` is set to `d3_bj`
+- **Description**: This damping function parameter is only relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://github.com/dftd3/simple-dftd3/blob/main/assets/parameters.toml). If not set, will search in ABACUS built-in dataset based on the `dft_functional` keywords. User set value will overwrite the searched value.
 
 ### vdw_d
 
@@ -3925,7 +3921,7 @@ Currently supported: `RPA`, `LDA`, `PBE`, `HSE`, `HF`.
 
 - **Type**: String
 - **Description**: The method to solve the Casida equation $AX=\Omega X$ in LR-TDDFT under Tamm-Dancoff approximation (TDA), where $A_{ai,bj}=(\epsilon_a-\epsilon_i)\delta_{ij}\delta_{ab}+(ai|f_{Hxc}|bj)+\alpha_{EX}(ab|ij)$ is the particle-hole excitation matrix and $X$ is the transition amplitude.
-  - `dav`: Construct $AX$ and diagonalize the Hamiltonian matrix iteratively with Davidson algorithm.
+  - `dav`/`dav_subspace`/ `cg`: Construct $AX$ and diagonalize the Hamiltonian matrix iteratively with Davidson/Non-ortho-Davidson/CG algorithm.
   - `lapack`: Construct the full $A$ matrix and directly diagonalize with LAPACK.
   - `spectrum`: Calculate absorption spectrum only without solving Casida equation. The `OUT.${suffix}/` directory should contain the
   files for LR-TDDFT eigenstates and eigenvalues, i.e. `Excitation_Energy.dat` and `Excitation_Amplitude_${processor_rank}.dat`

python/pyabacus/examples/diago_matrix.py

@@ -10,10 +10,13 @@ def load_mat(mat_file):
     return h_mat, nbasis, nband
 
 def calc_eig_pyabacus(mat_file, method):
-    algo = {
+    dav = {
         'dav_subspace': hsolver.dav_subspace,
         'davidson': hsolver.davidson
     }
+    cg = {
+        'cg': hsolver.cg
+    }
 
     h_mat, nbasis, nband = load_mat(mat_file)
 
@@ -22,25 +25,27 @@ def calc_eig_pyabacus(mat_file, method):
     diag_elem = np.where(np.abs(diag_elem) < 1e-8, 1e-8, diag_elem)
     precond = 1.0 / np.abs(diag_elem)
 
-    def mm_op(x):
+    def mvv_op(x):
         return h_mat.dot(x)
 
-    e, _ = algo[method](
-        mm_op,
-        v0,
-        nbasis,
-        nband,
-        precond,
-        dav_ndim=8,
-        tol=1e-8,
-        max_iter=1000
-    )
+    if method in dav:
+        algo = dav[method]
+        # args: mvvop, init_v, dim, num_eigs, precondition, dav_ndim, tol, max_iter
+        args = (mvv_op, v0, nbasis, nband, precond, 8, 1e-12, 5000)
+    elif method in cg:
+        algo = cg[method]
+        # args: mvvop, init_v, dim, num_eigs, precondition, tol, max_iter
+        args = (mvv_op, v0, nbasis, nband, precond, 1e-12, 5000)
+    else:
+        raise ValueError(f"Method {method} not available")
+    
+    e, _ = algo(*args)
 
     print(f'eigenvalues calculated by pyabacus-{method} is: \n', e)
 
     return e
 
-def calc_eig_scipy(mat_file):
+def calc_eigsh(mat_file):
     h_mat, _, nband = load_mat(mat_file)
     e, _ = scipy.sparse.linalg.eigsh(h_mat, k=nband, which='SA', maxiter=1000)
     e = np.sort(e)
@@ -50,13 +55,12 @@ def calc_eig_scipy(mat_file):
 
 if __name__ == '__main__':
     mat_file = './Si2.mat'
-    method = ['dav_subspace', 'davidson']
+    method = ['dav_subspace', 'davidson', 'cg']
 
     for m in method:
         print(f'\n====== Calculating eigenvalues using {m} method... ======')
         e_pyabacus = calc_eig_pyabacus(mat_file, m)
-        e_scipy = calc_eig_scipy(mat_file)
+        e_scipy = calc_eigsh(mat_file)
 
         print('eigenvalues difference: \n', e_pyabacus - e_scipy)
-    
-    
+        

python/pyabacus/src/hsolver/CMakeLists.txt

@@ -2,6 +2,7 @@
 list(APPEND _diago
     ${HSOLVER_PATH}/diago_dav_subspace.cpp
     ${HSOLVER_PATH}/diago_david.cpp
+    ${HSOLVER_PATH}/diago_cg.cpp
     ${HSOLVER_PATH}/diag_const_nums.cpp
     ${HSOLVER_PATH}/diago_iter_assist.cpp
 

python/pyabacus/src/hsolver/py_diago_cg.hpp

@@ -0,0 +1,192 @@
+#ifndef PYTHON_PYABACUS_SRC_PY_DIAGO_CG_HPP
+#define PYTHON_PYABACUS_SRC_PY_DIAGO_CG_HPP
+
+#include <complex>
+#include <functional>
+
+#include <pybind11/pybind11.h>
+#include <pybind11/complex.h>
+#include <pybind11/functional.h>
+#include <pybind11/numpy.h>
+#include <pybind11/stl.h>
+
+#include <ATen/core/tensor.h>
+#include <ATen/core/tensor_map.h>
+#include <ATen/core/tensor_types.h>
+
+#include "module_hsolver/diago_cg.h"
+#include "module_base/module_device/memory_op.h"
+
+namespace py = pybind11;
+
+namespace py_hsolver
+{
+
+class PyDiagoCG
+{
+public:
+    PyDiagoCG(int dim, int num_eigs) : dim{dim}, num_eigs{num_eigs} { }
+    PyDiagoCG(const PyDiagoCG&) = delete;
+    PyDiagoCG& operator=(const PyDiagoCG&) = delete;
+    PyDiagoCG(PyDiagoCG&& other)
+    {
+        psi = other.psi;
+        other.psi = nullptr;
+
+        eig = other.eig;
+        other.eig = nullptr;
+    }
+
+    ~PyDiagoCG() 
+    {
+        if (psi != nullptr) 
+        {
+            delete psi;
+            psi = nullptr;
+        }
+
+        if (eig != nullptr)
+        {
+            delete eig;
+            eig = nullptr;
+        }
+    }
+
+    void init_eig() 
+    {
+        eig = new ct::Tensor(ct::DataType::DT_DOUBLE, {num_eigs});
+        eig->zero();
+    }
+
+    py::array_t<double> get_eig() 
+    {
+        py::array_t<double> eig_out(eig->NumElements());
+        py::buffer_info eig_buf = eig_out.request();
+        double* eig_out_ptr = static_cast<double*>(eig_buf.ptr);
+
+        if (eig == nullptr) {
+            throw std::runtime_error("eig is not initialized");
+        }
+        double* eig_ptr = eig->data<double>();
+
+        std::copy(eig_ptr, eig_ptr + eig->NumElements(), eig_out_ptr);
+        return eig_out;
+    }
+
+    void set_psi(py::array_t<std::complex<double>> psi_in)
+    {
+        py::buffer_info psi_buf = psi_in.request();
+        std::complex<double>* psi_ptr = static_cast<std::complex<double>*>(psi_buf.ptr);
+
+        psi = new ct::TensorMap(
+            psi_ptr,
+            ct::DataType::DT_COMPLEX_DOUBLE,
+            ct::DeviceType::CpuDevice,
+            ct::TensorShape({num_eigs, dim})
+        );
+    }
+
+    py::array_t<std::complex<double>> get_psi()
+    {
+        py::array_t<std::complex<double>> psi_out({num_eigs, dim});
+        py::buffer_info psi_buf = psi_out.request();
+        std::complex<double>* psi_out_ptr = static_cast<std::complex<double>*>(psi_buf.ptr);
+
+        if (psi == nullptr) {
+            throw std::runtime_error("psi is not initialized");
+        }
+        std::complex<double>* psi_ptr = psi->data<std::complex<double>>();
+
+        std::copy(psi_ptr, psi_ptr + psi->NumElements(), psi_out_ptr);
+        return psi_out;
+    }
+
+    void set_prec(py::array_t<double> prec_in)
+    {
+        py::buffer_info prec_buf = prec_in.request();
+        double* prec_ptr = static_cast<double*>(prec_buf.ptr);
+
+        prec = new ct::TensorMap(
+            prec_ptr,
+            ct::DataType::DT_DOUBLE,
+            ct::DeviceType::CpuDevice,
+            ct::TensorShape({dim})
+        );
+    }
+
+    void diag(
+        std::function<py::array_t<std::complex<double>>(py::array_t<std::complex<double>>)> mm_op,
+        int diag_ndim,
+        double tol,
+        bool need_subspace,
+        bool scf_type,
+        int nproc_in_pool = 1
+    ) {
+        const std::string basis_type = "pw";
+        const std::string calculation = scf_type ? "scf" : "nscf";
+
+        auto hpsi_func = [mm_op] (const ct::Tensor& psi_in, ct::Tensor& hpsi_out) {
+            const auto ndim = psi_in.shape().ndim();
+            REQUIRES_OK(ndim <= 2, "dims of psi_in should be less than or equal to 2");
+            const int nvec   = ndim == 1 ? 1 : psi_in.shape().dim_size(0);
+            const int ld_psi = ndim == 1 ? psi_in.NumElements() : psi_in.shape().dim_size(1);
+
+            // Note: numpy's py::array_t is row-major, and
+            //       our tensor-array is row-major
+            py::array_t<std::complex<double>> psi({ld_psi, nvec});
+            py::buffer_info psi_buf = psi.request();
+            std::complex<double>* psi_ptr = static_cast<std::complex<double>*>(psi_buf.ptr);
+            std::copy(psi_in.data<std::complex<double>>(), psi_in.data<std::complex<double>>() + nvec * ld_psi, psi_ptr);
+
+            py::array_t<std::complex<double>> hpsi = mm_op(psi);
+
+            py::buffer_info hpsi_buf = hpsi.request();
+            std::complex<double>* hpsi_ptr = static_cast<std::complex<double>*>(hpsi_buf.ptr);
+            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 spsi_func = [this] (const ct::Tensor& psi_in, ct::Tensor& spsi_out) {
+            const auto ndim = psi_in.shape().ndim();
+            REQUIRES_OK(ndim <= 2, "dims of psi_in should be less than or equal to 2");
+            const int nrow   = ndim == 1 ? psi_in.NumElements() : psi_in.shape().dim_size(1);
+            const int nbands = ndim == 1 ? 1 : psi_in.shape().dim_size(0);
+            syncmem_z2z_h2h_op()(
+                this->ctx,
+                this->ctx,
+                spsi_out.data<std::complex<double>>(), 
+                psi_in.data<std::complex<double>>(), 
+                static_cast<size_t>(nrow * nbands)
+            );
+        };
+
+        cg = std::make_unique<hsolver::DiagoCG<std::complex<double>, base_device::DEVICE_CPU>>(
+            basis_type,
+            calculation,
+            need_subspace,
+            subspace_func,
+            tol,
+            diag_ndim,
+            nproc_in_pool
+        );
+
+        cg->diag(hpsi_func, spsi_func, *psi, *eig, *prec);
+    }
+
+private:
+    base_device::DEVICE_CPU* ctx = {};
+
+    int dim;
+    int num_eigs;
+
+    ct::Tensor* psi = nullptr;
+    ct::Tensor* eig = nullptr;
+    ct::Tensor* prec = nullptr;
+
+    std::unique_ptr<hsolver::DiagoCG<std::complex<double>, base_device::DEVICE_CPU>> cg;
+};
+
+} // namespace py_hsolver
+
+#endif