Ultrasoft pseudopotentials (#3122)

* Feature: PW_Basis_Sup class

* Refactor: class PW_Basis_Sup

* Feature: enable energy calculation of uspp

* Feature: cal force in uspp

* Featrue: stress calculation of uspp in pw

* Refactor: update GPU interface of some functions

* Test: add unittests

* Docs: add docs and examples

Author: YuLiu98

docs/advanced/input_files/input-main.md

@@ -39,7 +39,7 @@
     - [ecutwfc](#ecutwfc)
     - [ecutrho](#ecutrho)
     - [nx, ny, nz](#nx-ny-nz)
-    - [nsx, nsy, nsz](#nsx-nsy-nsz)
+    - [ndx, ndy, ndz](#ndx-ndy-ndz)
     - [pw\_seed](#pw_seed)
     - [pw\_diag\_thr](#pw_diag_thr)
     - [pw\_diag\_nmax](#pw_diag_nmax)
@@ -699,13 +699,19 @@ These variables are used to control the plane wave related parameters.
 ### nx, ny, nz
 
 - **Type**: Integer
-- **Description**: If set to a positive number, then the three variables specify the numbers of FFT grid points in x, y, z directions, respectively. If set to 0, the number will be calculated from ecutrho. Note: you must specify all three dimensions for this setting to be used.
+- **Description**: If set to a positive number, then the three variables specify the numbers of FFT grid points in x, y, z directions, respectively. If set to 0, the number will be calculated from ecutrho. 
+
+    Note: You must specify all three dimensions for this setting to be used.
 - **Default**: 0
 
-### nsx, nsy, nsz
+### ndx, ndy, ndz
 
 - **Type**: Integer
-- **Description**: If set to a positive number, then the three variables specify the numbers of FFT grid (for the smooth part of charge density in ultrasoft pseudopotential) points in x, y, z directions, respectively. If set to 0, the number will be calculated from ecutwfc. Note: you must specify all three dimensions for this setting to be used.
+- **Description**: If set to a positive number, then the three variables specify the numbers of FFT grid (for the dense part of charge density in ultrasoft pseudopotential) points in x, y, z directions, respectively. If set to 0, the number will be calculated from ecutwfc. 
+
+    Note: You must specify all three dimensions for this setting to be used.
+
+    Note: These parameters must be used combined with [nx,ny,nz](#nx-ny-nz). If [nx,ny,nz](#nx-ny-nz) are unset, ndx,ndy,ndz are used as [nx,ny,nz](#nx-ny-nz).
 - **Default**: 0
 
 ### pw_seed

docs/advanced/pp_orb.md

@@ -45,7 +45,11 @@ $$
 
 ## Pseudopotentials
 
-In ABACUS, we only support norm-conserving pseudopotentials. We support four different formats of the pseudopotential files: UPF, UPF2, VWR, and BLPS. For more information, check the `ATOMIC_SPECIES` section in the specification of the [STRU file](./input_files/stru.md).
+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:
 

examples/ultrasoft/README

@@ -0,0 +1,21 @@
+/*******************************************************************************/
+/
+/ These are the examples of ABACUS program.
+/
+/*******************************************************************************/
+
+ABACUS supports density functional theory calculations using ultrasoft pseudopotentials
+
+First one should prepare corresponding ultrasoft pseudopotentials and other input files.
+
+`ecutrho` should be no less than 4 times of `ecutwfc`. (`ecutrho` = 8 to 12 times `ecutwfc`, typically)
+
+```
+INPUT_PARAMETERS
+ecutwfc                 40
+ecutrho                 240
+```
+
+Detailed explaination of these keywords can be found [here](https://github.com/deepmodeling/abacus-develop/blob/develop/docs/input-main.md#plane-wave-related-variables).
+
+

examples/ultrasoft/bcc_Fe/INPUT

@@ -0,0 +1,31 @@
+INPUT_PARAMETERS
+#Parameters (1.General)
+calculation             scf
+nbands                  10
+symmetry                0
+latname                 bcc
+nspin                   2
+pseudo_dir              ../../../tests/PP_ORB
+
+#Parameters (2.Iteration)
+ecutwfc                 40
+ecutrho                 240
+scf_thr                 1e-7
+scf_nmax                100
+
+#Parameters (3.Basis)
+basis_type             pw
+
+#Parameters (4.Smearing)
+smearing_method        gaussian
+smearing_sigma         0.02
+
+#parameters (5.Mixing)
+mixing_type            pulay
+mixing_beta            0.4
+
+pseudo_mesh            1
+pseudo_rcut            10
+
+cal_force              1
+cal_stress             1

examples/ultrasoft/bcc_Fe/KPT

@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+8 8 8 0 0 0

examples/ultrasoft/bcc_Fe/STRU

@@ -0,0 +1,14 @@
+ATOMIC_SPECIES
+Fe 55.845 Fe.pbe-nd-rrkjus.UPF
+
+LATTICE_CONSTANT
+5.6
+
+ATOMIC_POSITIONS
+Direct
+
+Fe
+8
+1 
+0 0 0
+

examples/ultrasoft/bcc_Fe/run.sh

@@ -0,0 +1,18 @@
+#!/bin/bash
+
+ABACUS_PATH=$(awk -F "=" '$1=="ABACUS_PATH"{print $2}' ../../SETENV)
+ABACUS_NPROCS=$(awk -F "=" '$1=="ABACUS_NPROCS"{print $2}' ../../SETENV)
+ABACUS_THREADS=$(awk -F "=" '$1=="ABACUS_THREADS"{print $2}' ../../SETENV)
+
+OMP_NUM_THREADS=${ABACUS_THREADS} mpirun -np ${ABACUS_NPROCS} ${ABACUS_PATH} | tee scf.output
+
+if [[ ! -f scf.output ]] || 
+   [[ ! -f OUT.ABACUS/running_scf.log ]] ||
+   [[ ! ( "$(tail -1 OUT.ABACUS/running_scf.log)" == " Total  Time  :"* ) ]] 
+then
+	echo "job is failed!"
+	exit 1
+else
+	echo "job is successed!"
+	exit 0
+fi

examples/ultrasoft/runall.sh

@@ -0,0 +1,11 @@
+#!/bin/bash
+
+for ii in *
+do
+	if [ -d $ii ];then
+		cd ${ii}
+		echo "RUN: ${ii}"
+		bash run.sh
+		cd ..
+	fi
+done

source/Makefile.Objects

@@ -340,6 +340,7 @@ OBJS_PSI_INITIALIZER=psi_initializer.o\
 OBJS_PW=fft.o\
     pw_basis.o\
     pw_basis_k.o\
+    pw_basis_sup.o\
     pw_distributeg.o\
     pw_distributeg_method1.o\
     pw_distributeg_method2.o\
@@ -507,6 +508,7 @@ OBJS_SRCPW=H_Ewald_pw.o\
     charge_mixing.o\
     fp_energy.o\
     forces.o\
+    forces_us.o\
     force_op.o\
     stress_op.o\
     wf_op.o\
@@ -532,6 +534,7 @@ OBJS_SRCPW=H_Ewald_pw.o\
     stress_func_kin.o\
     stress_func_loc.o\
     stress_func_nl.o\
+    stress_func_us.o\
     stress_pw.o\
     of_stress_pw.o\
     symmetry_rho.o\

source/module_base/blas_connector.h

@@ -54,10 +54,26 @@ extern "C"
 		const double *beta, double *y, const int *incy);
 
     // A := alpha x * y.T + A
-	void dger_(int *m, int *n, double *alpha, double *x, int *incx, double *y, int *incy, double *a, int *lda);
-	void zgerc_(int *m, int *n, std::complex<double> *alpha,std::complex<double> *x, int *incx, std::complex<double> *y, int *incy,std::complex<double> *a, int *lda);
+    void dger_(const int* m,
+               const int* n,
+               const double* alpha,
+               const double* x,
+               const int* incx,
+               const double* y,
+               const int* incy,
+               double* a,
+               const int* lda);
+    void zgerc_(const int* m,
+                const int* n,
+                const std::complex<double>* alpha,
+                const std::complex<double>* x,
+                const int* incx,
+                const std::complex<double>* y,
+                const int* incy,
+                std::complex<double>* a,
+                const int* lda);
 
-	// level 3: matrix-matrix operations, O(n^2) data and O(n^3) work.
+    // level 3: matrix-matrix operations, O(n^2) data and O(n^3) work.
 
 	// Peize Lin add ?gemm 2017-10-27, to compute C = a * A.? * B.? + b * C
 	// A is general