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 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\ diff --git a/source/module_esolver/lcao_after_scf.cpp b/source/module_esolver/lcao_after_scf.cpp index 73fb309a059..7e8aa9093c3 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,29 @@ 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::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(/*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/CMakeLists.txt b/source/module_io/CMakeLists.txt index 677223d8444..49ccc527b4e 100644 --- a/source/module_io/CMakeLists.txt +++ b/source/module_io/CMakeLists.txt @@ -68,6 +68,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 new file mode 100644 index 00000000000..095e6f48a7c --- /dev/null +++ b/source/module_io/cal_pLpR.cpp @@ -0,0 +1,291 @@ +#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" +#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 + * + */ + +// 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)); +} + +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_); + return std::complex(mi) * 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); + 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_; + } + return std::complex(0, -0.5) * (val_plus - val_minus); +} + +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); + 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_; + } + return std::complex(0.5) * (val_plus + val_minus); +} + +ModuleIO::AngularMomentumCalculator::AngularMomentumCalculator( + 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, + const int rank) +{ + + // 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) + { + 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(); + + // 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, + temp, + tatom); +} + +void ModuleIO::AngularMomentumCalculator::kernel( + std::ofstream* ofs, + const UnitCell& ucell, + const char dir, + const int precision) +{ + if (!ofs->is_open()) + { + 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++) + { + 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); + 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()); + } + } + } + } + } + } + } + } + } +} + +void ModuleIO::AngularMomentumCalculator::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 new file mode 100644 index 00000000000..a5368d885b3 --- /dev/null +++ b/source/module_io/cal_pLpR.h @@ -0,0 +1,224 @@ +/** + * 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 + * ----------- + * + * 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 + * + * 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" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_cell/module_neighbor/sltk_atom_arrange.h" + +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, + 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 AngularMomentumCalculator + { + public: + // the default constructor is meaningless + AngularMomentumCalculator() = 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 + */ + AngularMomentumCalculator( + 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, + const int rank = 0); + ~AngularMomentumCalculator() = default; + + void calculate(const std::string& prefix, + const std::string& outdir, + const UnitCell& ucell, + const int precision = 10, + const int rank = 0); + + 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_; + + /** + * @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/read_input_item_output.cpp b/source/module_io/read_input_item_output.cpp index 0614ca1d08e..125ace481d9 100644 --- a/source/module_io/read_input_item_output.cpp +++ b/source/module_io/read_input_item_output.cpp @@ -280,6 +280,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_io/test/CMakeLists.txt b/source/module_io/test/CMakeLists.txt index e837d083cef..db7b727f384 100644 --- a/source/module_io/test/CMakeLists.txt +++ b/source/module_io/test/CMakeLists.txt @@ -260,4 +260,27 @@ add_test( NAME io_read_wfc_nao_test_parallel COMMAND mpirun -np 2 ./io_read_wfc_nao_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 neighbor + 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/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() 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 diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index 77ee42a2ed2..29ee8457946 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -371,6 +371,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;