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

Author: Qianruipku

docs/advanced/elec_properties/Berry_phase.md

@@ -23,7 +23,7 @@ pseudo_dir      ../../../tests/PP_ORB  //the path to locate the pesudopotential
 orbital_dir     ../../../tests/PP_ORB  //the path to locate the numerical orbital files
 ntype         3
 ecutwfc       50 // Ry
-symmetry      0 // turn off symmetry
+symmetry      -1 // turn off symmetry
 calculation   nscf // non-self-consistent calculation
 basis_type    lcao // atomic basis
 init_chg  file // read charge from files
@@ -70,4 +70,6 @@ The results are shown as follows:
  P =    0.8906925  (mod    2.1748536)  (   0.0000000,   0.0000000,   0.8906925) C/m^2
 ```
 
-The electric polarization **P** is multivalued, which modulo a quantum e**R**/V~cell~. Note: the values in parentheses are the components of the **P** along the c axis in the x, y, z Cartesian coordinates when set gdir = 3 in INPUT file.
+The electric polarization **P** is multivalued, which modulo a quantum e**R**/V~cell~. 
+
+Note: The vectors R1, R2, and R3 refer to the three lattice vectors of the unit cell. When gdir=3, the calculated polarization is along the R3 direction. The three values in parentheses represent the re-projection of the polarization along the R3 direction onto the Cartesian coordinate system (i.e., the xyz coordinate system). To obtain the full polarization components in the Cartesian system, you need to calculate the polarization for R1, R2, and R3 separately, and then sum their respective x, y, and z components.

docs/advanced/pp_orb.md

@@ -19,7 +19,20 @@ Inside ABACUS, orbitals in LCAO basis are arranged lexicographically by species-
 
 ## Generating atomic orbital bases
 
-Users may also choose to generate their own atomic obitals. In ABACUS, the atomic orbital bases are generated using a scheme developed in the [paper](https://iopscience.iop.org/article/10.1088/0953-8984/22/44/445501). A detailed description of the procedure for generating orbitals will be provided later.
+Users may also generate ABACUS numerical atomic obitals based on their own flavor. The theoretical background of orbital generation can be found in following works:
+
+- Spillage: [Chen M, Guo G C, He L. Systematically improvable optimized atomic basis sets for ab initio calculations[J]. Journal of Physics: Condensed Matter, 2010, 22(44): 445501.](https://iopscience.iop.org/article/10.1088/0953-8984/22/44/445501)
+- PTG DPSI: [Lin P, Ren X, He L. Strategy for constructing compact numerical atomic orbital basis sets by incorporating the gradients of reference wavefunctions[J]. Physical Review B, 2021, 103(23): 235131.](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.103.235131)
+
+Guidelines for generating atomic orbital bases are as follows:
+
+- [Numerical Atomic Orbitals 1: the nomenclature and usage of numerical atomic orbitals in ABACUS](https://mcresearch.github.io/abacus-user-guide/abacus-nac1.html) (Chinese)
+- [Numerical Atomic Orbitals 2: generate numerical atomic orbitals based on given norm-conserving pseudopotential](https://mcresearch.github.io/abacus-user-guide/abacus-nac1.html) (Chinese)
+- [Numerical Atomic Orbitals 3: generate high-precision numerical atomic orbitals](https://mcresearch.github.io/abacus-user-guide/abacus-nac1.html) (Chinese)
+
+Stable orbital generation programs can be found in guidelines above, there is also another developing version of orbital generation program, in which algorithms are consecutively improved: [Github repository of ABACUS ORBGEN project](https://github.com/kirk0830/ABACUS-ORBGEN), the usage of which can be found in README (in English) file.
+
+*NOTE*: users are encouraged to cite the above works when numerical atomic orbitals and its generation codes are used in their research.
 
 ## BSSE Correction
 
@@ -51,16 +64,46 @@ $$
 $$
 
 ## Pseudopotentials
-
-In ABACUS, we support norm-conserving and ultrasoft pseudopotentials. 
-For norm-conserving pseudopotentials, we support four different formats of the pseudopotential files: UPF, UPF2, VWR, and BLPS. 
-For ultrasoft pseudopotentials, currently we support only one format of the pseudopotential files: UPF2.
-
-For more information, check the `ATOMIC_SPECIES` section in the specification of the [STRU file](./input_files/stru.md).
-
-Here we list some common sources of the pseudopotential files:
-
-1. [Quantum ESPRESSO](http://www.quantum-espresso.org/pseudopotentials/).
-2. [SG15-ONCV](http://quantum-simulation.org/potentials/sg15_oncv/upf/).
-3. [DOJO](http://www.pseudo-dojo.org/).
-4. [BLPS](https://github.com/PrincetonUniversity/BLPSLibrary).
+### Supported formats
+ABACUS supports both norm-conserving and ultrasoft pseudopotentials. For norm-conserving pseudopotentials, UPF, UPF2, VWR, and BLPS formats are supported. For ultrasoft pseudopotentials, UPF and UPF2 formats are supported. 
+
+### Usage
+For more information about pseudopotential usage, check the `ATOMIC_SPECIES` section in the specification of the [STRU file](./input_files/stru.md).
+
+### Download
+Users can find pseudopotentials in the following links:
+
+**Website**
+- [Quantum ESPRESSO](https://www.quantum-espresso.org/pseudopotentials): the official website of Quantum ESPRESSO, where you can find a large number of pseudopotential files.
+- [Stantard Solid State Pseudopotential library](https://www.materialscloud.org/sssp): a library of **high-quality** pseudopotentials for solid-state calculations, with **a large number of tests on efficiency and precison**.
+- [PWmat](http://www.pwmat.com/potential-download): a website that provides a large number of pseudopotential files, various kinds of semi-core constructed pseudopotentials are included. **Several sets (with or without f-electrons/noncolinear core correction) of Lanthanide pseudopotentials are also available**.
+- [THEOS](http://theossrv1.epfl.ch/Main/Pseudopotentials): PSlibrary 0.3.1, a library of pseudopotentials for DFT calculations, including ultrasoft, paw, norm-conserving both full-relativistic and scalar-relativistic pseudopotentials.
+- [ABACUS@USTC](https://abacus.ustc.edu.cn/pseudo/list.htm): **ABACUS official website** where you can find a large number of pseudopotential files and numerical atomic orbital files.
+- [BLPS](https://github.com/PrincetonUniversity/BLPSLibrary): BLPS format pseudopotential library
+
+**Norm-conserving pseudopotentials**
+- [SG15](http://www.quantum-simulation.org/potentials/sg15_oncv/): **vastly used in ABACUS** DFT calculation and numerical atomic orbital generation.
+- [PseudoDOJO](http://www.pseudo-dojo.org/): another widely used pseudopotential database, developed by Abinit group, **including Lanthanide pseudopotentials (f-electrons frozen)**.
+- [The Rappe group](https://www.sas.upenn.edu/rappegroup/research/pseudo-potential-gga.html): a collection of GGA pseudopotentials which are generated with Opium code, several tests proves that are out-performing in alloy systems.
+- [Matteo Giantomassi's Github repo](https://github.com/gmatteo/pseudos_ac_she): a Github repository that contains norm-conserving pseudopotentials for **Actinides and superheavy elements to 120-th element**.
+
+**Ultrasoft pseudopotentials**
+- [Vanderbilt](http://www.physics.rutgers.edu/~dhv/uspp/): a collection of ultrasoft pseudopotentials generated by Vanderbilt group.
+- [GBRV](https://www.physics.rutgers.edu/gbrv/) by Kevin F. Garrity, Joseph W. Bennett, Karin M. Rabe, and David Vanderbilt: presently the most popular ultrasoft pseudpotentials in Quantum ESPRESSO user community.
+
+### Pseudopotential Generation
+For pseudopotential generation, please refer to the following links for more information:
+- [Quantum ESPRESSO](http://www.quantum-espresso.org/pseudopotentials/)
+- [ONCVPSP](http://www.mat-simresearch.com/)
+- [Opium](https://opium.sourceforge.net/)
+
+A Chinese guideline is also available here: [A brief introduction of norm-conserving pseudopotential generation](https://mcresearch.github.io/abacus-user-guide/abacus-upf.html)
+
+# ABACUS Pseudopotential-Numerical atomic orbital Square (APNS) project
+For the purpose of providing high-quality pseudopotentials and numerical atomic orbitals, we have initiated the APNS project. The project is aimed at providing a large number of high-quality pseudopotentials and numerical atomic orbitals, along with diverse test data for the ABACUS user community, reduce the cost of generating and testing pseudopotentials and numerical atomic orbitals by users, and promote the development of ABACUS software. The project is currently in the development stage, and we welcome contributions from the community. For more information, please refer to the following links:
+- [APNS website: test data and results](https://kirk0830.github.io/ABACUS-Pseudopot-Nao-Square/)
+- [APNS workflow (Github repository): high-throughput test of pseudopotentials and numerical atomic orbitals](https://github.com/kirk0830/ABACUS-Pseudopot-Nao-Square)
+
+There are also other excellent projects that provide high-quality pseudopotentials along with test data:
+- [Solid State Pseudopotential library](https://www.materialscloud.org/sssp)
+- [Verification of the precision of DFT implementation via AiiDA common workflows](https://acwf-verification.materialscloud.org/)

docs/conf.py

@@ -18,7 +18,7 @@
 # -- Project information -----------------------------------------------------
 
 project = 'ABACUS'
-copyright = '2022, ABACUS'
+copyright = '2024, ABACUS'
 author = 'ABACUS'
 
 # The full version, including alpha/beta/rc tags

examples/berryphase/lcao_PbTiO3/INPUT-nscf-a

@@ -3,7 +3,7 @@ pseudo_dir              ../../../tests/PP_ORB
 orbital_dir             ../../../tests/PP_ORB
 ntype       		3
 ecutwfc     		50
-symmetry    		0
+symmetry    		-1
 calculation 		nscf
 basis_type  		lcao
 init_chg  		file

examples/berryphase/lcao_PbTiO3/INPUT-nscf-b

@@ -3,7 +3,7 @@ pseudo_dir              ../../../tests/PP_ORB
 orbital_dir             ../../../tests/PP_ORB
 ntype       		3
 ecutwfc     		50
-symmetry    		0
+symmetry    		-1
 calculation 		nscf
 basis_type  		lcao
 init_chg  		file

examples/berryphase/lcao_PbTiO3/INPUT-nscf-c

@@ -3,7 +3,7 @@ pseudo_dir              ../../../tests/PP_ORB
 orbital_dir             ../../../tests/PP_ORB
 ntype       		3
 ecutwfc     		50
-symmetry    		0
+symmetry    		-1
 calculation 		nscf
 basis_type  		lcao
 init_chg  		file

examples/berryphase/pw_PbTiO3/INPUT-nscf-a

@@ -2,7 +2,7 @@ INPUT_PARAMETERS
 pseudo_dir              ../../../tests/PP_ORB
 ntype                   3
 ecutwfc                 50
-symmetry                0
+symmetry                -1
 pw_diag_thr             1e-10
 calculation             nscf
 basis_type              pw

examples/berryphase/pw_PbTiO3/INPUT-nscf-b

@@ -2,7 +2,7 @@ INPUT_PARAMETERS
 pseudo_dir              ../../../tests/PP_ORB
 ntype                   3
 ecutwfc                 50
-symmetry                0
+symmetry                -1
 pw_diag_thr             1e-10
 calculation             nscf
 basis_type              pw

examples/berryphase/pw_PbTiO3/INPUT-nscf-c

@@ -2,7 +2,7 @@ INPUT_PARAMETERS
 pseudo_dir              ../../../tests/PP_ORB
 ntype                   3
 ecutwfc                 50
-symmetry                0
+symmetry                -1
 pw_diag_thr             1e-10
 calculation             nscf
 basis_type              pw

source/module_base/grid/partition.cpp

@@ -0,0 +1,142 @@
+#include "module_base/grid/partition.h"
+#include "module_base/constants.h"
+
+#include <cmath>
+#include <functional>
+#include <numeric>
+#include <algorithm>
+#include <vector>
+#include <cassert>
+
+namespace Grid {
+namespace Partition {
+
+const double stratmann_a = 0.64;
+
+double w_becke(
+    int nR0,
+    const double* drR,
+    const double* dRR,
+    int nR,
+    const int* iR,
+    int c
+) {
+    assert(nR > 0 && nR0 >= nR);
+    std::vector<double> P(nR, 1.0);
+    for (int i = 0; i < nR; ++i) {
+        int I = iR[i];
+        for (int j = i + 1; j < nR; ++j) {
+            int J = iR[j];
+            double mu = (drR[I] - drR[J]) / dRR[I*nR0 + J];
+            double s = s_becke(mu);
+            P[I] *= s;
+            P[J] *= (1.0 - s); // s(-mu) = 1 - s(mu)
+        }
+    }
+    return P[c] / std::accumulate(P.begin(), P.end(), 0.0);
+}
+
+
+double s_becke(double mu) {
+    /* 
+     * Becke's iterated polynomials (3rd order)
+     *
+     * s(mu) = 0.5 * (1 - p(p(p(mu))))
+     *
+     * p(x) = 0.5 * x * (3 - x^2)
+     *
+     */
+    double p = 0.5 * mu * (3.0 - mu*mu);
+    p = 0.5 * p * (3.0 - p*p);
+    p = 0.5 * p * (3.0 - p*p);
+    return 0.5 * (1.0 - p);
+}
+
+
+double w_stratmann(
+    int nR0,
+    const double* drR,
+    const double* dRR,
+    const double* drR_thr,
+    int nR,
+    int* iR,
+    int c
+) {
+    assert(nR > 0 && nR0 >= nR);
+    int I = iR[c], J = 0;
+
+    // If r falls within the exclusive zone of a center, return immediately.
+    for (int j = 0; j < nR; ++j) {
+        J = iR[j];
+        if (drR[J] <= drR_thr[J]) {
+            return static_cast<double>(I == J);
+        }
+    }
+
+    // Even if the grid point does not fall within the exclusive zone of any
+    // center, the normalized weight could still be 0 or 1, and this can be
+    // figured out by examining the unnormalized weight alone.
+
+    // Swap the grid center to the first position in iteration for convenience.
+    // Restore the original order before return.
+    std::swap(iR[0], iR[c]);
+
+    std::vector<double> P(nR);
+    for (int j = 1; j < nR; ++j) {
+        J = iR[j];
+        double mu = (drR[I] - drR[J]) / dRR[I*nR0 + J];
+        P[j] = s_stratmann(mu);
+    }
+    P[0] = std::accumulate(P.begin() + 1, P.end(), 1.0,
+                           std::multiplies<double>());
+
+    if (P[0] == 0.0 || P[0] == 1.0) {
+        std::swap(iR[0], iR[c]); // restore the original order
+        return P[0];
+    }
+
+    // If it passes all the screening above, all unnormalized weights
+    // have to be calculated in order to get the normalized weight.
+
+    std::for_each(P.begin() + 1, P.end(), [](double& s) { s = 1.0 - s; });
+    for (int i = 1; i < nR; ++i) {
+        I = iR[i];
+        for (int j = i + 1; j < nR; ++j) {
+            J = iR[j];
+            double mu = (drR[I] - drR[J]) / dRR[I*nR0 + J];
+            double s = s_stratmann(mu);
+            P[i] *= s;
+            P[j] *= (1.0 - s); // s(-mu) = 1 - s(mu)
+        }
+    }
+
+    std::swap(iR[0], iR[c]); // restore the original order
+    return P[0] / std::accumulate(P.begin(), P.end(), 0.0);
+}
+
+
+double s_stratmann(double mu) {
+    /*
+     * Stratmann's piecewise cell function
+     *
+     * s(mu) = 0.5 * (1 - g(mu/a))
+     *
+     *        /             -1                          x <= -1
+     *        |
+     * g(x) = | (35x - 35x^3 + 21x^5 - 5x^7) / 16       |x| < 1
+     *        |
+     *        \             +1                          x >= +1
+     *
+     */
+    double x = mu / stratmann_a;
+    double x2 = x * x;
+    double h = 0.0625 * x * (35 + x2 * (-35 + x2 * (21 - 5 * x2)));
+
+    bool mid = std::abs(x) < 1;
+    double g = !mid * (1 - 2 * std::signbit(x)) + mid * h;
+    return 0.5 * (1.0 - g);
+}
+
+
+} // end of namespace Partition
+} // end of namespace Grid