Refactor atom_input and unify the fake atom data structure of grid and atom_input.

Author: goodchong

source/module_cell/module_neighbor/sltk_atom.h

@@ -20,6 +20,9 @@ class FAtom
 	int type;
 	int natom;
 
+	int cell_x;
+	int cell_y;
+	int cell_z;
 public:
 //==========================================================
 // Default Constructor and deconstructor
@@ -59,6 +62,10 @@ class FAtom
 	const double& z() const { return d_z; }
 	const int& getType() const { return type;}
 	const int& getNatom() const { return natom;}
+	const int& getCellX() const { return cell_x; }
+	const int& getCellY() const { return cell_y; }
+	const int& getCellZ() const { return cell_z; }
+
 
 //==========================================================
 // MEMBER FUNCTION :
@@ -69,6 +76,10 @@ class FAtom
 	void setZ(const double& r) { d_z = r; }
 	void setType(const int ntype) {type = ntype;}
 	void setNatom(const int atom) {natom = atom;}
+	void setCellX(const int nx) {cell_x = nx;}
+	void setCellY(const int ny) {cell_y = ny;}
+	void setCellZ(const int nz) {cell_z = nz;}
+	
 };
 
 #endif

source/module_cell/module_neighbor/sltk_atom_input.cpp

@@ -188,14 +188,6 @@ Atom_input::Atom_input
 
 Atom_input::~Atom_input()
 {
-	if (expand_flag)
-	{
-		delete [] store_x;
-		delete [] store_y;
-		delete [] store_z;
-		delete [] store_type;
-		delete [] store_natom;
-	}
 }
 
 //============================================
@@ -409,27 +401,6 @@ void Atom_input::Expand_Grid(const UnitCell &ucell, const int ntype)
 		GlobalV::ofs_running << " Good luck! " << std::endl;
 	}
 
-	double *x_old = new double[d_amount];
-	double *y_old = new double[d_amount];
-	double *z_old = new double[d_amount];
-
-	int *type_old = new int[d_amount];
-	int *natom_old = new int[d_amount];
-
-	int ia = 0;
-	for (int i = 0;i < ntype;i++)
-	{
-		for (int j = 0;j < ucell.atoms[i].na;j++)
-		{
-			x_old[ia] = ucell.atoms[i].tau[j].x;
-			y_old[ia] = ucell.atoms[i].tau[j].y;
-			z_old[ia] = ucell.atoms[i].tau[j].z;
-			type_old[ia] = i;
-			natom_old[ia] = j;
-			ia++;
-		}
-	}
-
 	// how many copys we need now.
 	const int gcopy =
 	    (glayerX + glayerX_minus) *
@@ -443,72 +414,37 @@ void Atom_input::Expand_Grid(const UnitCell &ucell, const int ntype)
 	if(test_atom_input)ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"Atom_number_now",d_amount_expand);
 
 	// store new atom positions.
-	this->store_x = new double[d_amount_expand];
-	this->store_y = new double[d_amount_expand];
-	this->store_z = new double[d_amount_expand];
-
-	ModuleBase::Memory::record("SLTK::Epd_Atom",sizeof(double) * d_amount_expand*3);
-
-	// store which grid the atom is in.
-	store_cell_x = new int[d_amount_expand];
-	store_cell_y = new int[d_amount_expand];
-	store_cell_z = new int[d_amount_expand];
-
-	this->store_type = new int[d_amount_expand];
-	this->store_natom = new int[d_amount_expand];
 
-	ModuleBase::Memory::record("SLTK::Epd_atom_info",sizeof(int) * d_amount_expand*5);
-
-	int ia_all = 0;
 
 	for (int ix = -glayerX_minus; ix < glayerX; ix++)
 	{
 		for (int iy = -glayerY_minus; iy < glayerY; iy++)
 		{
 			for (int iz = -glayerZ_minus; iz < glayerZ; iz++)
 			{
-				for (int ia = 0; ia < d_amount; ia++)
+				for (int i = 0;i < ntype;i++)
 				{
-					store_x[ia_all] = x_old[ia] + vec1[0] * ix + vec2[0] * iy + vec3[0] * iz;
-					store_y[ia_all] = y_old[ia] + vec1[1] * ix + vec2[1] * iy + vec3[1] * iz;
-					store_z[ia_all] = z_old[ia] + vec1[2] * ix + vec2[2] * iy + vec3[2] * iz;
-					store_type[ia_all] = type_old[ia];
-					store_natom[ia_all] = natom_old[ia];
-
-					// notice '0' is not the origin unitcell.
-					store_cell_x[ia_all] = ix + glayerX_minus;
-					store_cell_y[ia_all] = iy + glayerY_minus;
-					store_cell_z[ia_all] = iz + glayerZ_minus;
-
-					if (test_atom_input > 1)
+					for (int j = 0;j < ucell.atoms[i].na;j++)
 					{
-						if (d_amount_expand < 1000)
-						{
-							GlobalV::ofs_running << "\n" << std::setw(6) << ia_all
-							<< std::setw(10) << x_old[ia]
-							<< std::setw(10) << y_old[ia]
-							<< std::setw(10) << z_old[ia]
-							<< std::setw(10) << store_x[ia_all]
-							<< std::setw(10) << store_y[ia_all]
-							<< std::setw(10) << store_z[ia_all]
-							<< std::setw(6) << store_cell_x[ia_all]
-							<< std::setw(6) << store_cell_y[ia_all]
-							<< std::setw(6) << store_cell_z[ia_all];
-						}
+						FAtom fake_atom;
+						fake_atom.setX(ucell.atoms[i].tau[j].x + vec1[0] * ix + vec2[0] * iy + vec3[0] * iz);
+						fake_atom.setY(ucell.atoms[i].tau[j].y + vec1[1] * ix + vec2[1] * iy + vec3[1] * iz);
+						fake_atom.setZ(ucell.atoms[i].tau[j].z + vec1[2] * ix + vec2[2] * iy + vec3[2] * iz);
+						fake_atom.setType(i);
+						fake_atom.setNatom(j);
+						fake_atom.setCellX(ix);
+						fake_atom.setCellY(iy);
+						fake_atom.setCellZ(iz);
+
+						this->fake_atoms.push_back(fake_atom);
 					}
-
-					ia_all++;
 				}
 			}
 		}
 	}
 
-	// mohan fix bug 2012-04-02 by valgrind.
-	delete[] store_cell_x;
-	delete[] store_cell_y;
-	delete[] store_cell_z;
 
-	assert(ia_all == d_amount_expand);
+	assert(this->fake_atoms.size() == d_amount_expand);
 
 	// becareful! now the cell is not the origin cell,
 	// it's unitcell in expand case! so don't add
@@ -554,12 +490,6 @@ void Atom_input::Expand_Grid(const UnitCell &ucell, const int ntype)
 			<< " Ymax=" << y_max_expand
 			<< " Zmax=" << z_max_expand << std::endl;
 	}
-
-	delete[] x_old;
-	delete[] y_old;
-	delete[] z_old;
-	delete[] type_old;
-	delete[] natom_old;
 	return;
 }
 
@@ -614,11 +544,7 @@ void Atom_input::set_FAtom(const UnitCell &ucell, FAtom &a)const
 //----------------------------------------------------------
 	if (expand_flag)
 	{
-		a.setX(store_x[d_current]);
-		a.setY(store_y[d_current]);
-		a.setZ(store_z[d_current]);
-		a.setType(store_type[d_current]);
-		a.setNatom(store_natom[d_current]);
+		a = fake_atoms[d_current];
 		++ d_current;
 	}
 
@@ -630,9 +556,7 @@ void Atom_input::set_FAtom(const UnitCell &ucell, FAtom &a)const
 //----------------------------------------------------------
 	else
 	{
-		Load_atom(ucell
-
-		);
+		Load_atom(ucell);
 		a.setX(x);
 		a.setY(y);
 		a.setZ(z);

source/module_cell/module_neighbor/sltk_atom_input.h

@@ -118,6 +118,8 @@ class Atom_input
 
 	int getGrid_layerZ_minus(void) const { return glayerZ_minus;}
 
+	FAtom getFakeAtom(const int index) const { return fake_atoms[index];}
+
 private:
 	int test_atom_input;	//caoyu reconst 2021-05-24
 	int d_amount;//number of atoms.
@@ -154,14 +156,9 @@ class Atom_input
 // NAME : Expand_Grid
 //==========================================================
 	void Expand_Grid(const UnitCell& ucell, const int ntype);
-	double* store_x;
-	double* store_y;
-	double* store_z;
-	int* store_cell_x;
-	int* store_cell_y;
-	int* store_cell_z;
-	int* store_type;
-	int* store_natom;
+
+	std::vector<FAtom> fake_atoms;
+
 	double x_min_expand;
 	double y_min_expand;
 	double z_min_expand;

source/module_cell/module_neighbor/sltk_grid.cpp

@@ -401,49 +401,9 @@ void Grid::In_Which_Cell(const UnitCell &ucell, int &a, int &b, int &c, const FA
 // A bug remain!
 // d_minX, d_minY, d_minZ must be the cell origin
 //----------------------------------------------------------
-		double directx, directy, directz;
-		ModuleBase::Mathzone::Cartesian_to_Direct(
-		    atom.x(), atom.y(), atom.z(),
-		    vec1[0], vec1[1], vec1[2],
-		    vec2[0], vec2[1], vec2[2],
-		    vec3[0], vec3[1], vec3[2],
-		    directx, directy, directz);
-	
-		// cell_x_length
-		// mohan add the 'save_add' 2011-03-12
-		// for example, if a is 2.999999999999999999, the int command
-		// will make a = 2 (in fact we want 3)
-		// if a is 3.00000000000000001, it's save.
-		//static double save_add = 1.0e-15; // from -8 to -15, mohan add 2012-03-22
-		//a = static_cast<int>( (directx - this->d_minX) + save_add );
-		//b = static_cast<int>( (directy - this->d_minY) + save_add );
-		//c = static_cast<int>( (directz - this->d_minZ) + save_add );
-		int now_type = atom.getType();
-		int now_number = atom.getNatom();
-		double now_x_d = ucell.atoms[now_type].taud[now_number].x;
-		double now_y_d = ucell.atoms[now_type].taud[now_number].y;
-		double now_z_d = ucell.atoms[now_type].taud[now_number].z;
-		a = static_cast<int>(directx - now_x_d - this->d_minX + 0.5 );
-                b = static_cast<int>(directy - now_y_d - this->d_minY + 0.5 );
-                c = static_cast<int>(directz - now_z_d - this->d_minZ + 0.5 );
-
-		//ofs_running << std::setw(8) << atom.x() << std::setw(8) << atom.y() << std::setw(8) << atom.z()
-		//<< std::setw(12) << directx << std::setw(12) << directy << std::setw(12) << directz
-		//<< std::setw(6) << a << std::setw(6) << b << std::setw(6) << c << std::endl;
-	
-	/*	
-		if(a==0 && b==4 && c==0)
-		{
-			std::cout << std::setprecision(25) << std::endl;
-			std::cout << " save_add=" << save_add << std::endl;
-			std::cout << atom.x() << " " << atom.y() << " " << atom.z() << std::endl;
-			std::cout << " directy=" << directy << " d_minX=" << d_minY << " b=" << b << std::endl;
-			std::cout << " (int)directy=" << (int)directy << " (int)d_minY=" << (int)d_minY << " static_cast<int>( (directx - this->d_minX) )=" << static_cast<int>( (directy - this->d_minY) ) << std::endl;
-			std::cout << std::endl;
-			int ok;
-			cin >> ok;
-		}
-		*/
+		a = atom.getCellX() - this->d_minX;
+		b = atom.getCellY() - this->d_minY;
+		c = atom.getCellZ() - this->d_minZ;
 	}
 	else
 	{
@@ -500,7 +460,7 @@ void Grid::Build_Hash_Table(const UnitCell &ucell, AtomLink* const pointCache)
 	Hash_one_hit = 0; // mohan add 2010-06-25
 	for (; current < end; ++ current)
 	{
-		AtomLink* const hashTarget = this->getHashCode(ucell, current->fatom);
+				AtomLink* const hashTarget = this->getHashCode(ucell, current->fatom);
 
 		//================================================
 		// Find a new position

source/module_cell/module_neighbor/sltk_grid_driver.cpp

@@ -75,14 +75,17 @@ int Grid_Driver::Locate_offset(
 	temp.fatom.setX(cartesian_pos.x);
 	temp.fatom.setY(cartesian_pos.y);
 	temp.fatom.setZ(cartesian_pos.z);
+	temp.fatom.setCellX(0);
+	temp.fatom.setCellY(0);
+	temp.fatom.setCellZ(0);
 	temp.fatom.setType(ntype);
 	temp.fatom.setNatom(nnumber);
 
 //----------------------------------------------------------
 // EXPLAIN : Find the Hash number of this atom position
 //----------------------------------------------------------
 	AtomLink* Search = this->getHashCode(ucell, temp.fatom);
-
+	
 //----------------------------------------------------------
 // EXPLAIN : If we don't get the index for one Hash try,
 // we need to search in array.
@@ -196,28 +199,6 @@ void Grid_Driver::Find_adjacent_atom(const int offset, std::shared_ptr<AdjacentS
 	return;
 }
 
-//==========================================================
-// For expand case
-//==========================================================
-double Grid_Driver::Distance(const AtomLink& a1, const AtomLink& a2)const
-{
-	const double dx = a1.fatom.x() - a2.fatom.x();
-	const double dy = a1.fatom.y() - a2.fatom.y();
-	const double dz = a1.fatom.z() - a2.fatom.z();
-	return sqrt(dx*dx + dy*dy + dz*dz);
-}
-
-//==========================================================
-// For not_expand case
-//==========================================================
-double Grid_Driver::Distance(const AtomLink& a1, const ModuleBase::Vector3<double> &adjacent_site)const
-{
-	const double dx = a1.fatom.x() - adjacent_site.x;
-	const double dy = a1.fatom.y() - adjacent_site.y;
-	const double dz = a1.fatom.z() - adjacent_site.z;
-	return sqrt(dx*dx + dy*dy + dz*dz);
-}
-
 
 ModuleBase::Vector3<double> Grid_Driver::Calculate_adjacent_site
 (
@@ -241,30 +222,6 @@ ModuleBase::Vector3<double> Grid_Driver::Calculate_adjacent_site
 	return adjacent_site;
 }
 
-AdjacentAtomInfo Grid_Driver::get_adjs(const UnitCell& ucell_in, const size_t &iat)
-{
-    const int it = ucell_in.iat2it[iat];
-    const int ia = ucell_in.iat2ia[iat];
-    const ModuleBase::Vector3<double> &tau = ucell_in.atoms[it].tau[ia];
-
-    AdjacentAtomInfo adjs;
-    this->Find_atom(ucell_in, tau, it, ia, &adjs);
-    return adjs;
-}
-
-std::vector<AdjacentAtomInfo> Grid_Driver::get_adjs(const UnitCell& ucell_in)
-{
-    std::vector<AdjacentAtomInfo> adjs(ucell_in.nat);
-#ifdef _OPENMP
-#pragma omp parallel for
-#endif
-    for(size_t iat = 0; iat < ucell_in.nat; iat++)
-    {
-        adjs[iat] = Grid_Driver::get_adjs(ucell_in, iat);
-    }
-    return adjs;
-}
-
 // filter_adjs delete not adjacent atoms in adjs
 void filter_adjs(const std::vector<bool>& is_adj, AdjacentAtomInfo& adjs)
 {

source/module_cell/module_neighbor/sltk_grid_driver.h

@@ -77,13 +77,6 @@ class Grid_Driver : public Grid
 	const ModuleBase::Vector3<double>& getAdjacentTau(const int i) const { return adj_info.adjacent_tau[i]; } 
 	const ModuleBase::Vector3<int>& getBox(const int i) const {return adj_info.box[i];}
 
-	//==========================================================
-	// get_adjs will not store results in Grid_Driver::adj_info,
-	// but return the AdjacentAtomInfo object instead
-	//==========================================================
-    AdjacentAtomInfo get_adjs(const UnitCell& ucell_in, const size_t &iat);
-    std::vector<AdjacentAtomInfo> get_adjs(const UnitCell& ucell_in);
-
 private:
 
 	mutable AdjacentAtomInfo adj_info;
@@ -107,9 +100,6 @@ class Grid_Driver : public Grid
 		std::shared_ptr<AdjacentSet> as,
 		AdjacentAtomInfo &adjs)const;
 
-	double Distance(const AtomLink& a1, const ModuleBase::Vector3<double> &adjacent_site)const;
-	double Distance(const AtomLink& a1, const AtomLink& a2)const;
-
 //==========================================================
 // MEMBER FUNCTIONS :
 // NAME : Calculate_adjacent_site