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

… pw_soc_force

Dockerfile.cuda

@@ -1,9 +1,9 @@
 FROM nvidia/cuda:12.2.0-devel-ubuntu22.04
 
 RUN apt update && apt install -y --no-install-recommends \
-    libopenblas-openmp-dev liblapack-dev libscalapack-mpi-dev libelpa-dev libfftw3-dev libcereal-dev \
+    libopenblas-openmp-dev liblapack-dev libscalapack-mpi-dev libfftw3-dev libcereal-dev \
     libxc-dev libgtest-dev libgmock-dev libbenchmark-dev python3-numpy \
-    bc cmake git g++ make bc time sudo unzip vim wget
+    bc cmake git g++ make bc time sudo unzip vim wget libopenmpi-dev gfortran libtool-bin
 
 ENV GIT_SSL_NO_VERIFY=true TERM=xterm-256color \
     OMPI_ALLOW_RUN_AS_ROOT=1 OMPI_ALLOW_RUN_AS_ROOT_CONFIRM=1 \
@@ -13,6 +13,17 @@ RUN git clone https://github.com/llohse/libnpy.git && \
     cp libnpy/include/npy.hpp /usr/local/include && \
     rm -r libnpy
 
+RUN cd /tmp && \
+    ELPA_VER=2024.05.001 && \
+    wget -q https://elpa.mpcdf.mpg.de/software/tarball-archive/Releases/$ELPA_VER/elpa-$ELPA_VER.tar.gz && \
+    tar xzf elpa-$ELPA_VER.tar.gz  && rm elpa-$ELPA_VER.tar.gz && \
+    cd elpa-$ELPA_VER && \
+    ./configure CXX=mpic++ CFLAGS="-O3 -march=native" FCFLAGS="-O3" LDFLAGS="-L/usr/local/cuda/lib64 -lstdc++" NVCCFLAGS="-arch sm_75 -arch sm_80" --enable-openmp --enable-nvidia-gpu --with-NVIDIA-GPU-compute-capability="sm_70" --with-cuda-path=/usr/local/cuda/ && \
+    make -j`nproc` && \
+    make PREFIX=/usr/local install && \
+    ln -s /usr/local/include/elpa_openmp-$ELPA_VER/elpa /usr/local/include/ && \
+    cd /tmp && rm -rf elpa-$ELPA_VER
+
 ADD https://api.github.com/repos/deepmodeling/abacus-develop/git/refs/heads/develop /dev/null
 
 RUN git clone https://github.com/deepmodeling/abacus-develop.git --depth 1 && \

docs/advanced/elec_properties/charge.md

@@ -10,7 +10,7 @@ After finishing the calculation, the information of the charge density is stroed
 The SPIN${spin}_CHG.cube file looks like:
 
 ```
-Cubefile created from ABACUS SCF calculation
+STEP: 0  Cubefile created from ABACUS. Inner loop is z, followed by y and x
 2 (nspin) 0.914047 (fermi energy, in Ry)
 2 0.0 0.0 0.0 
 27 0.222222 0 0
@@ -31,7 +31,7 @@ Cubefile created from ABACUS SCF calculation
  ...
 ```
 
-The first line is a brief description.\
+The first line contains the current ion step.\
 The second line contains NSPIN and Fermi energy.\
 The following 4 lines are the informations of lattice, in order:\
  total number of atoms, the coordinate of original point.\

docs/advanced/input_files/input-main.md

@@ -87,6 +87,7 @@
     - [printe](#printe)
     - [scf\_nmax](#scf_nmax)
     - [scf\_thr](#scf_thr)
+    - [scf\_ene\_thr](#scf_ene_thr)
     - [scf\_thr\_type](#scf_thr_type)
     - [chg\_extrap](#chg_extrap)
     - [lspinorb](#lspinorb)
@@ -142,6 +143,7 @@
     - [out\_level](#out_level)
     - [out\_alllog](#out_alllog)
     - [out\_mat\_hs](#out_mat_hs)
+    - [out\_mat\_tk](#out_mat_tk)
     - [out\_mat\_r](#out_mat_r)
     - [out\_mat\_hs2](#out_mat_hs2)
     - [out\_mat\_t](#out_mat_t)
@@ -258,7 +260,6 @@
     - [md\_tfreq](#md_tfreq)
     - [md\_tchain](#md_tchain)
     - [md\_pmode](#md_pmode)
-    - [md\_prec\_level](#md_prec_level)
     - [ref\_cell\_factor](#ref_cell_factor)
     - [md\_pcouple](#md_pcouple)
     - [md\_pfirst, md\_plast](#md_pfirst-md_plast)
@@ -422,8 +423,9 @@
     - [lr\_nstates](#lr_nstates)
     - [abs\_wavelen\_range](#abs_wavelen_range)
     - [out\_wfc\_lr](#out_wfc_lr)
-[back to top](#full-list-of-input-keywords)
+    - [abs\_broadening](#abs_broadening)
 
+[back to top](#full-list-of-input-keywords)
 ## System variables
 
 These variables are used to control general system parameters.
@@ -1171,9 +1173,16 @@ Note: In new angle mixing, you should set `mixing_beta_mag >> mixing_beta`. The
 ### scf_thr
 
 - **Type**: Real
-- **Description**: It's the threshold for electronic iteration. It represents the charge density error between two sequential densities from electronic iterations. Usually for local orbitals, usually 1e-6 may be accurate enough.
+- **Description**: It's the density threshold for electronic iteration. It represents the charge density error between two sequential densities from electronic iterations. Usually for local orbitals, usually 1e-6 may be accurate enough.
 - **Default**: 1.0e-9 (plane-wave basis), or 1.0e-7 (localized atomic orbital basis).
 
+### scf_ene_thr
+
+- **Type**: Real
+- **Description**: It's the energy threshold for electronic iteration. It represents the total energy error between two sequential densities from electronic iterations.
+- **Default**: -1.0. If the user does not set this parameter, it will not take effect.
+- **Unit**: eV
+
 ### scf_thr_type
 
 - **Type**: Integer
@@ -1504,8 +1513,10 @@ These variables are used to control the output of properties.
 
 ### out_chg
 
-- **Type**: Integer
+- **Type**: Integer \[Integer\](optional)
 - **Description**: 
+
+  The first integer controls whether to output the charge density on real space grids:
   - 1. Output the charge density (in Bohr^-3) on real space grids into the density files in the folder `OUT.${suffix}`. The files are named as:
     - nspin = 1: SPIN1_CHG.cube;
     - nspin = 2: SPIN1_CHG.cube, and SPIN2_CHG.cube;
@@ -1514,11 +1525,16 @@ These variables are used to control the output of properties.
     - nspin = 1: SPIN1_CHG_INI.cube
     - nspin = 2: SPIN1_CHG_INI.cube, and SPIN2_CHG_INI.cube;
     - nspin = 4: SPIN1_CHG_INI.cube, SPIN2_CHG_INI.cube, SPIN3_CHG_INI.cube, and SPIN4_CHG_INI.cube.
+
+  The second integer controls the precision of the charge density output, if not given, will use `3` as default. For purpose restarting from this file and other high-precision involved calculation, recommend to use `10`.
 
+  ---
   The circle order of the charge density on real space grids is: x is the outer loop, then y and finally z (z is moving fastest).
 
   If EXX(exact exchange) is calculated, (i.e. *[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0/opt_orb* or *[rpa](#rpa)==True*), the Hexx(R) files will be output in the folder `OUT.${suffix}` too, which can be read in NSCF calculation.
-- **Default**: 0
+
+  In molecular dynamics calculations, the output frequency is controlled by [out_interval](#out_interval).
+- **Default**: 0 3
 
 ### out_pot
 
@@ -1535,6 +1551,8 @@ These variables are used to control the output of properties.
     - nspin = 1: SPIN1_POT_INI.cube;
     - nspin = 2: SPIN1_POT_INI.cube, and SPIN2_POT_INI.cube;
     - nspin = 4: SPIN1_POT_INI.cube, SPIN2_POT_INI.cube, SPIN3_POT_INI.cube, and SPIN4_POT_INI.cube.
+
+  In molecular dynamics calculations, the output frequency is controlled by [out_interval](#out_interval).
 - **Default**: 0
 
 ### out_dm
@@ -1601,7 +1619,7 @@ These variables are used to control the output of properties.
 
 ### out_band
 
-- **Type**: Boolean Integer(optional)
+- **Type**: Boolean \[Integer\](optional)
 - **Description**: Whether to output the band structure (in eV), optionally output precision can be set by a second parameter, default is 8. For more information, refer to the [band.md](../elec_properties/band.md)
 - **Default**: False
 
@@ -1645,11 +1663,18 @@ These variables are used to control the output of properties.
 
 ### out_mat_hs
 
-- **Type**: Boolean Integer(optional)
+- **Type**: Boolean \[Integer\](optional)
 - **Availability**: Numerical atomic orbital basis
 - **Description**: Whether to print the upper triangular part of the Hamiltonian matrices (in Ry) and overlap matrices for each k point into files in the directory `OUT.${suffix}`. The second number controls precision. For more information, please refer to [hs_matrix.md](../elec_properties/hs_matrix.md#out_mat_hs). Also controled by [out_interval](#out_interval) and [out_app_flag](#out_app_flag).
 - **Default**: False 8
 
+### out_mat_tk
+
+- **Type**: Boolean \[Integer\](optional)
+- **Availability**: Numerical atomic orbital basis
+- **Description**: Whether to print the upper triangular part of the kinetic matrices (in Ry) for each k point into `OUT.${suffix}/data-i-T`, where i is the index of k points (see `OUT.${suffix}/kpoints`). One may optionally provide a second parameter to specify the precision. 
+- **Default**: False \[8\]
+
 ### out_mat_r
 
 - **Type**: Boolean
@@ -1724,8 +1749,7 @@ The band (KS orbital) energy for each (k-point, spin, band) will be printed in t
 ### out_interval
 
 - **Type**: Integer
-- **Availability**: Numerical atomic orbital basis
-- **Description**: Control the interval for printing Mulliken population analysis, $r(R)$, $H(R)$, $S(R)$, $T(R)$, $dH(R)$, $H(k)$, $S(k)$ and $wfc(k)$ matrices during molecular dynamics calculations. Check input parameters [out_mul](#out_mul), [out_mat_r](#out_mat_r), [out_mat_hs2](#out_mat_hs2), [out_mat_t](#out_mat_t), [out_mat_dh](#out_mat_dh), [out_mat_hs](#out_mat_hs) and [out_wfc_lcao](#out_wfc_lcao) for more information, respectively.
+- **Description**: Control the interval for printing charge density, local potential, electrostatic potential, Mulliken population analysis, $r(R)$, $H(R)$, $S(R)$, $T(R)$, $dH(R)$, $H(k)$, $S(k)$ and $\psi(k)$ matrices during molecular dynamics calculations. Check input parameters [out_chg](#out_chg), [out_pot](#out_pot), [out_mul](#out_mul), [out_mat_r](#out_mat_r), [out_mat_hs2](#out_mat_hs2), [out_mat_t](#out_mat_t), [out_mat_dh](#out_mat_dh), [out_mat_hs](#out_mat_hs) and [out_wfc_lcao](#out_wfc_lcao) for more information, respectively.
 - **Default**: 1
 
 ### out_element_info
@@ -1869,7 +1893,7 @@ Warning: this function is not robust enough for the current version. Please try
 
 - **Type**: Boolean
 - **Availability**: numerical atomic orbital basis
-- **Description**: print energy and force labels and descriptors for DeePKS training
+- **Description**: Print labels and descriptors for DeePKS training in OUT.${suffix}. The names of these files start with "deepks".
 - **Note**: In `LCAO` calculation, the path of a numerical descriptor (an `orb` file) is needed to be specified under the `NUMERICAL_DESCRIPTOR` tag in the `STRU` file. For example:
 
   ```text
@@ -2533,14 +2557,14 @@ These variables are used to control molecular dynamics calculations. For more in
 - **Default**: iso
 - **Relavent**: [md_tfreq](#md_tfreq), [md_tchain](#md_tchain), [md_pcouple](#md_pcouple), [md_pfreq](#md_pfreq), and [md_pchain](#md_pchain).
 
-### md_prec_level
+<!-- ### md_prec_level
 
 - **Type**: Integer
 - **Description**: Determine the precision level of variable-cell molecular dynamics calculations.
   - 0: FFT grids do not change, only G vectors and K vectors are changed due to the change of lattice vector. This level is suitable for cases where the variation of the volume and shape is not large, and the efficiency is relatively higher.
   - 2: FFT grids change per step. This level is suitable for cases where the variation of the volume and shape is large, such as the MSST method. However, accuracy comes at the cost of efficiency.
 
-- **Default**: 0
+- **Default**: 0 -->
 
 ### ref_cell_factor
 

docs/advanced/pp_orb.md

@@ -8,7 +8,14 @@ The default value of basis_type is pw. The size of pw basis set is controlled by
 
 When choosing lcao basis set, users need to prepare a set of atomic orbitals. Such files may be downloaded from the [official website](http://abacus.ustc.edu.cn/pseudo/list.htm). For more information, also check the `NUMERICAL_ORBITAL` section in the specification of the [STRU file](./input_files/stru.md).
 
-The sequence of orbitals in lcao basis set is as follows. First, all the orbitals belonging to one particular atom are put together. These atom orbitals are arranged as the atom order specified in the STRU file. Then, the orbitals of each atom are arranged according to the orbital files. If the orbital file says that the number of s、p、d...orbitals is $n_s$、$n_p$、$n_d$...then the orbitals are aligned as first $n_s$ s orbitals, then $n_p$ p orbitals, and then $n_d$ d orbitals...Last, the angular part of each orbital is real spherical harmonic function. They are aligned as Y00, Y10, Y11, Y1-1, Y20, Y21, Y2-1, Y22, Y2-2, which is s,$p_z$,$p_x$,$p_y$,$d_{z^2}$,$d_{xz}$,$d_{yz}$,$d_{x^2-y^2}$,$d_{xy}$. The corresponding formula can be seen in [Table of spherical harmonics - Wikipedia](https://en.wikipedia.org/wiki/Table_of_spherical_harmonics#Real_spherical_harmonics). Note that these formula lack of the Condon–Shortley phase $(-1)^m$, which is presented in the lcao orbitals of ABACUS.
+The angular part of orbitals are real spherical harmonics defined (in terms of conventional spherical harmonics in quantum mechanical literature) as
+
+$$Y_{lm} = \left\{\begin{matrix}\sqrt{2}~\text{Im} Y_l^{|m|} & m \lt 0 \\[6pt] Y_l^0 & m = 0 \\[6pt] \sqrt{2}~\text{Re}Y_{l}^{|m|} & m \gt 0\end{matrix}\right. $$
+
+Note that real spherical harmonics adopted by ABACUS differ from some other definition, e.g. [Table of spherical harmonics - Wikipedia](https://en.wikipedia.org/wiki/Table_of_spherical_harmonics#Real_spherical_harmonics), by a factor of $(-1)^m$.
+
+Inside ABACUS, orbitals in LCAO basis are arranged lexicographically by species-position-l-zeta-m **except for the intra-m ordering**. Specifically, orbitals are first ordered by their atomic species in accordance with the `ATOMIC_SPECIES` section of the STRU file. For orbitals of the same species, orbitals belonging to each atom are put together, with their overall order following the `ATOMIC_POSITIONS` section of the STRU file. Orbitals on each atom are further ascendingly ordered by their angular momentum (s,p,d,f,...), followed by an order based on their their zeta number. Finally, m is ordered as 0, 1, -1, 2, 2, $\ldots$, l, -l, which is the only exception of the lexicographic order.
+
 
 ## Generating atomic orbital bases
 

python/pyabacus/CMakeLists.txt

@@ -11,6 +11,9 @@ find_package(pybind11 CONFIG REQUIRED)
 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")
 
@@ -100,6 +103,29 @@ list(APPEND _naos
 add_library(naopack SHARED 
     ${_naos}
     )
+# add diago shared library
+list(APPEND _diago
+    ${HSOLVER_PATH}/diago_dav_subspace.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
@@ -125,9 +151,10 @@ 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_diago_dav_subspace.cpp
     )
 pybind11_add_module(_core MODULE ${_sources})
-target_link_libraries(_core PRIVATE pybind11::headers naopack)
+target_link_libraries(_core PRIVATE pybind11::headers naopack diagopack)
 target_compile_definitions(_core PRIVATE VERSION_INFO=${PROJECT_VERSION})
 # set RPATH
 execute_process(
@@ -141,5 +168,7 @@ 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})
 

python/pyabacus/README.md

@@ -1,33 +1,33 @@
 Build Example: TwoCenterIntegral Section in ABACUS
-==============
+==================================================
 
 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).
 
-
 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
 -----------
 
 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,
-using [`cibuildwheel`][].
+using .
 
 Use `pytest -v` to run all the unit tests for pyabacus in the local machine.
+
 ```shell
 $ cd tests/
 $ pytest -v
 ```
 
 Run `python vis_nao.py` to visualize the numerical orbital.
+
 ```shell
 $ cd examples/
 $ python vis_nao.py
@@ -41,6 +41,16 @@ $ python ex_s_rotate.py
 norm(S_e3 - S_numer) =  3.341208104032616e-15
 ```
 
+Run `python diago_matrix.py` in `examples` to check the diagonalization of a matrix.
+
+```shell
+$ cd examples/
+$ python diago_matrix.py
+eigenvalues calculated by pyabacus: [-0.38440611 0.24221155 0.31593272 0.53144616 0.85155108 1.06950155 1.11142051 1.12462152]
+eigenvalues calculated by scipy: [-0.38440611 0.24221155 0.31593272 0.53144616 0.85155108 1.06950154 1.11142051 1.12462151]
+error: [9.26164700e-12 2.42959514e-10 2.96529468e-11 7.77933273e-12 7.53686002e-12 2.95628810e-09 1.04678111e-09 7.79106313e-09]
+```
+
 License
 -------
 
@@ -58,4 +68,4 @@ s.sphbesj(1, 0.0)
 0.0
 ```
 
-[`cibuildwheel`]:          https://cibuildwheel.readthedocs.io
+[`cibuildwheel`]: https://cibuildwheel.readthedocs.io

python/pyabacus/examples/diago_matrix.py

@@ -0,0 +1,38 @@
+from pyabacus import hsolver
+import numpy as np
+import scipy
+
+h_mat = scipy.io.loadmat('./Si2.mat')['Problem']['A'][0, 0]
+
+nbasis = h_mat.shape[0]
+nband = 8
+
+v0 = np.random.rand(nbasis, nband)
+
+diag_elem = h_mat.diagonal()
+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):
+    return h_mat.dot(x)
+
+e, v = hsolver.dav_subspace(
+    mm_op,
+    v0,
+    nbasis,
+    nband,
+    precond,
+    dav_ndim=8,
+    tol=1e-8,
+    max_iter=1000,
+    scf_type=False
+)
+
+print('eigenvalues calculated by pyabacus: ', e)
+
+e_scipy, v_scipy = scipy.sparse.linalg.eigsh(h_mat, k=nband, which='SA', maxiter=1000)
+e_scipy = np.sort(e_scipy)
+print('eigenvalues calculated by scipy: ', e_scipy)
+
+print('error:', e - e_scipy)

python/pyabacus/pyproject.toml

@@ -10,6 +10,7 @@ description="A minimal two_center_integral package (with pybind11)"
 readme = "README.md"
 authors = [
   { name = "Jie Li", email = "[email protected]" },
+  { name = "Chenxu Bai", email = "[email protected]" },
 ]
 requires-python = ">=3.7"
 classifiers = [