Feature : support PAW in PW calculations (in progress) (#2891)

* Feature : add paw_vnl_psi and test

* Feature : added (S+I)|psi>

* add input variable "use_paw"

* added global class paw_cell

* fix some bugs

* fix macro and cmake

* fix cmake

* fix cmake

* add set_paw_k and init_paw_cell in main program

* put get_vkb into main program

* added accumulate_rhoij to main program

* added set_rhoij

* added prepare_libpaw

* add set_rho_core_paw

* added init_rho

* fix test_paw3

* modify preparation of paw

* add get_nhat

* change ecut

* bug fix

* fix bug

* fix some errors

* modify v_hartree

* fix bug

* fix bug

---------

Co-authored-by: wenfei-li <liwenfei@gmail.com>

CMakeLists.txt

@@ -612,6 +612,12 @@ target_link_libraries(${ABACUS_BIN_NAME}
     container
 )
 
+if(ENABLE_PAW)
+  target_link_libraries(${ABACUS_BIN_NAME}
+  paw
+  )
+endif()
+
 if(ENABLE_LCAO)
   target_link_libraries(${ABACUS_BIN_NAME}
       hamilt_lcao

source/Makefile.Objects

@@ -176,6 +176,7 @@ OBJS_ELECSTAT=elecstate.o\
     potential_new.o\
     potential_types.o\
     pot_local.o\
+    pot_local_paw.o\
     H_Hartree_pw.o\
     H_TDDFT_pw.o\
     pot_xc.o\

source/module_base/global_variable.cpp

@@ -60,6 +60,8 @@ int RELAX_NMAX = 20;
 int md_prec_level = 0;
 int SCF_NMAX = 100;
 
+bool use_paw = false;
+
 std::string BASIS_TYPE = "pw"; // xiaohui add 2013-09-01
 std::string KS_SOLVER = "cg"; // xiaohui add 2013-09-01
 double SEARCH_RADIUS = -1.0;

source/module_base/global_variable.h

@@ -58,6 +58,8 @@ extern int OUT_FREQ_ION;
 extern double relax_scale_force;
 extern bool relax_new;
 
+extern bool use_paw;
+
 extern bool fixed_atoms;
 
 extern int RELAX_NMAX; // 8.3

source/module_basis/module_pw/pw_basis_k.cpp

@@ -330,6 +330,55 @@ void PW_Basis_K::get_ig2ixyz_k()
 #endif
 }
 
+std::vector<int> PW_Basis_K::get_ig2ix(const int ik) const
+{
+    std::vector<int> ig_to_ix;
+    ig_to_ix.resize(npwk[ik]);
+
+    for(int ig = 0; ig < npwk[ik]; ig++)
+    {
+        int isz = this->igl2isz_k[ig + ik * npwk_max];
+        int is = isz / this->nz;
+        int ixy = this->is2fftixy[is];
+        int ix = ixy / this->ny;
+        if (ix < (nx / 2) + 1) ix += nx;
+        ig_to_ix[ig] = ix;
+    }
+    return ig_to_ix;
+}
+
+std::vector<int> PW_Basis_K::get_ig2iy(const int ik) const
+{
+    std::vector<int> ig_to_iy;
+    ig_to_iy.resize(npwk[ik]);
+
+    for(int ig = 0; ig < npwk[ik]; ig++)
+    {
+        int isz = this->igl2isz_k[ig + ik * npwk_max];
+        int is = isz / this->nz;
+        int ixy = this->is2fftixy[is];
+        int iy = ixy % this->ny;
+        if (iy < (ny / 2) + 1) iy += ny;
+        ig_to_iy[ig] = iy;
+    }
+    return ig_to_iy;
+}
+
+std::vector<int> PW_Basis_K::get_ig2iz(const int ik) const
+{
+    std::vector<int> ig_to_iz;
+    ig_to_iz.resize(npwk[ik]);
+
+    for(int ig = 0; ig < npwk[ik]; ig++)
+    {
+        int isz = this->igl2isz_k[ig + ik * npwk_max];
+        int iz = isz % this->nz;
+        if (iz < (nz / 2) + 1) iz += nz;
+        ig_to_iz[ig] = iz;
+    }
+    return ig_to_iz;
+}
+
 template <>
 float * PW_Basis_K::get_kvec_c_data() const {
     return this->s_kvec_c;

source/module_basis/module_pw/pw_basis_k.h

@@ -156,6 +156,13 @@ class PW_Basis_K : public PW_Basis
     //get igl2ig_k or igk(ik,ig) in older ABACUS
     int& getigl2ig(const int ik, const int igl) const;
 
+    //get ig_to_ix
+    std::vector<int> get_ig2ix(const int ik) const;
+    //get ig_to_iy
+    std::vector<int> get_ig2iy(const int ik) const;
+    //get ig_to_iz
+    std::vector<int> get_ig2iz(const int ik) const;
+
     template <typename FPTYPE> FPTYPE * get_gk2_data() const;
     template <typename FPTYPE> FPTYPE * get_gcar_data() const;
     template <typename FPTYPE> FPTYPE * get_kvec_c_data() const;

source/module_cell/CMakeLists.txt

@@ -1,6 +1,8 @@
 add_subdirectory(module_symmetry)
 add_subdirectory(module_neighbor)
-add_subdirectory(module_paw)
+if(ENABLE_PAW)
+  add_subdirectory(module_paw)
+endif()
 
 add_library(
     cell

source/module_cell/module_paw/CMakeLists.txt

@@ -5,8 +5,12 @@ add_library(
     paw_sphbes.cpp
     paw_cell.cpp
     paw_cell_libpaw.cpp
+    paw_atom.cpp
 )
 
+target_link_libraries(paw libpaw_interface -lifcore)
+add_dependencies(paw libpaw_interface)
+
 if(ENABLE_COVERAGE)
   add_coverage(paw)
 endif()

source/module_cell/module_paw/paw_atom.cpp

@@ -59,7 +59,7 @@ void Paw_Atom::convert_rhoij()
     {
         if(std::abs(rhoij[i]) > 1e-10)
         {
-            rhoijselect[nrhoijsel] = i;
+            rhoijselect[nrhoijsel] = i+1; //index in fortran
             rhoijp[nrhoijsel] = rhoij[i];
             nrhoijsel ++;
         }

source/module_cell/module_paw/paw_cell.cpp

@@ -1,23 +1,26 @@
 #include "paw_cell.h"
 #include "module_base/tool_title.h"
 #include "module_base/tool_quit.h"
+#ifdef __MPI
+#include "module_base/parallel_reduce.h"
+#endif
+
+namespace GlobalC
+{
+    Paw_Cell paw_cell;
+}
 
 void Paw_Cell::init_paw_cell(
     const double ecutwfc_in, const double cell_factor_in,
     const double omega_in,
     const int nat_in, const int ntyp_in,
     const int * atom_type_in, const double ** atom_coord_in,
-    const std::vector<std::string> & filename_list_in,
-    const int nx_in, const int ny_in, const int nz_in,
-    const std::complex<double> * eigts1_in, const std::complex<double> * eigts2_in, const std::complex<double> * eigts3_in)
+    const std::vector<std::string> & filename_list_in)
 {
     ModuleBase::TITLE("Paw_Element","init_paw_cell");
 
     this -> nat = nat_in;
     this -> ntyp = ntyp_in;
-    this -> nx = nx_in;
-    this -> ny = ny_in;
-    this -> nz = nz_in;
     this -> omega = omega_in;
 
     atom_coord.resize(nat);
@@ -71,6 +74,16 @@ void Paw_Cell::init_paw_cell(
         int nproj = paw_element_list[it].get_mstates();
         paw_atom_list[iat].init_paw_atom(nproj);
     }
+}
+
+void Paw_Cell::set_eigts(const int nx_in, const int ny_in, const int nz_in,
+        const std::complex<double> * eigts1_in,
+        const std::complex<double> * eigts2_in,
+        const std::complex<double> * eigts3_in)
+{
+    this -> nx = nx_in;
+    this -> ny = ny_in;
+    this -> nz = nz_in;
 
     eigts1.resize(nat);
     eigts2.resize(nat);
@@ -422,6 +435,14 @@ void Paw_Cell::get_vkb()
     }
 }
 
+void Paw_Cell::reset_rhoij()
+{
+    for(int iat = 0; iat < nat; iat ++)
+    {
+        paw_atom_list[iat].reset_rhoij();
+    }
+}
+
 void Paw_Cell::accumulate_rhoij(const std::complex<double> * psi, const double weight)
 {
     ModuleBase::TITLE("Paw_Cell","accumulate_rhoij");
@@ -451,9 +472,9 @@ void Paw_Cell::accumulate_rhoij(const std::complex<double> * psi, const double w
             }
         }
 
-        // ca should be summed over MPI ranks since planewave basis is distributed
-        // but not for now (I'll make sure serial version works first)
-        // Parallel_Reduce::reduce_complex_double_pool(ca.data(), nproj);
+#ifdef __MPI
+         Parallel_Reduce::reduce_complex_double_pool(ca.data(), nproj);
+#endif
 
         paw_atom_list[iat].set_ca(ca, weight);
         paw_atom_list[iat].accumulate_rhoij();
@@ -530,9 +551,9 @@ void Paw_Cell::paw_nl_psi(const int mode, const std::complex<double> * psi, std:
             }
         }
 
-        // ca should be summed over MPI ranks since planewave basis is distributed
-        // but not for now (I'll make sure serial version works first)
-        // Parallel_Reduce::reduce_complex_double_pool(ca.data(), nproj);
+#ifdef __MPI
+        Parallel_Reduce::reduce_complex_double_pool(ca.data(), nproj);
+#endif
 
         // sum_ij D_ij ca_j
         std::vector<std::complex<double>> v_ca;

source/module_cell/module_paw/paw_cell.h

@@ -28,9 +28,12 @@ class Paw_Cell
         const double omega_in,
         const int nat_in, const int ntyp_in,
         const int * atom_type_in, const double ** atom_coord_in,
-        const std::vector<std::string> & filename_list_in,
-        const int nx_in, const int ny_in, const int nz_in,
-        const std::complex<double> * eigts1_in, const std::complex<double> * eigts2_in, const std::complex<double> * eigts3_in);
+        const std::vector<std::string> & filename_list_in);
+
+    void set_eigts(const int nx_in, const int ny_in, const int nz_in,
+        const std::complex<double> * eigts1_in,
+        const std::complex<double> * eigts2_in,
+        const std::complex<double> * eigts3_in);
 
     // Given a list of k points, calculate the structure factors
     // exp(-i(k+G)R_I) = exp(-ikR_I) exp(-iG_xR_Ix) exp(-iG_yR_Iy) exp(-iG_zR_Iz)
@@ -45,6 +48,7 @@ class Paw_Cell
     // then accumulates the contribution of this wavefunction to rhoij
     // Note k-point information is not passed here, but prepared in set_paw_k
     void accumulate_rhoij(const std::complex<double> * psi, const double weight);
+    void reset_rhoij();
 
     // returns rhoij for each atom
     std::vector<std::vector<double>> get_rhoij();
@@ -204,6 +208,10 @@ class Paw_Cell
     int get_libpaw_ixc() {return ixc;}
     int get_libpaw_xclevel() {return xclevel;}
     int get_nspin() {return nspden;}
+
+    int get_nrxx() {return nx*ny*nz;}
+    int get_val(const int it) {return paw_element_list[it].get_zval();}
+    int get_zat(const int it) {return paw_element_list[it].get_zat();}
 
     private:
 // Info to be passed to libpaw_interface:
@@ -234,14 +242,18 @@ class Paw_Cell
 
 // Part IV. Calling Fortran subroutines from libpaw_interface
     public:
-#ifdef USE_PAW
     void prepare_paw();
     void get_vloc_ncoret(double* vloc, double* ncoret);
+    void init_rho(double** rho);
     void set_rhoij(int iat, int nrhoijsel, int size_rhoij, int* rhoijselect, double* rhoijp);
-    void get_nhat(double* nhat, double* nhatgr);
+    void get_nhat(double** nhat, double* nhatgr);
     void calculate_dij(double* vks, double* vxc);
     void get_dij(int iat, int size_dij, double* dij);
-#endif
 };
 
+namespace GlobalC
+{
+    extern Paw_Cell paw_cell;
+}
+
 #endif