diff --git a/source/module_base/intarray.h b/source/module_base/intarray.h
index f8db1250753..64cacf014b7 100644
--- a/source/module_base/intarray.h
+++ b/source/module_base/intarray.h
@@ -5,10 +5,10 @@
 #ifndef INTARRAY_H
 #define INTARRAY_H
 
-#include <iostream>
+#include <cassert>
 #include <fstream>
 #include <iomanip>
-#include <cassert>
+#include <iostream>
 
 #ifdef _MCD_CHECK
 //#include "./src_parallel/mcd.h"
@@ -16,63 +16,141 @@
 
 namespace ModuleBase
 {
-
+/**
+ * @brief Integer array
+ *
+ */
 class IntArray
 {
-public:
-	int *ptr;
-
-	// Constructors for different dimesnions
-	IntArray(const int d1 = 1, const int d2 = 1);
-	IntArray(const int d1, const int d2,const int d3);
-	IntArray(const int d1, const int d2,const int d3,const int d4);
-	IntArray(const int d1, const int d2,const int d3,const int d4,const int d5);
-	IntArray(const int d1, const int d2,const int d3,const int d4,const int d5,const int d6);
-
-	~IntArray();
-
-	void create(const int d1, const int d2);
-	void create(const int d1, const int d2, const int d3);
-	void create(const int d1, const int d2, const int d3, const int d4);
-	void create(const int d1, const int d2, const int d3, const int d4, const int d5);
-	void create(const int d1, const int d2, const int d3, const int d4, const int d5, const int d6);
-
-	const IntArray &operator=(const IntArray &right);
-	const IntArray &operator=(const int &right);
-
-	int &operator()(const int d1, const int d2);
-	int &operator()(const int d1, const int d2, const int d3);
-	int &operator()(const int d1, const int d2, const int d3,const int d4);
-	int &operator()(const int d1, const int d2, const int d3, const int d4, const int d5);
-	int &operator()(const int d1, const int d2, const int d3, const int d4, const int d5, const int d6);
-
-	const int &operator()(const int d1,const int d2)const;
-	const int &operator()(const int d1,const int d2,const int d3)const;
-	const int &operator()(const int d1,const int d2,const int d3,const int d4)const;
-	const int &operator()(const int d1,const int d2,const int d3,const int d4, const int d5)const;
-	const int &operator()(const int d1,const int d2,const int d3,const int d4, const int d5, const int d6)const;
-
-	void zero_out(void);
-
-	int getSize() const{ return size;}
-	int getDim() const{ return dim;}
-	int getBound1() const{ return bound1;}
-	int getBound2() const{ return bound2;}
-	int getBound3() const{ return bound3;}
-	int getBound4() const { return bound4;}
-	int getBound5() const { return bound5;}
-	int getBound6() const { return bound6;}
-
-	static int getArrayCount(void)
-	{ return arrayCount;}
-
-private:
-	int size;
-	int dim;
-	int bound1, bound2, bound3, bound4, bound5, bound6;
-	static int arrayCount;
-	void freemem();
+  public:
+    int *ptr;
+
+    /**
+     * @brief Construct a new Int Array object
+     *
+     * @param d1 The first dimension size
+     * @param d2 The second dimension size
+     */
+    IntArray(const int d1 = 1, const int d2 = 1);
+    IntArray(const int d1, const int d2, const int d3);
+    IntArray(const int d1, const int d2, const int d3, const int d4);
+    IntArray(const int d1, const int d2, const int d3, const int d4, const int d5);
+    IntArray(const int d1, const int d2, const int d3, const int d4, const int d5, const int d6);
+
+    ~IntArray();
+
+    /**
+     * @brief Create integer arrays
+     *
+     * @param[in] d1
+     * @param[in] d2
+     */
+    void create(const int d1, const int d2);
+    void create(const int d1, const int d2, const int d3);
+    void create(const int d1, const int d2, const int d3, const int d4);
+    void create(const int d1, const int d2, const int d3, const int d4, const int d5);
+    void create(const int d1, const int d2, const int d3, const int d4, const int d5, const int d6);
+
+    /**
+     * @brief Equal an IntArray to another one
+     *
+     * @param right
+     * @return const IntArray&
+     */
+    const IntArray &operator=(const IntArray &right);
+
+    /**
+     * @brief Equal all elements of an IntArray to an
+     * integer
+     *
+     * @param right
+     * @return const IntArray&
+     */
+    const IntArray &operator=(const int &right);
+
+    /**
+     * @brief Access elements by using operator "()"
+     *
+     * @param d1
+     * @param d2
+     * @return int&
+     */
+    int &operator()(const int d1, const int d2);
+    int &operator()(const int d1, const int d2, const int d3);
+    int &operator()(const int d1, const int d2, const int d3, const int d4);
+    int &operator()(const int d1, const int d2, const int d3, const int d4, const int d5);
+    int &operator()(const int d1, const int d2, const int d3, const int d4, const int d5, const int d6);
+
+    /**
+     * @brief Access elements by using "()" through pointer
+     * without changing its elements
+     *
+     * @param d1
+     * @param d2
+     * @return const int&
+     */
+    const int &operator()(const int d1, const int d2) const;
+    const int &operator()(const int d1, const int d2, const int d3) const;
+    const int &operator()(const int d1, const int d2, const int d3, const int d4) const;
+    const int &operator()(const int d1, const int d2, const int d3, const int d4, const int d5) const;
+    const int &operator()(const int d1, const int d2, const int d3, const int d4, const int d5, const int d6) const;
+
+    /**
+     * @brief Set all elements of an IntArray to zero
+     *
+     */
+    void zero_out(void);
+
+    int getSize() const
+    {
+        return size;
+    }
+    int getDim() const
+    {
+        return dim;
+    }
+    int getBound1() const
+    {
+        return bound1;
+    }
+    int getBound2() const
+    {
+        return bound2;
+    }
+    int getBound3() const
+    {
+        return bound3;
+    }
+    int getBound4() const
+    {
+        return bound4;
+    }
+    int getBound5() const
+    {
+        return bound5;
+    }
+    int getBound6() const
+    {
+        return bound6;
+    }
+
+    /**
+     * @brief Get the Array Count object
+     *
+     * @return int
+     */
+    static int getArrayCount(void)
+    {
+        return arrayCount;
+    }
+
+  private:
+    int size;
+    int dim;
+    int bound1, bound2, bound3, bound4, bound5, bound6;
+    static int arrayCount;
+    void freemem();
 };
-}
+} // namespace ModuleBase
 
-#endif	// IntArray class
+#endif // IntArray class
diff --git a/source/module_base/mathzone.h b/source/module_base/mathzone.h
index 99622ac29a1..6e5180afd8c 100644
--- a/source/module_base/mathzone.h
+++ b/source/module_base/mathzone.h
@@ -1,69 +1,82 @@
 #ifndef MATHZONE_H
 #define MATHZONE_H
 
-#include "realarray.h"
-#include "matrix3.h"
 #include "global_function.h"
-#include <vector>
-#include <map>
+#include "matrix3.h"
+#include "realarray.h"
+
 #include <cassert>
 #include <complex>
+#include <map>
+#include <vector>
 namespace ModuleBase
 {
 
+/**
+ * @brief atomic coordinates conversion functions
+ *
+ */
 class Mathzone
 {
-	public:
-
+  public:
     Mathzone();
     ~Mathzone();
 
-    template<class T>
-    static T Max3(const T &a,const T &b,const T &c)
+  public:
+    /**
+     * @brief Pointwise product of two vectors with same size
+     *
+     * @tparam Type
+     * @param[in] f1
+     * @param[in] f2
+     * @return std::vector<Type>
+     * @author Peize Lin (2016-08-03)
+     */
+    template <typename Type>
+    static std::vector<Type> Pointwise_Product(const std::vector<Type> &f1, const std::vector<Type> &f2)
     {
-        if (a>=b && a>=c) return a;
-        else if (b>=a && b>=c) return b;
-        else if (c>=a && c>=b) return c;
-		else throw std::runtime_error(ModuleBase::GlobalFunc::TO_STRING(__FILE__)+" line "+ModuleBase::GlobalFunc::TO_STRING(__LINE__));
+        assert(f1.size() == f2.size());
+        std::vector<Type> f(f1.size());
+        for (int ir = 0; ir != f.size(); ++ir)
+            f[ir] = f1[ir] * f2[ir];
+        return f;
     }
 
-    // be careful, this can only be used for plane wave
-    // during parallel calculation
-
-public:
-
-
-
-	// Peize Lin add 2016-08-03
-	template< typename Type >
-	static std::vector<Type> Pointwise_Product( const std::vector<Type> &f1, const std::vector<Type> &f2 )
-	{
-		assert(f1.size()==f2.size());
-		std::vector<Type> f(f1.size());
-		for( int ir=0; ir!=f.size(); ++ir )
-			f[ir] = f1[ir] * f2[ir];
-		return f;
-	}
-
-//==========================================================
-// MEMBER FUNCTION :
-// NAME : Direct_to_Cartesian
-// use lattice vector matrix R
-// change the direct std::vector (dx,dy,dz) to cartesuab vectir
-// (cx,cy,cz)
-// (dx,dy,dz) = (cx,cy,cz) * R
-//
-// NAME : Cartesian_to_Direct
-// the same as above
-// (cx,cy,cz) = (dx,dy,dz) * R^(-1)
-//==========================================================
-    static inline void Direct_to_Cartesian
-    (
-        const double &dx,const double &dy,const double &dz,
-        const double &R11,const double &R12,const double &R13,
-        const double &R21,const double &R22,const double &R23,
-        const double &R31,const double &R32,const double &R33,
-        double &cx,double &cy,double &cz)
+    /**
+     * @brief change direct coordinate (dx,dy,dz) to
+     * Cartesian coordinate (cx,cy,cz), (dx,dy,dz) = (cx,cy,cz) * R
+     *
+     * @param[in] dx Direct coordinats
+     * @param[in] dy
+     * @param[in] dz
+     * @param[in] R11 Lattice vector matrix R_ij: i_row, j_column
+     * @param[in] R12
+     * @param[in] R13
+     * @param[in] R21
+     * @param[in] R22
+     * @param[in] R23
+     * @param[in] R31
+     * @param[in] R32
+     * @param[in] R33
+     * @param[out] cx Cartesian coordinats
+     * @param[out] cy
+     * @param[out] cz 
+     */
+    static inline void Direct_to_Cartesian(const double &dx,
+                                           const double &dy,
+                                           const double &dz,
+                                           const double &R11,
+                                           const double &R12,
+                                           const double &R13,
+                                           const double &R21,
+                                           const double &R22,
+                                           const double &R23,
+                                           const double &R31,
+                                           const double &R32,
+                                           const double &R33,
+                                           double &cx,
+                                           double &cy,
+                                           double &cz)
     {
         static ModuleBase::Matrix3 lattice_vector;
         static ModuleBase::Vector3<double> direct_vec, cartesian_vec;
@@ -88,13 +101,41 @@ class Mathzone
         return;
     }
 
-    static inline void Cartesian_to_Direct
-    (
-        const double &cx,const double &cy,const double &cz,
-        const double &R11,const double &R12,const double &R13,
-        const double &R21,const double &R22,const double &R23,
-        const double &R31,const double &R32,const double &R33,
-        double &dx,double &dy,double &dz)
+    /**
+     * @brief Change Cartesian coordinate (cx,cy,cz) to
+     * direct coordinate (dx,dy,dz), (cx,cy,cz) = (dx,dy,dz) * R^(-1)
+     *
+     * @param[in] cx Cartesian coordinats
+     * @param[in] cy
+     * @param[in] cz
+     * @param[in] R11 Lattice vector matrix R_ij: i_row, j_column
+     * @param[in] R12
+     * @param[in] R13
+     * @param[in] R21
+     * @param[in] R22
+     * @param[in] R23
+     * @param[in] R31
+     * @param[in] R32
+     * @param[in] R33
+     * @param[out] dx Direct coordinats
+     * @param[out] dy
+     * @param[out] dz
+     */
+    static inline void Cartesian_to_Direct(const double &cx,
+                                           const double &cy,
+                                           const double &cz,
+                                           const double &R11,
+                                           const double &R12,
+                                           const double &R13,
+                                           const double &R21,
+                                           const double &R22,
+                                           const double &R23,
+                                           const double &R31,
+                                           const double &R32,
+                                           const double &R33,
+                                           double &dx,
+                                           double &dy,
+                                           double &dz)
     {
         static ModuleBase::Matrix3 lattice_vector, inv_lat;
         lattice_vector.e11 = R11;
@@ -120,45 +161,17 @@ class Mathzone
         return;
     }
 
-
-    static void To_Polar_Coordinate
+    void To_Polar_Coordinate
     (
-        const double &x_cartesian,
-        const double &y_cartesian,
-        const double &z_cartesian,
-        double &r,
-        double &theta,
-        double &phi
-    );
-	
-	
-	// coeff1 * x1 + (1-coeff1) * x2
-	// Peize Lin add 2017-08-09
-	template< typename T, typename T_coeff >
-	static T Linear_Mixing( const T & x1, const T & x2, const T_coeff & coeff1 )
-	{
-		return coeff1 * x1 + (1-coeff1) * x2;
-	}
-	template< typename T, typename T_coeff >
-	static std::vector<T> Linear_Mixing( const std::vector<T> & x1, const std::vector<T> & x2, const T_coeff & coeff1 )
-	{
-		assert(x1.size()==x2.size());
-		std::vector<T> x;
-		for( size_t i=0; i!=x1.size(); ++i )
-			x.push_back( Linear_Mixing( x1[i], x2[i], coeff1 ) );
-		return x;
-	}
-	template< typename T1, typename T2, typename T_coeff >
-	static std::map<T1,T2> Linear_Mixing( const std::map<T1,T2> & x1, const std::map<T1,T2> & x2, const T_coeff & coeff1 )
-	{
-		std::map<T1,T2> x;
-		for( const auto & x1i : x1 )
-			x.insert( make_pair( x1i.first, Linear_Mixing( x1i.second, x2.at(x1i.first), coeff1 ) ) );
-		return x;
-	}	
+    const double &x_cartesian,
+    const double &y_cartesian,
+    const double &z_cartesian,
+    double &r,
+    double &theta,
+    double &phi);
 
 };
 
-}
+} // namespace ModuleBase
 
 #endif
diff --git a/source/module_base/realarray.cpp b/source/module_base/realarray.cpp
index e6a0e60ddbe..757bb9ee386 100644
--- a/source/module_base/realarray.cpp
+++ b/source/module_base/realarray.cpp
@@ -53,6 +53,22 @@ realArray::realArray(const int d1,const int d2,const int d3,const int d4)
 	++arrayCount;
 }
 
+realArray::realArray(const realArray &cd)
+{
+	this->size = cd.getSize();
+	this->ptr = new double[size];
+	for (int i = 0; i < size; i++)
+		this->ptr[i] = cd.ptr[i];
+	this->dim = cd.dim;
+	this->bound1 = cd.bound1;
+	this->bound2 = cd.bound2;
+	this->bound3 = cd.bound3;
+	this->bound4 = cd.bound4;
+
+	++arrayCount;
+}
+
+
 //********************************
 //
 // Destructor for class realArray
@@ -65,8 +81,8 @@ realArray ::~realArray()
 
 void realArray::freemem()
 {
-	delete [] ptr;
-	ptr = NULL;
+		delete [] ptr;
+		ptr = NULL;
 }
 
 void realArray::create(const int d1,const int d2,const int d3,const int d4)
@@ -173,4 +189,4 @@ void realArray::zero_out(void)
 	return;
 }
 
-}
\ No newline at end of file
+}
diff --git a/source/module_base/realarray.h b/source/module_base/realarray.h
index 45cc5fc84c4..1dc2d7deeff 100644
--- a/source/module_base/realarray.h
+++ b/source/module_base/realarray.h
@@ -5,82 +5,162 @@
 #ifndef REALARRAY_H
 #define REALARRAY_H
 
-#include <iostream>
+#include <cassert>
 #include <fstream>
 #include <iomanip>
-#include <cassert>
+#include <iostream>
 
 #ifdef _MCD_CHECK
 //#include "./src_parallel/mcd.h"
 #endif
 
-
 namespace ModuleBase
 {
-
+/**
+ * @brief double float array
+ *
+ */
 class realArray
 {
-public:
-	double *ptr;
-
-	realArray(const int d1 = 1 ,const int d2 = 1,const int d3 = 1);
-	realArray(const int d1, const int d2,const int d3,const int d4);
-	~realArray();
-
-	void create(const int d1,const int d2,const int d3);
-	void create(const int d1,const int d2,const int d3,const int d4);
-
-	const realArray &operator=(const realArray &right);
-	const realArray &operator=(const double &right);
-
-	double &operator()(const int d1,const int d2,const int d3);
-	double &operator()(const int d1,const int d2,const int d3,const int d4);
-
-	const double &operator()(const int d1,const int d2,const int d3)const;
-	const double &operator()(const int d1,const int d2,const int d3,const int d4)const;
-
-	void zero_out(void);
-
-	int getSize() const
-	{ return size;}
-
-	int getDim() const
-	{ return dim;}
-
-	int getBound1() const
-	{ return bound1;}
-
-	int getBound2() const
-	{ return bound2;}
-
-	int getBound3() const
-	{ return bound3;}
-
-	int getBound4() const
-	{ return bound4;}
-
-	static int getArrayCount(void)
-	{ return arrayCount;}
-
-private:
-	int size;
-	int dim;
-	int bound1, bound2, bound3, bound4;
-	static int arrayCount;
-
-	void freemem();
+  public:
+    double *ptr;
+
+    realArray(const int d1 = 1, const int d2 = 1, const int d3 = 1);
+    realArray(const int d1, const int d2, const int d3, const int d4);
+    ~realArray();
+
+    /**
+     * @brief create 3 dimensional real array
+     *
+     * @param[in] d1 The first dimension size
+     * @param[in] d2 The second dimension size
+     * @param[in] d3 The third dimension size
+     */
+    void create(const int d1, const int d2, const int d3);
+    void create(const int d1, const int d2, const int d3, const int d4);
+
+    realArray(const realArray &cd);
+
+    /**
+     * @brief Equal a realArray to another one
+     *
+     * @param right
+     * @return const realArray&
+     */
+    const realArray &operator=(const realArray &right);
+    /**
+     * @brief Set all value of an array to a double float number
+     *
+     * @param right
+     * @return const realArray&
+     */
+    const realArray &operator=(const double &right);
+
+    /**
+     * @brief Access elements by using operator "()"
+     *
+     * @param d1
+     * @param d2
+     * @param d3
+     * @return double&
+     */
+    double &operator()(const int d1, const int d2, const int d3);
+    double &operator()(const int d1, const int d2, const int d3, const int d4);
+
+    /**
+     * @brief Access elements by using "()" through pointer
+     * without changing its elements
+     *
+     * @param d1
+     * @param d2
+     * @param d3
+     * @return const double&
+     */
+    const double &operator()(const int d1, const int d2, const int d3) const;
+    const double &operator()(const int d1, const int d2, const int d3, const int d4) const;
+
+    /**
+     * @brief Set all elements of an IntArray to zero
+     *
+     */
+    void zero_out(void);
+
+    /**
+     * @brief Get the Size object
+     *
+     * @return int
+     */
+    int getSize() const
+    {
+        return size;
+    }
+
+    /**
+     * @brief Get the Dim object
+     * i.e. the dimension of a real array
+     *
+     * @return int
+     */
+    int getDim() const
+    {
+        return dim;
+    }
+
+    /**
+     * @brief Get the Bound1 object
+     * i.e. the first dimension size
+     *
+     * @return int
+     */
+    int getBound1() const
+    {
+        return bound1;
+    }
+
+    int getBound2() const
+    {
+        return bound2;
+    }
+
+    int getBound3() const
+    {
+        return bound3;
+    }
+
+    int getBound4() const
+    {
+        return bound4;
+    }
+
+    /**
+     * @brief Get the Array Count object
+     *
+     * @return int
+     */
+    static int getArrayCount(void)
+    {
+        return arrayCount;
+    }
+
+  private:
+    int size;
+    int dim;
+    int bound1, bound2, bound3, bound4;
+    static int arrayCount;
+
+    void freemem();
 };
 
 //**************************************************
 // set elements of a as zeros which a is 1_d array.
 //**************************************************
-template<class T>
-void zeros(T *u,const int n)
+template <class T> void zeros(T *u, const int n)
 {
-	assert(n>0);
-	for (int i = 0;i < n;i++) u[i] = 0;
+    assert(n > 0);
+    for (int i = 0; i < n; i++)
+        u[i] = 0;
 }
 
-}
+} // namespace ModuleBase
 
-#endif	// realArray class
+#endif // realArray class
diff --git a/source/module_base/test/CMakeLists.txt b/source/module_base/test/CMakeLists.txt
index b7e9a9b7281..51f65a49720 100644
--- a/source/module_base/test/CMakeLists.txt
+++ b/source/module_base/test/CMakeLists.txt
@@ -39,7 +39,24 @@ AddTest(
   SOURCES complexmatrix_test.cpp ../complexmatrix.cpp ../matrix.cpp 
 )
 AddTest(
-  TARGET base_matrix
+ TARGET base_matrix
   LIBS ${math_libs}
   SOURCES matrix_test.cpp ../matrix.cpp
+)
+AddTest(
+  TARGET base_realarray
+  SOURCES realarray_test.cpp ../realarray.cpp
+)
+AddTest(
+  TARGET base_intarray
+  SOURCES intarray_test.cpp ../intarray.cpp
+)
+AddTest(
+  TARGET base_vector3
+  SOURCES vector3_test.cpp 
+)
+AddTest(
+  TARGET base_mathzone
+  LIBS ${math_libs}
+  SOURCES mathzone_test.cpp ../mathzone.cpp ../matrix3.cpp ../matrix.cpp ../tool_quit.cpp ../global_variable.cpp ../global_file.cpp ../memory.cpp ../timer.cpp
 )
\ No newline at end of file
diff --git a/source/module_base/test/intarray_test.cpp b/source/module_base/test/intarray_test.cpp
new file mode 100644
index 00000000000..d6694abea38
--- /dev/null
+++ b/source/module_base/test/intarray_test.cpp
@@ -0,0 +1,308 @@
+#include "../intarray.h"
+#include "gtest/gtest.h"
+
+/************************************************
+ *  unit test of class IntArray
+ ***********************************************/
+
+/**
+ * - Tested Functions:
+ *   - Construct
+ *     - construct an int array (2 to 6 dimensions)
+ *   - Creat 
+ *     - create an int array (2 to 6 dimensions)
+ *   - GetArrayCount
+ *     - get the total number of int array created
+ *   - GetSize
+ *     - get the total size of an int array
+ *   - GetDim
+ *     - get the dimension of an int array
+ *   - ZeroOut
+ *     - set all elements of an int array to zero
+ *   - GetBound
+ *     - get the size of each dimension of an int array
+ *   - ArrayEqReal
+ *     - set all value of an array to an int number
+ *   - ArrayEqArray
+ *     - equal an intarray to another intarray
+ *   - Parentheses
+ *     - access element by using operator"()"
+ *   - ConstParentheses
+ *     - access element by using "()" through pointer
+ *     - without changing its elements
+ */
+
+class IntArrayTest : public testing::Test
+{
+protected:
+	ModuleBase::IntArray a2, a3, a4, a5, a6;
+	int aa = 11;
+	int bb = 1;
+	int count0;
+	int count1;
+	const int zero = 0;
+};
+
+TEST_F(IntArrayTest,GetArrayCount)
+{
+	count0 = ModuleBase::IntArray::getArrayCount();
+	ModuleBase::IntArray c3, c4;
+	count1 = ModuleBase::IntArray::getArrayCount();
+	EXPECT_EQ((count1-count0),2);
+}
+
+TEST_F(IntArrayTest,Construct)
+{
+	ModuleBase::IntArray x2(1,5);
+	ModuleBase::IntArray x3(1,5,3);
+	ModuleBase::IntArray x4(1,7,3,4);
+	ModuleBase::IntArray x5(1,5,3,8,2);
+	ModuleBase::IntArray x6(1,7,3,4,3,2);
+	EXPECT_EQ(x2.getSize(),5);
+	EXPECT_EQ(x3.getSize(),15);
+	EXPECT_EQ(x4.getSize(),84);
+	EXPECT_EQ(x5.getSize(),240);
+	EXPECT_EQ(x6.getSize(),504);
+}
+
+TEST_F(IntArrayTest,Create)
+{
+	a2.create(2,1);
+	a3.create(3,2,1);
+	a4.create(4,3,2,1);
+	a5.create(5,4,3,2,1);
+	a6.create(6,5,4,3,2,1);
+	EXPECT_EQ(a2.getSize(),2);
+	EXPECT_EQ(a3.getSize(),6);
+	EXPECT_EQ(a4.getSize(),24);
+	EXPECT_EQ(a5.getSize(),120);
+	EXPECT_EQ(a6.getSize(),720);
+}
+
+
+TEST_F(IntArrayTest,GetSize)
+{
+	ModuleBase::IntArray x3(1,5,3);
+	ModuleBase::IntArray x4(1,7,3,4);
+	EXPECT_EQ(x3.getSize(),15);
+	EXPECT_EQ(x4.getSize(),84);
+}
+
+TEST_F(IntArrayTest,GetDim)
+{
+	a2.create(2,3);
+	a3.create(3,5,1);
+	a4.create(4,3,7,1);
+	a5.create(5,4,1,2,1);
+	a6.create(6,5,9,3,2,1);
+	EXPECT_EQ(a2.getDim(),2);
+	EXPECT_EQ(a3.getDim(),3);
+	EXPECT_EQ(a4.getDim(),4);
+	EXPECT_EQ(a5.getDim(),5);
+	EXPECT_EQ(a6.getDim(),6);
+}
+
+TEST_F(IntArrayTest,ZeroOut)
+{
+	a2.create(2,3);
+	a3.create(3,5,1);
+	a4.create(4,3,7,1);
+	a5.create(5,4,1,2,1);
+	a6.create(6,5,9,3,2,1);
+	a2.zero_out();
+	a3.zero_out();
+	a4.zero_out();
+	a5.zero_out();
+	a6.zero_out();
+	for (int i=0;i<a2.getSize();++i)
+	{
+		EXPECT_EQ(a2.ptr[i],zero);
+    	}
+	for (int i=0;i<a3.getSize();++i)
+	{
+		EXPECT_EQ(a3.ptr[i],zero);
+    	}
+	for (int i=0;i<a4.getSize();++i)
+	{
+		EXPECT_EQ(a4.ptr[i],zero);
+    	}
+	for (int i=0;i<a5.getSize();++i)
+	{
+		EXPECT_EQ(a5.ptr[i],zero);
+    	}
+	for (int i=0;i<a6.getSize();++i)
+	{
+		EXPECT_EQ(a6.ptr[i],zero);
+    	}
+}
+
+TEST_F(IntArrayTest,GetBound)
+{
+	// 2d
+	a2.create(2,3);
+	EXPECT_EQ(a2.getBound1(),2);
+	EXPECT_EQ(a2.getBound2(),3);
+	// 3d
+	a3.create(3,5,1);
+	EXPECT_EQ(a3.getBound1(),3);
+	EXPECT_EQ(a3.getBound2(),5);
+	EXPECT_EQ(a3.getBound3(),1);
+	EXPECT_EQ(a3.getBound4(),1);
+
+	ModuleBase::IntArray x3(1,5,3);
+	EXPECT_EQ(x3.getBound1(),1);
+	EXPECT_EQ(x3.getBound2(),5);
+	EXPECT_EQ(x3.getBound3(),3);
+	EXPECT_EQ(x3.getBound4(),0);
+	EXPECT_EQ(x3.getBound5(),0);
+	EXPECT_EQ(x3.getBound6(),0);
+	// 4d
+	a4.create(4,3,7,1);
+	EXPECT_EQ(a4.getBound1(),4);
+	EXPECT_EQ(a4.getBound2(),3);
+	EXPECT_EQ(a4.getBound3(),7);
+	EXPECT_EQ(a4.getBound4(),1);
+	// 5d
+	a5.create(5,4,1,2,1);
+	EXPECT_EQ(a5.getBound1(),5);
+	EXPECT_EQ(a5.getBound2(),4);
+	EXPECT_EQ(a5.getBound3(),1);
+	EXPECT_EQ(a5.getBound4(),2);
+	EXPECT_EQ(a5.getBound5(),1);
+	// 6d
+	a6.create(6,5,9,3,2,1);
+	EXPECT_EQ(a6.getBound1(),6);
+	EXPECT_EQ(a6.getBound2(),5);
+	EXPECT_EQ(a6.getBound3(),9);
+	EXPECT_EQ(a6.getBound4(),3);
+	EXPECT_EQ(a6.getBound5(),2);
+	EXPECT_EQ(a6.getBound6(),1);
+}
+
+TEST_F(IntArrayTest,ArrayEqInt)
+{
+	a2.create(2,3);
+	a3.create(3,5,1);
+	a4.create(4,3,7,1);
+	a5.create(5,4,1,2,1);
+	a6.create(6,5,9,3,2,1);
+	a2 = aa;
+	a3 = bb;
+	a4 = aa+1;
+	a5 = bb -1;
+	a6 = 10;
+	for (int i=0;i<a2.getSize();++i)
+	{
+		EXPECT_EQ(a2.ptr[i],aa);
+    	}
+	for (int i=0;i<a3.getSize();++i)
+	{
+		EXPECT_EQ(a3.ptr[i],bb);
+    	}
+	for (int i=0;i<a4.getSize();++i)
+	{
+		EXPECT_EQ(a4.ptr[i],aa+1);
+    	}
+	for (int i=0;i<a5.getSize();++i)
+	{
+		EXPECT_EQ(a5.ptr[i],bb-1);
+    	}
+	for (int i=0;i<a6.getSize();++i)
+	{
+		EXPECT_EQ(a6.ptr[i],10);
+    	}
+}
+
+TEST_F(IntArrayTest,ArrayEqArray)
+{
+	a2.create(2,3);
+	a3.create(3,5,1);
+	a4.create(4,3,7,1);
+	a5.create(5,4,1,2,1);
+	a6.create(6,5,9,3,2,1);
+	a2 = aa;
+	a3 = bb;
+	a4 = aa+1;
+	a5 = bb -1;
+	a6 = 10;
+	ModuleBase::IntArray x2(2,3);
+	ModuleBase::IntArray x3(3,5,1);
+	ModuleBase::IntArray x4(4,3,7,1);
+	ModuleBase::IntArray x5(5,4,1,2,1);
+	ModuleBase::IntArray x6(6,5,9,3,2,1);
+	x2 = a2; x3 = a3; x4 = a4; x5 = a5; x6 = a6;
+	for (int i=0;i<x2.getSize();++i)
+	{
+		EXPECT_EQ(x2.ptr[i],aa);
+    	}
+	for (int i=0;i<x3.getSize();++i)
+	{
+		EXPECT_EQ(x3.ptr[i],bb);
+    	}
+	for (int i=0;i<x4.getSize();++i)
+	{
+		EXPECT_EQ(x4.ptr[i],aa+1);
+    	}
+	for (int i=0;i<x5.getSize();++i)
+	{
+		EXPECT_EQ(x5.ptr[i],bb-1);
+    	}
+	for (int i=0;i<x6.getSize();++i)
+	{
+		EXPECT_EQ(x6.ptr[i],10);
+    	}
+}
+
+TEST_F(IntArrayTest,Parentheses)
+{
+	a2.create(2,3);
+	a3.create(3,5,1);
+	a4.create(4,3,7,1);
+	a5.create(5,4,1,2,1);
+	a6.create(6,5,9,3,2,1);
+	a2 = aa;
+	a3 = bb;
+	a4 = aa+1;
+	a5 = bb -1;
+	a6 = 10;
+	EXPECT_EQ(a2(0,0),aa);
+	EXPECT_EQ(a2(1,2),aa);
+	EXPECT_EQ(a3(0,0,0),bb);
+	EXPECT_EQ(a3(2,4,0),bb);
+	EXPECT_EQ(a4(0,0,0,0),aa+1);
+	EXPECT_EQ(a4(3,2,6,0),aa+1);
+	EXPECT_EQ(a5(0,0,0,0,0),bb-1);
+	EXPECT_EQ(a5(4,3,0,1,0),bb-1);
+	EXPECT_EQ(a6(0,0,0,0,0,0),10);
+	EXPECT_EQ(a6(5,4,8,2,1,0),10);
+}
+
+TEST_F(IntArrayTest,ConstParentheses)
+{
+	a2.create(2,3);
+	a3.create(3,5,1);
+	a4.create(4,3,7,1);
+	a5.create(5,4,1,2,1);
+	a6.create(6,5,9,3,2,1);
+	a2 = aa;
+	a3 = bb;
+	a4 = aa+1;
+	a5 = bb -1;
+	a6 = 10;
+	const ModuleBase::IntArray *x2(&a2);
+	const ModuleBase::IntArray *x3(&a3);
+	const ModuleBase::IntArray *x4(&a4);
+	const ModuleBase::IntArray *x5(&a5);
+	const ModuleBase::IntArray *x6(&a6);
+	EXPECT_EQ((*x2)(0,0),aa);
+	EXPECT_EQ((*x2)(1,2),aa);
+	EXPECT_EQ((*x3)(0,0,0),bb);
+	EXPECT_EQ((*x3)(2,4,0),bb);
+	EXPECT_EQ((*x4)(0,0,0,0),aa+1);
+	EXPECT_EQ((*x4)(3,2,6,0),aa+1);
+	EXPECT_EQ((*x5)(0,0,0,0,0),bb-1);
+	EXPECT_EQ((*x5)(4,3,0,1,0),bb-1);
+	EXPECT_EQ((*x6)(0,0,0,0,0,0),10);
+	EXPECT_EQ((*x6)(5,4,8,2,1,0),10);
+}
+
diff --git a/source/module_base/test/mathzone_test.cpp b/source/module_base/test/mathzone_test.cpp
new file mode 100644
index 00000000000..1ad11dd16ff
--- /dev/null
+++ b/source/module_base/test/mathzone_test.cpp
@@ -0,0 +1,86 @@
+#include "../mathzone.h"
+#include "../matrix3.h"
+#include "../vector3.h"
+#include "gtest/gtest.h"
+#include "gmock/gmock.h"
+
+/************************************************
+ *  unit test of class Mathzone
+ ***********************************************/
+
+/**
+ * - Tested Functions:
+ *   - PointwiseProduct
+ *     - return a vector, which is the pointwise
+ *     - product of another two vectors of the same
+ *     - length
+ *   - Direct2Cartesian
+ *     - change atomic coordinates from direct
+ *     - to Cartesian
+ *   - Cartesian2Direct
+ *     - change atomic coordinates from Cartesian
+ *     - to Direct
+ */
+
+class MathzoneTest : public testing::Test
+{
+protected:
+	double R11 = 3.68; double R12 = 0.00; double R13 = 0.00;
+	double R21 = 0.00; double R22 = 10.1; double R23 = 0.00;
+	double R31 = 0.00; double R32 = 0.00; double R33 = 26.7;
+	ModuleBase::Matrix3 lattice;
+	ModuleBase::Vector3<double> direct, cartesian;
+};
+
+TEST_F(MathzoneTest, PointwiseProduct)
+{
+	std::vector<double> aa, bb, cc;
+	for(int i=0;i<10;i++)
+	{
+		aa.push_back(i*i);
+		bb.push_back(i*2);
+	}
+	cc = ModuleBase::Mathzone::Pointwise_Product(aa,bb);
+	for(int i=0;i<10;i++)
+	{
+		EXPECT_EQ(cc[i],i*i*i*2);
+	}
+}
+
+TEST_F(MathzoneTest, Direct2Cartesian)
+{
+	direct.set(0.1,0.2,0.4);
+	cartesian.set(0.368,2.02,10.68);
+	ModuleBase::Vector3<double> cartnew;
+	ModuleBase::Mathzone::Direct_to_Cartesian(direct.x,
+			direct.y,
+			direct.z,
+			R11, R12, R13,
+			R21, R22, R23,
+			R31, R32, R33,
+			cartnew.x,
+			cartnew.y,
+			cartnew.z);
+	EXPECT_NEAR(cartnew.x,cartesian.x, 1e-15);
+	EXPECT_NEAR(cartnew.y,cartesian.y, 1e-15);
+	EXPECT_NEAR(cartnew.z,cartesian.z, 1e-15);
+}
+
+TEST_F(MathzoneTest, Cartesian2Direct)
+{
+	direct.set(0.1,0.2,0.4);
+	cartesian.set(0.368,2.02,10.68);
+	ModuleBase::Vector3<double> directnew;
+	ModuleBase::Mathzone::Cartesian_to_Direct(cartesian.x,
+			cartesian.y,
+			cartesian.z,
+			R11, R12, R13,
+			R21, R22, R23,
+			R31, R32, R33,
+			directnew.x,
+			directnew.y,
+			directnew.z);
+	EXPECT_NEAR(directnew.x,direct.x, 1e-15);
+	EXPECT_NEAR(directnew.y,direct.y, 1e-15);
+	EXPECT_NEAR(directnew.z,direct.z, 1e-15);
+}
diff --git a/source/module_base/test/realarray_test.cpp b/source/module_base/test/realarray_test.cpp
new file mode 100644
index 00000000000..a6c861c2fb4
--- /dev/null
+++ b/source/module_base/test/realarray_test.cpp
@@ -0,0 +1,184 @@
+#include "../realarray.h"
+#include "gtest/gtest.h"
+
+/************************************************
+ *  unit test of class realArray
+ ***********************************************/
+
+/**
+ * - Tested Functions:
+ *   - GetArrayCount
+ *     - get the total number of real array created
+ *   - Construct
+ *     - construct real arrays (3 or 4 dimensions)
+ *   - Creat 
+ *     - create a real array (3 or 4 dimensions)
+ *   - GetSize
+ *     - get the total size of a real array
+ *   - GetDim
+ *     - get the dimension of a real array
+ *   - ZeroOut
+ *     - set all elements of a real array to zero
+ *   - GetBound
+ *     - get the size of each dimension of a real array
+ *   - ArrayEqReal
+ *     - set all value of an array to a double float 
+ *   - ArrayEqArray
+ *     - equal a realarray to another one
+ *   - Parentheses
+ *     - access element by using operator"()"
+ *   - ConstParentheses
+ *     - access element by using "()" through pointer
+ *     - without changing its elements
+ */
+
+class realArrayTest : public testing::Test
+{
+protected:
+	ModuleBase::realArray a3, a4, b3, b4;
+	double aa = 11.0;
+	double bb = 1.0;
+	int count0;
+	int count1;
+	const double zero = 0.0;
+};
+
+TEST_F(realArrayTest,GetArrayCount)
+{
+	count0 = ModuleBase::realArray::getArrayCount();
+	ModuleBase::realArray c3, c4;
+	count1 = ModuleBase::realArray::getArrayCount();
+	EXPECT_EQ((count1-count0),2);
+}
+
+TEST_F(realArrayTest,Construct)
+{
+	ModuleBase::realArray x3(1,5,3);
+	ModuleBase::realArray xp3(x3);
+	ModuleBase::realArray x4(1,7,3,4);
+	ModuleBase::realArray xp4(x4);
+	EXPECT_EQ(x3.getSize(),15);
+	EXPECT_EQ(xp3.getSize(),15);
+	EXPECT_EQ(x4.getSize(),84);
+	EXPECT_EQ(xp4.getSize(),84);
+}
+
+TEST_F(realArrayTest,Create)
+{
+	a3.create(1,2,3);
+	a4.create(1,2,3,4);
+	EXPECT_EQ(a3.getSize(),6);
+	EXPECT_EQ(a4.getSize(),24);
+}
+
+TEST_F(realArrayTest,GetSize)
+{
+	ModuleBase::realArray a3(1,5,3);
+	//std::cout<< &a3 << &(this->a3) <<std::endl;
+	ModuleBase::realArray a4(1,7,3,4);
+	EXPECT_EQ(a3.getSize(),15);
+	EXPECT_EQ(a4.getSize(),84);
+}
+
+TEST_F(realArrayTest,GetDim)
+{
+	ModuleBase::realArray a3(1,5,3);
+	ModuleBase::realArray a4(1,7,3,4);
+	EXPECT_EQ(a3.getDim(),3);
+	EXPECT_EQ(a4.getDim(),4);
+}
+
+TEST_F(realArrayTest,ZeroOut)
+{
+	a3.create(1,2,3);
+	a4.create(1,2,3,4);
+	a3.zero_out();
+	a4.zero_out();
+	for (int i=0;i<a3.getSize();++i)
+	{
+		EXPECT_EQ(a3.ptr[i],zero);
+    	}
+	for (int i=0;i<a4.getSize();++i)
+	{
+		EXPECT_EQ(a4.ptr[i],zero);
+    	}
+}
+
+TEST_F(realArrayTest,GetBound)
+{
+	ModuleBase::realArray a3(1,5,3);
+	ModuleBase::realArray a4(1,7,3,4);
+	EXPECT_EQ(a3.getBound1(),1);
+	EXPECT_EQ(a3.getBound2(),5);
+	EXPECT_EQ(a3.getBound3(),3);
+	EXPECT_EQ(a4.getBound1(),1);
+	EXPECT_EQ(a4.getBound2(),7);
+	EXPECT_EQ(a4.getBound3(),3);
+	EXPECT_EQ(a4.getBound4(),4);
+}
+
+TEST_F(realArrayTest,ArrayEqReal)
+{
+	ModuleBase::realArray a3(1,5,3);
+	ModuleBase::realArray a4(1,7,3,4);
+	a3 = aa;
+	a4 = bb;
+	for (int i=0;i<a3.getSize();++i)
+	{
+		EXPECT_EQ(a3.ptr[i],aa);
+    	}
+	for (int i=0;i<a4.getSize();++i)
+	{
+		EXPECT_EQ(a4.ptr[i],bb);
+    	}
+}
+
+TEST_F(realArrayTest,ArrayEqArray)
+{
+	ModuleBase::realArray a3(1,5,3);
+	ModuleBase::realArray a4(1,7,3,4);
+	a3 = aa;
+	a4 = bb;
+	b3.create(1,5,3);
+	b4.create(1,7,3,4);
+	b3 = a3;
+	b4 = a4;
+	for (int i=0;i<b3.getSize();++i)
+	{
+		EXPECT_EQ(b3.ptr[i],aa);
+    	}
+	for (int i=0;i<b4.getSize();++i)
+	{
+		EXPECT_EQ(b4.ptr[i],bb);
+    	}
+}
+
+TEST_F(realArrayTest,Parentheses)
+{
+	ModuleBase::realArray a3(1,5,3);
+	ModuleBase::realArray a4(1,7,3,4);
+	a3 = aa;
+	a4 = bb;
+	EXPECT_EQ(a3(0,0,0),aa);
+	EXPECT_EQ(a3(0,4,2),aa);
+	EXPECT_NE(a3(0,0,0),bb);
+	EXPECT_NE(a3(0,4,2),bb);
+	EXPECT_EQ(a4(0,0,0,0),bb);
+	EXPECT_EQ(a4(0,6,2,3),bb);
+	EXPECT_NE(a4(0,0,0,0),aa);
+	EXPECT_NE(a4(0,6,2,3),aa);
+}
+
+TEST_F(realArrayTest,ConstParenthesesExpert)
+{
+	ModuleBase::realArray a3(1,5,3);
+	a3 = aa;
+	const ModuleBase::realArray b3 (a3);
+	EXPECT_EQ(b3(0,0,0),aa);
+	EXPECT_EQ(b3(0,4,2),aa);
+	ModuleBase::realArray a4(1,5,3,2);
+	a4 = bb;
+	const ModuleBase::realArray b4 (a4);
+	EXPECT_EQ(b4(0,1,2,1),bb);
+	EXPECT_EQ(b4(0,0,0,0),bb);
+}
diff --git a/source/module_base/test/vector3_test.cpp b/source/module_base/test/vector3_test.cpp
new file mode 100644
index 00000000000..e806802c702
--- /dev/null
+++ b/source/module_base/test/vector3_test.cpp
@@ -0,0 +1,722 @@
+#include "../vector3.h"
+#include "gtest/gtest.h"
+#include "gmock/gmock.h"
+
+/************************************************
+ *  unit test of class Vector3
+ ***********************************************/
+
+/**
+ * - Tested Functions:
+ *   - Construct
+ *     - two ways of constructing a 3d vector
+ *   - Set
+ *     - set a 3d vector
+ *   - Equal
+ *     - overload operator "=" for 3d vector
+ *   - PlusEqual
+ *     - overload operator "+=" for 3d vector
+ *   - MinusEqual
+ *     - overload operator "-=" for 3d vector
+ *   - MultiplyEqual
+ *     - overload operator "*=" for (3d vector) * scalar
+ *   - OverEqual
+ *     - overload operator "/=" for (3d vector) / scalar
+ *   - Negative
+ *     - overload operator "-" to get - Vector3
+ *   - Reverse
+ *     - same as negative
+ *   - Access
+ *     - access elements by using "[]"
+ *   - ConstAccess
+ *     - access elements by using "[]" through pinters
+ *     - withough chaning element values
+ *   - VectorPlus
+ *     - overload operator "+"  for two 3d vectors
+ *   - VectorMinus
+ *     - overload operator "-"  for two 3d vectors
+ *   - Norm2
+ *     - get the square of norm of a 3d vector
+ *   - Norm
+ *     - get the norm of a 3d vector
+ *   - Normalize
+ *     - normalize a 3d vector
+ *   - VmultiplyV
+ *     - overload operator "*" to calculate
+ *     - the dot product of two 3d vectors
+ *   - VdotV
+ *     - dot product of two 3d vectors
+ *   - VmultiplyNum
+ *     - overload operator "*" to calculate
+ *     - the product of a 3d vector with a scalar
+ *     - of the product of a scalar with a 3d vector
+ *   - VoverNum
+ *     - overload operator "/" to calculate
+ *     -  a 3d vector over a scalar
+ *   - OperatorCaret
+ *     - overload operator "^" to calculate
+ *     - the cross product of two 3d vectors
+ *   - VeqV
+ *     - overload operator "==" to assert
+ *     - the equality between two 3d vectors
+ *   - VneV
+ *     - overload operator "!=" to assert
+ *     - the inequality between two 3d vectors
+ *   - StdOutV
+ *     - overload operator "<<" to print out
+ *     - a 3d vectors on standard output
+ *   - PrintV
+ *     - print a 3d vectors on standard output
+ *     - with formats
+ */
+
+class Vector3Test : public testing::Test
+{
+protected:
+	double da = 3.0;
+	double db = 4.0;
+	double dc = 5.0;
+	int    ia = 3;
+	int    ib = 4;
+	int    ic = 5;
+	float  fa = 3.0;
+	float  fb = 4.0;
+	float  fc = 5.0;
+	// for capturing stdout
+	std::string output;
+};
+
+TEST_F(Vector3Test,Construct)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u (da,db,dc);
+	ModuleBase::Vector3<double> up (u);
+	EXPECT_EQ(u.x,3.0);
+	EXPECT_EQ(u.y,4.0);
+	EXPECT_EQ(u.z,5.0);
+	EXPECT_EQ(up.x,3.0);
+	EXPECT_EQ(up.y,4.0);
+	EXPECT_EQ(up.z,5.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v (fa,fb,fc);
+	ModuleBase::Vector3<float> vp (v);
+	EXPECT_EQ(v.x,3.0);
+	EXPECT_EQ(v.y,4.0);
+	EXPECT_EQ(v.z,5.0);
+	EXPECT_EQ(vp.x,3.0);
+	EXPECT_EQ(vp.y,4.0);
+	EXPECT_EQ(vp.z,5.0);
+	// int Vector3
+	ModuleBase::Vector3<int> w (ia,ib,ic);
+	ModuleBase::Vector3<int> wp (w);
+	EXPECT_EQ(w.x,3);
+	EXPECT_EQ(w.y,4);
+	EXPECT_EQ(w.z,5);
+	EXPECT_EQ(wp.x,3);
+	EXPECT_EQ(wp.y,4);
+	EXPECT_EQ(wp.z,5);
+}
+
+TEST_F(Vector3Test,Set)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u;
+	u.set(da,db,dc);
+	EXPECT_EQ(u.x,3.0);
+	EXPECT_EQ(u.y,4.0);
+	EXPECT_EQ(u.z,5.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v;
+	v.set(fa,fb,fc);
+	EXPECT_EQ(v.x,3.0);
+	EXPECT_EQ(v.y,4.0);
+	EXPECT_EQ(v.z,5.0);
+	// int Vector3
+	ModuleBase::Vector3<int> w;
+	w.set(ia,ib,ic);
+	EXPECT_EQ(w.x,3);
+	EXPECT_EQ(w.y,4);
+	EXPECT_EQ(w.z,5);
+}
+
+TEST_F(Vector3Test,Equal)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u, up;
+	u.set(da,db,dc);
+	up = u;
+	EXPECT_EQ(up.x,3.0);
+	EXPECT_EQ(up.y,4.0);
+	EXPECT_EQ(up.z,5.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v, vp;
+	v.set(fa,fb,fc);
+	vp = v;
+	EXPECT_EQ(vp.x,3.0);
+	EXPECT_EQ(vp.y,4.0);
+	EXPECT_EQ(vp.z,5.0);
+	// int Vector3
+	ModuleBase::Vector3<int> w, wp;
+	w.set(ia,ib,ic);
+	wp = w;
+	EXPECT_EQ(wp.x,3);
+	EXPECT_EQ(wp.y,4);
+	EXPECT_EQ(wp.z,5);
+}
+
+TEST_F(Vector3Test,PlusEqual)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u, up;
+	u.set(da,db,dc);
+	up.set(da,db,dc);
+	up += u;
+	EXPECT_EQ(up.x,6.0);
+	EXPECT_EQ(up.y,8.0);
+	EXPECT_EQ(up.z,10.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v, vp;
+	v.set(fa,fb,fc);
+	vp.set(fa,fb,fc);
+	vp += v;
+	EXPECT_EQ(vp.x,6.0);
+	EXPECT_EQ(vp.y,8.0);
+	EXPECT_EQ(vp.z,10.0);
+	// int Vector3
+	ModuleBase::Vector3<int> w, wp;
+	w.set(ia,ib,ic);
+	wp.set(ia,ib,ic);
+	wp += w;
+	EXPECT_EQ(wp.x,6);
+	EXPECT_EQ(wp.y,8);
+	EXPECT_EQ(wp.z,10);
+}
+
+TEST_F(Vector3Test,MinusEqual)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u, up;
+	u.set(da,db,dc);
+	up.set(3*da,3*db,3*dc);
+	up -= u;
+	EXPECT_EQ(up.x,6.0);
+	EXPECT_EQ(up.y,8.0);
+	EXPECT_EQ(up.z,10.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v, vp;
+	v.set(fa,fb,fc);
+	vp.set(3*fa,3*fb,3*fc);
+	vp -= v;
+	EXPECT_EQ(vp.x,6.0);
+	EXPECT_EQ(vp.y,8.0);
+	EXPECT_EQ(vp.z,10.0);
+	// int Vector3
+	ModuleBase::Vector3<int> w, wp;
+	w.set(ia,ib,ic);
+	wp.set(3*ia,3*ib,3*ic);
+	wp -= w;
+	EXPECT_EQ(wp.x,6);
+	EXPECT_EQ(wp.y,8);
+	EXPECT_EQ(wp.z,10);
+}
+
+TEST_F(Vector3Test,MultiplyEqual)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u;
+	u.set(da,db,dc);
+	u *= 2;
+	EXPECT_EQ(u.x,6.0);
+	EXPECT_EQ(u.y,8.0);
+	EXPECT_EQ(u.z,10.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v;
+	v.set(fa,fb,fc);
+	v *= 2;
+	EXPECT_EQ(v.x,6.0);
+	EXPECT_EQ(v.y,8.0);
+	EXPECT_EQ(v.z,10.0);
+	// int Vector3
+	ModuleBase::Vector3<int> w;
+	w.set(ia,ib,ic);
+	w *= 2;
+	EXPECT_EQ(w.x,6);
+	EXPECT_EQ(w.y,8);
+	EXPECT_EQ(w.z,10);
+}
+
+TEST_F(Vector3Test,OverEqual)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u;
+	u.set(4*da,4*db,4*dc);
+	u /= 2;
+	EXPECT_EQ(u.x,6.0);
+	EXPECT_EQ(u.y,8.0);
+	EXPECT_EQ(u.z,10.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v;
+	v.set(4*fa,4*fb,4*fc);
+	v /= 2;
+	EXPECT_EQ(v.x,6.0);
+	EXPECT_EQ(v.y,8.0);
+	EXPECT_EQ(v.z,10.0);
+	// int Vector3
+	ModuleBase::Vector3<int> w;
+	w.set(4*ia,4*ib,4*ic);
+	w /= 2;
+	EXPECT_EQ(w.x,6);
+	EXPECT_EQ(w.y,8);
+	EXPECT_EQ(w.z,10);
+}
+
+TEST_F(Vector3Test,Negative)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u, up;
+	u.set(da,db,dc);
+	up = -u;
+	EXPECT_EQ(up.x,-3.0);
+	EXPECT_EQ(up.y,-4.0);
+	EXPECT_EQ(up.z,-5.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v, vp;
+	v.set(fa,fb,fc);
+	vp = -v;
+	EXPECT_EQ(vp.x,-3.0);
+	EXPECT_EQ(vp.y,-4.0);
+	EXPECT_EQ(vp.z,-5.0);
+	// int Vector3
+	ModuleBase::Vector3<int> w, wp;
+	w.set(ia,ib,ic);
+	wp = -w;
+	EXPECT_EQ(wp.x,-3);
+	EXPECT_EQ(wp.y,-4);
+	EXPECT_EQ(wp.z,-5);
+}
+
+TEST_F(Vector3Test,Access)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u;
+	u.set(da,db,dc);
+	EXPECT_EQ(u[0],3.0);
+	EXPECT_EQ(u[1],4.0);
+	EXPECT_EQ(u[2],5.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v;
+	v.set(fa,fb,fc);
+	EXPECT_EQ(v.x,3.0);
+	EXPECT_EQ(v.y,4.0);
+	EXPECT_EQ(v.z,5.0);
+	// int Vector3
+	ModuleBase::Vector3<int> w;
+	w.set(ia,ib,ic);
+	EXPECT_EQ(w.x,3);
+	EXPECT_EQ(w.y,4);
+	EXPECT_EQ(w.z,5);
+}
+
+
+TEST_F(Vector3Test,ConstAccess)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u;
+	u.set(da,db,dc);
+	const ModuleBase::Vector3<double> *up(&u);
+	EXPECT_EQ((*up)[0],3.0);
+	EXPECT_EQ((*up)[1],4.0);
+	EXPECT_EQ((*up)[2],5.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v;
+	const ModuleBase::Vector3<float> *vp(&v);
+	v.set(fa,fb,fc);
+	EXPECT_EQ((*vp).x,3.0);
+	EXPECT_EQ((*vp).y,4.0);
+	EXPECT_EQ((*vp).z,5.0);
+	// int Vector3
+	//ModuleBase::Vector3<int> w;
+	//w.set(ia,ib,ic);
+	//EXPECT_EQ(w.x,3);
+	//EXPECT_EQ(w.y,4);
+	//EXPECT_EQ(w.z,5);
+}
+
+
+TEST_F(Vector3Test,Reverse)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u;
+	u.set(da,db,dc);
+	u.reverse();
+	EXPECT_EQ(u.x,-3.0);
+	EXPECT_EQ(u.y,-4.0);
+	EXPECT_EQ(u.z,-5.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v;
+	v.set(fa,fb,fc);
+	v.reverse();
+	EXPECT_EQ(v.x,-3.0);
+	EXPECT_EQ(v.y,-4.0);
+	EXPECT_EQ(v.z,-5.0);
+	// int Vector3
+	ModuleBase::Vector3<int> w;
+	w.set(ia,ib,ic);
+	w.reverse();
+	EXPECT_EQ(w.x,-3);
+	EXPECT_EQ(w.y,-4);
+	EXPECT_EQ(w.z,-5);
+}
+
+TEST_F(Vector3Test,VectorPlus)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u,up,upp;
+	u.set(da,db,dc);
+	up.set(da,db,dc);
+	upp = u + up;
+	EXPECT_EQ(upp[0],6.0);
+	EXPECT_EQ(upp[1],8.0);
+	EXPECT_EQ(upp[2],10.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v,vp,vpp;
+	v.set(fa,fb,fc);
+	vp.set(fa,fb,fc);
+	vpp = v + vp;
+	EXPECT_EQ(vpp.x,6.0);
+	EXPECT_EQ(vpp.y,8.0);
+	EXPECT_EQ(vpp.z,10.0);
+	// int Vector3
+	ModuleBase::Vector3<int> w,wp,wpp;
+	w.set(ia,ib,ic);
+	wp.set(ia,ib,ic);
+	wpp = w + wp;
+	EXPECT_EQ(wpp.x,6);
+	EXPECT_EQ(wpp.y,8);
+	EXPECT_EQ(wpp.z,10);
+}
+
+TEST_F(Vector3Test,VectorMinus)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u,up,upp;
+	u.set(da,db,dc);
+	up.set(2*da,2*db,2*dc);
+	upp = u - up;
+	EXPECT_EQ(upp[0],-3.0);
+	EXPECT_EQ(upp[1],-4.0);
+	EXPECT_EQ(upp[2],-5.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v,vp,vpp;
+	v.set(fa,fb,fc);
+	vp.set(3*fa,3*fb,3*fc);
+	vpp = v - vp;
+	EXPECT_EQ(vpp.x,-6.0);
+	EXPECT_EQ(vpp.y,-8.0);
+	EXPECT_EQ(vpp.z,-10.0);
+	// int Vector3
+	ModuleBase::Vector3<int> w,wp,wpp;
+	w.set(3*ia,3*ib,3*ic);
+	wp.set(ia,ib,ic);
+	wpp = w - wp;
+	EXPECT_EQ(wpp.x,6);
+	EXPECT_EQ(wpp.y,8);
+	EXPECT_EQ(wpp.z,10);
+}
+
+TEST_F(Vector3Test,Norm2)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u;
+	u.set(da,db,dc);
+	EXPECT_EQ(u.norm2(),50.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v;
+	v.set(fa,fb,fc);
+	EXPECT_EQ(v.norm2(),50.0);
+	// int Vector3
+	ModuleBase::Vector3<int> w;
+	w.set(ia,ib,ic);
+	EXPECT_EQ(w.norm2(),50);
+}
+
+
+TEST_F(Vector3Test,Norm)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u;
+	u.set(da,db,dc);
+	double nm = u.norm();
+	double nm2= sqrt(50.0);
+	EXPECT_DOUBLE_EQ(nm,nm2);
+	EXPECT_FLOAT_EQ(nm,sqrt(50.0));
+	// float Vector3
+	ModuleBase::Vector3<float> v;
+	v.set(fa,fb,fc);
+	float nmp = v.norm();
+	float nmp2= sqrt(50.0);
+	EXPECT_FLOAT_EQ(nmp,sqrt(50.0));
+}
+
+
+TEST_F(Vector3Test,Normalize)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u;
+	u.set(da,db,dc);
+	u.normalize();
+	EXPECT_DOUBLE_EQ(u.norm(),1.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v;
+	v.set(fa,fb,fc);
+	v.normalize();
+	EXPECT_FLOAT_EQ(v.norm(),1.0);
+}
+
+TEST_F(Vector3Test,VmultiplyV)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u,up;
+	u.set(da,db,dc);
+	up.set(da,db,dc);
+	double mpd = u * up;
+	EXPECT_EQ(mpd,50.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v,vp;
+	v.set(fa,fb,fc);
+	vp.set(fa,fb,fc);
+	float mpf = v*vp;
+	EXPECT_EQ(mpf,50.0);
+	// int Vector3
+	ModuleBase::Vector3<int> w,wp;
+	w.set(ia,ib,ic);
+	wp.set(ia,ib,ic);
+	int mpi = w*wp;
+	EXPECT_EQ(mpf,50);
+}
+
+TEST_F(Vector3Test,VdotV)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u,up;
+	u.set(da,db,dc);
+	up.set(da,db,dc);
+	double mpd = dot(u,up);
+	EXPECT_EQ(mpd,50.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v,vp;
+	v.set(fa,fb,fc);
+	vp.set(fa,fb,fc);
+	float mpf = dot(v,vp);
+	EXPECT_EQ(mpf,50.0);
+	// int Vector3
+	ModuleBase::Vector3<int> w,wp;
+	w.set(ia,ib,ic);
+	wp.set(ia,ib,ic);
+	int mpi = dot(w,wp);
+	EXPECT_EQ(mpf,50);
+}
+
+TEST_F(Vector3Test,VmultiplyNum)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u,up,upp;
+	u.set(da,db,dc);
+	double s = 3.0;
+	up = s*u; upp = u*s;
+	EXPECT_EQ(upp[0],up[0]);
+	EXPECT_EQ(upp[1],up[1]);
+	EXPECT_EQ(upp[2],up[2]);
+	EXPECT_EQ(upp[0],9.0);
+	EXPECT_EQ(upp[1],12.0);
+	EXPECT_EQ(upp[2],15.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v,vp,vpp;
+	v.set(fa,fb,fc);
+	float t = 3.0;
+	vp = t*v; vpp = v*t;
+	EXPECT_EQ(vpp[0],vp[0]);
+	EXPECT_EQ(vpp[1],vp[1]);
+	EXPECT_EQ(vpp[2],vp[2]);
+	EXPECT_EQ(vpp[0],9.0);
+	EXPECT_EQ(vpp[1],12.0);
+	EXPECT_EQ(vpp[2],15.0);
+	// int Vector3
+	ModuleBase::Vector3<int> w,wp,wpp;
+	w.set(ia,ib,ic);
+	int q = 3;
+	wp = q*w; wpp = w*q;
+	EXPECT_EQ(wpp[0],wp[0]);
+	EXPECT_EQ(wpp[1],wp[1]);
+	EXPECT_EQ(wpp[2],wp[2]);
+	EXPECT_EQ(wpp[0],9.0);
+	EXPECT_EQ(wpp[1],12.0);
+	EXPECT_EQ(wpp[2],15.0);
+}
+
+TEST_F(Vector3Test,VoverNum)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u,up;
+	u.set(2*da,2*db,2*dc);
+	double s = 2.0;
+	up = u/s;
+	EXPECT_EQ(up.x,3.0);
+	EXPECT_EQ(up.y,4.0);
+	EXPECT_EQ(up.z,5.0);
+	// float Vector3
+	ModuleBase::Vector3<float> v,vp;
+	v.set(2*fa,2*fb,2*fc);
+	float t = 2.0;
+	vp = v/t;
+	EXPECT_EQ(vp.x,3.0);
+	EXPECT_EQ(vp.y,4.0);
+	EXPECT_EQ(vp.z,5.0);
+	// int Vector3
+	ModuleBase::Vector3<int> w,wp;
+	w.set(2*ia,2*ib,2*ic);
+	int q = 2;
+	wp = w/q;
+	EXPECT_EQ(wp.x,3);
+	EXPECT_EQ(wp.y,4);
+	EXPECT_EQ(wp.z,5);
+}
+
+TEST_F(Vector3Test,OperatorCaret)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u,up,upp;
+	u.set(da,db,dc);
+	up.set(da,db,dc);
+	upp = u^up;
+	EXPECT_EQ(upp.x,u.y*up.z - u.z*up.y);
+	EXPECT_EQ(upp.y,u.z*up.x - u.x*up.z);
+	EXPECT_EQ(upp.z,u.x*up.y - u.y*up.x);
+	// float Vector3
+	ModuleBase::Vector3<float> v,vp,vpp;
+	v.set(2*fa,2*fb,2*fc);
+	vp.set(fa,fb,fc);
+	vpp = v^vp;
+	EXPECT_EQ(vpp.x,v.y*vp.z - v.z*vp.y);
+	EXPECT_EQ(vpp.y,v.z*vp.x - v.x*vp.z);
+	EXPECT_EQ(vpp.z,v.x*vp.y - v.y*vp.x);
+	// int Vector3
+	ModuleBase::Vector3<int> w,wp,wpp;
+	w.set(2*ia,2*ib,2*ic);
+	wp.set(ia,ib,ic);
+	wpp = w^wp;
+	EXPECT_EQ(wpp.x,w.y*wp.z - w.z*wp.y);
+	EXPECT_EQ(wpp.y,w.z*wp.x - w.x*wp.z);
+	EXPECT_EQ(wpp.z,w.x*wp.y - w.y*wp.x);
+}
+
+TEST_F(Vector3Test,Cross)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u,up,upp;
+	u.set(da,db,dc);
+	up.set(da,db,dc);
+	upp = cross(u,up);
+	EXPECT_EQ(upp.x,u.y*up.z - u.z*up.y);
+	EXPECT_EQ(upp.y,u.z*up.x - u.x*up.z);
+	EXPECT_EQ(upp.z,u.x*up.y - u.y*up.x);
+	// float Vector3
+	ModuleBase::Vector3<float> v,vp,vpp;
+	v.set(2*fa,2*fb,2*fc);
+	vp.set(fa,fb,fc);
+	vpp = cross(v,vp);
+	EXPECT_EQ(vpp.x,v.y*vp.z - v.z*vp.y);
+	EXPECT_EQ(vpp.y,v.z*vp.x - v.x*vp.z);
+	EXPECT_EQ(vpp.z,v.x*vp.y - v.y*vp.x);
+	// int Vector3
+	ModuleBase::Vector3<int> w,wp,wpp;
+	w.set(2*ia,2*ib,2*ic);
+	wp.set(ia,ib,ic);
+	wpp = cross(w,wp);
+	EXPECT_EQ(wpp.x,w.y*wp.z - w.z*wp.y);
+	EXPECT_EQ(wpp.y,w.z*wp.x - w.x*wp.z);
+	EXPECT_EQ(wpp.z,w.x*wp.y - w.y*wp.x);
+}
+
+TEST_F(Vector3Test,VeqV)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u,up;
+	u.set(da,db,dc);
+	up.set(da,db,dc);
+	EXPECT_TRUE(up == u);
+	// float Vector3
+	ModuleBase::Vector3<float> v,vp;
+	v.set(fa,fb,fc);
+	vp.set(fa,fb,fc);
+	EXPECT_TRUE(vp == v);
+	// int Vector3
+	ModuleBase::Vector3<int> w,wp;
+	w.set(ia,ib,ic);
+	wp.set(ia,ib,ic);
+	EXPECT_TRUE(wp == w);
+}
+
+TEST_F(Vector3Test,VneV)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u,up;
+	u.set(da,db,dc);
+	up.set(da,db,2*dc);
+	EXPECT_TRUE(up != u);
+	// float Vector3
+	ModuleBase::Vector3<float> v,vp;
+	v.set(fa,fb,2*fc);
+	vp.set(fa,fb,fc);
+	EXPECT_TRUE(vp != v);
+	// int Vector3
+	ModuleBase::Vector3<int> w,wp;
+	w.set(ia,ib,2*ic);
+	wp.set(ia,ib,ic);
+	EXPECT_TRUE(wp != w);
+}
+
+TEST_F(Vector3Test,StdOutV)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u(da,db,dc);
+	testing::internal::CaptureStdout();
+	std::cout << u << std::endl;
+	output = testing::internal::GetCapturedStdout();
+	EXPECT_THAT(output,testing::HasSubstr("("));
+	// float Vector3
+	ModuleBase::Vector3<float> v(fa,fb,fc);
+	testing::internal::CaptureStdout();
+	std::cout << v << std::endl;
+	output = testing::internal::GetCapturedStdout();
+	EXPECT_THAT(output,testing::HasSubstr(","));
+	// int Vector3
+	ModuleBase::Vector3<int> w(ia,ib,ic);
+	testing::internal::CaptureStdout();
+	std::cout << w << std::endl;
+	output = testing::internal::GetCapturedStdout();
+	EXPECT_THAT(output,testing::HasSubstr(")"));
+}
+
+TEST_F(Vector3Test,PrintV)
+{
+	// double Vector3
+	ModuleBase::Vector3<double> u(3.1415926,db,dc);
+	testing::internal::CaptureStdout();
+	u.print();
+	output = testing::internal::GetCapturedStdout();
+	EXPECT_THAT(output,testing::HasSubstr("3.1416"));
+	// float Vector3
+	ModuleBase::Vector3<float> v(fa,fb,3.14);
+	testing::internal::CaptureStdout();
+	v.print();
+	output = testing::internal::GetCapturedStdout();
+	EXPECT_THAT(output,testing::HasSubstr("3.14"));
+	// int Vector3
+	ModuleBase::Vector3<int> w(ia,101,ic);
+	testing::internal::CaptureStdout();
+	w.print();
+	output = testing::internal::GetCapturedStdout();
+	EXPECT_THAT(output,testing::HasSubstr("101"));
+}
+
diff --git a/source/module_base/vector3.h b/source/module_base/vector3.h
index f230f616f3c..94a09afcfdb 100644
--- a/source/module_base/vector3.h
+++ b/source/module_base/vector3.h
@@ -6,103 +6,372 @@
 #endif
 
 #include <cmath>
-#include <iostream>
 #include <iomanip>
+#include <iostream>
 namespace ModuleBase
 {
 
-template <class T>
-class Vector3
+/**
+ * @brief 3 elements vector
+ *
+ * @tparam T
+ */
+template <class T> class Vector3
 {
-public:
-	T x;
-	T y;
-	T z;
-	
-	Vector3(const T &x1 = 0,const T &y1 = 0,const T &z1 = 0) :x(x1),y(y1),z(z1){};
-	Vector3(const Vector3<T> &v) :x(v.x),y(v.y),z(v.z){};		// Peize Lin add 2018-07-16
-	void set(const T &x1, const T &y1,const T &z1) { x = x1; y = y1; z = z1; }
-
-	Vector3<T>& operator =(const Vector3<T> &u) { x=u.x; y=u.y; z=u.z;     return *this; }
-	Vector3<T>& operator+=(const Vector3<T> &u) { x+=u.x; y+=u.y; z+=u.z;  return *this; }
-	Vector3<T>& operator-=(const Vector3<T> &u) { x-=u.x; y-=u.y; z-=u.z;  return *this; }
-	Vector3<T>& operator*=(const Vector3<T> &u); 
-	Vector3<T>& operator*=(const T &s)          { x*=s; y*=s; z*=s;        return *this; }
-	Vector3<T>& operator/=(const Vector3<T> &u);
-	Vector3<T>& operator/=(const T &s)          { x/=s; y/=s; z/=s;        return *this; }
-	Vector3<T>  operator -() const              { return Vector3<T>(-x,-y,-z); }			// Peize Lin add 2017-01-10
-
-	T operator[](int index)const { return (&x)[index]; }
-	T& operator[](int index) { return (&x)[index]; }
-
-	T norm2(void) const	{ return x*x + y*y + z*z; }
-	T norm(void) const	{ return sqrt(norm2()); }
-	Vector3<T>& normalize(void){ const T m=norm(); x/=m; y/=m; z/=m; return *this; }	// Peize Lin update return 2019-09-08
-	Vector3<T>& reverse(void){ x=-x; y=-y; z=-z; return *this; }						// Peize Lin update return 2019-09-08
-
-	void print(void)const ;		// mohan add 2009-11-29
+  public:
+    T x;
+    T y;
+    T z;
+
+    /**
+     * @brief Construct a new Vector 3 object
+     *
+     * @param x1
+     * @param y1
+     * @param z1
+     */
+    Vector3(const T &x1 = 0, const T &y1 = 0, const T &z1 = 0) : x(x1), y(y1), z(z1){};
+    Vector3(const Vector3<T> &v) : x(v.x), y(v.y), z(v.z){}; // Peize Lin add 2018-07-16
+
+    /**
+     * @brief set a 3d vector
+     *
+     * @param x1
+     * @param y1
+     * @param z1
+     */
+    void set(const T &x1, const T &y1, const T &z1)
+    {
+        x = x1;
+        y = y1;
+        z = z1;
+    }
+
+    /**
+     * @brief Overload operator "=" for Vector3
+     *
+     * @param u
+     * @return Vector3<T>&
+     */
+    Vector3<T> &operator=(const Vector3<T> &u)
+    {
+        x = u.x;
+        y = u.y;
+        z = u.z;
+        return *this;
+    }
+
+    /**
+     * @brief Overload operator "+=" for Vector3
+     *
+     * @param u
+     * @return Vector3<T>&
+     */
+    Vector3<T> &operator+=(const Vector3<T> &u)
+    {
+        x += u.x;
+        y += u.y;
+        z += u.z;
+        return *this;
+    }
+
+    /**
+     * @brief Overload operator "-=" for Vector3
+     *
+     * @param u
+     * @return Vector3<T>&
+     */
+    Vector3<T> &operator-=(const Vector3<T> &u)
+    {
+        x -= u.x;
+        y -= u.y;
+        z -= u.z;
+        return *this;
+    }
+
+    /**
+     * @brief Overload operator "*=" for (Vector3)*scalar
+     *
+     * @param s
+     * @return Vector3<T>&
+     */
+    Vector3<T> &operator*=(const T &s)
+    {
+        x *= s;
+        y *= s;
+        z *= s;
+        return *this;
+    }
+
+    /**
+     * @brief Overload operator "/=" for (Vector3)/scalar
+     *
+     * @param s
+     * @return Vector3<T>&
+     */
+    Vector3<T> &operator/=(const T &s)
+    {
+        x /= s;
+        y /= s;
+        z /= s;
+        return *this;
+    }
+
+    /**
+     * @brief Overload operator "-" to get (-Vector3)
+     *
+     * @return Vector3<T>
+     */
+    Vector3<T> operator-() const
+    {
+        return Vector3<T>(-x, -y, -z);
+    } // Peize Lin add 2017-01-10
+
+    /**
+     * @brief Over load "[]" for accessing elements with pointers
+     *
+     * @param index
+     * @return T
+     */
+    T operator[](int index) const
+    {
+        return (&x)[index];
+    }
+
+    /**
+     * @brief Overload operator "[]" for accesing elements
+     *
+     * @param index
+     * @return T&
+     */
+    T &operator[](int index)
+    {
+        return (&x)[index];
+    }
+
+    /**
+     * @brief Get the square of nomr of a Vector3
+     *
+     * @return T
+     */
+    T norm2(void) const
+    {
+        return x * x + y * y + z * z;
+    }
+
+    /**
+     * @brief Get the norm of a Vector3
+     *
+     * @return T
+     */
+    T norm(void) const
+    {
+        return sqrt(norm2());
+    }
+
+    /**
+     * @brief Normalize a Vector3
+     *
+     * @return Vector3<T>&
+     */
+    Vector3<T> &normalize(void)
+    {
+        const T m = norm();
+        x /= m;
+        y /= m;
+        z /= m;
+        return *this;
+    } // Peize Lin update return 2019-09-08
+
+    /**
+     * @brief Get (-Vector3)
+     *
+     * @return Vector3<T>&
+     */
+    Vector3<T> &reverse(void)
+    {
+        x = -x;
+        y = -y;
+        z = -z;
+        return *this;
+    } // Peize Lin update return 2019-09-08
+
+    /**
+     * @brief Print a Vector3 on standard output
+     * with formats
+     *
+     */
+    void print(void) const; // mohan add 2009-11-29
 };
 
-template <class T> inline Vector3<T> operator+( const Vector3<T> &u, const Vector3<T> &v ) { return Vector3<T>( u.x+v.x, u.y+v.y, u.z+v.z ); }
-template <class T> inline Vector3<T> operator-( const Vector3<T> &u, const Vector3<T> &v ) { return Vector3<T>( u.x-v.x, u.y-v.y, u.z-v.z ); }
-//u.v=(ux*vx)+(uy*vy)+(uz*vz)                                                     
-template <class T> inline T          operator*( const Vector3<T> &u, const Vector3<T> &v ) { return ( u.x*v.x + u.y*v.y + u.z*v.z ); }
-template <class T> inline Vector3<T> operator*( const T &s,          const Vector3<T> &u ) { return Vector3<T>( u.x*s, u.y*s, u.z*s ); }
-template <class T> inline Vector3<T> operator*( const Vector3<T> &u, const T &s          ) { return Vector3<T>( u.x*s, u.y*s, u.z*s ); } // mohan add 2009-5-10
-template <class T> inline Vector3<T> operator/( const Vector3<T> &u, const T &s          ) { return Vector3<T>( u.x/s, u.y/s, u.z/s ); }
-//u.v=(ux*vx)+(uy*vy)+(uz*vz)
-template <class T> inline T          dot      ( const Vector3<T> &u, const Vector3<T> &v ) { return ( u.x*v.x + u.y*v.y + u.z*v.z ); }
-// | i  j  k  |
-// | ux uy uz |
-// | vx vy vz |
-// u.v=(uy*vz-uz*vy)i+(-ux*vz+uz*vx)j+(ux*vy-uy*vx)k
-template <class T> inline Vector3<T> operator^(const Vector3<T> &u,const Vector3<T> &v)
-{	
-	return Vector3<T> ( u.y * v.z - u.z * v.y,
-	                   -u.x * v.z + u.z * v.x,
-	                    u.x * v.y - u.y * v.x);
+/**
+ * @brief Overload "+" for two Vector3
+ *
+ * @param[in] u
+ * @param[in] v
+ * @return Vector3<T>
+ */
+template <class T> inline Vector3<T> operator+(const Vector3<T> &u, const Vector3<T> &v)
+{
+    return Vector3<T>(u.x + v.x, u.y + v.y, u.z + v.z);
 }
-// | i  j  k  |
-// | ux uy uz |
-// | vx vy vz |
-// u.v=(uy*vz-uz*vy)i+(-ux*vz+uz*vx)j+(ux*vy-uy*vzx)k
-template <class T> inline Vector3<T> cross(const Vector3<T> &u,const Vector3<T> &v)
+
+/**
+ * @brief Overload "-" for two Vector3
+ *
+ * @param[in] u
+ * @param[in] v
+ * @return Vector3<T>
+ */
+template <class T> inline Vector3<T> operator-(const Vector3<T> &u, const Vector3<T> &v)
+{
+    return Vector3<T>(u.x - v.x, u.y - v.y, u.z - v.z);
+}
+
+/**
+ * @brief Overload "*" to calculate the dot product
+ * of two Vector3
+ *
+ * @param u
+ * @param v
+ * @return template <class T>
+ */
+template <class T> inline T operator*(const Vector3<T> &u, const Vector3<T> &v)
 {
-	return Vector3<T> ( u.y * v.z - u.z * v.y,
-	                   -u.x * v.z + u.z * v.x,
-	                    u.x * v.y - u.y * v.x); 
+    return (u.x * v.x + u.y * v.y + u.z * v.z);
 }
-//s = u.(v x w)
-//template <class T> T TripleScalarProduct(Vector3<T> u, Vector3<T> v, Vector3<T> w)
+
+/**
+ * @brief Overload "*" to calculate (Vector3)*scalar
+ *
+ * @param[in] s
+ * @param[in] u
+ * @return Vector3<T>
+ */
+template <class T> inline Vector3<T> operator*(const T &s, const Vector3<T> &u)
+{
+    return Vector3<T>(u.x * s, u.y * s, u.z * s);
+}
+
+/**
+ * @brief Overload "*" to calculate scalar*(Vector3)
+ *
+ * @param u
+ * @param s
+ * @return Vector3<T>
+ */
+template <class T> inline Vector3<T> operator*(const Vector3<T> &u, const T &s)
+{
+    return Vector3<T>(u.x * s, u.y * s, u.z * s);
+} // mohan add 2009-5-10
+
+/**
+ * @brief Overload "/" to calculate Vector3/scalar
+ *
+ * @tparam T
+ * @param u
+ * @param s
+ * @return Vector3<T>
+ */
+template <class T> inline Vector3<T> operator/(const Vector3<T> &u, const T &s)
+{
+    return Vector3<T>(u.x / s, u.y / s, u.z / s);
+}
+
+/**
+ * @brief Dot productor of two Vector3
+ *
+ * @param u
+ * @param v
+ * @return T
+ * @note u.v=(ux*vx)+(uy*vy)+(uz*vz)
+ */
+template <class T> inline T dot(const Vector3<T> &u, const Vector3<T> &v)
+{
+    return (u.x * v.x + u.y * v.y + u.z * v.z);
+}
+
+/**
+ * @brief Overload "^" for cross product of two Vector3
+ *
+ * @param u
+ * @param v
+ * @return template <class T>
+ * @note
+ * | i  j  k  |
+ * | ux uy uz |
+ * | vx vy vz |
+ * u.v=(uy*vz-uz*vy)i+(-ux*vz+uz*vx)j+(ux*vy-uy*vx)k
+ */
+template <class T> inline Vector3<T> operator^(const Vector3<T> &u, const Vector3<T> &v)
+{
+    return Vector3<T>(u.y * v.z - u.z * v.y, -u.x * v.z + u.z * v.x, u.x * v.y - u.y * v.x);
+}
+
+/**
+ * @brief Cross product of two Vector3
+ *
+ * @param u
+ * @param v
+ * @return template <class T>
+ * @note
+ * | i  j  k  |
+ * | ux uy uz |
+ * | vx vy vz |
+ * u.v=(uy*vz-uz*vy)i+(-ux*vz+uz*vx)j+(ux*vy-uy*vx)k
+ */
+template <class T> inline Vector3<T> cross(const Vector3<T> &u, const Vector3<T> &v)
+{
+    return Vector3<T>(u.y * v.z - u.z * v.y, -u.x * v.z + u.z * v.x, u.x * v.y - u.y * v.x);
+}
+// s = u.(v x w)
+// template <class T> T TripleScalarProduct(Vector3<T> u, Vector3<T> v, Vector3<T> w)
 //{
 //	return T((u.x * (v.y * w.z - v.z * w.y)) +
 //	         (u.y * (-v.x * w.z + v.z * w.x)) +
 //	         (u.z * (v.x * w.y - v.y * w.x)));
-//}
+// }
 
-//whether m1 != m2
-template <class T> inline bool operator !=(const Vector3<T> &u, const Vector3<T> &v){ return !(u == v); }
-//whether u == v
-template <class T> inline bool operator ==(const Vector3<T> &u, const Vector3<T> &v)
+// whether m1 != m2
+template <class T> inline bool operator!=(const Vector3<T> &u, const Vector3<T> &v)
 {
-	if(u.x == v.x && u.y == v.y && u.z == v.z)
-		return true;
-	return false;
+    return !(u == v);
 }
-
-template <class T> void Vector3<T>::print(void)const
+// whether u == v
+template <class T> inline bool operator==(const Vector3<T> &u, const Vector3<T> &v)
 {
-	std::cout.precision(5) ;
-	std::cout << "(" << std::setw(10) << x << "," << std::setw(10) << y << ","
-	<< std::setw(10) << z  << ")"  << std::endl ;
-	return ;
+    if (u.x == v.x && u.y == v.y && u.z == v.z)
+        return true;
+    return false;
 }
-template<class T> static std::ostream & operator << ( std::ostream &os, const Vector3<T> &u )
+
+/**
+ * @brief Print a Vector3 on standard output with formats
+ *
+ */
+template <class T> void Vector3<T>::print(void) const
 {
-	os << "(" << std::setw(10) << u.x << "," << std::setw(10) << u.y << "," << std::setw(10) << u.z  << ")";
-	return os;
+    std::cout.precision(5);
+    std::cout << "(" << std::setw(10) << x << "," << std::setw(10) << y << "," << std::setw(10) << z << ")"
+              << std::endl;
+    return;
 }
 
+/**
+ * @brief Overload "<<" tor print out a
+ * Vector3 on standard output
+ *
+ * @tparam T
+ * @param[in] os
+ * @param[in] u
+ * @return std::ostream&
+ */
+template <class T> static std::ostream &operator<<(std::ostream &os, const Vector3<T> &u)
+{
+    os << "(" << std::setw(10) << u.x << "," << std::setw(10) << u.y << "," << std::setw(10) << u.z << ")";
+    return os;
 }
 
+} // namespace ModuleBase
+
 #endif