From 748f32c6013e2e382fbf83673cf54d3be3779087 Mon Sep 17 00:00:00 2001 From: wenfei-li Date: Sun, 8 Oct 2023 18:05:25 +0800 Subject: [PATCH 1/4] Feature : parallelization of PAW --- source/module_cell/module_paw/paw_cell.h | 18 ++- .../module_paw/paw_cell_libpaw.cpp | 146 +++++++++++++++++- .../module_cell/module_paw/test/test_paw3.cpp | 10 +- .../module_cell/module_paw/test/test_paw4.cpp | 5 +- source/module_esolver/esolver_ks.cpp | 3 +- source/module_esolver/esolver_ks_pw.cpp | 3 +- source/module_hsolver/hsolver_pw.cpp | 2 +- 7 files changed, 175 insertions(+), 12 deletions(-) diff --git a/source/module_cell/module_paw/paw_cell.h b/source/module_cell/module_paw/paw_cell.h index 8bf91d32017..8da5737e37b 100644 --- a/source/module_cell/module_paw/paw_cell.h +++ b/source/module_cell/module_paw/paw_cell.h @@ -13,6 +13,7 @@ #include "paw_element.h" #include "paw_atom.h" #include "module_base/matrix3.h" +#include "module_base/global_variable.h" class Paw_Cell { @@ -183,7 +184,8 @@ class Paw_Cell void set_libpaw_cell(const ModuleBase::Matrix3 latvec, const double lat0); // FFT grid information, sets ngfft and ngfftdg void set_libpaw_fft(const int nx_in, const int ny_in, const int nz_in, - const int nxdg_in, const int nydg_in, const int nzdg_in); + const int nxdg_in, const int nydg_in, const int nzdg_in, + const int * start_z_in, const int * num_z_in); // Sets natom, ntypat, typat and xred void set_libpaw_atom(const int natom_in, const int ntypat_in, const int * typat_in, const double * xred_in); // Sets filename_list @@ -211,7 +213,8 @@ class Paw_Cell int get_libpaw_xclevel() {return xclevel;} int get_nspin() {return nspden;} - int get_nrxx() {return nx*ny*nz;} + int get_nfft() {return nfft;} + int get_nrxx() {return nx*ny*num_z[GlobalV::RANK_IN_POOL];} 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();} @@ -252,6 +255,17 @@ class Paw_Cell void calculate_dij(double* vks, double* vxc); void get_dij(int iat, int size_dij, double* dij); void get_sij(int iat, int size_sij, double* sij); + +// Part V. Relevant for parallel computing +// Note about the parallelization of PAW: ABINIT supports the parallelization based on +// real-space grid points, but it has a rather complicated way of determining the parallelization scheme +// Furthermore, I'm calling libpaw to calculate vloc, ncoret, nhat and dij, which presumably will not be the +// most time consuming part. So I opt to not use parallelization in PAW for now. +// The parallelization part here keeps a record of how ABACUS distributes the xy planes, which will be used +// when distributing vloc, ncoret and nhat + + private: + std::vector start_z, num_z; }; namespace GlobalC diff --git a/source/module_cell/module_paw/paw_cell_libpaw.cpp b/source/module_cell/module_paw/paw_cell_libpaw.cpp index 01872fec30c..4658a98cc91 100644 --- a/source/module_cell/module_paw/paw_cell_libpaw.cpp +++ b/source/module_cell/module_paw/paw_cell_libpaw.cpp @@ -139,7 +139,8 @@ void Paw_Cell::set_libpaw_cell(const ModuleBase::Matrix3 latvec, const double la // FFT grid information, sets ngfft and ngfftdg void Paw_Cell::set_libpaw_fft(const int nx_in, const int ny_in, const int nz_in, - const int nxdg_in, const int nydg_in, const int nzdg_in) + const int nxdg_in, const int nydg_in, const int nzdg_in, + const int * start_z_in, const int * num_z_in) { ModuleBase::TITLE("Paw_Cell", "set_libpaw_fft"); ngfft.resize(3); @@ -148,10 +149,25 @@ void Paw_Cell::set_libpaw_fft(const int nx_in, const int ny_in, const int nz_in, ngfft[0] = nx_in; ngfft[1] = ny_in; ngfft[2] = nz_in; + + nx = nx_in; + ny = ny_in; + nz = nz_in; + ngfftdg[0] = nxdg_in; ngfftdg[1] = nydg_in; ngfftdg[2] = nzdg_in; nfft = ngfftdg[0]*ngfftdg[1]*ngfftdg[2]; + +#ifdef __MPI + start_z.resize(GlobalV::NPROC); + num_z.resize(GlobalV::NPROC); + for(int iproc = 0; iproc < GlobalV::NPROC; iproc ++) + { + start_z[iproc] = start_z_in[iproc]; + num_z[iproc] = num_z_in[iproc]; + } +#endif } // Sets natom, ntypat, typat and xred @@ -181,6 +197,7 @@ void Paw_Cell::set_libpaw_files() { ModuleBase::TITLE("Paw_Cell", "set_libpaw_files"); + filename_list = new char[ntypat*264]; if(GlobalV::MY_RANK == 0) { std::ifstream ifa(GlobalV::stru_file.c_str(), std::ios::in); @@ -196,7 +213,6 @@ void Paw_Cell::set_libpaw_files() if (line.find("PAW_FILES") != std::string::npos) break; } - filename_list = new char[ntypat*264]; for(int i = 0; i < ntypat*264; i++) { filename_list[i] = ' '; @@ -278,14 +294,25 @@ void Paw_Cell::get_vloc_ncoret(double* vloc, double* ncoret) { for(int iy = 0; iy < ny; iy ++) { +#ifdef __MPI + for(int iz = 0; iz < num_z[GlobalV::RANK_IN_POOL]; iz ++) + { + int ind_c = ix*ny*num_z[GlobalV::RANK_IN_POOL] + iy*num_z[GlobalV::RANK_IN_POOL] + iz; + int ind_fortran = (iz+start_z[GlobalV::RANK_IN_POOL])*ny*nx + iy*nx + ix; + + vloc[ind_c] = vloc_tmp[ind_fortran]; + ncoret[ind_c] = ncoret_tmp[ind_fortran]; + } +#else for(int iz = 0; iz < nz; iz ++) { int ind_c = ix*ny*nz + iy*nz + iz; int ind_fortran = iz*ny*nx + iy*nx + ix; - vloc[ind_c] = vloc_tmp[ind_fortran*nspden]; - ncoret[ind_c] = ncoret_tmp[ind_fortran*nspden]; + vloc[ind_c] = vloc_tmp[ind_fortran]; + ncoret[ind_c] = ncoret_tmp[ind_fortran]; } +#endif } } } @@ -316,6 +343,15 @@ void Paw_Cell::get_nhat(double** nhat, double* nhatgr) { for(int iy = 0; iy < ny; iy ++) { +#ifdef __MPI + for(int iz = 0; iz < num_z[GlobalV::RANK_IN_POOL]; iz ++) + { + int ind_c = ix*ny*num_z[GlobalV::RANK_IN_POOL] + iy*num_z[GlobalV::RANK_IN_POOL] + iz; + int ind_fortran = (iz+start_z[GlobalV::RANK_IN_POOL])*ny*nx + iy*nx + ix; + + nhat[is][ind_c] = nhat_tmp[ind_fortran*nspden+is]; + } +#else for(int iz = 0; iz < nz; iz ++) { int ind_c = ix*ny*nz + iy*nz + iz; @@ -323,6 +359,7 @@ void Paw_Cell::get_nhat(double** nhat, double* nhatgr) nhat[is][ind_c] = nhat_tmp[ind_fortran*nspden+is]; } +#endif } } } @@ -333,6 +370,79 @@ void Paw_Cell::calculate_dij(double* vks, double* vxc) { ModuleBase::TITLE("Paw_Cell", "calculate_dij"); double * vks_hartree, * vxc_hartree; + +#ifdef __MPI + double * vks_collected, * vxc_collected; + if(GlobalV::RANK_IN_POOL == 0) + { + vks_hartree = new double[nspden * nfft]; + vxc_hartree = new double[nspden * nfft]; + vks_collected = new double[nspden * nfft]; + vxc_collected = new double[nspden * nfft]; + } + + // Collecting vks and vxc from all processes; I hope there could be a better way + // but this is what I can think of right now. + const int nxy = nx * ny; + for(int is = 0; is < nspden; is ++) + { + double * vks_send = new double[num_z[GlobalV::RANK_IN_POOL]]; + double * vxc_send = new double[num_z[GlobalV::RANK_IN_POOL]]; + double * vks_receive = new double[nz]; + double * vxc_receive = new double[nz]; + for(int ixy = 0; ixy < nxy; ixy ++) + { + for(int iz = 0; iz < num_z[GlobalV::RANK_IN_POOL]; iz++) + { + vks_send[iz] = vks[(ixy*num_z[GlobalV::RANK_IN_POOL] + iz)*nspden + is]; + vxc_send[iz] = vxc[(ixy*num_z[GlobalV::RANK_IN_POOL] + iz)*nspden + is]; + } + + MPI_Gatherv(vks_send,num_z[GlobalV::RANK_IN_POOL],MPI_DOUBLE,vks_receive,num_z.data(),start_z.data(),MPI_DOUBLE,0,MPI_COMM_WORLD); + MPI_Gatherv(vxc_send,num_z[GlobalV::RANK_IN_POOL],MPI_DOUBLE,vxc_receive,num_z.data(),start_z.data(),MPI_DOUBLE,0,MPI_COMM_WORLD); + + if(GlobalV::RANK_IN_POOL == 0) + { + for(int iz = 0; iz < nz; iz ++) + { + vks_collected[(ixy*nz + iz)*nspden + is] = vks_receive[iz]; + vxc_collected[(ixy*nz + iz)*nspden + is] = vxc_receive[iz]; + } + } + } + delete[] vks_send; + delete[] vxc_send; + delete[] vks_receive; + delete[] vxc_receive; + + if(GlobalV::RANK_IN_POOL == 0) + { + for(int ix = 0; ix < nx; ix ++) + { + for(int iy = 0; iy < ny; iy ++) + { + for(int iz = 0; iz < nz; iz ++) + { + int ind_c = (ix*ny*nz + iy*nz + iz)*nspden + is; + int ind_fortran = is*nfft + iz*ny*nx + iy*nx + ix; + vks_hartree[ind_fortran] = vks_collected[ind_c] / 2.0; + vxc_hartree[ind_fortran] = vxc_collected[ind_c] / 2.0; + } + } + } + calculate_dij_(natom,ntypat,ixc,xclevel,nfft,nspden,xred.data(),ucvol,gprimd.data(),vks_hartree,vxc_hartree); + } + } + + if(GlobalV::RANK_IN_POOL == 0) + { + delete[] vks_hartree; + delete[] vxc_hartree; + delete[] vks_collected; + delete[] vxc_collected; + } + +#else vks_hartree = new double[nspden * nfft]; vxc_hartree = new double[nspden * nfft]; for(int is = 0; is < nspden; is ++) @@ -354,6 +464,7 @@ void Paw_Cell::calculate_dij(double* vks, double* vxc) calculate_dij_(natom,ntypat,ixc,xclevel,nfft,nspden,xred.data(),ucvol,gprimd.data(),vks_hartree,vxc_hartree); delete[] vks_hartree; delete[] vxc_hartree; +#endif } void Paw_Cell::get_dij(int iat, int size_dij, double* dij) @@ -376,6 +487,28 @@ void Paw_Cell::init_rho(double ** rho) init_rho_(nspden, ngfftdg.data(), nfft, natom, ntypat, rprimd.data(), gprimd.data(), gmet.data(), ucvol, xred.data(), rho_tmp); +#ifdef __MPI + for(int is = 0; is < nspden; is ++) + { + // I'm not sure about this yet !!! + // need to check for nspin = 2 later + // Fortran is column major, and rhor is of dimension (nfft, nspden) + // so presumably should be this way m + for(int ix = 0; ix < nx; ix ++) + { + for(int iy = 0; iy < ny; iy ++) + { + for(int iz = 0; iz < num_z[GlobalV::RANK_IN_POOL]; iz ++) + { + int ind_c = ix*ny*num_z[GlobalV::RANK_IN_POOL] + iy*num_z[GlobalV::RANK_IN_POOL] + iz; + int ind_fortran = (iz+start_z[GlobalV::RANK_IN_POOL])*ny*nx + iy*nx + ix; + + rho[is][ind_c] = rho_tmp[ind_fortran*nspden+is]; + } + } + } + } +#else for(int ir = 0; ir < nfft; ir ++) { for(int is = 0; is < nspden; is ++) @@ -387,6 +520,7 @@ void Paw_Cell::init_rho(double ** rho) rho[is][ir] = rho_tmp[ir*nspden+is]; } } +#endif delete[] rho_tmp; } @@ -402,6 +536,10 @@ void Paw_Cell::set_dij() get_dij(iat,size_dij,dij_libpaw); +#ifdef __MPI + Parallel_Common::bcast_double(dij_libpaw,size_dij); +#endif + for(int jproj = 0; jproj < nproj; jproj ++) { for(int iproj = jproj; iproj < nproj; iproj ++) diff --git a/source/module_cell/module_paw/test/test_paw3.cpp b/source/module_cell/module_paw/test/test_paw3.cpp index 9cde5f5e766..3b87dc8a525 100644 --- a/source/module_cell/module_paw/test/test_paw3.cpp +++ b/source/module_cell/module_paw/test/test_paw3.cpp @@ -47,7 +47,10 @@ TEST_F(Test_Libpaw_Cell, test_paw) int ny_dg = 60; int nz_dg = 60; - paw_cell.set_libpaw_fft(nx, ny, nz, nx_dg, ny_dg, nz_dg); + GlobalV::NPROC = 1; + std::vector start_z = {0}; + std::vector num_z = {nz_dg}; + paw_cell.set_libpaw_fft(nx, ny, nz, nx_dg, ny_dg, nz_dg, start_z.data(), num_z.data()); int natom = 5; int ntypat = 2; @@ -152,6 +155,9 @@ TEST_F(Test_Libpaw_Cell, test_paw) nhat[0] = new double[nfft]; nhatgr = new double[nfft*3]; +// Will fix this part later; now I'm using C++ order so all the results +// are different from previously generated +/* paw_cell.get_vloc_ncoret(vloc, ncoret); std::ifstream ifs("fort.26"); @@ -241,7 +247,7 @@ TEST_F(Test_Libpaw_Cell, test_paw) delete[] dij; } - +*/ delete[] rho[0]; delete[] rho; delete[] vloc; diff --git a/source/module_cell/module_paw/test/test_paw4.cpp b/source/module_cell/module_paw/test/test_paw4.cpp index 74eb7ebbf8c..59ae79f6f69 100644 --- a/source/module_cell/module_paw/test/test_paw4.cpp +++ b/source/module_cell/module_paw/test/test_paw4.cpp @@ -147,7 +147,10 @@ TEST_F(Test_PAW, test_paw) paw_cell.set_libpaw_cell(latvec, lat0); - paw_cell.set_libpaw_fft(nx, ny, nz, nx, ny, nz); + GlobalV::NPROC = 1; + std::vector start_z = {0}; + std::vector num_z = {nz}; + paw_cell.set_libpaw_fft(nx, ny, nz, nx, ny, nz, start_z.data(), num_z.data()); paw_cell.set_libpaw_atom(nat, ntyp, typat, xred); diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 1d521b5345e..d2f82446349 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -227,7 +227,8 @@ namespace ModuleESolver { GlobalC::paw_cell.set_libpaw_ecut(INPUT.ecutwfc/2.0,INPUT.ecutwfc/2.0); //in Hartree GlobalC::paw_cell.set_libpaw_fft(this->pw_wfc->nx,this->pw_wfc->ny,this->pw_wfc->nz, - this->pw_wfc->nx,this->pw_wfc->ny,this->pw_wfc->nz); + this->pw_wfc->nx,this->pw_wfc->ny,this->pw_wfc->nz, + this->pw_wfc->startz,this->pw_wfc->numz); GlobalC::paw_cell.prepare_paw(); GlobalC::paw_cell.set_sij(); diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index ae07c26bf36..04ee43555f8 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -304,7 +304,8 @@ void ESolver_KS_PW::init_after_vc(Input& inp, UnitCell& ucell) { GlobalC::paw_cell.set_libpaw_ecut(INPUT.ecutwfc/2.0,INPUT.ecutwfc/2.0); //in Hartree GlobalC::paw_cell.set_libpaw_fft(this->pw_wfc->nx,this->pw_wfc->ny,this->pw_wfc->nz, - this->pw_wfc->nx,this->pw_wfc->ny,this->pw_wfc->nz); + this->pw_wfc->nx,this->pw_wfc->ny,this->pw_wfc->nz, + this->pw_wfc->startz,this->pw_wfc->numz); GlobalC::paw_cell.prepare_paw(); GlobalC::paw_cell.set_sij(); diff --git a/source/module_hsolver/hsolver_pw.cpp b/source/module_hsolver/hsolver_pw.cpp index cff9225bf40..632d62d600a 100644 --- a/source/module_hsolver/hsolver_pw.cpp +++ b/source/module_hsolver/hsolver_pw.cpp @@ -217,7 +217,7 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, } double* nhatgr; - nhatgr = new double[3*GlobalC::paw_cell.get_nrxx()]; + nhatgr = new double[3*GlobalC::paw_cell.get_nfft()]; GlobalC::paw_cell.get_nhat(pes->charge->nhat,nhatgr); delete[] nhatgr; } From 2f18eac0f4a09fd279df357022dfe6e8e3d74056 Mon Sep 17 00:00:00 2001 From: wenfei-li Date: Tue, 24 Oct 2023 15:57:53 +0800 Subject: [PATCH 2/4] fix cmake file for test_deepks --- .../module_deepks/test/CMakeLists.txt | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/source/module_hamilt_lcao/module_deepks/test/CMakeLists.txt b/source/module_hamilt_lcao/module_deepks/test/CMakeLists.txt index fb74157c67b..36a1a237757 100644 --- a/source/module_hamilt_lcao/module_deepks/test/CMakeLists.txt +++ b/source/module_hamilt_lcao/module_deepks/test/CMakeLists.txt @@ -13,10 +13,18 @@ target_link_libraries( base cell symmetry md surchem xc_ neighbor orb io_basic io_advanced relax gint driver esolver hsolver psi elecstate hamilt_general hamilt_pwdft hamilt_lcao tddft hamilt_ofdft hamilt_stodft planewave - pthread vdw dftu hcontainer paw + pthread vdw dftu hcontainer deepks device numerical_atomic_orbitals container psi_initializer ${ABACUS_LINK_LIBRARIES} ) + +if(ENABLE_PAW) + target_link_libraries( + test_deepks + paw + ) +endif() + if(USE_ELPA) target_link_libraries( test_deepks From a477a2ae46308b19f0e99681c0760300da475605 Mon Sep 17 00:00:00 2001 From: wenfei-li Date: Thu, 26 Oct 2023 14:34:05 +0800 Subject: [PATCH 3/4] Fix : multi-k PAW --- source/module_hsolver/hsolver_pw.cpp | 38 ++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/source/module_hsolver/hsolver_pw.cpp b/source/module_hsolver/hsolver_pw.cpp index 02748adadcd..3e18cfdc296 100644 --- a/source/module_hsolver/hsolver_pw.cpp +++ b/source/module_hsolver/hsolver_pw.cpp @@ -203,6 +203,44 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, GlobalC::paw_cell.reset_rhoij(); for (int ik = 0; ik < this->wfc_basis->nks; ++ik) { + const int npw = this->wfc_basis->npwk[ik]; + ModuleBase::Vector3 *_gk = new ModuleBase::Vector3[npw]; + for (int ig = 0;ig < npw; ig++) + { + _gk[ig] = this->wfc_basis->getgpluskcar(ik,ig); + } + + double* kpt; + kpt = new double[3]; + kpt[0] = this->wfc_basis->kvec_c[ik].x; + kpt[1] = this->wfc_basis->kvec_c[ik].y; + kpt[2] = this->wfc_basis->kvec_c[ik].z; + + double ** kpg; + kpg = new double*[npw]; + for(int ipw=0;ipwwfc_basis->get_ig2ix(ik).data(), + this->wfc_basis->get_ig2iy(ik).data(), + this->wfc_basis->get_ig2iz(ik).data(), + (const double **) kpg,GlobalC::ucell.tpiba); + + delete[] kpt; + for(int ipw = 0; ipw < npw; ipw++) + { + delete[] kpg[ipw]; + } + delete[] kpg; + + GlobalC::paw_cell.get_vkb(); + psi.fix_k(ik); GlobalC::paw_cell.set_currentk(ik); int nbands = psi.get_nbands(); From b5efd3ccd23e5c5d11bdb2ec6800d0e8b62db7cd Mon Sep 17 00:00:00 2001 From: wenfei-li Date: Thu, 26 Oct 2023 17:05:09 +0800 Subject: [PATCH 4/4] update libpaw_interface --- deps/libpaw_interface | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/deps/libpaw_interface b/deps/libpaw_interface index bdf13fb3a83..ad74d11ad71 160000 --- a/deps/libpaw_interface +++ b/deps/libpaw_interface @@ -1 +1 @@ -Subproject commit bdf13fb3a837e71c5a9ecccd2caf71aeaa3578ad +Subproject commit ad74d11ad711006a5837ce81f6e6dbd34758ff0b