tmp debug

Author: maki49

source/source_cell/module_neighbor/sltk_grid_driver.cpp

@@ -34,7 +34,7 @@ void Grid_Driver::Find_atom(const UnitCell& ucell,
     AdjacentAtomInfo* local_adjs = adjs == nullptr ? &this->adj_info : adjs;
     local_adjs->clear();
     const std::vector<FAtom*>& all_atom = all_adj_info[ntype][nnumber];
-
+    std::cout << "all_atom size: " << all_atom.size() << std::endl;
     for (const FAtom* atom: all_atom)
     {
         local_adjs->ntype.push_back(atom->type);

source/source_esolver/lcao_before_scf.cpp

@@ -50,7 +50,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
                                                    orb_.get_rcutmax_Phi(),
                                                    ucell.infoNL.get_rcutmax_Beta(),
                                                    PARAM.globalv.gamma_only_local);
-
+    std::cout << "search radius for neighbor list: " << search_radius << std::endl;
     //! 3) use search_radius to search adj atoms
     atom_arrange::search(PARAM.globalv.search_pbc,
                          GlobalV::ofs_running,

source/source_io/cal_pLpR.cpp

@@ -156,13 +156,15 @@ ModuleIO::AngularMomentumCalculator::AngularMomentumCalculator(
                                    ucell.infoNL.get_rcutmax_Beta(),
                                    PARAM.globalv.gamma_only_local);
     temp = std::max(temp, search_radius);
+    std::cout << "search radius for neighbor list: " << temp << std::endl;
     this->neighbor_searcher_ = std::unique_ptr<Grid_Driver>(new Grid_Driver(tdestructor, tgrid));
     atom_arrange::search(searchpbc,
                          *ofs_,
                          *neighbor_searcher_,
                          ucell,
                          temp,
                          tatom);
+    std::cout << "neighbored atom num:" << neighbor_searcher_->getAdjacentNum() << std::endl;
 }
 
 void ModuleIO::AngularMomentumCalculator::kernel(
@@ -171,10 +173,12 @@ void ModuleIO::AngularMomentumCalculator::kernel(
     const char dir,
     const int precision)
 {
+    std::cout << "begin kernel" << std::endl;
     if (!ofs->is_open())
     {
         return;
     }
+    std::cout << "file check ok" << std::endl;
     // an easy sanity check
     assert(dir == 'x' || dir == 'y' || dir == 'z');
 
@@ -197,8 +201,11 @@ void ModuleIO::AngularMomentumCalculator::kernel(
         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++)
+            AdjacentAtomInfo adjs;
+            neighbor_searcher_->Find_atom(ucell, ri, it, ia, &adjs);
+            std::cout << "adjacent num: " << adjs.adj_num << std::endl;
+            // bug: adj num is 0 here?!
+            for (int ia_adj = 0; ia_adj < adjs.adj_num; ia_adj++)
             {
                 rj = neighbor_searcher_->getAdjacentTau(ia_adj);
                 int jt = neighbor_searcher_->getType(ia_adj);
@@ -207,11 +214,14 @@ void ModuleIO::AngularMomentumCalculator::kernel(
                 dr = (ri - rj) * ucell.lat0;
                 const ModuleBase::Vector3<int> iR = neighbor_searcher_->getBox(ia_adj);
                 // the two-center-integral
-
+                std::cout << "ia=" << ia << " ja=" << ja << " R=(" << iR.x << "," << iR.y << "," << iR.z << ")" << std::endl;
+                std::cout << "nwl_i=" << atyp_i.nwl << " nw_i=" << atyp_i.nw << std::endl;
                 for (int li = 0; li < atyp_i.nwl + 1; li++)
                 {
+                    std::cout << "nchi=" << atyp_i.l_nchi[li] << std::endl;
                     for (int iz = 0; iz < atyp_i.l_nchi[li]; iz++)
                     {
+                        std::cout << "iz=" << iz << " li=" << li << std::endl;
                         for (int mi = -li; mi <= li; mi++)
                         {
                             for (int lj = 0; lj < atyp_j.nwl + 1; lj++)
@@ -236,7 +246,11 @@ void ModuleIO::AngularMomentumCalculator::kernel(
                                             val = cal_LzijR(calculator_, 
                                                 it, ia, li, iz, mi, jt, ja, lj, jz, mj, dr);
                                         }
-
+                                        std::cout << fmt.format(
+                                            it, ia, li, iz, mi,
+                                            iR.x, iR.y, iR.z,
+                                            jt, ja, lj, jz, mj,
+                                            val.real(), val.imag());
                                         *ofs << fmt.format(
                                             it, ia, li, iz, mi,
                                             iR.x, iR.y, iR.z,

source/source_io/ctrl_output_lcao.cpp

@@ -258,7 +258,7 @@ void ctrl_output_lcao(UnitCell& ucell,
         ModuleIO::AngularMomentumCalculator mylcalculator(
             inp.orbital_dir,
             ucell,
-            inp.search_radius,
+            std::max(inp.search_radius, orb.get_rcutmax_Phi()),
             inp.test_deconstructor,
             inp.test_grid,
             inp.test_atom_input,