rewrite module neighbour to parition atom into boxes

Author: goodchong

source/module_cell/module_neighbor/sltk_atom_arrange.cpp

@@ -70,7 +70,10 @@ void atom_arrange::search(const bool pbc_flag,
     ModuleBase::TITLE("atom_arrange", "search");
     ModuleBase::timer::tick("atom_arrange", "search");
 
-    assert(search_radius_bohr > 0.0);
+    if (search_radius_bohr < 0.0)
+    {
+        ModuleBase::WARNING_QUIT("atom_arrange::search", " search_radius_bohr < 0,forbidden");
+    }
 
     ModuleBase::GlobalFunc::OUT(ofs_in, "searching radius is (Bohr))", search_radius_bohr);
     ModuleBase::GlobalFunc::OUT(ofs_in, "searching radius unit is (Bohr))", ucell.lat0);
@@ -82,7 +85,7 @@ void atom_arrange::search(const bool pbc_flag,
 
     const double radius_lat0unit = search_radius_bohr / ucell.lat0;
 
-	
+
     Atom_input at(ofs_in, ucell, pbc_flag, radius_lat0unit, test_atom_in);
 
     grid_d.init(ofs_in, ucell, at);
@@ -138,5 +141,5 @@ void atom_arrange::delete_vector(std::ofstream& ofs_in,
 
     grid_d.delete_vector(at2.getGrid_layerX_minus(), at2.getGrid_layerY_minus(), at2.getGrid_layerZ_minus());
 
-    grid_d.delete_Cell();
+    grid_d.clear_atoms();
 }

source/module_cell/module_neighbor/sltk_atom_input.cpp

@@ -12,8 +12,10 @@ Atom_input::Atom_input(std::ofstream& ofs_in,
                        const bool boundary_in,
                        const double radius_in,
                        const int& test_atom_in)
-    : periodic_boundary(boundary_in), radius(radius_in), expand_flag(false), glayerX(1), glayerX_minus(0), glayerY(1),
-      glayerY_minus(0), glayerZ(1), glayerZ_minus(0), test_atom_input(test_atom_in)
+    : periodic_boundary(boundary_in), radius(radius_in),
+      glayerX(1), glayerX_minus(0), glayerY(1),
+      glayerY_minus(0), glayerZ(1), glayerZ_minus(0),
+      test_atom_input(test_atom_in)
 {
     ModuleBase::TITLE("Atom_input", "Atom_input");
 
@@ -23,70 +25,13 @@ Atom_input::Atom_input(std::ofstream& ofs_in,
         ModuleBase::GlobalFunc::OUT(ofs_in, "Searching radius(lat0)", radius);
     }
 
-    if (radius < 0.0)
-    {
-        ModuleBase::WARNING_QUIT("atom_arrange::init", " search radius < 0,forbidden");
-    }
-    // random selection, in order to estimate again.
-    this->x_min = ucell.atoms[0].tau[0].x;
-    this->y_min = ucell.atoms[0].tau[0].y;
-    this->z_min = ucell.atoms[0].tau[0].z;
-    this->x_max = ucell.atoms[0].tau[0].x;
-    this->y_max = ucell.atoms[0].tau[0].y;
-    this->z_max = ucell.atoms[0].tau[0].z;
-
-    // calculate min & max value
-    for (int i = 0; i < ucell.ntype; i++)
-    {
-        for (int j = 0; j < ucell.atoms[i].na; j++)
-        {
-            x_min = std::min(x_min, ucell.atoms[i].tau[j].x);
-            x_max = std::max(x_max, ucell.atoms[i].tau[j].x);
-            y_min = std::min(y_min, ucell.atoms[i].tau[j].y);
-            y_max = std::max(y_max, ucell.atoms[i].tau[j].y);
-            z_min = std::min(z_min, ucell.atoms[i].tau[j].z);
-            z_max = std::max(z_max, ucell.atoms[i].tau[j].z);
-        }
-    }
-
-    if (test_atom_input)
-    {
-        ModuleBase::GlobalFunc::OUT(ofs_in, "Find the coordinate range of the input atom(unit:lat0).");
-        ModuleBase::GlobalFunc::OUT(ofs_in, "min_tau", x_min, y_min, z_min);
-        ModuleBase::GlobalFunc::OUT(ofs_in, "max_tau", x_max, y_max, z_max);
-    }
-
-    //----------------------------------------------------------
-    // CALL MEMBER FUNCTION :
-    // NAME : Check_Expand_Condition(check if swe need to
-    // expand grid,and generate 6 MEMBER VARIABLE(number of
-    // layers for 6 dimension)
-    // initial value for "glayerX,Y,Z" : 1
-    // (if > 2 ,expand flag = 1)
-    // initial value for "glayerX,Y,Z_minus" : 0
-    // ( if > 1 ,expand flag = 1)
-    //----------------------------------------------------------
-
     this->Check_Expand_Condition(ucell);
 
     if (test_atom_input)
     {
         ModuleBase::GlobalFunc::OUT(ofs_in, "glayer+", glayerX, glayerY, glayerZ);
         ModuleBase::GlobalFunc::OUT(ofs_in, "glayer-", glayerX_minus, glayerY_minus, glayerZ_minus);
-        ModuleBase::GlobalFunc::OUT(ofs_in, "expand_flag", expand_flag);
-    }
-
-    //----------------------------------------------------------
-    // CALL MEMBER FUNCTION :
-    // NAME : calculate_cells
-    // Calculate how many cells we need in each direction.
-    //----------------------------------------------------------
-    this->calculate_cells();
-    if (test_atom_input)
-    {
-        ModuleBase::GlobalFunc::OUT(ofs_in, "CellDim", cell_nx, cell_ny, cell_nz);
     }
-    return;
 }
 
 Atom_input::~Atom_input()
@@ -180,143 +125,28 @@ void Atom_input::Check_Expand_Condition(const UnitCell& ucell)
     double a23_norm = sqrt(a23_1 * a23_1 + a23_2 * a23_2 + a23_3 * a23_3);
     double extend_v = a23_norm * radius;
     double extend_d1 = extend_v / ucell.omega * ucell.lat0 * ucell.lat0 * ucell.lat0;
-    int extend_d11 = static_cast<int>(extend_d1);
-    // 2016-09-05, LiuXh
-    if (extend_d1 - extend_d11 > 0.0)
-    {
-        extend_d11 += 1;
-    }
+    int extend_d11 = std::ceil(extend_d1);
 
     double a31_1 = ucell.latvec.e32 * ucell.latvec.e13 - ucell.latvec.e33 * ucell.latvec.e12;
     double a31_2 = ucell.latvec.e31 * ucell.latvec.e13 - ucell.latvec.e33 * ucell.latvec.e11;
     double a31_3 = ucell.latvec.e31 * ucell.latvec.e12 - ucell.latvec.e32 * ucell.latvec.e11;
     double a31_norm = sqrt(a31_1 * a31_1 + a31_2 * a31_2 + a31_3 * a31_3);
     double extend_d2 = a31_norm * radius / ucell.omega * ucell.lat0 * ucell.lat0 * ucell.lat0;
-    int extend_d22 = static_cast<int>(extend_d2);
-    // 2016-09-05, LiuXh
-    if (extend_d2 - extend_d22 > 0.0)
-    {
-        extend_d22 += 1;
-    }
+    int extend_d22 = std::ceil(extend_d2);
 
     double a12_1 = ucell.latvec.e12 * ucell.latvec.e23 - ucell.latvec.e13 * ucell.latvec.e22;
     double a12_2 = ucell.latvec.e11 * ucell.latvec.e23 - ucell.latvec.e13 * ucell.latvec.e21;
     double a12_3 = ucell.latvec.e11 * ucell.latvec.e22 - ucell.latvec.e12 * ucell.latvec.e21;
     double a12_norm = sqrt(a12_1 * a12_1 + a12_2 * a12_2 + a12_3 * a12_3);
     double extend_d3 = a12_norm * radius / ucell.omega * ucell.lat0 * ucell.lat0 * ucell.lat0;
-    int extend_d33 = static_cast<int>(extend_d3);
+    int extend_d33 = std::ceil(extend_d3);
     // 2016-09-05, LiuXh
-    if (extend_d3 - extend_d33 > 0.0)
-    {
-        extend_d33 += 1;
-    }
 
     glayerX = extend_d11 + 1;
     glayerY = extend_d22 + 1;
     glayerZ = extend_d33 + 1;
-    // Begin, 2016-09-05, LiuXh
-    // glayerX_minus = extend_d11 +1;
-    // glayerY_minus = extend_d22 +1;
-    // glayerZ_minus = extend_d33 +1;
     glayerX_minus = extend_d11;
     glayerY_minus = extend_d22;
     glayerZ_minus = extend_d33;
     // End, 2016-09-05, LiuXh
-
-    if (glayerX == 1)
-    {
-        glayerX++;
-    }
-    if (glayerY == 1)
-    {
-        glayerY++;
-    }
-    if (glayerZ == 1)
-    {
-        glayerZ++;
-    }
-    if (glayerX_minus == 1)
-    {
-        glayerX_minus++;
-    }
-    if (glayerY_minus == 1)
-    {
-        glayerY_minus++;
-    }
-    if (glayerZ_minus == 1)
-    {
-        glayerZ_minus++;
-    }
-    // End, 2016-07-19, LiuXh
-    /*
-        if(test_atom_input)
-        {
-            GlobalV::ofs_running << " Extend distance from the (minX,minY,minZ) direct position in this unitcell: " <<
-       std::endl;
-        }
-
-        if(test_atom_input)OUT(GlobalV::ofs_running,"ExtentDim-",extent_1DX_minus,extent_1DY_minus,extent_1DZ_minus);
-    */
-    //----------------------------------------------------------
-    // EXPLAIN : if extent don't satisfty the searching
-    // requiment, we must expand one more layer
-    //----------------------------------------------------------
-
-    if (glayerX > 2 || glayerY > 2 || glayerZ > 2)
-    {
-        this->expand_flag = true;
-    }
-    else if (glayerX_minus > 1 || glayerX_minus > 1 || glayerX_minus > 1)
-    {
-        this->expand_flag = true;
-    }
-    else
-    {
-        this->expand_flag = false;
-    }
-    return;
-}
-
-void Atom_input::calculate_cells()
-{
-    ModuleBase::TITLE("Atom_input", "calculate_cells");
-    //----------------------------------------------------------
-    // EXPLAIN :
-    // Expand_Case : Simple , we already know the cell numbres,
-    // all the trouble is only to construct adjacentset using all
-    // the cells.
-    // Not_Expand_Case : Using searching radius to construct
-    // the cells ,  trouble here,but is the convenience of searching
-    // time , we then only need to search 27-adjacent cell for each cell.
-    //----------------------------------------------------------
-    if (expand_flag)
-    {
-        cell_nx = glayerX + glayerX_minus;
-        cell_ny = glayerY + glayerY_minus;
-        cell_nz = glayerZ + glayerZ_minus;
-    }
-    else
-    {
-        // maybe a bug, if we don't use direct
-        // coordinates, mohan note 2011-04-14
-        double real_nx, real_ny, real_nz;
-        real_nx = (x_max - x_min) / radius;
-        real_ny = (y_max - y_min) / radius;
-        real_nz = (z_max - z_min) / radius;
-        cell_nx = static_cast<int>(real_nx) + 1;
-        cell_ny = static_cast<int>(real_ny) + 1;
-        cell_nz = static_cast<int>(real_nz) + 1;
-    }
-
-    //================
-    // Wrong !
-    //================
-    //	if(int_nx != real_nx) this->cell_nx++;
-    //	if(int_ny != real_ny) this->cell_ny++;
-    //	if(int_nz != real_nz) this->cell_nz++;
-    //=======================================
-    // Not need because if int_nx = real_nx,
-    // the position belong to the next cell
-    //=======================================
-    return;
 }

source/module_cell/module_neighbor/sltk_atom_input.h

@@ -19,11 +19,6 @@ class Atom_input
     ~Atom_input();
 
   public:
-    bool getExpandFlag() const
-    {
-        return expand_flag;
-    }
-
     int getBoundary() const
     {
         return periodic_boundary;
@@ -34,45 +29,6 @@ class Atom_input
         return radius;
     }
 
-    //==========================================================
-    //
-    //==========================================================
-    double minX() const
-    {
-		return x_min;
-    }
-
-    double minY() const
-    {
-		return y_min;
-    }
-
-    double minZ() const
-    {
-		return z_min;
-    }
-
-    //==========================================================
-    //
-    //==========================================================
-    int getCell_nX() const
-    {
-        return cell_nx;
-    }
-
-    int getCell_nY() const
-    {
-        return cell_ny;
-    }
-
-    int getCell_nZ() const
-    {
-        return cell_nz;
-    }
-
-    //==========================================================
-    //
-    //==========================================================
     int getGrid_layerX() const
     {
         return glayerX;
@@ -106,36 +62,15 @@ class Atom_input
   private:
     int test_atom_input; // caoyu reconst 2021-05-24
     bool periodic_boundary;
-
     double radius;
 
-    double x_min;
-    double y_min;
-    double z_min;
-    double x_max;
-    double y_max;
-    double z_max;
-    //==========================================================
-    // MEMBRE FUNCTION :
-    // NAME : Check_Expand_Condition
-    //==========================================================
     void Check_Expand_Condition(const UnitCell& ucell);
-    bool expand_flag;
     int glayerX;
     int glayerX_minus;
     int glayerY;
     int glayerY_minus;
     int glayerZ;
     int glayerZ_minus;
-
-    //==========================================================
-    // MEMBRE FUNCTION :
-    // NAME : Expand_Grid
-    //==========================================================
-    void calculate_cells();
-    int cell_nx;
-    int cell_ny;
-    int cell_nz;
 };
 
 #endif