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 <[email protected]>
deepmodeling · Sep 11, 2023 · 4f4b1d8 · 4f4b1d8
1 parent 1158c14
commit 4f4b1d8
Show file tree

Hide file tree

Showing 53 changed files with 217,194 additions and 523 deletions.
diff --git a/CMakeLists.txt b/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

diff --git a/source/Makefile.Objects b/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\

diff --git a/source/module_base/global_variable.cpp b/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;

diff --git a/source/module_base/global_variable.h b/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

diff --git a/source/module_basis/module_pw/pw_basis_k.cpp b/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;

diff --git a/source/module_basis/module_pw/pw_basis_k.h b/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;

diff --git a/source/module_cell/CMakeLists.txt b/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

diff --git a/source/module_cell/module_paw/CMakeLists.txt b/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()

diff --git a/source/module_cell/module_paw/paw_atom.cpp b/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 ++;
         }

diff --git a/source/module_cell/module_paw/paw_cell.cpp b/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;

diff --git a/source/module_cell/module_paw/paw_cell.h b/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