diff --git a/source/module_cell/module_paw/CMakeLists.txt b/source/module_cell/module_paw/CMakeLists.txt index 694cd9e7d3..0c0af45197 100644 --- a/source/module_cell/module_paw/CMakeLists.txt +++ b/source/module_cell/module_paw/CMakeLists.txt @@ -4,6 +4,7 @@ add_library( paw_element.cpp paw_sphbes.cpp paw_cell.cpp + paw_cell_libpaw.cpp ) if(ENABLE_COVERAGE) diff --git a/source/module_cell/module_paw/paw_cell.h b/source/module_cell/module_paw/paw_cell.h index 0afe9adb00..8ecd33cd5b 100644 --- a/source/module_cell/module_paw/paw_cell.h +++ b/source/module_cell/module_paw/paw_cell.h @@ -12,6 +12,7 @@ #include "paw_element.h" #include "paw_atom.h" +#include "module_base/matrix3.h" class Paw_Cell { @@ -20,6 +21,8 @@ class Paw_Cell Paw_Cell(){}; ~Paw_Cell(){}; +// PART I. Operations in ABACUS + void init_paw_cell( const double ecutwfc_in, const double cell_factor_in, const double omega_in, @@ -151,6 +154,63 @@ class Paw_Cell void set_ylm(const int npw_in, const double ** kpg); +// Part II. Passing infor for the initialization of PAW + + public: + // The following gathers information needed by LibPAW, they will be inserted + // to proper places in the main program + + // Cutoff energies, sets ecut and ecutpaw + void set_libpaw_ecut(const double ecut_in, const double ecutpaw_in); + // Sets rprimd, gprimd, gmet and ucvol + // Input is latvec and lat0, will calculate the required info + 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); + // 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 + void set_libpaw_files(); + + // Extract the information + double get_libpaw_ecut() {return ecut;} + double get_libpaw_ecutpaw() {return ecutpaw;} + std::vector get_libpaw_rprimd() {return rprimd;} + std::vector get_libpaw_gprimd() {return gprimd;} + std::vector get_libpaw_gmet() {return gmet;} + double get_libpaw_ucvol() {return ucvol;} + std::vector get_libpaw_ngfft() {return ngfft;} + std::vector get_libpaw_ngfftdg() {return ngfftdg;} + int get_libpaw_natom() {return natom;} + int get_libpaw_ntypat() {return ntypat;} + std::vector get_libpaw_typat() {return typat;} + std::vector get_libpaw_xred() {return xred;} + char* get_libpaw_filename_list() {return filename_list;} + + private: +// Info to be passed to libpaw_interface: +// 1. ecut, ecutpaw : kinetic energy cutoff of the planewave basis set +// there will be one coarse grid for density/potential, and a fine grid for PAW +// the unit is in Hartree +// 2. rprimd, gprimd : real and reciprocal space lattice vectors, respectively +// unit for rprimd is in Bohr, and for gprimd is in Bohr^-1 +// 3. gmet : reciprocal space metric (bohr^-2) +// 4. ucvol : volume of unit cell (Bohr^3) +// 5. ngfft, ngfftdg : dimension of FFT grids of the corase and fine grids +// 6. natom, ntypat, typat: #. atoms, #. element types +// and typat records the type of each atom +// 7. xred : coordinate of each atom, in terms of rprimd (namely, direct coordinate) +// 8. filename_list : filename of the PAW xml files for each element + + double ecut, ecutpaw; + std::vector rprimd, gprimd, gmet; + double ucvol; + std::vector ngfft, ngfftdg; + int natom, ntypat; + std::vector typat; + std::vector xred; + char* filename_list; }; #endif \ No newline at end of file diff --git a/source/module_cell/module_paw/paw_cell_libpaw.cpp b/source/module_cell/module_paw/paw_cell_libpaw.cpp new file mode 100644 index 0000000000..308e2cc5c6 --- /dev/null +++ b/source/module_cell/module_paw/paw_cell_libpaw.cpp @@ -0,0 +1,206 @@ +#include "module_base/tool_quit.h" +#include "module_base/tool_title.h" +#include "paw_cell.h" +#include "module_base/global_variable.h" + +// The subroutines here are used to gather information from the main ABACUS program +// 1. ecut, ecutpaw : kinetic energy cutoff of the planewave basis set +// there will be one coarse grid for density/potential, and a fine grid for PAW +// the unit is in Hartree +// 2. rprimd, gprimd : real and reciprocal space lattice vectors, respectively +// unit for rprimd is in Bohr, and for gprimd is in Bohr^-1 +// 3. gmet : reciprocal space metric (bohr^-2) +// 4. ucvol : volume of unit cell (Bohr^3) +// 5. ngfft, ngfftdg : dimension of FFT grids of the corase and fine grids +// 6. natom, ntypat, typat: #. atoms, #. element types +// and typat records the type of each atom +// 7. xred : coordinate of each atom, in terms of rprimd (namely, direct coordinate) +// 8. filename_list : filename of the PAW xml files for each element + +// Cutoff energies, sets ecut and ecutpaw +void Paw_Cell::set_libpaw_ecut(const double ecut_in, const double ecutpaw_in) +{ + ModuleBase::TITLE("Paw_Cell", "set_libpaw_ecut"); + ecut = ecut_in; + ecutpaw = ecutpaw_in; +} + +// inverse of 3 by 3 matrix, needed by set_libpaw_cell to calculate gprimd +// adapted from m_symtk/matr3inv of ABINIT + +void matr3inv(std::vector& mat_in, std::vector& mat_out) +{ + + assert(mat_in.size() == 9); + assert(mat_out.size() == 9); + + double t1 = mat_in[4] * mat_in[8] - mat_in[7] * mat_in[5]; + double t2 = mat_in[7] * mat_in[2] - mat_in[1] * mat_in[8]; + double t3 = mat_in[1] * mat_in[5] - mat_in[4] * mat_in[2]; + double det = mat_in[0] * t1 + mat_in[3] * t2 + mat_in[6] * t3; + + double dd; + if (std::abs(det) > 1e-16) + { + dd = 1.0 / det; + } + else + { + ModuleBase::WARNING_QUIT("matr3inv", "matrix is singular!"); + } + + mat_out[0] = t1 * dd; + mat_out[3] = t2 * dd; + mat_out[6] = t3 * dd; + mat_out[1] = (mat_in[6] * mat_in[5] - mat_in[3] * mat_in[8]) * dd; + mat_out[4] = (mat_in[0] * mat_in[8] - mat_in[6] * mat_in[2]) * dd; + mat_out[7] = (mat_in[3] * mat_in[2] - mat_in[0] * mat_in[5]) * dd; + mat_out[2] = (mat_in[3] * mat_in[7] - mat_in[6] * mat_in[4]) * dd; + mat_out[5] = (mat_in[6] * mat_in[1] - mat_in[0] * mat_in[7]) * dd; + mat_out[8] = (mat_in[0] * mat_in[4] - mat_in[3] * mat_in[1]) * dd; +} + +// calculates G = A^T A for 3 by 3 matrix +// G_ij = sum_k A_ki A_kj + +void mattmat(std::vector& mat_in, std::vector& mat_out) +{ + mat_out[0] = mat_in[0] * mat_in[0] + mat_in[1] * mat_in[1] + mat_in[2] * mat_in[2]; + mat_out[1] = mat_in[0] * mat_in[3] + mat_in[1] * mat_in[4] + mat_in[2] * mat_in[5]; + mat_out[2] = mat_in[0] * mat_in[6] + mat_in[1] * mat_in[7] + mat_in[2] * mat_in[8]; + mat_out[3] = mat_in[3] * mat_in[0] + mat_in[4] * mat_in[1] + mat_in[5] * mat_in[2]; + mat_out[4] = mat_in[3] * mat_in[3] + mat_in[4] * mat_in[4] + mat_in[5] * mat_in[5]; + mat_out[5] = mat_in[3] * mat_in[6] + mat_in[4] * mat_in[7] + mat_in[5] * mat_in[8]; + mat_out[6] = mat_in[6] * mat_in[0] + mat_in[7] * mat_in[1] + mat_in[8] * mat_in[2]; + mat_out[7] = mat_in[6] * mat_in[3] + mat_in[7] * mat_in[4] + mat_in[8] * mat_in[5]; + mat_out[8] = mat_in[6] * mat_in[6] + mat_in[7] * mat_in[7] + mat_in[8] * mat_in[8]; +} + +// Sets rprimd, gprimd, gmet and ucvol +// Only real space lattice vector needed, others are to be calculated +void Paw_Cell::set_libpaw_cell(const ModuleBase::Matrix3 latvec, const double lat0) +{ + ModuleBase::TITLE("Paw_Cell", "set_libpaw_cell"); + + rprimd.resize(9); + gprimd.resize(9); + gmet.resize(9); + + rprimd[0] = latvec.e11 * lat0; + rprimd[1] = latvec.e12 * lat0; + rprimd[2] = latvec.e13 * lat0; + rprimd[3] = latvec.e21 * lat0; + rprimd[4] = latvec.e22 * lat0; + rprimd[5] = latvec.e23 * lat0; + rprimd[6] = latvec.e31 * lat0; + rprimd[7] = latvec.e32 * lat0; + rprimd[8] = latvec.e33 * lat0; + + // calculating gprimd, gmet and ucvol, adapted from m_geometry/metric of ABINIT + + // Compute first dimensional primitive translation vectors in reciprocal space + // gprimd from rprimd + // Then, computes metrics for real and recip space rmet and gmet using length + // dimensional primitive translation vectors in columns of rprimd(3,3) and gprimd(3,3). + // gprimd is the inverse transpose of rprimd. + // i.e. $ rmet_{i,j}= \sum_k ( rprimd_{k,i}*rprimd_{k,j} ) $ + // $ gmet_{i,j}= \sum_k ( gprimd_{k,i}*gprimd_{k,j} ) $ + // Also computes unit cell volume ucvol in $\textrm{bohr}^3$ + + // Compute unit cell volume + // ucvol=rprimd(1,1)*(rprimd(2,2)*rprimd(3,3)-rprimd(3,2)*rprimd(2,3))+& + // rprimd(2,1)*(rprimd(3,2)*rprimd(1,3)-rprimd(1,2)*rprimd(3,3))+& + // rprimd(3,1)*(rprimd(1,2)*rprimd(2,3)-rprimd(2,2)*rprimd(1,3)) + ucvol = rprimd[0] * (rprimd[4] * rprimd[8] - rprimd[7] * rprimd[5]) + + rprimd[3] * (rprimd[7] * rprimd[2] - rprimd[1] * rprimd[8]) + + rprimd[6] * (rprimd[1] * rprimd[5] - rprimd[4] * rprimd[2]); + // ucvol = det3r(rprimd) + + if (std::abs(ucvol) < 1e-12) + { + ModuleBase::WARNING_QUIT("set_libpaw_cell", "ucvol vanishingly small"); + } + + // ABACUS allows negative volume, but it seems ABINIT do not + // I am not quite sure why, but to avoid possible complications I will follow ABINIT here + // as the user can always just exchange two axis to make it positive + if (ucvol < 0) + { + ModuleBase::WARNING_QUIT("set_libpaw_cell", "ucvol negative, one way to solve is to exchange two cell axis"); + } + + // Generate gprimd + matr3inv(rprimd, gprimd); + + // Compute reciprocal space metric. + mattmat(gprimd, gmet); +} + +// 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) +{ + ngfft.resize(3); + ngfftdg.resize(3); + + ngfft[0] = nx_in; + ngfft[1] = ny_in; + ngfft[2] = nz_in; + ngfftdg[0] = nxdg_in; + ngfftdg[1] = nydg_in; + ngfftdg[2] = nzdg_in; +} + +// Sets natom, ntypat, typat and xred +// !!!!!!!Note : the index stored in typat here will start from 1, not 0 !!!!!! +void Paw_Cell::set_libpaw_atom(const int natom_in, const int ntypat_in, const int* typat_in, const double* xred_in) +{ + natom = natom_in; + ntypat = ntypat_in; + + typat.resize(natom); + xred.resize(3 * natom); + + for(int iat = 0; iat < natom; iat ++) + { + typat[iat] = typat_in[iat]; + for(int j = 0; j < 3; j ++) + { + xred[3 * iat + j] = xred_in[3 * iat + j]; + } + } +} + +// Sets filename_list +// I'm going to read directly from STRU file +void Paw_Cell::set_libpaw_files() +{ + ModuleBase::TITLE("Paw_Cell", "set_libpaw_files"); + + if(GlobalV::MY_RANK == 0) + { + std::ifstream ifa(GlobalV::stru_file.c_str(), std::ios::in); + if (!ifa) + { + ModuleBase::WARNING_QUIT("set_libpaw_files", "can not open stru file"); + } + + std::string line; + while(!ifa.eof()) + { + getline(ifa,line); + if (line.find("PAW_FILES") != std::string::npos) break; + } + + filename_list = new char[ntypat*264]; + for(int i = 0; i < ntypat; i++) + { + std::string filename; + ifa >> filename; + for(int j = 0; j < filename.size(); j++) + { + filename_list[264*i+j] = filename[j]; + } + } + } +} diff --git a/source/module_cell/module_paw/test/CMakeLists.txt b/source/module_cell/module_paw/test/CMakeLists.txt index bee39a448d..501c99354d 100644 --- a/source/module_cell/module_paw/test/CMakeLists.txt +++ b/source/module_cell/module_paw/test/CMakeLists.txt @@ -6,6 +6,7 @@ install(FILES sphbes_ref.dat func.dat qlist.dat fq_ref.dat eigts.dat ca.dat rhoij.dat vkb_ref.dat Si_test.xml gnorm.dat ptilde_ref.dat psi.dat rhoij1.dat H.LDA_PW-JTH.xml Fe.GGA_PBE-JTH.xml O.GGA_PBE-JTH.xml Si.GGA_PBE-JTH.xml + STRU DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) AddTest( @@ -24,4 +25,10 @@ AddTest( TARGET Test_Paw2 LIBS ${math_libs} base device SOURCES test_paw2.cpp ../paw_atom.cpp +) + +AddTest( + TARGET Test_Paw3 + LIBS ${math_libs} base device + SOURCES test_paw3.cpp ../paw_atom.cpp ../paw_cell_libpaw.cpp ) \ No newline at end of file diff --git a/source/module_cell/module_paw/test/STRU b/source/module_cell/module_paw/test/STRU new file mode 100644 index 0000000000..48299f50ae --- /dev/null +++ b/source/module_cell/module_paw/test/STRU @@ -0,0 +1,3 @@ +PAW_FILES +C.xml +H.xml diff --git a/source/module_cell/module_paw/test/test_paw.cpp b/source/module_cell/module_paw/test/test_paw.cpp index 8439900c61..530d73736f 100644 --- a/source/module_cell/module_paw/test/test_paw.cpp +++ b/source/module_cell/module_paw/test/test_paw.cpp @@ -3,6 +3,19 @@ #include "../paw_element.h" +/* +Subroutines related to the processing of PAW xml files: + +1. init_paw_element, which prepares some related values +2. read_paw_xml, which reads the xml file and stores the information + +plus subroutines related to the processing of projector functions: +1. spherical_bessel_function: helper functions, returns values of spherical bessel functions +2. spherical_bessel_transform: carries out spherical bessel transform +3. splint and spline: performs cubic spline interpolation + +*/ + class Test_Read_Paw : public testing::Test { protected: diff --git a/source/module_cell/module_paw/test/test_paw1.cpp b/source/module_cell/module_paw/test/test_paw1.cpp index d7d5f51213..9814de3d54 100644 --- a/source/module_cell/module_paw/test/test_paw1.cpp +++ b/source/module_cell/module_paw/test/test_paw1.cpp @@ -3,6 +3,19 @@ #include "../paw_cell.h" +/* +Unit Test for subroutines related to the calculation of rhoij, +where reciprocal space wavefunctions are passed from outside, and rhoij is calculated: + +1. init_paw_cell and set_paw_k, which collects necessary information +2. accumulate_rhoij, which calculates sum_n + +plus two mathematical subroutines: +3. calc_ylm, which gives values of spherical harmonics +4. ass_leg_pol, which gives values of legendre polynomials + +*/ + class Test_Paw_Cell : public testing::Test { protected: diff --git a/source/module_cell/module_paw/test/test_paw2.cpp b/source/module_cell/module_paw/test/test_paw2.cpp index 4491484311..c63212c47a 100644 --- a/source/module_cell/module_paw/test/test_paw2.cpp +++ b/source/module_cell/module_paw/test/test_paw2.cpp @@ -3,6 +3,15 @@ #include #include "../paw_atom.h" +/* + +Unit Test for subroutines related to the processing of rhoij, +the on-site density matrix to be used for PAW: +1. set_ca, which passes from outside and saves it +2. accumulate_rhoij, which accumulates the contribution of one band +3. convert_rhoij, which converts rhoij to format required by PAW + +*/ class Test_Paw_Atom : public testing::Test { diff --git a/source/module_cell/module_paw/test/test_paw3.cpp b/source/module_cell/module_paw/test/test_paw3.cpp new file mode 100644 index 0000000000..cc5e258d8a --- /dev/null +++ b/source/module_cell/module_paw/test/test_paw3.cpp @@ -0,0 +1,115 @@ +#include "gtest/gtest.h" +#include +#include + +#include "../paw_cell.h" + +/* + +Unit Test for the following subroutines, which are used to pass information +from main ABACUS program to LibPAW: + +1. set_libpaw_ecut, which sets kinetic energy cutoff +2. set_libpaw_cell, which sets quantities related to cell parameters +3. set_libpaw_fft, which sets the real-space FFT grid +4. set_libpaw_atom, which sets information of atoms in the unit cell +5. set_libpaw_files, which sets the names of PAW xml files + +*/ + +class Test_Libpaw_Cell : public testing::Test +{ + protected: + + Paw_Cell paw_cell; +}; + +TEST_F(Test_Libpaw_Cell, test_paw) +{ + ModuleBase::Matrix3 latvec; + + double ecut = 30.0; + paw_cell.set_libpaw_ecut(ecut,ecut); + + latvec.e11 = 0.5; latvec.e12 = 0.5; latvec.e13 = 0.0; + latvec.e21 = 0.0; latvec.e22 = 0.5; latvec.e23 = 0.5; + latvec.e31 = 0.5; latvec.e32 = 0.0; latvec.e33 = 0.5; + + double lat0 = 10.2; + + paw_cell.set_libpaw_cell(latvec, lat0); + + int nx = 30; + int ny = 30; + int nz = 30; + + paw_cell.set_libpaw_fft(nx, ny, nz, nx, ny, nz); + + int natom = 2; + int ntypat = 2; + int typat[2] = {1,2}; + double xred[6] = {0.0, 0.0, 0.0, 0.5, 0.5, 0.5}; + + paw_cell.set_libpaw_atom(natom, ntypat, typat, xred); + + paw_cell.set_libpaw_files(); + + EXPECT_NEAR(paw_cell.get_libpaw_ecut(),30.0,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_ecutpaw(),30.0,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_ucvol(),265.302,1e-10); + + EXPECT_NEAR(paw_cell.get_libpaw_rprimd()[0],5.1,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_rprimd()[1],5.1,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_rprimd()[2],0.0,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_rprimd()[3],0.0,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_rprimd()[4],5.1,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_rprimd()[5],5.1,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_rprimd()[6],5.1,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_rprimd()[7],0.0,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_rprimd()[8],5.1,1e-10); + + EXPECT_NEAR(paw_cell.get_libpaw_gprimd()[0],0.09803921568,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_gprimd()[1],0.09803921568,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_gprimd()[2],-0.09803921568,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_gprimd()[3],-0.09803921568,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_gprimd()[4],0.09803921568,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_gprimd()[5],0.09803921568,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_gprimd()[6],0.09803921568,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_gprimd()[7],-0.09803921568,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_gprimd()[8],0.09803921568,1e-10); + + EXPECT_NEAR(paw_cell.get_libpaw_gmet()[0],0.02883506343,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_gmet()[1],-0.00961168781,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_gmet()[2],-0.00961168781,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_gmet()[3],-0.00961168781,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_gmet()[4],0.02883506343,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_gmet()[5],-0.00961168781,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_gmet()[6],-0.00961168781,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_gmet()[7],-0.00961168781,1e-10); + EXPECT_NEAR(paw_cell.get_libpaw_gmet()[8],0.02883506343,1e-10); + + EXPECT_EQ(paw_cell.get_libpaw_ngfft()[0],30); + EXPECT_EQ(paw_cell.get_libpaw_ngfft()[1],30); + EXPECT_EQ(paw_cell.get_libpaw_ngfft()[2],30); + + EXPECT_EQ(paw_cell.get_libpaw_ngfftdg()[0],30); + EXPECT_EQ(paw_cell.get_libpaw_ngfftdg()[1],30); + EXPECT_EQ(paw_cell.get_libpaw_ngfftdg()[2],30); + + EXPECT_EQ(paw_cell.get_libpaw_natom(),2); + EXPECT_EQ(paw_cell.get_libpaw_ntypat(),2); + + EXPECT_EQ(paw_cell.get_libpaw_typat()[0],1); + EXPECT_EQ(paw_cell.get_libpaw_typat()[1],2); + + EXPECT_EQ(paw_cell.get_libpaw_filename_list()[0],'C'); + EXPECT_EQ(paw_cell.get_libpaw_filename_list()[1],'.'); + EXPECT_EQ(paw_cell.get_libpaw_filename_list()[2],'x'); + EXPECT_EQ(paw_cell.get_libpaw_filename_list()[3],'m'); + EXPECT_EQ(paw_cell.get_libpaw_filename_list()[4],'l'); + EXPECT_EQ(paw_cell.get_libpaw_filename_list()[264],'H'); + EXPECT_EQ(paw_cell.get_libpaw_filename_list()[265],'.'); + EXPECT_EQ(paw_cell.get_libpaw_filename_list()[266],'x'); + EXPECT_EQ(paw_cell.get_libpaw_filename_list()[267],'m'); + EXPECT_EQ(paw_cell.get_libpaw_filename_list()[268],'l'); +} \ No newline at end of file