Ultrasoft pseudopotentials

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

source/module_base/global_variable.cpp

@@ -62,6 +62,7 @@ int SCF_NMAX = 100;
 
 bool use_paw = false;
 bool use_uspp = false;
+bool double_grid = false;
 
 std::string BASIS_TYPE = "pw"; // xiaohui add 2013-09-01
 std::string KS_SOLVER = "cg"; // xiaohui add 2013-09-01

source/module_base/global_variable.h

@@ -60,6 +60,7 @@ extern bool relax_new;
 
 extern bool use_paw;
 extern bool use_uspp;
+extern bool double_grid;
 
 extern bool fixed_atoms;
 

source/module_base/module_container/base/third_party/blas.h

@@ -65,8 +65,24 @@ void dsymv_(const char *uplo, const int *n,
             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.
 

source/module_basis/module_pw/CMakeLists.txt

@@ -2,6 +2,7 @@ list(APPEND objects
     fft.cpp
     pw_basis.cpp
     pw_basis_k.cpp
+    pw_basis_sup.cpp
     pw_distributeg.cpp
     pw_distributeg_method1.cpp
     pw_distributeg_method2.cpp

source/module_basis/module_pw/pw_basis.h

@@ -293,4 +293,5 @@ class PW_Basis
 }
 #endif // PWBASIS_H
 
-#include "./pw_basis_big.h" //temporary it will be removed
+#include "./pw_basis_big.h" //temporary it will be removed
+#include "pw_basis_sup.h"

source/module_basis/module_pw/pw_basis_big.h

@@ -15,15 +15,12 @@ namespace ModulePW
 class PW_Basis_Big: public PW_Basis
 {
 public:
-
-    // combine [bx,by,bz] FFT grids into a big one
-	// typical values are bx=2, by=2, bz=2
-	// nbx=nx/bx, nby=ny/by, nbz=nz/bz, 
-    PW_Basis_Big(){
-        bx = 1;
-        by = 1;
-        bz = 1;
-    }
+  // combine [bx,by,bz] FFT grids into a big one
+  // typical values are bx=2, by=2, bz=2
+  // nbx=nx/bx, nby=ny/by, nbz=nz/bz,
+  PW_Basis_Big()
+  {
+  }
     PW_Basis_Big(std::string device_, std::string precision_) : PW_Basis(device_, precision_) {}
 
     ~PW_Basis_Big(){};
@@ -34,7 +31,7 @@ class PW_Basis_Big: public PW_Basis
         bz = bz_in;
         bxyz = bx * by * bz;
     }
-    int bx,by,bz,bxyz;
+    int bx = 1, by = 1, bz = 1, bxyz = 1;
     int nbx, nby, nbz;
     int nbzp;
     int nbxx;