From afcfe07a096074596588bed4b0ef8863e376212b Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Thu, 21 Nov 2024 17:58:55 +0800 Subject: [PATCH 01/10] initial commit --- source/module_io/cal_pLpR.cpp | 196 ++++++++++++++++++++++++++++++++++ source/module_io/cal_pLpR.h | 67 ++++++++++++ 2 files changed, 263 insertions(+) create mode 100644 source/module_io/cal_pLpR.cpp create mode 100644 source/module_io/cal_pLpR.h diff --git a/source/module_io/cal_pLpR.cpp b/source/module_io/cal_pLpR.cpp new file mode 100644 index 00000000000..9e8f194194a --- /dev/null +++ b/source/module_io/cal_pLpR.cpp @@ -0,0 +1,196 @@ +#include +#include +#include +#include +#include +#include +#include +#include "module_cell/unitcell.h" +#include "module_base/spherical_bessel_transformer.h" +#include "module_basis/module_nao/two_center_integrator.h" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_cell/module_neighbor/sltk_atom_arrange.h" +#include "module_parameter/parameter.h" + +void _radial_indexing(const UnitCell& ucell, + std::vector>& lin2comp, + std::map, int>& comp2lin) +{ + int i = 0; + for (int it = 0; it < ucell.ntype; it++) + { + const Atom& atom = ucell.atoms[it]; + const int lmax = atom.nwl; + for (int ia = 0; ia < atom.na; ia++) + { + for (int l = 0; l < lmax; l++) + { + const int nzeta = atom.l_nchi[l]; + for (int iz = 0; iz < nzeta; iz++) + { + for (int m = -l; m <= l; m++) + { + lin2comp.push_back(std::make_tuple(it, ia, l, iz, m)); + comp2lin[std::make_tuple(it, ia, l, iz, m)] = i; + i++; + } + } + } + } + } +} + +/** + * + * FIXME: the following part will be transfered to TwoCenterIntegrator soon + * + */ + +// L+|l, m> = sqrt((l-m)(l+m+1))|l, m+1>, return the sqrt((l-m)(l+m+1)) +double _lplus_on_ylm(const int l, const int m) +{ + return std::sqrt((l - m) * (l + m + 1)); +} + +// L-|l, m> = sqrt((l+m)(l-m+1))|l, m-1>, return the sqrt((l+m)(l-m+1)) +double _lminus_on_ylm(const int l, const int m) +{ + return std::sqrt((l + m) * (l - m + 1)); +} + +void _cal_pLzpR(const std::unique_ptr& calculator, + const int it, const int ia, const int il, const int iz, const int mi, + const int jt, const int ja, const int jl, const int jz, const int mj, + const ModuleBase::Vector3& vR, + std::complex& val) +{ + double val_ = 0; + calculator->calculate(it, il, iz, mi, jt, jl, jz, mj, vR, &val_); + val = std::complex(mi) * val_; +} + +void _cal_pLypR(const std::unique_ptr& calculator, + const int it, const int ia, const int il, const int iz, const int im, + const int jt, const int ja, const int jl, const int jz, const int jm, + const ModuleBase::Vector3& vR, + std::complex& val) +{ + // Ly = -i/2 * (L+ - L-) + const double plus_ = _lplus_on_ylm(jl, jm); + const double minus_ = _lminus_on_ylm(jl, jm); + double val_plus = 0, val_minus = 0; + if (plus_ != 0) + { + calculator->calculate(it, il, iz, im, jt, jl, jz, jm + 1, vR, &val_plus); + val_plus *= plus_; + } + if (minus_ != 0) + { + calculator->calculate(it, il, iz, im, jt, jl, jz, jm - 1, vR, &val_minus); + val_minus *= minus_; + } + val = std::complex(0, -0.5) * (val_plus - val_minus); +} + +void _cal_pLxpR(const std::unique_ptr& calculator, + const int it, const int ia, const int il, const int iz, const int im, + const int jt, const int ja, const int jl, const int jz, const int jm, + const ModuleBase::Vector3& vR, + std::complex& val) +{ + // Lx = 1/2 * (L+ + L-) + const double plus_ = _lplus_on_ylm(jl, jm); + const double minus_ = _lminus_on_ylm(jl, jm); + double val_plus = 0, val_minus = 0; + if (plus_ != 0) + { + calculator->calculate(it, il, iz, im, jt, jl, jz, jm + 1, vR, &val_plus); + val_plus *= plus_; + } + if (minus_ != 0) + { + calculator->calculate(it, il, iz, im, jt, jl, jz, jm - 1, vR, &val_minus); + val_minus *= minus_; + } + val = std::complex(0.5) * (val_plus + val_minus); +} + +void _cal_pLpR(const std::unique_ptr& calculator, + const int it, const int ia, const int il, const int iz, const int im, + const int jt, const int ja, const int jl, const int jz, const int jm, + const ModuleBase::Vector3& vR, + std::complex& val, + const char dir = 'z') +{ + switch (dir) + { + case 'z': + _cal_pLzpR(calculator, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR, val); + break; + case 'y': + _cal_pLypR(calculator, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR, val); + break; + case 'x': + _cal_pLxpR(calculator, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR, val); + break; + default: + break; + } +} + +// calculate the matrix of +void calculate(const std::unique_ptr& calculator, + const UnitCell& ucell, + const double rcut_max, + std::vector& out, + std::ofstream* ptrlog = nullptr) +{ + // I follow the void ESolver_LJ::runner to find the neighboring atoms + const int t_atom = PARAM.inp.test_atom_input; + const int t_deconstructor = PARAM.inp.test_deconstructor; + const int t_grid = PARAM.inp.test_grid; + Grid_Driver neighbor_searcher(t_deconstructor, t_grid); + atom_arrange::search(PARAM.inp.search_pbc, + *(ptrlog), + neighbor_searcher, + ucell, + rcut_max, + t_atom); + + // calculate the matrix elements + +} + +/** + * @brief Write a 2D matrix to a plain text file. + * + * @tparam T the type of the matrix elements + * @param matrix the matrix that will be written + * @param ncols the number of columns of the matrix + * @param ofs the output file stream + * @param oflog the log file stream + * @param rank the identifier to distinguish different MPI processes + */ +template +void _write_plain_matrix(const std::vector& matrix, + const int ncols, + std::ofstream& ofs, + std::ofstream& oflog, + const int precision = 8, + const int rank = 0) +{ + if (!ofs.good()) { + return; + } + + if (rank == 0) { + for (int i = 0; i < matrix.size(); i++) { + ofs << std::setw(precision + 6) << std::right << std::scientific << matrix[i]; + if ((i + 1) % ncols == 0) { + ofs << std::endl; + } else { + ofs << " "; + } + } + } +} \ No newline at end of file diff --git a/source/module_io/cal_pLpR.h b/source/module_io/cal_pLpR.h new file mode 100644 index 00000000000..e158651f5c9 --- /dev/null +++ b/source/module_io/cal_pLpR.h @@ -0,0 +1,67 @@ +/** + * calculate the matrix elements + * + * Formulation + * ----------- + * + * Calculate the with ladder operator L+ and L-. + * + * The relation between Lx, Ly and L+, L- are: + * + * Lx = (L+ + L-) / 2 + * Ly = (L+ - L-) / 2i + * + * With L+, the spherical harmonic function Ylm (denoted as |l, m> in the following) + * can be raised: + * + * L+|l, m> = sqrt((l-m)(l+m+1))|l, m+1> + * + * Likely, with L-, the spherical harmonic function Ylm can be lowered: + * + * L-|l, m> = sqrt((l+m)(l-m+1))|l, m-1> + * + * Therefore the Lx matrix element can be calculated as: + * + * = sqrt((l-m)(l+m+1)) * delta(m, m'+1) / 2 + * + sqrt((l+m)(l-m+1)) * delta(m, m'-1) / 2 + * + * The Ly matrix element can be calculated as: + * + * = sqrt((l-m)(l+m+1)) * delta(m, m'+1) / 2i + * - sqrt((l+m)(l-m+1)) * delta(m, m'-1) / 2i + * + * The Lz matrix element can be calculated as: + * + * = m * delta(m, m') + * + * However, things will change when there are more than one centers. + * + * Technical Details + * ----------------- + * + * 0. The calculation of matrix elements involves the two-center-integral calculation, + * this is supported by ABACUS built-in class TwoCenterIntegrator. + * see: source/module_basis/module_nao/two_center_integrator.h. + * + * 1. The interface of it is RadialCollection, which is a collection of radial functions. + * see: source/module_basis/module_nao/radial_collection.h + * + * 2. The radial functions are stored in AtomicRadials class, + * see: source/module_basis/module_nao/atomic_radials.h + * + * 3. The construction of AtomicRadials involves the filename of orbital, it is stored + * in the UnitCell instance + * + * + */ +#include +#include +#include +#include +#include +#include "module_cell/unitcell.h" +#include "module_basis/module_nao/two_center_integrator.h" +namespace ModuleIO +{ + +} // namespace ModuleIO \ No newline at end of file From 7221861d8e2e9b1db48077a330ab554c7b9831f4 Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Tue, 1 Apr 2025 19:17:25 +0800 Subject: [PATCH 02/10] temp save --- source/module_io/CMakeLists.txt | 1 + source/module_io/cal_pLpR.cpp | 191 ++++++++++-------------- source/module_io/cal_pLpR.h | 70 ++++++++- source/module_io/test/CMakeLists.txt | 26 +++- source/module_io/test/cal_pLpR_test.cpp | 168 +++++++++++++++++++++ 5 files changed, 341 insertions(+), 115 deletions(-) create mode 100644 source/module_io/test/cal_pLpR_test.cpp diff --git a/source/module_io/CMakeLists.txt b/source/module_io/CMakeLists.txt index e9f7d0b9822..8fee206c7e0 100644 --- a/source/module_io/CMakeLists.txt +++ b/source/module_io/CMakeLists.txt @@ -67,6 +67,7 @@ if(ENABLE_LCAO) output_dmk.cpp output_mulliken.cpp io_npz.cpp + cal_pLpR.cpp ) list(APPEND objects_advanced unk_overlap_lcao.cpp diff --git a/source/module_io/cal_pLpR.cpp b/source/module_io/cal_pLpR.cpp index 9e8f194194a..8423078af22 100644 --- a/source/module_io/cal_pLpR.cpp +++ b/source/module_io/cal_pLpR.cpp @@ -11,34 +11,7 @@ #include "module_cell/module_neighbor/sltk_grid_driver.h" #include "module_cell/module_neighbor/sltk_atom_arrange.h" #include "module_parameter/parameter.h" - -void _radial_indexing(const UnitCell& ucell, - std::vector>& lin2comp, - std::map, int>& comp2lin) -{ - int i = 0; - for (int it = 0; it < ucell.ntype; it++) - { - const Atom& atom = ucell.atoms[it]; - const int lmax = atom.nwl; - for (int ia = 0; ia < atom.na; ia++) - { - for (int l = 0; l < lmax; l++) - { - const int nzeta = atom.l_nchi[l]; - for (int iz = 0; iz < nzeta; iz++) - { - for (int m = -l; m <= l; m++) - { - lin2comp.push_back(std::make_tuple(it, ia, l, iz, m)); - comp2lin[std::make_tuple(it, ia, l, iz, m)] = i; - i++; - } - } - } - } - } -} +#include "module_io/cal_pLpR.h" /** * @@ -58,22 +31,22 @@ double _lminus_on_ylm(const int l, const int m) return std::sqrt((l + m) * (l - m + 1)); } -void _cal_pLzpR(const std::unique_ptr& calculator, - const int it, const int ia, const int il, const int iz, const int mi, - const int jt, const int ja, const int jl, const int jz, const int mj, - const ModuleBase::Vector3& vR, - std::complex& val) +std::complex ModuleIO::cal_LzijR( + const std::unique_ptr& calculator, + const int it, const int ia, const int il, const int iz, const int mi, + const int jt, const int ja, const int jl, const int jz, const int mj, + const ModuleBase::Vector3& vR) { double val_ = 0; calculator->calculate(it, il, iz, mi, jt, jl, jz, mj, vR, &val_); - val = std::complex(mi) * val_; + return std::complex(mi) * val_; } -void _cal_pLypR(const std::unique_ptr& calculator, - const int it, const int ia, const int il, const int iz, const int im, - const int jt, const int ja, const int jl, const int jz, const int jm, - const ModuleBase::Vector3& vR, - std::complex& val) +std::complex ModuleIO::cal_LyijR( + const std::unique_ptr& calculator, + const int it, const int ia, const int il, const int iz, const int im, + const int jt, const int ja, const int jl, const int jz, const int jm, + const ModuleBase::Vector3& vR) { // Ly = -i/2 * (L+ - L-) const double plus_ = _lplus_on_ylm(jl, jm); @@ -89,14 +62,14 @@ void _cal_pLypR(const std::unique_ptr& calculator, calculator->calculate(it, il, iz, im, jt, jl, jz, jm - 1, vR, &val_minus); val_minus *= minus_; } - val = std::complex(0, -0.5) * (val_plus - val_minus); + return std::complex(0, -0.5) * (val_plus - val_minus); } -void _cal_pLxpR(const std::unique_ptr& calculator, - const int it, const int ia, const int il, const int iz, const int im, - const int jt, const int ja, const int jl, const int jz, const int jm, - const ModuleBase::Vector3& vR, - std::complex& val) +std::complex ModuleIO::cal_LxijR( + const std::unique_ptr& calculator, + const int it, const int ia, const int il, const int iz, const int im, + const int jt, const int ja, const int jl, const int jz, const int jm, + const ModuleBase::Vector3& vR) { // Lx = 1/2 * (L+ + L-) const double plus_ = _lplus_on_ylm(jl, jm); @@ -112,84 +85,78 @@ void _cal_pLxpR(const std::unique_ptr& calculator, calculator->calculate(it, il, iz, im, jt, jl, jz, jm - 1, vR, &val_minus); val_minus *= minus_; } - val = std::complex(0.5) * (val_plus + val_minus); + return std::complex(0.5) * (val_plus + val_minus); } -void _cal_pLpR(const std::unique_ptr& calculator, - const int it, const int ia, const int il, const int iz, const int im, - const int jt, const int ja, const int jl, const int jz, const int jm, - const ModuleBase::Vector3& vR, - std::complex& val, - const char dir = 'z') +ModuleIO::AngularMomentumExpectationCalculator::AngularMomentumExpectationCalculator( + const std::string& orbital_dir, + const UnitCell& ucell, + const double& search_radius, + const int tdestructor, + const int tgrid, + const int tatom, + const bool searchpbc, + std::ofstream* ptr_log) { - switch (dir) + std::vector forb(ucell.ntype); + for (int i = 0; i < ucell.ntype; ++i) { - case 'z': - _cal_pLzpR(calculator, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR, val); - break; - case 'y': - _cal_pLypR(calculator, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR, val); - break; - case 'x': - _cal_pLxpR(calculator, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR, val); - break; - default: - break; + forb[i] = orbital_dir + ucell.orbital_fn[i]; } -} + this->orb_ = std::unique_ptr(new RadialCollection); + this->orb_->build(ucell.ntype, forb.data(), 'o'); -// calculate the matrix of -void calculate(const std::unique_ptr& calculator, - const UnitCell& ucell, - const double rcut_max, - std::vector& out, - std::ofstream* ptrlog = nullptr) -{ - // I follow the void ESolver_LJ::runner to find the neighboring atoms - const int t_atom = PARAM.inp.test_atom_input; - const int t_deconstructor = PARAM.inp.test_deconstructor; - const int t_grid = PARAM.inp.test_grid; - Grid_Driver neighbor_searcher(t_deconstructor, t_grid); - atom_arrange::search(PARAM.inp.search_pbc, - *(ptrlog), - neighbor_searcher, - ucell, - rcut_max, - t_atom); - - // calculate the matrix elements - + ModuleBase::SphericalBesselTransformer sbt(true); + this->orb_->set_transformer(sbt); + + const double rcut_max = orb_->rcut_max(); + const int ngrid = int(rcut_max / 0.01) + 1; + const double cutoff = 2.0 * rcut_max; + this->orb_->set_uniform_grid(true, ngrid, cutoff, 'i', true); + + this->calculator_ = std::unique_ptr(new TwoCenterIntegrator); + this->calculator_->tabulate(*orb_, *orb_, 'S', ngrid, cutoff); + + // Initialize Ylm coefficients + ModuleBase::Ylm::set_coefficients(); + + // ofs_running + this->ofs_ = ptr_log; + + // for neighbor list search + this->neighbor_searcher_ = std::unique_ptr(new Grid_Driver(tdestructor, tgrid)); + atom_arrange::search( + searchpbc, + *ofs_, + *neighbor_searcher_, + ucell, + search_radius, + tatom); } -/** - * @brief Write a 2D matrix to a plain text file. - * - * @tparam T the type of the matrix elements - * @param matrix the matrix that will be written - * @param ncols the number of columns of the matrix - * @param ofs the output file stream - * @param oflog the log file stream - * @param rank the identifier to distinguish different MPI processes - */ -template -void _write_plain_matrix(const std::vector& matrix, - const int ncols, - std::ofstream& ofs, - std::ofstream& oflog, - const int precision = 8, - const int rank = 0) +void ModuleIO::AngularMomentumExpectationCalculator::calculate( + std::ofstream* ofs, + const UnitCell& ucell, + const char dir) { - if (!ofs.good()) { + if (ofs == nullptr) + { return; } - - if (rank == 0) { - for (int i = 0; i < matrix.size(); i++) { - ofs << std::setw(precision + 6) << std::right << std::scientific << matrix[i]; - if ((i + 1) % ncols == 0) { - ofs << std::endl; - } else { - ofs << " "; + ModuleBase::Vector3 ri, rj, dr; + for (int it = 0; it < ucell.ntype; it++) + { + const Atom& atyp_i = ucell.atoms[it]; + for (int ia = 0; ia < atyp_i.na; ia++) + { + ri = atyp_i.tau[ia]; + neighbor_searcher_->Find_atom(ucell, ri, it, ia); + for (int ia_adj = 0; ia_adj < neighbor_searcher_->getAdjacentNum(); ia_adj++) + { + rj = neighbor_searcher_->getAdjacentTau(ia_adj); + int jt = neighbor_searcher_->getType(ia_adj); + dr = (ri - rj) * ucell.lat0; + const ModuleBase::Vector3 iR = neighbor_searcher_->getBox(ia_adj); } } } diff --git a/source/module_io/cal_pLpR.h b/source/module_io/cal_pLpR.h index e158651f5c9..654f016a076 100644 --- a/source/module_io/cal_pLpR.h +++ b/source/module_io/cal_pLpR.h @@ -1,5 +1,7 @@ /** - * calculate the matrix elements + * calculate the matrix elements, in which the Lx/Ly/Lz + * are the angular momentum operators, |phi_i> and |phi_j> are the numerical + * atomic orbitals (NAOs). * * Formulation * ----------- @@ -52,16 +54,80 @@ * 3. The construction of AtomicRadials involves the filename of orbital, it is stored * in the UnitCell instance * - * + * 4. The calculation will run over all supercells in which the two-center-integral + * is not zero. This is done by the class SltkGridDriver, which is a driver for + * searching the neighboring cells. */ #include #include #include #include #include +#include #include "module_cell/unitcell.h" #include "module_basis/module_nao/two_center_integrator.h" + namespace ModuleIO { + std::complex cal_LzijR( + const std::unique_ptr& calculator, + const int it, const int ia, const int il, const int iz, const int mi, + const int jt, const int ja, const int jl, const int jz, const int mj, + const ModuleBase::Vector3& vR); + + std::complex cal_LyijR( + const std::unique_ptr& calculator, + const int it, const int ia, const int il, const int iz, const int im, + const int jt, const int ja, const int jl, const int jz, const int jm, + const ModuleBase::Vector3& vR); + + std::complex cal_LxijR( + const std::unique_ptr& calculator, + const int it, const int ia, const int il, const int iz, const int im, + const int jt, const int ja, const int jl, const int jz, const int jm, + const ModuleBase::Vector3& vR); + + // the calculation of matrix elements will be outputted + // in the way that indexed by: + // it, ia, il, iz, im, iRx, iRy, iRz, jt, ja, jl, jz, jm, in which the + // iRx, iRy, iRz are the indices of the supercell in which the two-center-integral + // it and jt are indexes of atomtypes, + // ia and ja are indexes of atoms within the atomtypes, + // il and jl are indexes of the angular momentum, + // iz and jz are indexes of the zeta functions + // im and jm are indexes of the magnetic quantum numbers. + // The output is a complex number, which is the value of the matrix element. + // Always the matrix is quite large, so direct print to file. + class AngularMomentumExpectationCalculator + { + public: + AngularMomentumExpectationCalculator() = delete; // the default constructor is meaningless + AngularMomentumExpectationCalculator( + const std::string& orbital_dir, + const UnitCell& ucell, + const double& search_radius, + const int tdestructor, + const int tgrid, + const int tatom, + const bool searchpbc, + std::ofstream* ptr_log = nullptr); + ~AngularMomentumExpectationCalculator() = default; + + void calculate(std::ofstream* ofs, + const UnitCell& ucell, + const char dir = 'x'); + + private: + // ofsrunning + std::ofstream* ofs_; + // the two-center-integrator + std::unique_ptr calculator_; + // the spherical bessel transformer + ModuleBase::SphericalBesselTransformer sbt_; + // the radial collection + std::unique_ptr orb_; + // neighboring searcher + std::unique_ptr neighbor_searcher_; + }; } // namespace ModuleIO \ No newline at end of file diff --git a/source/module_io/test/CMakeLists.txt b/source/module_io/test/CMakeLists.txt index 543c4f363c2..19068397e07 100644 --- a/source/module_io/test/CMakeLists.txt +++ b/source/module_io/test/CMakeLists.txt @@ -237,4 +237,28 @@ AddTest( add_test(NAME orb_io_test_parallel COMMAND mpirun -np 4 ./orb_io_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} -) \ No newline at end of file +) + +if(ENABLE_LCAO) +AddTest( + TARGET cal_pLpR_test + LIBS parameter base ${math_libs} device + SOURCES + cal_pLpR_test.cpp + ../cal_pLpR.cpp + ../../module_basis/module_ao/ORB_atomic_lm.cpp + ../../module_basis/module_ao/ORB_atomic.cpp + ../../module_basis/module_nao/radial_set.cpp + ../../module_basis/module_nao/numerical_radial.cpp + ../../module_basis/module_nao/beta_radials.cpp + ../../module_basis/module_nao/hydrogen_radials.cpp + ../../module_basis/module_nao/pswfc_radials.cpp + ../../module_basis/module_nao/atomic_radials.cpp + ../../module_basis/module_nao/sphbes_radials.cpp + ../../module_basis/module_nao/projgen.cpp + ../../module_basis/module_nao/two_center_integrator.cpp + ../../module_basis/module_nao/two_center_table.cpp + ../../module_basis/module_nao/real_gaunt_table.cpp + ../../module_basis/module_nao/radial_collection.cpp +) +endif() \ No newline at end of file diff --git a/source/module_io/test/cal_pLpR_test.cpp b/source/module_io/test/cal_pLpR_test.cpp new file mode 100644 index 00000000000..d9f9b92085e --- /dev/null +++ b/source/module_io/test/cal_pLpR_test.cpp @@ -0,0 +1,168 @@ +#include +#include +#include +#include +#include + +#include "module_io/cal_pLpR.h" +#include "module_basis/module_nao/two_center_integrator.h" +#include "module_basis/module_nao/radial_collection.h" +#include "module_base/spherical_bessel_transformer.h" +#include "module_base/ylm.h" + +#define DOUBLETHRESHOLD 1e-12 + +class CalpLpRTest : public ::testing::Test +{ +protected: + void SetUp() override + { + ModuleBase::SphericalBesselTransformer sbt(true); + + this->orb_ = std::unique_ptr(new RadialCollection); + this->orb_->build(1, &forb_, 'o'); + this->orb_->set_transformer(sbt); + + const double rcut = 8.0; + const int ngrid = int(rcut / 0.01) + 1; + const double cutoff = 2.0 * rcut; + this->orb_->set_uniform_grid(true, ngrid, cutoff, 'i', true); + + this->calculator_ = std::unique_ptr(new TwoCenterIntegrator); + this->calculator_->tabulate(*orb_, *orb_, 'S', ngrid, cutoff); + + // Initialize Ylm coefficients + ModuleBase::Ylm::set_coefficients(); + } + + void TearDown() override + { + // Cleanup code here, if needed + } + + const std::string forb_ = "../../../../tests/PP_ORB/Si_gga_8au_100Ry_2s2p1d.orb"; + + std::unique_ptr calculator_; + std::unique_ptr orb_; +}; + +TEST_F(CalpLpRTest, CalLzijRTest) +{ + std::complex out; + + ModuleBase::Vector3 vR(0., 0., 0.); // home-cell + int it = 0, ia = 0, il = 0, iz = 0, mi = 0; + int jt = 0, ja = 0, jl = 0, jz = 0, mj = 0; // self, the first s + + // = 0: no magnetic moment + out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); + EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); + // = 0: orthogonal + jl = 1; + out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); + EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); + // = 0: orthogonal + il = 1; mi = -1; jl = 1; mj = 0; + out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); + EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); + // = 1: same + il = 1; mi = 1; jl = 1; mj = 1; + out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); + EXPECT_NEAR(out.real(), 1.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); + // = 0: orthogonal + il = 2; mi = -1; jl = 2; mj = 0; + out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); + EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); + // = 1: same + il = 2; mi = 1; jl = 2; mj = 1; + out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); + EXPECT_NEAR(out.real(), 1.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); +} + +TEST_F(CalpLpRTest, CalLxijRTest) +{ + std::complex out; + + ModuleBase::Vector3 vR(0., 0., 0.); // home-cell + int it = 0, ia = 0, il = 0, iz = 0, im = 0; + int jt = 0, ja = 0, jl = 0, jz = 0, jm = 0; // self, the first s + + // = 0: no anisotropy + out = ModuleIO::cal_LxijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); + EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); + // = 0: orthogonal + jl = 1; + out = ModuleIO::cal_LxijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); + EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); + // = 0.5: Lx|p(m=0)> = 1/2 (|p(m=-1)> + |p(m=1)>) + il = 1; im = -1; jl = 1; jm = 0; + out = ModuleIO::cal_LxijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); + EXPECT_NEAR(out.real(), 0.5*sqrt(2.0), DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); + // = 0: expectation value is 0 + il = 1; im = 1; jl = 1; jm = 1; + out = ModuleIO::cal_LxijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); + EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); + // = 0.5: Lx|d(m=0)> = 1/2 (|d(m=-1)> + |d(m=1)>) + il = 2; im = -1; jl = 2; jm = 0; + out = ModuleIO::cal_LxijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); + EXPECT_NEAR(out.real(), 0.5*sqrt(6.0), DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); +} + +TEST_F(CalpLpRTest, CalLyijRTest) +{ + std::complex out; + + ModuleBase::Vector3 vR(0., 0., 0.); // home-cell + int it = 0, ia = 0, il = 0, iz = 0, im = 0; + int jt = 0, ja = 0, jl = 0, jz = 0, jm = 0; // self, the first s + + // = 0: no anisotropy + out = ModuleIO::cal_LyijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); + EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); + // = 0: orthogonal + jl = 1; + out = ModuleIO::cal_LyijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); + EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); + // = -i/2: Ly|p(m=0)> = -i/2 (|p(m=1)> - |p(m=-1)>) + il = 1; im = -1; jl = 1; jm = 0; + out = ModuleIO::cal_LyijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); + EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 0.5*sqrt(2.0), DOUBLETHRESHOLD); + // = 0: expectation value is 0 + il = 1; im = 1; jl = 1; jm = 1; + out = ModuleIO::cal_LyijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); + EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); + // = -i/2: Ly|d(m=0)> = -i/2 (|d(m=1)> - |d(m=-1)>) + il = 2; im = -1; jl = 2; jm = 0; + out = ModuleIO::cal_LyijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); + EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 0.5*sqrt(6.0), DOUBLETHRESHOLD); +} + +int main(int argc, char **argv) +{ +#ifdef __MPI + MPI_Init(&argc, &argv); +#endif + ::testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); + +#ifdef __MPI + MPI_Finalize(); +#endif + return 0; +} \ No newline at end of file From 43d1c0e12a6feebb57bfa78ade935251820ed6bd Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Tue, 1 Apr 2025 19:58:34 +0800 Subject: [PATCH 03/10] nearly complete --- source/module_esolver/esolver_ks_lcao.cpp | 16 ++++++ source/module_io/cal_pLpR.cpp | 63 ++++++++++++++++++++- source/module_io/cal_pLpR.h | 5 +- source/module_io/read_input_item_output.cpp | 15 +++++ source/module_parameter/input_parameter.h | 1 + 5 files changed, 97 insertions(+), 3 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index b088ebab02f..afe4fddedb3 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -28,6 +28,7 @@ #include "module_io/write_proj_band_lcao.h" #include "module_io/write_vxc.hpp" #include "module_io/write_wfc_nao.h" +#include "module_io/cal_pLpR.h" #include "module_parameter/parameter.h" //--------------temporary---------------------------- @@ -1242,6 +1243,21 @@ void ESolver_KS_LCAO::after_scf(const int istep) delete ekinetic; } + if (PARAM.inp.out_mat_l[0]) + { + ModuleIO::AngularMomentumExpectationCalculator mylcalculator( + PARAM.inp.orbital_dir, + GlobalC::ucell, + PARAM.inp.search_radius, + PARAM.inp.test_deconstructor, + PARAM.inp.test_grid, + PARAM.inp.test_atom_input, + PARAM.inp.search_pbc, + &GlobalV::ofs_running + ); + + } + // add by jingan in 2018.11.7 if (PARAM.inp.calculation == "nscf" && PARAM.inp.towannier90) { diff --git a/source/module_io/cal_pLpR.cpp b/source/module_io/cal_pLpR.cpp index 8423078af22..2a952c5abfb 100644 --- a/source/module_io/cal_pLpR.cpp +++ b/source/module_io/cal_pLpR.cpp @@ -12,7 +12,7 @@ #include "module_cell/module_neighbor/sltk_atom_arrange.h" #include "module_parameter/parameter.h" #include "module_io/cal_pLpR.h" - +#include "module_base/formatter.h" /** * * FIXME: the following part will be transfered to TwoCenterIntegrator soon @@ -137,12 +137,28 @@ ModuleIO::AngularMomentumExpectationCalculator::AngularMomentumExpectationCalcul void ModuleIO::AngularMomentumExpectationCalculator::calculate( std::ofstream* ofs, const UnitCell& ucell, - const char dir) + const char dir, + const int precision) { if (ofs == nullptr) { return; } + // an easy sanity check + assert(dir == 'x' || dir == 'y' || dir == 'z'); + + // it, ia, il, iz, im, iRx, iRy, iRz, jt, ja, jl, jz, jm + // the iRx, iRy, iRz are the indices of the supercell in which the two-center-integral + // it and jt are indexes of atomtypes, + // ia and ja are indexes of atoms within the atomtypes, + // il and jl are indexes of the angular momentum, + // iz and jz are indexes of the zeta functions + // im and jm are indexes of the magnetic quantum numbers. + std::string fmtstr = "%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d"; + fmtstr += "%" + std::to_string(precision*2) + "." + std::to_string(precision) + "e"; + fmtstr += "%" + std::to_string(precision*2) + "." + std::to_string(precision) + "e\n"; + FmtCore fmt(fmtstr); + ModuleBase::Vector3 ri, rj, dr; for (int it = 0; it < ucell.ntype; it++) { @@ -155,8 +171,51 @@ void ModuleIO::AngularMomentumExpectationCalculator::calculate( { rj = neighbor_searcher_->getAdjacentTau(ia_adj); int jt = neighbor_searcher_->getType(ia_adj); + const Atom& atyp_j = ucell.atoms[jt]; + int ja = neighbor_searcher_->getNatom(ia_adj); dr = (ri - rj) * ucell.lat0; const ModuleBase::Vector3 iR = neighbor_searcher_->getBox(ia_adj); + // the two-center-integral + + for (int li = 0; li < atyp_i.nwl + 1; li++) + { + for (int iz = 0; iz < atyp_i.l_nchi[li]; iz++) + { + for (int mi = -li; mi <= li; mi++) + { + for (int lj = 0; lj < atyp_j.nwl + 1; lj++) + { + for (int jz = 0; jz < atyp_j.l_nchi[lj]; jz++) + { + for (int mj = -lj; mj <= lj; mj++) + { + std::complex val = 0; + if (dir == 'x') + { + val = cal_LxijR(calculator_, + it, ia, li, iz, mi, jt, ja, lj, jz, mj, dr); + } + else if (dir == 'y') + { + val = cal_LyijR(calculator_, + it, ia, li, iz, mi, jt, ja, lj, jz, mj, dr); + } + else if (dir == 'z') + { + val = cal_LzijR(calculator_, + it, ia, li, iz, mi, jt, ja, lj, jz, mj, dr); + } + *ofs_ << fmt.format( + it, ia, li, iz, mi, + iR.x, iR.y, iR.z, + jt, ja, lj, jz, mj, + val.real(), val.imag()); + } + } + } + } + } + } } } } diff --git a/source/module_io/cal_pLpR.h b/source/module_io/cal_pLpR.h index 654f016a076..005671c1ae0 100644 --- a/source/module_io/cal_pLpR.h +++ b/source/module_io/cal_pLpR.h @@ -66,6 +66,8 @@ #include #include "module_cell/unitcell.h" #include "module_basis/module_nao/two_center_integrator.h" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_cell/module_neighbor/sltk_atom_arrange.h" namespace ModuleIO { @@ -115,7 +117,8 @@ namespace ModuleIO void calculate(std::ofstream* ofs, const UnitCell& ucell, - const char dir = 'x'); + const char dir = 'x', + const int precision = 10); private: // ofsrunning diff --git a/source/module_io/read_input_item_output.cpp b/source/module_io/read_input_item_output.cpp index 03742d24d04..94b07e84be6 100644 --- a/source/module_io/read_input_item_output.cpp +++ b/source/module_io/read_input_item_output.cpp @@ -275,6 +275,21 @@ void ReadInput::item_output() }; this->add_item(item); } + { + Input_Item item("out_mat_l"); + item.annotation = "output the expectation values of angular momentum operators"; + item.read_value = [](const Input_Item& item, Parameter& para) { + const size_t count = item.get_size(); + if (count != 1 && count != 2) + { + ModuleBase::WARNING_QUIT("ReadInput", "out_mat_l should have 1 or 2 values"); + } + para.input.out_mat_l[0] = assume_as_boolean(item.str_values[0]); + para.input.out_mat_l[1] = (count == 2) ? std::stoi(item.str_values[1]) : 8; + }; + sync_intvec(input.out_mat_l, 2, 0); + this->add_item(item); + } { Input_Item item("out_mat_dh"); item.annotation = "output of derivative of H(R) matrix"; diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index bfd9ae91cbc..936dc220e29 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -331,6 +331,7 @@ struct Input_para bool out_bandgap = false; ///< QO added for bandgap printing std::vector out_mat_hs = {0, 8}; ///< output H matrix and S matrix in local basis. std::vector out_mat_tk = {0, 8}; ///< output T(k) matrix in local basis. + std::vector out_mat_l = {0, 8}; ///< output L matrix in local basis. bool out_mat_hs2 = false; ///< LiuXh add 2019-07-16, output H(R) matrix and ///< S(R) matrix in local basis. bool out_mat_dh = false; From 4c6e7c8b2cad12e047fd2afab2728d189a4c2379 Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Wed, 2 Apr 2025 12:48:32 +0800 Subject: [PATCH 04/10] complete --- source/module_esolver/esolver_ks_lcao.cpp | 9 +- source/module_io/cal_pLpR.cpp | 113 +++++++++++++++++----- source/module_io/cal_pLpR.h | 100 +++++++++++++++++-- source/module_io/test/CMakeLists.txt | 2 +- 4 files changed, 193 insertions(+), 31 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index afe4fddedb3..0e26a770d2b 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -1253,9 +1253,14 @@ void ESolver_KS_LCAO::after_scf(const int istep) PARAM.inp.test_grid, PARAM.inp.test_atom_input, PARAM.inp.search_pbc, - &GlobalV::ofs_running + &GlobalV::ofs_running, + GlobalV::MY_RANK ); - + mylcalculator.calculate(PARAM.inp.suffix, + PARAM.inp.read_file_dir, + GlobalC::ucell, + PARAM.inp.out_mat_l[1], + GlobalV::MY_RANK); } // add by jingan in 2018.11.7 diff --git a/source/module_io/cal_pLpR.cpp b/source/module_io/cal_pLpR.cpp index 2a952c5abfb..01b15499ba3 100644 --- a/source/module_io/cal_pLpR.cpp +++ b/source/module_io/cal_pLpR.cpp @@ -13,6 +13,7 @@ #include "module_parameter/parameter.h" #include "module_io/cal_pLpR.h" #include "module_base/formatter.h" +#include "module_base/parallel_common.h" /** * * FIXME: the following part will be transfered to TwoCenterIntegrator soon @@ -96,51 +97,81 @@ ModuleIO::AngularMomentumExpectationCalculator::AngularMomentumExpectationCalcul const int tgrid, const int tatom, const bool searchpbc, - std::ofstream* ptr_log) + std::ofstream* ptr_log, + const int rank) { - std::vector forb(ucell.ntype); - for (int i = 0; i < ucell.ntype; ++i) + + // ofs_running + this->ofs_ = ptr_log; + *ofs_ << "\n\n\n\n"; + *ofs_ << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; + *ofs_ << " | |" << std::endl; + *ofs_ << " | Angular momentum expectation value calculation: |" << std::endl; + *ofs_ << " | This is a post-processing step. The expectation value of operator |" << std::endl; + *ofs_ << " | Lx, Ly, Lz (, in which a and b are ABACUS numerical atomic |" << std::endl; + *ofs_ << " | orbitals) will be calculated. |" << std::endl; + *ofs_ << " | The result will be printed to file with name ${suffix}_Lx/y/z.dat |" << std::endl; + *ofs_ << " | |" << std::endl; + *ofs_ << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; + *ofs_ << "\n\n\n\n"; + + int ntype_ = ucell.ntype; +#ifdef __MPI + Parallel_Common::bcast_int(ntype_); +#endif + std::vector forb(ntype_); + if (rank == 0) { - forb[i] = orbital_dir + ucell.orbital_fn[i]; + for (int i = 0; i < ucell.ntype; ++i) + { + forb[i] = orbital_dir + ucell.orbital_fn[i]; + } } +#ifdef __MPI + Parallel_Common::bcast_string(forb.data(), ntype_); +#endif + this->orb_ = std::unique_ptr(new RadialCollection); this->orb_->build(ucell.ntype, forb.data(), 'o'); - + ModuleBase::SphericalBesselTransformer sbt(true); this->orb_->set_transformer(sbt); - + const double rcut_max = orb_->rcut_max(); const int ngrid = int(rcut_max / 0.01) + 1; const double cutoff = 2.0 * rcut_max; this->orb_->set_uniform_grid(true, ngrid, cutoff, 'i', true); - + this->calculator_ = std::unique_ptr(new TwoCenterIntegrator); this->calculator_->tabulate(*orb_, *orb_, 'S', ngrid, cutoff); - + // Initialize Ylm coefficients ModuleBase::Ylm::set_coefficients(); - - // ofs_running - this->ofs_ = ptr_log; - + // for neighbor list search + double temp = -1.0; + temp = atom_arrange::set_sr_NL(*ofs_, + PARAM.inp.out_level, + search_radius, + ucell.infoNL.get_rcutmax_Beta(), + PARAM.globalv.gamma_only_local); + temp = std::max(temp, search_radius); this->neighbor_searcher_ = std::unique_ptr(new Grid_Driver(tdestructor, tgrid)); - atom_arrange::search( - searchpbc, - *ofs_, - *neighbor_searcher_, - ucell, - search_radius, - tatom); + atom_arrange::search(searchpbc, + *ofs_, + *neighbor_searcher_, + ucell, + temp, + tatom); } -void ModuleIO::AngularMomentumExpectationCalculator::calculate( +void ModuleIO::AngularMomentumExpectationCalculator::kernel( std::ofstream* ofs, const UnitCell& ucell, const char dir, const int precision) { - if (ofs == nullptr) + if (!ofs->is_open()) { return; } @@ -205,7 +236,8 @@ void ModuleIO::AngularMomentumExpectationCalculator::calculate( val = cal_LzijR(calculator_, it, ia, li, iz, mi, jt, ja, lj, jz, mj, dr); } - *ofs_ << fmt.format( + + *ofs << fmt.format( it, ia, li, iz, mi, iR.x, iR.y, iR.z, jt, ja, lj, jz, mj, @@ -219,4 +251,41 @@ void ModuleIO::AngularMomentumExpectationCalculator::calculate( } } } +} + +void ModuleIO::AngularMomentumExpectationCalculator::calculate( + const std::string& prefix, + const std::string& outdir, + const UnitCell& ucell, + const int precision, + const int rank) +{ + if (rank != 0) + { + return; + } + std::ofstream ofout; + const std::string dir = "xyz"; + const std::string title = "# it ia il iz im iRx iRy iRz jt ja jl jz jm \n" + "# it: atomtype index of the first atom\n" + "# ia: atomic index of the first atom within the atomtype\n" + "# il: angular momentum index of the first atom\n" + "# iz: zeta function index of the first atom\n" + "# im: magnetic quantum number of the first atom\n" + "# iRx, iRy, iRz: the indices of the supercell\n" + "# jt: atomtype index of the second atom\n" + "# ja: atomic index of the second atom within the atomtype\n" + "# jl: angular momentum index of the second atom\n" + "# jz: zeta function index of the second atom\n" + "# jm: magnetic quantum number of the second atom\n" + "# : the value of the matrix element\n"; + + for (char d : dir) + { + std::string fn = outdir + prefix + "_L" + d + ".dat"; + ofout.open(fn, std::ios::out); + ofout << title; + this->kernel(&ofout, ucell, d, precision); + ofout.close(); + } } \ No newline at end of file diff --git a/source/module_io/cal_pLpR.h b/source/module_io/cal_pLpR.h index 005671c1ae0..e55e328c061 100644 --- a/source/module_io/cal_pLpR.h +++ b/source/module_io/cal_pLpR.h @@ -71,18 +71,75 @@ namespace ModuleIO { + /** + * @brief calculate the matrix elements, in which the Lz + * are the angular momentum operators, |phi_i> and |phi_j> are the numerical + * atomic orbitals (NAOs). + * + * @param calculator the std::unique_ptr instance + * @param it atomtype index of the first atom + * @param ia atomic index of the first atom within the atomtype + * @param il angular momentum index of the first atom + * @param iz zeta function index of the first atom + * @param mi magnetic quantum number of the first atom + * @param jt atomtype index of the second atom + * @param ja atomic index of the second atom within the atomtype + * @param jl angular momentum index of the second atom + * @param jz zeta function index of the second atom + * @param mj magnetic quantum number of the second atom + * @param vR the vector from the first atom to the second atom + * @return std::complex + */ std::complex cal_LzijR( const std::unique_ptr& calculator, const int it, const int ia, const int il, const int iz, const int mi, const int jt, const int ja, const int jl, const int jz, const int mj, const ModuleBase::Vector3& vR); + /** + * @brief calculate the matrix elements, in which the Lz + * are the angular momentum operators, |phi_i> and |phi_j> are the numerical + * atomic orbitals (NAOs). + * + * @param calculator the std::unique_ptr instance + * @param it atomtype index of the first atom + * @param ia atomic index of the first atom within the atomtype + * @param il angular momentum index of the first atom + * @param iz zeta function index of the first atom + * @param mi magnetic quantum number of the first atom + * @param jt atomtype index of the second atom + * @param ja atomic index of the second atom within the atomtype + * @param jl angular momentum index of the second atom + * @param jz zeta function index of the second atom + * @param mj magnetic quantum number of the second atom + * @param vR the vector from the first atom to the second atom + * @return std::complex + */ std::complex cal_LyijR( const std::unique_ptr& calculator, const int it, const int ia, const int il, const int iz, const int im, const int jt, const int ja, const int jl, const int jz, const int jm, const ModuleBase::Vector3& vR); + /** + * @brief calculate the matrix elements, in which the Lz + * are the angular momentum operators, |phi_i> and |phi_j> are the numerical + * atomic orbitals (NAOs). + * + * @param calculator the std::unique_ptr instance + * @param it atomtype index of the first atom + * @param ia atomic index of the first atom within the atomtype + * @param il angular momentum index of the first atom + * @param iz zeta function index of the first atom + * @param mi magnetic quantum number of the first atom + * @param jt atomtype index of the second atom + * @param ja atomic index of the second atom within the atomtype + * @param jl angular momentum index of the second atom + * @param jz zeta function index of the second atom + * @param mj magnetic quantum number of the second atom + * @param vR the vector from the first atom to the second atom + * @return std::complex + */ std::complex cal_LxijR( const std::unique_ptr& calculator, const int it, const int ia, const int il, const int iz, const int im, @@ -103,7 +160,20 @@ namespace ModuleIO class AngularMomentumExpectationCalculator { public: - AngularMomentumExpectationCalculator() = delete; // the default constructor is meaningless + // the default constructor is meaningless + AngularMomentumExpectationCalculator() = delete; + /** + * @brief Construct a new Angular Momentum Expectation Calculator object + * + * @param orbital_dir the directory of the orbital file + * @param ucell the unit cell object + * @param search_radius the search radius for the neighboring atoms + * @param tdestructor test flag, for destructor + * @param tgrid test flag, for grid + * @param tatom test flag, for atom input + * @param searchpbc + * @param ptr_log pointer to the ofstream object for logging + */ AngularMomentumExpectationCalculator( const std::string& orbital_dir, const UnitCell& ucell, @@ -112,13 +182,15 @@ namespace ModuleIO const int tgrid, const int tatom, const bool searchpbc, - std::ofstream* ptr_log = nullptr); + std::ofstream* ptr_log = nullptr, + const int rank = 0); ~AngularMomentumExpectationCalculator() = default; - void calculate(std::ofstream* ofs, - const UnitCell& ucell, - const char dir = 'x', - const int precision = 10); + void calculate(const std::string& prefix, + const std::string& outdir, + const UnitCell& ucell, + const int precision = 10, + const int rank = 0); private: // ofsrunning @@ -132,5 +204,21 @@ namespace ModuleIO // neighboring searcher std::unique_ptr neighbor_searcher_; + + /** + * @brief calculate the matrix elements. Due to + * the large size of the matrix, the result will be printed to file + * directly. + * + * @param ofs pointer to the ofstream object for printing, if nullptr, + * the result will not be printed + * @param ucell the unit cell object + * @param dir the direction of the angular momentum operator, 'x', 'y' or 'z' + * @param precision the precision of the output, default is 10 + */ + void kernel(std::ofstream* ofs, + const UnitCell& ucell, + const char dir = 'x', + const int precision = 10); }; } // namespace ModuleIO \ No newline at end of file diff --git a/source/module_io/test/CMakeLists.txt b/source/module_io/test/CMakeLists.txt index 19068397e07..7b7a7916c39 100644 --- a/source/module_io/test/CMakeLists.txt +++ b/source/module_io/test/CMakeLists.txt @@ -242,7 +242,7 @@ add_test(NAME orb_io_test_parallel if(ENABLE_LCAO) AddTest( TARGET cal_pLpR_test - LIBS parameter base ${math_libs} device + LIBS parameter base ${math_libs} device neighbor SOURCES cal_pLpR_test.cpp ../cal_pLpR.cpp From daad6e6081b28419a4b93eaecd1c2ac9087e9fe4 Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Wed, 2 Apr 2025 13:03:56 +0800 Subject: [PATCH 05/10] update the doc --- docs/advanced/input_files/input-main.md | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index f8eec32d73d..9dae471ef06 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -154,6 +154,7 @@ - [out\_mat\_dh](#out_mat_dh) - [out\_mat\_xc](#out_mat_xc) - [out\_mat\_xc2](#out_mat_xc2) + - [out\_mat\_l](#out_mat_l) - [out\_eband\_terms](#out_eband_terms) - [out\_hr\_npz/out\_dm\_npz](#out_hr_npzout_dm_npz) - [dm\_to\_rho](#dm_to_rho) @@ -1807,6 +1808,13 @@ The band (KS orbital) energy for each (k-point, spin, band) will be printed in t - **Description**: Whether to print the exchange-correlation matrices in **numerical orbital representation** (unit: Ry): $\braket{\phi_i|V_\text{xc}^\text{(semi-)local}+V_\text{exx}+V_\text{DFTU}|\phi_j}(\mathbf{R})$ in CSR format (the same format as [out_mat_hs2](../elec_properties/hs_matrix.md#out_mat_hs2)) in the directory `OUT.${suffix}`. (Note that currently DeePKS term is not included. ) The files are named `Vxc_R_spin$s`. - **Default**: False +### out_mat_l + +- **Type**: Boolean [Integer\](optional) +- **Availability**: Numerical atomic orbital (NAO) basis +- **Description**: Whether to print the expectation value of the angular momentum operator $\hat{L}_x$, $\hat{L}_y$, and $\hat{L}_z$ in the basis of the localized atomic orbitals. The files are named `OUT.${suffix}/${suffix}_Lx.dat`, `OUT.${suffix}/${suffix}_Ly.dat`, and `OUT.${suffix}/${suffix}_Lz.dat`. The second integer controls the precision of the output. +- **Default**: False 8 + ### out_eband_terms - **Type**: Boolean From f4086519b54ef21f1d0daf9a0e2166fc8d7c825c Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Wed, 2 Apr 2025 13:14:46 +0800 Subject: [PATCH 06/10] update makefile --- source/Makefile.Objects | 1 + 1 file changed, 1 insertion(+) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 9b690c5fdf8..9348b870d2c 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -554,6 +554,7 @@ OBJS_IO=input_conv.o\ read_input_item_output.o\ read_set_globalv.o\ orb_io.o\ + cal_pLpR.o\ OBJS_IO_LCAO=cal_r_overlap_R.o\ write_orb_info.o\ From e7a0bcbb8d3bc91a86c923c908368f308945b5d6 Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Wed, 2 Apr 2025 13:39:04 +0800 Subject: [PATCH 07/10] move functions, correct CMakeLists --- source/module_esolver/esolver_ks_lcao.cpp | 20 -------------------- source/module_esolver/lcao_after_scf.cpp | 21 +++++++++++++++++++++ source/module_io/test/CMakeLists.txt | 1 - 3 files changed, 21 insertions(+), 21 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index a75d71be5f9..e4d8e56082b 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -563,26 +563,6 @@ void ESolver_KS_LCAO::after_all_runners(UnitCell& ucell) #endif ); } - - if (PARAM.inp.out_mat_l[0]) - { - ModuleIO::AngularMomentumExpectationCalculator mylcalculator( - PARAM.inp.orbital_dir, - GlobalC::ucell, - PARAM.inp.search_radius, - PARAM.inp.test_deconstructor, - PARAM.inp.test_grid, - PARAM.inp.test_atom_input, - PARAM.inp.search_pbc, - &GlobalV::ofs_running, - GlobalV::MY_RANK - ); - mylcalculator.calculate(PARAM.inp.suffix, - PARAM.inp.read_file_dir, - GlobalC::ucell, - PARAM.inp.out_mat_l[1], - GlobalV::MY_RANK); - } ModuleBase::timer::tick("ESolver_KS_LCAO", "after_all_runners"); } diff --git a/source/module_esolver/lcao_after_scf.cpp b/source/module_esolver/lcao_after_scf.cpp index 73fb309a059..c93346af4b9 100644 --- a/source/module_esolver/lcao_after_scf.cpp +++ b/source/module_esolver/lcao_after_scf.cpp @@ -25,6 +25,7 @@ #include "module_io/write_istate_info.h" #include "module_io/write_proj_band_lcao.h" #include "module_io/write_wfc_nao.h" +#include "module_io/cal_pLpR.h" #include "module_parameter/parameter.h" //--------------temporary---------------------------- @@ -481,6 +482,26 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep, const RA.delete_grid(); } + if (PARAM.inp.out_mat_l[0]) + { + ModuleIO::AngularMomentumExpectationCalculator mylcalculator( + PARAM.inp.orbital_dir, + ucell, + PARAM.inp.search_radius, + PARAM.inp.test_deconstructor, + PARAM.inp.test_grid, + PARAM.inp.test_atom_input, + PARAM.inp.search_pbc, + &GlobalV::ofs_running, + GlobalV::MY_RANK + ); + mylcalculator.calculate(PARAM.inp.suffix, + PARAM.inp.read_file_dir, + ucell, + PARAM.inp.out_mat_l[1], + GlobalV::MY_RANK); + } + ModuleBase::timer::tick("ESolver_KS_LCAO", "after_scf"); } diff --git a/source/module_io/test/CMakeLists.txt b/source/module_io/test/CMakeLists.txt index d322e5cd4ac..db7b727f384 100644 --- a/source/module_io/test/CMakeLists.txt +++ b/source/module_io/test/CMakeLists.txt @@ -278,7 +278,6 @@ AddTest( ../../module_basis/module_nao/pswfc_radials.cpp ../../module_basis/module_nao/atomic_radials.cpp ../../module_basis/module_nao/sphbes_radials.cpp - ../../module_basis/module_nao/projgen.cpp ../../module_basis/module_nao/two_center_integrator.cpp ../../module_basis/module_nao/two_center_table.cpp ../../module_basis/module_nao/real_gaunt_table.cpp From 806f97bb5e9effbd99f92fbea59cfd9106948e47 Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Wed, 2 Apr 2025 15:17:18 +0800 Subject: [PATCH 08/10] update annotation and shorten the class name --- source/module_esolver/lcao_after_scf.cpp | 33 +++++++++++++----------- source/module_io/cal_pLpR.cpp | 6 ++--- source/module_io/cal_pLpR.h | 8 +++--- 3 files changed, 25 insertions(+), 22 deletions(-) diff --git a/source/module_esolver/lcao_after_scf.cpp b/source/module_esolver/lcao_after_scf.cpp index c93346af4b9..7e8aa9093c3 100644 --- a/source/module_esolver/lcao_after_scf.cpp +++ b/source/module_esolver/lcao_after_scf.cpp @@ -482,24 +482,27 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep, const RA.delete_grid(); } + //------------------------------------------------------------------ + //! 20) calculate expectation of angular momentum operator in LCAO basis + //------------------------------------------------------------------ if (PARAM.inp.out_mat_l[0]) { - ModuleIO::AngularMomentumExpectationCalculator mylcalculator( - PARAM.inp.orbital_dir, - ucell, - PARAM.inp.search_radius, - PARAM.inp.test_deconstructor, - PARAM.inp.test_grid, - PARAM.inp.test_atom_input, - PARAM.inp.search_pbc, - &GlobalV::ofs_running, - GlobalV::MY_RANK + ModuleIO::AngularMomentumCalculator mylcalculator( + /*orbital_dir=*/PARAM.inp.orbital_dir, + /*ucell=*/ucell, + /*search_radius=*/PARAM.inp.search_radius, + /*test_deconstructor=*/PARAM.inp.test_deconstructor, + /*test_grid=*/PARAM.inp.test_grid, + /*test_atom_input=*/PARAM.inp.test_atom_input, + /*search_pbc=*/PARAM.inp.search_pbc, + /*ofs=*/&GlobalV::ofs_running, + /*rank=*/GlobalV::MY_RANK ); - mylcalculator.calculate(PARAM.inp.suffix, - PARAM.inp.read_file_dir, - ucell, - PARAM.inp.out_mat_l[1], - GlobalV::MY_RANK); + mylcalculator.calculate(/*suffix=*/PARAM.inp.suffix, + /*outdir=*/PARAM.globalv.global_out_dir, + /*ucell=*/ucell, + /*precision=*/PARAM.inp.out_mat_l[1], + /*rank=*/GlobalV::MY_RANK); } ModuleBase::timer::tick("ESolver_KS_LCAO", "after_scf"); diff --git a/source/module_io/cal_pLpR.cpp b/source/module_io/cal_pLpR.cpp index 01b15499ba3..095e6f48a7c 100644 --- a/source/module_io/cal_pLpR.cpp +++ b/source/module_io/cal_pLpR.cpp @@ -89,7 +89,7 @@ std::complex ModuleIO::cal_LxijR( return std::complex(0.5) * (val_plus + val_minus); } -ModuleIO::AngularMomentumExpectationCalculator::AngularMomentumExpectationCalculator( +ModuleIO::AngularMomentumCalculator::AngularMomentumCalculator( const std::string& orbital_dir, const UnitCell& ucell, const double& search_radius, @@ -165,7 +165,7 @@ ModuleIO::AngularMomentumExpectationCalculator::AngularMomentumExpectationCalcul tatom); } -void ModuleIO::AngularMomentumExpectationCalculator::kernel( +void ModuleIO::AngularMomentumCalculator::kernel( std::ofstream* ofs, const UnitCell& ucell, const char dir, @@ -253,7 +253,7 @@ void ModuleIO::AngularMomentumExpectationCalculator::kernel( } } -void ModuleIO::AngularMomentumExpectationCalculator::calculate( +void ModuleIO::AngularMomentumCalculator::calculate( const std::string& prefix, const std::string& outdir, const UnitCell& ucell, diff --git a/source/module_io/cal_pLpR.h b/source/module_io/cal_pLpR.h index e55e328c061..a5368d885b3 100644 --- a/source/module_io/cal_pLpR.h +++ b/source/module_io/cal_pLpR.h @@ -157,11 +157,11 @@ namespace ModuleIO // im and jm are indexes of the magnetic quantum numbers. // The output is a complex number, which is the value of the matrix element. // Always the matrix is quite large, so direct print to file. - class AngularMomentumExpectationCalculator + class AngularMomentumCalculator { public: // the default constructor is meaningless - AngularMomentumExpectationCalculator() = delete; + AngularMomentumCalculator() = delete; /** * @brief Construct a new Angular Momentum Expectation Calculator object * @@ -174,7 +174,7 @@ namespace ModuleIO * @param searchpbc * @param ptr_log pointer to the ofstream object for logging */ - AngularMomentumExpectationCalculator( + AngularMomentumCalculator( const std::string& orbital_dir, const UnitCell& ucell, const double& search_radius, @@ -184,7 +184,7 @@ namespace ModuleIO const bool searchpbc, std::ofstream* ptr_log = nullptr, const int rank = 0); - ~AngularMomentumExpectationCalculator() = default; + ~AngularMomentumCalculator() = default; void calculate(const std::string& prefix, const std::string& outdir, From 119c54e6e4b89962c989475bfff2be321ab03909 Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Wed, 2 Apr 2025 16:57:45 +0800 Subject: [PATCH 09/10] delete the inclusion phrase in esolver_ks_lcao.cpp --- source/module_esolver/esolver_ks_lcao.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index e4d8e56082b..cacede3babe 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -28,7 +28,6 @@ #include "module_io/write_istate_info.h" #include "module_io/write_proj_band_lcao.h" #include "module_io/write_wfc_nao.h" -#include "module_io/cal_pLpR.h" #include "module_parameter/parameter.h" // be careful of hpp, there may be multiple definitions of functions, 20250302, mohan From 1bdf6e754bd0cfe0ff2f09f8ee2774612d754b86 Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Wed, 2 Apr 2025 16:59:04 +0800 Subject: [PATCH 10/10] further recover the esolver_ks_lcao.cpp --- source/module_esolver/esolver_ks_lcao.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index cacede3babe..941e4938be4 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -562,7 +562,7 @@ void ESolver_KS_LCAO::after_all_runners(UnitCell& ucell) #endif ); } - + ModuleBase::timer::tick("ESolver_KS_LCAO", "after_all_runners"); }