diff --git a/tools/01_NAO_generation/SIAB/src_tools/Gram_Schmidt_Orth.h b/tools/01_NAO_generation/SIAB/src_tools/Gram_Schmidt_Orth.h
index 272874c5b3c..7971d8e657c 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/Gram_Schmidt_Orth.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/Gram_Schmidt_Orth.h
@@ -1,156 +1,199 @@
 #ifndef GRAM_SCHMIDT_ORTH_H
 #define GRAM_SCHMIDT_ORTH_H
 
-#include<cstddef>
-#include<exception>
+#include <cstddef>
+#include <exception>
 using std::exception;
 
-class GS_Orth_Coordinate: public exception{};
+class GS_Orth_Coordinate : public exception
+{
+};
 
-template<typename Func_Type, typename R_Type=double>
+template <typename Func_Type, typename R_Type = double>
 class Gram_Schmidt_Orth
 {
-public:
-
-	~Gram_Schmidt_Orth();
-	Func_Type **orth();
-	
-private:
-	int coordinate;			// 0:cartesian, 1:sphere
-	size_t func_num;
-	size_t variable_num;
-	Func_Type **func{nullptr};
-	R_Type *r{nullptr};
-	R_Type *r_2{nullptr};	
-	R_Type *rab{nullptr};
-
-private:
-	bool finish_init{false};
-	
-public:
-
-	size_t get_coordinate() const { return coordinate; }
-	Gram_Schmidt_Orth &set_coordinate(size_t coordinate_in)
-	{
-		if( 0 != coordinate_in && 1 != coordinate_in) { throw GS_Orth_Coordinate();}
-		coordinate = coordinate_in; 
-		return *this; 
-	}
-
-	size_t get_func_num() const { return func_num; }
-	Gram_Schmidt_Orth &set_func_num(size_t func_num_in){ func_num = func_num_in; return *this; }
-
-	size_t get_variable_num() const { return variable_num; }
-	Gram_Schmidt_Orth &set_variable_num(size_t variable_num_in){ variable_num = variable_num_in; return *this; }
-
-	Func_Type **get_func() const { return func; }
-	Gram_Schmidt_Orth &set_func(Func_Type **func_in){ func = func_in; return *this; }
-
-	R_Type *get_r() const { return r; }
-	Gram_Schmidt_Orth &set_r(R_Type *r_in){ r = r_in; return *this; }
-
-	R_Type *get_r_2() const { return r_2; }
-
-	R_Type *get_rab() const { return rab; }
-	Gram_Schmidt_Orth &set_rab_element(const R_Type &rab_element_in)
-	{
-		if(NULL==rab)
-			rab = new R_Type [variable_num];
-		for( size_t i(0); i<variable_num; ++i) 
-			rab[i] = rab_element_in;
-		return *this; 
-	}
-		
-private:
-	void init();
-	Func_Type cal_norm(Func_Type *func_in);
-
+  public:
+    ~Gram_Schmidt_Orth();
+    Func_Type** orth();
+
+  private:
+    int coordinate; // 0:cartesian, 1:sphere
+    size_t func_num;
+    size_t variable_num;
+    Func_Type** func{nullptr};
+    R_Type* r{nullptr};
+    R_Type* r_2{nullptr};
+    R_Type* rab{nullptr};
+
+  private:
+    bool finish_init{false};
+
+  public:
+    size_t get_coordinate() const
+    {
+        return coordinate;
+    }
+    Gram_Schmidt_Orth& set_coordinate(size_t coordinate_in)
+    {
+        if (0 != coordinate_in && 1 != coordinate_in)
+        {
+            throw GS_Orth_Coordinate();
+        }
+        coordinate = coordinate_in;
+        return *this;
+    }
+
+    size_t get_func_num() const
+    {
+        return func_num;
+    }
+    Gram_Schmidt_Orth& set_func_num(size_t func_num_in)
+    {
+        func_num = func_num_in;
+        return *this;
+    }
+
+    size_t get_variable_num() const
+    {
+        return variable_num;
+    }
+    Gram_Schmidt_Orth& set_variable_num(size_t variable_num_in)
+    {
+        variable_num = variable_num_in;
+        return *this;
+    }
+
+    Func_Type** get_func() const
+    {
+        return func;
+    }
+    Gram_Schmidt_Orth& set_func(Func_Type** func_in)
+    {
+        func = func_in;
+        return *this;
+    }
+
+    R_Type* get_r() const
+    {
+        return r;
+    }
+    Gram_Schmidt_Orth& set_r(R_Type* r_in)
+    {
+        r = r_in;
+        return *this;
+    }
+
+    R_Type* get_r_2() const
+    {
+        return r_2;
+    }
+
+    R_Type* get_rab() const
+    {
+        return rab;
+    }
+    Gram_Schmidt_Orth& set_rab_element(const R_Type& rab_element_in)
+    {
+        if (NULL == rab)
+            rab = new R_Type[variable_num];
+        for (size_t i(0); i < variable_num; ++i)
+            rab[i] = rab_element_in;
+        return *this;
+    }
+
+  private:
+    void init();
+    Func_Type cal_norm(Func_Type* func_in);
 };
 
-#include "mathzone.h"
 #include "lapack_connector.h"
+#include "mathzone.h"
 
-template<typename Func_Type, typename R_Type>
-void Gram_Schmidt_Orth<Func_Type,R_Type>::init()
+template <typename Func_Type, typename R_Type>
+void Gram_Schmidt_Orth<Func_Type, R_Type>::init()
 {
-	if(this->finish_init) return;
-	if(nullptr == this->r_2)
-	{
-		this->r_2 = new R_Type[variable_num];
-		Mathzone::Pointwise_Product( this->variable_num, this->r, this->r, this->r_2);		
-	}
-	if(nullptr == rab)
-	{
-		this->rab = new R_Type[variable_num];
-		this->rab[0] = 0.0;
-		for( size_t ir(1); ir<variable_num; ++ir)
-		{
-			this->rab[ir] = this->r[ir] - this->r[ir-1];
-		}
-	}
-	this->finish_init = true;
-	return;
+    if (this->finish_init)
+        return;
+    if (nullptr == this->r_2)
+    {
+        this->r_2 = new R_Type[variable_num];
+        Mathzone::Pointwise_Product(this->variable_num, this->r, this->r, this->r_2);
+    }
+    if (nullptr == rab)
+    {
+        this->rab = new R_Type[variable_num];
+        this->rab[0] = 0.0;
+        for (size_t ir(1); ir < variable_num; ++ir)
+        {
+            this->rab[ir] = this->r[ir] - this->r[ir - 1];
+        }
+    }
+    this->finish_init = true;
+    return;
 }
 
-template<typename Func_Type, typename R_Type>
-Gram_Schmidt_Orth<Func_Type,R_Type>::~Gram_Schmidt_Orth()
+template <typename Func_Type, typename R_Type>
+Gram_Schmidt_Orth<Func_Type, R_Type>::~Gram_Schmidt_Orth()
 {
-	delete[]rab;
-	delete[]r_2;
+    delete[] rab;
+    delete[] r_2;
 }
 
-template<typename Func_Type, typename R_Type>
-Func_Type **Gram_Schmidt_Orth<Func_Type,R_Type>::orth()
+template <typename Func_Type, typename R_Type>
+Func_Type** Gram_Schmidt_Orth<Func_Type, R_Type>::orth()
 {
-	init();
-
-	//	Schmidt: hn to en
-	//		e1 = h1 / ||h1||
-	//		gn = hn - \sum{i=1 to n-1}(hn,ei)ei
-	//		en = gn / ||gn||
-
-	for( size_t if1(0); if1<func_num; ++if1)
-	{
-		for( size_t if2(0); if2<if1; ++if2)
-		{
-			// (hn,ei)
-			Func_Type *mul_func = new Func_Type[this->variable_num];
-			Mathzone::Pointwise_Product( this->variable_num, this->func[if1], this->func[if2], mul_func);		
-			const Func_Type in_product = cal_norm(mul_func);
-			delete[]mul_func; mul_func=nullptr;
-
-			// hn -(hn,ei)ei
-			LapackConnector::daxpy( this->variable_num, -in_product, this->func[if2], 1, this->func[if1], 1);
-		}
-		
-		// ||gn||
-		Func_Type *func_2 = new Func_Type[this->variable_num];
-		Mathzone::Pointwise_Product( this->variable_num, this->func[if1], this->func[if1], func_2);
-		const Func_Type norm = sqrt(cal_norm(func_2));
-		delete[]func_2; func_2=nullptr;
-		
-		// en = gn / ||gn||
-		LapackConnector::dscal( this->variable_num, 1.0/norm, this->func[if1], 1);
-	}
-	return func;
+    init();
+
+    //	Schmidt: hn to en
+    //		e1 = h1 / ||h1||
+    //		gn = hn - \sum{i=1 to n-1}(hn,ei)ei
+    //		en = gn / ||gn||
+
+    for (size_t if1(0); if1 < func_num; ++if1)
+    {
+        for (size_t if2(0); if2 < if1; ++if2)
+        {
+            // (hn,ei)
+            Func_Type* mul_func = new Func_Type[this->variable_num];
+            Mathzone::Pointwise_Product(this->variable_num, this->func[if1], this->func[if2], mul_func);
+            const Func_Type in_product = cal_norm(mul_func);
+            delete[] mul_func;
+            mul_func = nullptr;
+
+            // hn -(hn,ei)ei
+            LapackConnector::daxpy(this->variable_num, -in_product, this->func[if2], 1, this->func[if1], 1);
+        }
+
+        // ||gn||
+        Func_Type* func_2 = new Func_Type[this->variable_num];
+        Mathzone::Pointwise_Product(this->variable_num, this->func[if1], this->func[if1], func_2);
+        const Func_Type norm = sqrt(cal_norm(func_2));
+        delete[] func_2;
+        func_2 = nullptr;
+
+        // en = gn / ||gn||
+        LapackConnector::dscal(this->variable_num, 1.0 / norm, this->func[if1], 1);
+    }
+    return func;
 }
 
-template<typename Func_Type, typename R_Type>
-Func_Type Gram_Schmidt_Orth<Func_Type,R_Type>::cal_norm(Func_Type *func_in)
+template <typename Func_Type, typename R_Type>
+Func_Type Gram_Schmidt_Orth<Func_Type, R_Type>::cal_norm(Func_Type* func_in)
 {
-	Func_Type norm(0.0);
-	switch(coordinate)
-	{
-		case 0:
-			Mathzone::Simpson_Integral( this->variable_num, func_in, this->rab, norm);		
-			break;
-		case 1:	
-			Func_Type *tmp_func = new Func_Type[this->variable_num];
-			Mathzone::Pointwise_Product( this->variable_num, func_in, this->r_2, tmp_func);
-			Mathzone::Simpson_Integral( this->variable_num, tmp_func, this->rab, norm);	
-			delete[]tmp_func; tmp_func=nullptr;
-			break;
-	}
-	return norm;
+    Func_Type norm(0.0);
+    switch (coordinate)
+    {
+    case 0:
+        Mathzone::Simpson_Integral(this->variable_num, func_in, this->rab, norm);
+        break;
+    case 1:
+        Func_Type* tmp_func = new Func_Type[this->variable_num];
+        Mathzone::Pointwise_Product(this->variable_num, func_in, this->r_2, tmp_func);
+        Mathzone::Simpson_Integral(this->variable_num, tmp_func, this->rab, norm);
+        delete[] tmp_func;
+        tmp_func = nullptr;
+        break;
+    }
+    return norm;
 }
 #endif
diff --git a/tools/01_NAO_generation/SIAB/src_tools/Gram_Schmidt_Orth_bak.cpp b/tools/01_NAO_generation/SIAB/src_tools/Gram_Schmidt_Orth_bak.cpp
index a1a58a47f90..285d6c93077 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/Gram_Schmidt_Orth_bak.cpp
+++ b/tools/01_NAO_generation/SIAB/src_tools/Gram_Schmidt_Orth_bak.cpp
@@ -1,83 +1,87 @@
 #include "Gram_Schmidt_Orth.h"
-#include "mathzone.h"
 #include "lapack_connector.h"
+#include "mathzone.h"
 
 void Gram_Schmidt_Orth::init()
 {
-	if(this->finish_init) return;
-	if(nullptr == this->r_2)
-	{
-		this->r_2 = new double[variable_num];
-		Mathzone::Pointwise_Product( this->variable_num, this->r, this->r, this->r_2);		
-	}
-	if(nullptr == rab)
-	{
-		this->rab = new double[variable_num];
-		this->rab[0] = 0.0;
-		for( size_t ir(1); ir<variable_num; ++ir)
-		{
-			this->rab[ir] = this->r[ir] - this->r[ir-1];
-		}
-	}
-	this->finish_init = true;
-	return;
+    if (this->finish_init)
+        return;
+    if (nullptr == this->r_2)
+    {
+        this->r_2 = new double[variable_num];
+        Mathzone::Pointwise_Product(this->variable_num, this->r, this->r, this->r_2);
+    }
+    if (nullptr == rab)
+    {
+        this->rab = new double[variable_num];
+        this->rab[0] = 0.0;
+        for (size_t ir(1); ir < variable_num; ++ir)
+        {
+            this->rab[ir] = this->r[ir] - this->r[ir - 1];
+        }
+    }
+    this->finish_init = true;
+    return;
 }
 
 Gram_Schmidt_Orth::~Gram_Schmidt_Orth()
 {
-	delete[]rab;
-	delete[]r_2;
+    delete[] rab;
+    delete[] r_2;
 }
 
-double **Gram_Schmidt_Orth::orth()
+double** Gram_Schmidt_Orth::orth()
 {
-	init();
+    init();
+
+    //	Schmidt: hn to en
+    //		e1 = h1 / ||h1||
+    //		gn = hn - \sum{i=1 to n-1}(hn,ei)ei
+    //		en = gn / ||gn||
+
+    for (size_t if1(0); if1 < func_num; ++if1)
+    {
+        for (size_t if2(0); if2 < if1; ++if2)
+        {
+            // (hn,ei)
+            double* mul_func = new double[this->variable_num];
+            Mathzone::Pointwise_Product(this->variable_num, this->func[if1], this->func[if2], mul_func);
+            const double in_product = cal_norm(mul_func);
+            delete[] mul_func;
+            mul_func = nullptr;
 
-	//	Schmidt: hn to en
-	//		e1 = h1 / ||h1||
-	//		gn = hn - \sum{i=1 to n-1}(hn,ei)ei
-	//		en = gn / ||gn||
+            // hn -(hn,ei)ei
+            daxpy(this->variable_num, -in_product, this->func[if2], 1, this->func[if1], 1);
+        }
 
-	for( size_t if1(0); if1<func_num; ++if1)
-	{
-		for( size_t if2(0); if2<if1; ++if2)
-		{
-			// (hn,ei)
-			double *mul_func = new double[this->variable_num];
-			Mathzone::Pointwise_Product( this->variable_num, this->func[if1], this->func[if2], mul_func);		
-			const double in_product = cal_norm(mul_func);
-			delete[]mul_func; mul_func=nullptr;
+        // ||gn||
+        double* func_2 = new double[this->variable_num];
+        Mathzone::Pointwise_Product(this->variable_num, this->func[if1], this->func[if1], func_2);
+        const double norm = sqrt(cal_norm(func_2));
+        delete[] func_2;
+        func_2 = nullptr;
 
-			// hn -(hn,ei)ei
-			daxpy( this->variable_num, -in_product, this->func[if2], 1, this->func[if1], 1);
-		}
-		
-		// ||gn||
-		double *func_2 = new double[this->variable_num];
-		Mathzone::Pointwise_Product( this->variable_num, this->func[if1], this->func[if1], func_2);
-		const double norm = sqrt(cal_norm(func_2));
-		delete[]func_2; func_2=nullptr;
-		
-		// en = gn / ||gn||
-		dscal( this->variable_num, 1.0/norm, this->func[if1], 1);
-	}
-	return func;
+        // en = gn / ||gn||
+        dscal(this->variable_num, 1.0 / norm, this->func[if1], 1);
+    }
+    return func;
 }
 
-double Gram_Schmidt_Orth::cal_norm(double *func_in)
+double Gram_Schmidt_Orth::cal_norm(double* func_in)
 {
-	double norm(0.0);
-	switch(coordinate)
-	{
-		case 0:
-			Mathzone::Simpson_Integral( this->variable_num, func_in, this->rab, norm);		
-			break;
-		case 1:	
-			double *tmp_func = new double[this->variable_num];
-			Mathzone::Pointwise_Product( this->variable_num, func_in, this->r_2, tmp_func);
-			Mathzone::Simpson_Integral( this->variable_num, tmp_func, this->rab, norm);	
-			delete[]tmp_func; tmp_func=nullptr;
-			break;
-	}
-	return norm;
+    double norm(0.0);
+    switch (coordinate)
+    {
+    case 0:
+        Mathzone::Simpson_Integral(this->variable_num, func_in, this->rab, norm);
+        break;
+    case 1:
+        double* tmp_func = new double[this->variable_num];
+        Mathzone::Pointwise_Product(this->variable_num, func_in, this->r_2, tmp_func);
+        Mathzone::Simpson_Integral(this->variable_num, tmp_func, this->rab, norm);
+        delete[] tmp_func;
+        tmp_func = nullptr;
+        break;
+    }
+    return norm;
 }
\ No newline at end of file
diff --git a/tools/01_NAO_generation/SIAB/src_tools/Gram_Schmidt_Orth_bak.h b/tools/01_NAO_generation/SIAB/src_tools/Gram_Schmidt_Orth_bak.h
index 18551a8c82c..fbb157b0090 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/Gram_Schmidt_Orth_bak.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/Gram_Schmidt_Orth_bak.h
@@ -1,69 +1,108 @@
 #ifndef GRAM_SCHMIDT_ORTH_H
 #define GRAM_SCHMIDT_ORTH_H
 
-#include<cstddef>
-#include<exception>
+#include <cstddef>
+#include <exception>
 using std::exception;
 
-class GS_Orth_Coordinate: public exception{};
+class GS_Orth_Coordinate : public exception
+{
+};
 
 class Gram_Schmidt_Orth
 {
-public:
+  public:
+    ~Gram_Schmidt_Orth();
+    double** orth();
 
-	~Gram_Schmidt_Orth();
-	double **orth();
-	
-private:
-	int coordinate;			// 0:cartesian, 1:sphere
-	size_t func_num;
-	size_t variable_num;
-	double **func{nullptr};
-	double *r{nullptr};
-	double *r_2{nullptr};	
-	double *rab{nullptr};
+  private:
+    int coordinate; // 0:cartesian, 1:sphere
+    size_t func_num;
+    size_t variable_num;
+    double** func{nullptr};
+    double* r{nullptr};
+    double* r_2{nullptr};
+    double* rab{nullptr};
 
-private:
-	bool finish_init{false};
-	
-public:
+  private:
+    bool finish_init{false};
 
-	size_t get_coordinate() const { return coordinate; }
-	Gram_Schmidt_Orth &set_coordinate(size_t coordinate_in)
-	{
-		if( 0 != coordinate_in && 1 != coordinate_in) { throw GS_Orth_Coordinate();}
-		coordinate = coordinate_in; 
-		return *this; 
-	}
+  public:
+    size_t get_coordinate() const
+    {
+        return coordinate;
+    }
+    Gram_Schmidt_Orth& set_coordinate(size_t coordinate_in)
+    {
+        if (0 != coordinate_in && 1 != coordinate_in)
+        {
+            throw GS_Orth_Coordinate();
+        }
+        coordinate = coordinate_in;
+        return *this;
+    }
 
-	size_t get_func_num() const { return func_num; }
-	Gram_Schmidt_Orth &set_func_num(size_t func_num_in){ func_num = func_num_in; return *this; }
+    size_t get_func_num() const
+    {
+        return func_num;
+    }
+    Gram_Schmidt_Orth& set_func_num(size_t func_num_in)
+    {
+        func_num = func_num_in;
+        return *this;
+    }
 
-	size_t get_variable_num() const { return variable_num; }
-	Gram_Schmidt_Orth &set_variable_num(size_t variable_num_in){ variable_num = variable_num_in; return *this; }
+    size_t get_variable_num() const
+    {
+        return variable_num;
+    }
+    Gram_Schmidt_Orth& set_variable_num(size_t variable_num_in)
+    {
+        variable_num = variable_num_in;
+        return *this;
+    }
 
-	double **get_func() const { return func; }
-	Gram_Schmidt_Orth &set_func(double **func_in){ func = func_in; return *this; }
+    double** get_func() const
+    {
+        return func;
+    }
+    Gram_Schmidt_Orth& set_func(double** func_in)
+    {
+        func = func_in;
+        return *this;
+    }
 
-	double *get_r() const { return r; }
-	Gram_Schmidt_Orth &set_r(double *r_in){ r = r_in; return *this; }
+    double* get_r() const
+    {
+        return r;
+    }
+    Gram_Schmidt_Orth& set_r(double* r_in)
+    {
+        r = r_in;
+        return *this;
+    }
 
-	double *get_r_2() const { return r_2; }
+    double* get_r_2() const
+    {
+        return r_2;
+    }
 
-	double *get_rab() const { return rab; }
-	Gram_Schmidt_Orth &set_rab_element(const double &rab_element_in)
-	{
-		if(NULL==rab)
-			rab = new double [variable_num];
-		for( size_t i(0); i<variable_num; ++i) 
-			rab[i] = rab_element_in;
-		return *this; 
-	}
-		
-private:
-	void init();
-	double cal_norm(double *func_in);
+    double* get_rab() const
+    {
+        return rab;
+    }
+    Gram_Schmidt_Orth& set_rab_element(const double& rab_element_in)
+    {
+        if (NULL == rab)
+            rab = new double[variable_num];
+        for (size_t i(0); i < variable_num; ++i)
+            rab[i] = rab_element_in;
+        return *this;
+    }
 
+  private:
+    void init();
+    double cal_norm(double* func_in);
 };
 
 #endif
\ No newline at end of file
diff --git a/tools/01_NAO_generation/SIAB/src_tools/Random.h b/tools/01_NAO_generation/SIAB/src_tools/Random.h
index 32d43e4f637..37df2fde2e0 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/Random.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/Random.h
@@ -1,46 +1,47 @@
 #ifndef RANDOM_H
 #define RANDOM_H
-#include <cstdlib>
 #include <cassert>
+#include <cstdlib>
 
 class Random
 {
-	public:
-	Random();
-	~Random();
-
-	static void between0and1( double *v, const int &num )
-	{
-		assert( v!= NULL);
-		assert( num > 1);
-		for(int i=0; i<num; i++)
-		{
-			v[i] = static_cast<double>( std::rand() ) / RAND_MAX; 
-		}
-	}
-
-	static double betweenMinus2and2(void)
-	{
-		return 2.0*betweenMinus1and1();
-	}
-
-	static double betweenMinus1and1(void)
-	{
-		int a = std::rand() % 2;
-		if(a==0) return between0and1();
-		else if(a==1) return betweenMinus1and0();
-	}
-
-	static double between0and1(void)
-	{
-		return static_cast<double>( std::rand() )/RAND_MAX;
-	}
-
-	static double betweenMinus1and0(void)
-	{
-		return -static_cast<double>( std::rand() )/RAND_MAX;
-	}
-
+  public:
+    Random();
+    ~Random();
+
+    static void between0and1(double* v, const int& num)
+    {
+        assert(v != NULL);
+        assert(num > 1);
+        for (int i = 0; i < num; i++)
+        {
+            v[i] = static_cast<double>(std::rand()) / RAND_MAX;
+        }
+    }
+
+    static double betweenMinus2and2(void)
+    {
+        return 2.0 * betweenMinus1and1();
+    }
+
+    static double betweenMinus1and1(void)
+    {
+        int a = std::rand() % 2;
+        if (a == 0)
+            return between0and1();
+        else if (a == 1)
+            return betweenMinus1and0();
+    }
+
+    static double between0and1(void)
+    {
+        return static_cast<double>(std::rand()) / RAND_MAX;
+    }
+
+    static double betweenMinus1and0(void)
+    {
+        return -static_cast<double>(std::rand()) / RAND_MAX;
+    }
 };
 
 #endif
diff --git a/tools/01_NAO_generation/SIAB/src_tools/Simulated_Annealing.cpp b/tools/01_NAO_generation/SIAB/src_tools/Simulated_Annealing.cpp
index 14d48ea16cc..94671c769e1 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/Simulated_Annealing.cpp
+++ b/tools/01_NAO_generation/SIAB/src_tools/Simulated_Annealing.cpp
@@ -1,89 +1,137 @@
 #include "Simulated_Annealing.h"
-#include "mathzone.h"
+
 #include "Random.h"
+#include "mathzone.h"
 
 void Simulated_Annealing::cal_all_temperature(void)
 {
-	preprocess_all_temperature();
-	for( size_t itemperature(0); itemperature != this->temperature_num; ++itemperature)
-	{
-		cal_each_temperature(itemperature);
-		if(quit_loop_all_temperature()){ break; }
-		update_temperature(itemperature);
-	}
-	reprocess_all_temperature();
-	return;
+    preprocess_all_temperature();
+    for (size_t itemperature(0); itemperature != this->temperature_num; ++itemperature)
+    {
+        cal_each_temperature(itemperature);
+        if (quit_loop_all_temperature())
+        {
+            break;
+        }
+        update_temperature(itemperature);
+    }
+    reprocess_all_temperature();
+    return;
 }
 
-void Simulated_Annealing::cal_each_temperature( size_t itemperature)
+void Simulated_Annealing::cal_each_temperature(size_t itemperature)
 {
-	preprocess_each_temperature(itemperature);
-	for( size_t istep(0); istep != this->step_num; ++istep)
-	{
-		cal_each_step(itemperature, istep);
-		if(quit_loop_each_temperature(itemperature)){ break; }
-	}
-	reprocess_each_temperature(itemperature);
-	return;
+    preprocess_each_temperature(itemperature);
+    for (size_t istep(0); istep != this->step_num; ++istep)
+    {
+        cal_each_step(itemperature, istep);
+        if (quit_loop_each_temperature(itemperature))
+        {
+            break;
+        }
+    }
+    reprocess_each_temperature(itemperature);
+    return;
 }
 
-void Simulated_Annealing::cal_each_step( size_t itemperature, size_t istep)
+void Simulated_Annealing::cal_each_step(size_t itemperature, size_t istep)
 {
-	preprocess_each_step(itemperature, istep);
-	move_variable(itemperature, istep);
-	double new_function;
-	bool accept_flag(false), reject_flag(false);
-	judge_accept_reject(new_function, itemperature, istep, accept_flag, reject_flag);
-	process_accept_reject(new_function, itemperature, istep, accept_flag, reject_flag);
-	reprocess_each_step(itemperature, istep);
-	return;
+    preprocess_each_step(itemperature, istep);
+    move_variable(itemperature, istep);
+    double new_function;
+    bool accept_flag(false), reject_flag(false);
+    judge_accept_reject(new_function, itemperature, istep, accept_flag, reject_flag);
+    process_accept_reject(new_function, itemperature, istep, accept_flag, reject_flag);
+    reprocess_each_step(itemperature, istep);
+    return;
 }
 
-void Simulated_Annealing::judge_accept_reject( double &new_function, size_t itemperature, size_t istep, bool &accept_flag, bool &reject_flag)
+void Simulated_Annealing::judge_accept_reject(double& new_function,
+                                              size_t itemperature,
+                                              size_t istep,
+                                              bool& accept_flag,
+                                              bool& reject_flag)
 {
-	if(must_accept_before(new_function, itemperature, istep)) { accept_flag = true; }
-	else if(must_reject_before(new_function, itemperature, istep)) { reject_flag = true; }
-	if(accept_flag || reject_flag){ return; }
-	else
-	{
-		new_function = cal_new_function(itemperature, istep);
-		if(must_accept_after(new_function, itemperature, istep)) { accept_flag = true; }
-		else if(must_reject_after(new_function, itemperature, istep)) { reject_flag = true; }
-	}
+    if (must_accept_before(new_function, itemperature, istep))
+    {
+        accept_flag = true;
+    }
+    else if (must_reject_before(new_function, itemperature, istep))
+    {
+        reject_flag = true;
+    }
+    if (accept_flag || reject_flag)
+    {
+        return;
+    }
+    else
+    {
+        new_function = cal_new_function(itemperature, istep);
+        if (must_accept_after(new_function, itemperature, istep))
+        {
+            accept_flag = true;
+        }
+        else if (must_reject_after(new_function, itemperature, istep))
+        {
+            reject_flag = true;
+        }
+    }
 
-	if(accept_flag)
-	{
-		if(reject_flag) { throw Accept_Reject_Conflict(); }
-		else { return; }
-	}
-	else if(reject_flag) { return; }
-	else
-	{
-		if(new_function<=old_function) { accept_flag = true; }
-		else
-		{
-			const double p = Mathzone::Boltzmann_dist(new_function - old_function, this->temperature);
-			const double r = Random::between0and1();
-			if (r < p) { accept_flag = true; }
-			else { reject_flag = true; }		
-		}
-	}	
+    if (accept_flag)
+    {
+        if (reject_flag)
+        {
+            throw Accept_Reject_Conflict();
+        }
+        else
+        {
+            return;
+        }
+    }
+    else if (reject_flag)
+    {
+        return;
+    }
+    else
+    {
+        if (new_function <= old_function)
+        {
+            accept_flag = true;
+        }
+        else
+        {
+            const double p = Mathzone::Boltzmann_dist(new_function - old_function, this->temperature);
+            const double r = Random::between0and1();
+            if (r < p)
+            {
+                accept_flag = true;
+            }
+            else
+            {
+                reject_flag = true;
+            }
+        }
+    }
 }
 
-void Simulated_Annealing::process_accept_reject( double &new_function, size_t itemperature, size_t istep, bool &accept_flag, bool &reject_flag)
+void Simulated_Annealing::process_accept_reject(double& new_function,
+                                                size_t itemperature,
+                                                size_t istep,
+                                                bool& accept_flag,
+                                                bool& reject_flag)
 {
-	// accept and reject may conflict without checking before
-	if(accept_flag)
-	{
-		accept_process( new_function, itemperature, istep); 
-	}
-	if(reject_flag)
-	{
-		reject_process( new_function, itemperature, istep);		
-	}
+    // accept and reject may conflict without checking before
+    if (accept_flag)
+    {
+        accept_process(new_function, itemperature, istep);
+    }
+    if (reject_flag)
+    {
+        reject_process(new_function, itemperature, istep);
+    }
 }
 
-void Simulated_Annealing::accept_process( double new_function, size_t itemperature, size_t istep)
+void Simulated_Annealing::accept_process(double new_function, size_t itemperature, size_t istep)
 {
-	this->old_function = new_function;
+    this->old_function = new_function;
 }
diff --git a/tools/01_NAO_generation/SIAB/src_tools/Simulated_Annealing.h b/tools/01_NAO_generation/SIAB/src_tools/Simulated_Annealing.h
index 5e8c1f017a1..0e2bc6205dc 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/Simulated_Annealing.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/Simulated_Annealing.h
@@ -1,56 +1,115 @@
 #ifndef SIMULATED_ANNEALING_H
 #define SIMULATED_ANNEALING_H
 
-#include<cstddef>
-#include<exception>
+#include <cstddef>
+#include <exception>
 using std::exception;
 
-class Accept_Reject_Conflict: public exception{};
+class Accept_Reject_Conflict : public exception
+{
+};
 
 class Simulated_Annealing
 {
-	public:
-	
-	virtual void cal_all_temperature(void);
-	virtual void cal_each_temperature( size_t itemperature);
-	virtual void cal_each_step( size_t itemperature, size_t istep);
-	virtual void judge_accept_reject( double &new_function, size_t itemperature, size_t istep, bool &accept_flag, bool &reject_flag);
-	virtual void process_accept_reject( double &new_function, size_t itemperature, size_t istep, bool &accept_flag, bool &reject_flag);
-	virtual void accept_process( double new_function, size_t itemperature, size_t istep);
-
-	virtual void move_variable(size_t itemperature, size_t istep)=0;
-	virtual double cal_new_function(size_t itemperature, size_t istep)=0;
-	
-	virtual void preprocess_all_temperature(){}
-	virtual bool quit_loop_all_temperature(){return false;}
-	virtual void update_temperature(size_t itemperature){}
-	virtual void reprocess_all_temperature(){}
-	
-	virtual void preprocess_each_temperature(size_t itemperature){}
-	virtual bool quit_loop_each_temperature(size_t itemperature){return false;}
-	virtual void reprocess_each_temperature(size_t itemperature){}
-	
-	virtual void preprocess_each_step(size_t itemperature, size_t istep){}
-	virtual void reprocess_each_step(size_t itemperature, size_t iste){}
-	
-	virtual bool must_accept_before(double new_function, size_t itemperature, size_t istep){return false;}
-	virtual bool must_reject_before(double new_function, size_t itemperature, size_t istep){return false;}
-	virtual bool must_accept_after(double new_function, size_t itemperature, size_t istep){return false;}
-	virtual bool must_reject_after(double new_function, size_t itemperature, size_t istep){return false;}
-	
-	virtual void reject_process(double new_function, size_t itemperature, size_t istep){}
-	
-	protected:
-	size_t temperature_num;
-	size_t step_num;
-	double temperature;
-	double old_function;
-	
-	public:
-	Simulated_Annealing &set_temperature_num(size_t temperature_num_in){temperature_num = temperature_num_in; return *this;}
-	Simulated_Annealing &set_step_num(size_t step_num_in){step_num = step_num_in; return *this;}
-	Simulated_Annealing &set_temperature(double temperature_in){temperature = temperature_in; return *this;}
-	Simulated_Annealing &set_old_function(double old_function_in){old_function = old_function_in; return *this;}
+  public:
+    virtual void cal_all_temperature(void);
+    virtual void cal_each_temperature(size_t itemperature);
+    virtual void cal_each_step(size_t itemperature, size_t istep);
+    virtual void judge_accept_reject(double& new_function,
+                                     size_t itemperature,
+                                     size_t istep,
+                                     bool& accept_flag,
+                                     bool& reject_flag);
+    virtual void process_accept_reject(double& new_function,
+                                       size_t itemperature,
+                                       size_t istep,
+                                       bool& accept_flag,
+                                       bool& reject_flag);
+    virtual void accept_process(double new_function, size_t itemperature, size_t istep);
+
+    virtual void move_variable(size_t itemperature, size_t istep) = 0;
+    virtual double cal_new_function(size_t itemperature, size_t istep) = 0;
+
+    virtual void preprocess_all_temperature()
+    {
+    }
+    virtual bool quit_loop_all_temperature()
+    {
+        return false;
+    }
+    virtual void update_temperature(size_t itemperature)
+    {
+    }
+    virtual void reprocess_all_temperature()
+    {
+    }
+
+    virtual void preprocess_each_temperature(size_t itemperature)
+    {
+    }
+    virtual bool quit_loop_each_temperature(size_t itemperature)
+    {
+        return false;
+    }
+    virtual void reprocess_each_temperature(size_t itemperature)
+    {
+    }
+
+    virtual void preprocess_each_step(size_t itemperature, size_t istep)
+    {
+    }
+    virtual void reprocess_each_step(size_t itemperature, size_t iste)
+    {
+    }
+
+    virtual bool must_accept_before(double new_function, size_t itemperature, size_t istep)
+    {
+        return false;
+    }
+    virtual bool must_reject_before(double new_function, size_t itemperature, size_t istep)
+    {
+        return false;
+    }
+    virtual bool must_accept_after(double new_function, size_t itemperature, size_t istep)
+    {
+        return false;
+    }
+    virtual bool must_reject_after(double new_function, size_t itemperature, size_t istep)
+    {
+        return false;
+    }
+
+    virtual void reject_process(double new_function, size_t itemperature, size_t istep)
+    {
+    }
+
+  protected:
+    size_t temperature_num;
+    size_t step_num;
+    double temperature;
+    double old_function;
+
+  public:
+    Simulated_Annealing& set_temperature_num(size_t temperature_num_in)
+    {
+        temperature_num = temperature_num_in;
+        return *this;
+    }
+    Simulated_Annealing& set_step_num(size_t step_num_in)
+    {
+        step_num = step_num_in;
+        return *this;
+    }
+    Simulated_Annealing& set_temperature(double temperature_in)
+    {
+        temperature = temperature_in;
+        return *this;
+    }
+    Simulated_Annealing& set_old_function(double old_function_in)
+    {
+        old_function = old_function_in;
+        return *this;
+    }
 };
 
 #endif
\ No newline at end of file
diff --git a/tools/01_NAO_generation/SIAB/src_tools/complexarray.cpp b/tools/01_NAO_generation/SIAB/src_tools/complexarray.cpp
index 22519a92c40..e271ab0f784 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/complexarray.cpp
+++ b/tools/01_NAO_generation/SIAB/src_tools/complexarray.cpp
@@ -2,51 +2,52 @@
  * ESCP:Electro-Structure Calculate Package.
  ********************************************/
 
-#include <cstdlib>
 #include "complexarray.h"
 
+#include <cstdlib>
+
 int ComplexArray::arrayCount = 0;
 
 void ComplexArrayAlloc()
 {
-	cout << "\n Allocation error for ComplexArray " << endl;
-	exit(0);
+    cout << "\n Allocation error for ComplexArray " << endl;
+    exit(0);
 }
 
-ComplexArray::ComplexArray(const int d1,const int d2,const int d3)
+ComplexArray::ComplexArray(const int d1, const int d2, const int d3)
 {
-	dim = 3;
-	bound1 = (d1 <= 0) ? 1 : d1;
-	bound2 = (d2 <= 0) ? 1 : d2;
-	bound3 = (d3 <= 0) ? 1 : d3;
-	bound4 = 0;
+    dim = 3;
+    bound1 = (d1 <= 0) ? 1 : d1;
+    bound2 = (d2 <= 0) ? 1 : d2;
+    bound3 = (d3 <= 0) ? 1 : d3;
+    bound4 = 0;
 
-	set_new_handler(ComplexArrayAlloc);
+    set_new_handler(ComplexArrayAlloc);
 
-	size = bound1 * bound2 * bound3 ;	//* sizeof(float);
+    size = bound1 * bound2 * bound3; //* sizeof(float);
 
-	ptr = new std::complex<double>[size]();
-	assert(ptr != 0);
+    ptr = new std::complex<double>[size]();
+    assert(ptr != 0);
 
-	++arrayCount;
+    ++arrayCount;
 }
 
-ComplexArray::ComplexArray(const int d1,const int d2,const int d3,const int d4)
+ComplexArray::ComplexArray(const int d1, const int d2, const int d3, const int d4)
 {
-	dim = 4;
-	bound1 = (d1 <= 0) ? 1 : d1;
-	bound2 = (d2 <= 0) ? 1 : d2;
-	bound3 = (d3 <= 0) ? 1 : d3;
-	bound4 = (d4 <= 0) ? 1 : d4;
+    dim = 4;
+    bound1 = (d1 <= 0) ? 1 : d1;
+    bound2 = (d2 <= 0) ? 1 : d2;
+    bound3 = (d3 <= 0) ? 1 : d3;
+    bound4 = (d4 <= 0) ? 1 : d4;
 
-	set_new_handler(ComplexArrayAlloc);
+    set_new_handler(ComplexArrayAlloc);
 
-	size = bound1 * bound2 * bound3 * bound4 ;	//* sizeof(float);
+    size = bound1 * bound2 * bound3 * bound4; //* sizeof(float);
 
-	ptr = new std::complex<double>[size]();
-	assert(ptr != 0);
+    ptr = new std::complex<double>[size]();
+    assert(ptr != 0);
 
-	++arrayCount;
+    ++arrayCount;
 }
 
 //********************************
@@ -54,61 +55,63 @@ ComplexArray::ComplexArray(const int d1,const int d2,const int d3,const int d4)
 // Destructor for class ComplexArray
 //
 //********************************
-ComplexArray ::~ComplexArray()
+ComplexArray::~ComplexArray()
 {
     freemem();
 }
 
 void ComplexArray::freemem()
 {
-	delete [] ptr;
-	ptr = NULL;
+    delete[] ptr;
+    ptr = NULL;
 }
 
-void ComplexArray::create(const int d1,const int d2,const int d3,const int d4)
+void ComplexArray::create(const int d1, const int d2, const int d3, const int d4)
 {
-	size = d1 * d2 * d3 * d4;
-	assert(size>0);
+    size = d1 * d2 * d3 * d4;
+    assert(size > 0);
 
-	dim = 4;
+    dim = 4;
 
-	bound1 = d1;
-	bound2 = d2;
-	bound3 = d3;
-	bound4 = d4;
+    bound1 = d1;
+    bound2 = d2;
+    bound3 = d3;
+    bound4 = d4;
 
-	delete [] ptr;
-	ptr = new std::complex<double>[size]();
-	assert(ptr != 0);
+    delete[] ptr;
+    ptr = new std::complex<double>[size]();
+    assert(ptr != 0);
 }
 
-void ComplexArray::create(const int d1,const int d2,const int d3)
+void ComplexArray::create(const int d1, const int d2, const int d3)
 {
-	size = d1 * d2 * d3;
-	assert(size>0);
+    size = d1 * d2 * d3;
+    assert(size > 0);
 
-	dim = 3;
+    dim = 3;
 
-	bound1 = d1;
-	bound2 = d2;
-	bound3 = d3;
-	bound4 = 1;
+    bound1 = d1;
+    bound2 = d2;
+    bound3 = d3;
+    bound4 = 1;
 
-	delete [] ptr;
-	ptr = new std::complex<double>[size]();
-	assert(ptr != 0);
+    delete[] ptr;
+    ptr = new std::complex<double>[size]();
+    assert(ptr != 0);
 }
 
-const ComplexArray &ComplexArray::operator=(const ComplexArray &right)
+const ComplexArray& ComplexArray::operator=(const ComplexArray& right)
 {
-	for (int i = 0;i < size;i++) ptr[i] = right.ptr[i];
-	return *this;// enables x = y = z;
+    for (int i = 0; i < size; i++)
+        ptr[i] = right.ptr[i];
+    return *this; // enables x = y = z;
 }
 
-const ComplexArray &ComplexArray::operator=(const std::complex<double> &right)
+const ComplexArray& ComplexArray::operator=(const std::complex<double>& right)
 {
-	for (int i = 0;i < size;i++) ptr[i] = right;
-	return *this;// enables x = y = z;
+    for (int i = 0; i < size; i++)
+        ptr[i] = right;
+    return *this; // enables x = y = z;
 }
 
 //********************************************************
@@ -117,25 +120,25 @@ const ComplexArray &ComplexArray::operator=(const std::complex<double> &right)
 // const reference return create an cvakue
 //
 //********************************************************
-const std::complex<double> &ComplexArray::operator()
-(const int ind1,const int ind2,const int ind3)const
+const std::complex<double>& ComplexArray::operator()(const int ind1, const int ind2, const int ind3) const
 {
-	const int ind = (ind1 * bound2 + ind2) * bound3 + ind3 ;
-	return ptr[ind];
+    const int ind = (ind1 * bound2 + ind2) * bound3 + ind3;
+    return ptr[ind];
 }
 
-
 //********************************************************
 //
 // overloaded subscript operator for const real Array
 // const reference return creates an cvakue
 //
 //********************************************************
-const std::complex<double> &ComplexArray::operator()
-(const int ind1,const int ind2,const int ind3,const int ind4)const
+const std::complex<double>& ComplexArray::operator()(const int ind1,
+                                                     const int ind2,
+                                                     const int ind3,
+                                                     const int ind4) const
 {
-	const int ind = ((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4;
-	return ptr[ind];
+    const int ind = ((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4;
+    return ptr[ind];
 }
 
 //********************************************************
@@ -144,10 +147,10 @@ const std::complex<double> &ComplexArray::operator()
 // const reference return creates an lvakue
 //
 //********************************************************
-complex<double> &ComplexArray::operator()(const int ind1,const int ind2,const int ind3)
+complex<double>& ComplexArray::operator()(const int ind1, const int ind2, const int ind3)
 {
-	const int ind = (ind1 * bound2 + ind2) * bound3 + ind3;
-	return ptr[ind]; // reference return
+    const int ind = (ind1 * bound2 + ind2) * bound3 + ind3;
+    return ptr[ind]; // reference return
 }
 
 //********************************************************
@@ -156,10 +159,10 @@ complex<double> &ComplexArray::operator()(const int ind1,const int ind2,const in
 // const reference return creates an lvakue
 //
 //********************************************************
-complex<double> &ComplexArray::operator()(const int ind1,const int ind2,const int ind3,const int ind4)
+complex<double>& ComplexArray::operator()(const int ind1, const int ind2, const int ind3, const int ind4)
 {
-	const int ind = ((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4;
-	return ptr[ind];// reference return
+    const int ind = ((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4;
+    return ptr[ind]; // reference return
 }
 
 //****************************
@@ -169,7 +172,9 @@ complex<double> &ComplexArray::operator()(const int ind1,const int ind2,const in
 //****************************
 void ComplexArray::zero_out(void)
 {
-	if (size <= 0) return;
-	for (int i = 0;i < size; i++) ptr[i] = 0;
-	return;
+    if (size <= 0)
+        return;
+    for (int i = 0; i < size; i++)
+        ptr[i] = 0;
+    return;
 }
diff --git a/tools/01_NAO_generation/SIAB/src_tools/complexarray.h b/tools/01_NAO_generation/SIAB/src_tools/complexarray.h
index 6506fce530d..2f6820664e2 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/complexarray.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/complexarray.h
@@ -4,65 +4,79 @@
 #ifndef COMPLEXARRAY_H
 #define COMPLEXARRAY_H
 
-#include <iostream>
-#include <fstream>
-#include <iomanip>
 #include <cassert>
 #include <complex>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
 
 using namespace std;
 
 class ComplexArray
 {
-public:
-	complex<double> *ptr;
+  public:
+    complex<double>* ptr;
 
-	ComplexArray(const int d1 = 1 ,const int d2 = 1,const int d3 = 1);
-	ComplexArray(const int d1, const int d2,const int d3,const int d4);
-	~ComplexArray();
+    ComplexArray(const int d1 = 1, const int d2 = 1, const int d3 = 1);
+    ComplexArray(const int d1, const int d2, const int d3, const int d4);
+    ~ComplexArray();
 
-	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);
+    void create(const int d1, const int d2, const int d3, const int d4);
 
-	const ComplexArray &operator=(const ComplexArray &right);
-	const ComplexArray &operator=(const std::complex<double> &right);
+    const ComplexArray& operator=(const ComplexArray& right);
+    const ComplexArray& operator=(const std::complex<double>& right);
 
-	complex<double> &operator()(const int d1,const int d2,const int d3);
-	complex<double> &operator()(const int d1,const int d2,const int d3,const int d4);
+    complex<double>& operator()(const int d1, const int d2, const int d3);
+    complex<double>& operator()(const int d1, const int d2, const int d3, const int d4);
 
-	const std::complex<double> &operator()(const int d1,const int d2,const int d3)const;
-	const std::complex<double> &operator()(const int d1,const int d2,const int d3,const int d4)const;
+    const std::complex<double>& operator()(const int d1, const int d2, const int d3) const;
+    const std::complex<double>& operator()(const int d1, const int d2, const int d3, const int d4) const;
 
-	void zero_out(void);
+    void zero_out(void);
 
-	int getSize() const
-	{ return size;}
+    int getSize() const
+    {
+        return size;
+    }
 
-	int getDim() const
-	{ return dim;}
+    int getDim() const
+    {
+        return dim;
+    }
 
-	int getBound1() const
-	{ return bound1;}
+    int getBound1() const
+    {
+        return bound1;
+    }
 
-	int getBound2() const
-	{ return bound2;}
+    int getBound2() const
+    {
+        return bound2;
+    }
 
-	int getBound3() const
-	{ return bound3;}
+    int getBound3() const
+    {
+        return bound3;
+    }
 
-	int getBound4() const
-	{ return bound4;}
+    int getBound4() const
+    {
+        return bound4;
+    }
 
     int getArrayCount(void)
-	{ return arrayCount;}
+    {
+        return arrayCount;
+    }
 
-private:
-	int size;
-	int dim;
-	int bound1, bound2, bound3, bound4;
-	static int arrayCount;
+  private:
+    int size;
+    int dim;
+    int bound1, bound2, bound3, bound4;
+    static int arrayCount;
 
-	void freemem();
+    void freemem();
 };
 
-#endif	// ComplexArray class
+#endif // ComplexArray class
diff --git a/tools/01_NAO_generation/SIAB/src_tools/complexmatrix.h b/tools/01_NAO_generation/SIAB/src_tools/complexmatrix.h
index 9d75a910af2..d3225369979 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/complexmatrix.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/complexmatrix.h
@@ -11,110 +11,110 @@ using namespace std;
 class ComplexMatrix
 {
 
-public:
-
-	complex<double> *c;
-	int nr;
-	int nc;
-	int size;
-
-	//==============
-	// Constructors
-	//==============
-	inline ComplexMatrix(const int nrows = 1,const int ncols = 1);
-
-	//===================
-	// Copy constructor
-	//==================
-	inline ComplexMatrix(const ComplexMatrix &m1);
-	inline ~ComplexMatrix();
-
-	//============
-	// Operators
-	//============
-	complex<double> &operator()(const int i,const int j)
-	{
-		return c[nc * i + j];//mohan modify in-line 2007-10-1
-	}
-	const std::complex<double> &operator()(const int i,const int j)const
-	{
-		return c[nc * i + j];//mohan modify in-line 2007-10-13
-	}
-
-	friend ComplexMatrix operator+(const ComplexMatrix &m1,  const ComplexMatrix &m2);
-	friend ComplexMatrix operator-(const ComplexMatrix &m1,  const ComplexMatrix &m2);
-	friend ComplexMatrix operator*(const ComplexMatrix &m1,  const ComplexMatrix &m2);
-	friend ComplexMatrix multiply_special(const char Symmetry, const char Side, const char Uplo, const ComplexMatrix &m1, const ComplexMatrix &m2);
-	friend ComplexMatrix operator*(const std::complex<double> &s, const ComplexMatrix &m);
-	friend ComplexMatrix operator*(const ComplexMatrix &m,   const std::complex<double> &s);
-	friend ComplexMatrix operator*(const double &s,          const ComplexMatrix &m);
-	friend ComplexMatrix operator*(const ComplexMatrix &m,   const double &s);
-	inline ComplexMatrix& operator=(const ComplexMatrix &m);
-	inline ComplexMatrix& operator*=(const std::complex<double> &s);
-	inline ComplexMatrix& operator+=(const ComplexMatrix &m);
-	inline ComplexMatrix& operator-=(const ComplexMatrix &m);
-
-	inline ComplexMatrix &add(const ComplexMatrix &m1, const ComplexMatrix &m2);
-	inline ComplexMatrix &minus(const ComplexMatrix &m1, const ComplexMatrix &m2);
-	inline ComplexMatrix &multiply(const ComplexMatrix &m1, const ComplexMatrix &m2);
-	inline ComplexMatrix &multiply_special(const char Symmetry, const char Side, const char Uplo, const ComplexMatrix &m1, const ComplexMatrix &m2);
-	
-	//==================
-	// member function:
-	//==================
-	inline void create(const int nrow,const int ncol);
-	inline void zero_out(void);
-	static int& getMCount(void){return mCount;}
-	inline void set_as_identity_matrix(void);
-private:
-
-	static int mCount;
-	inline void freemem(void);//mohan add 2007-11-20
-	inline void init(const int nrows,const int ncols);
+  public:
+    complex<double>* c;
+    int nr;
+    int nc;
+    int size;
+
+    //==============
+    // Constructors
+    //==============
+    inline ComplexMatrix(const int nrows = 1, const int ncols = 1);
+
+    //===================
+    // Copy constructor
+    //==================
+    inline ComplexMatrix(const ComplexMatrix& m1);
+    inline ~ComplexMatrix();
+
+    //============
+    // Operators
+    //============
+    complex<double>& operator()(const int i, const int j)
+    {
+        return c[nc * i + j]; // mohan modify in-line 2007-10-1
+    }
+    const std::complex<double>& operator()(const int i, const int j) const
+    {
+        return c[nc * i + j]; // mohan modify in-line 2007-10-13
+    }
+
+    friend ComplexMatrix operator+(const ComplexMatrix& m1, const ComplexMatrix& m2);
+    friend ComplexMatrix operator-(const ComplexMatrix& m1, const ComplexMatrix& m2);
+    friend ComplexMatrix operator*(const ComplexMatrix& m1, const ComplexMatrix& m2);
+    friend ComplexMatrix multiply_special(const char Symmetry,
+                                          const char Side,
+                                          const char Uplo,
+                                          const ComplexMatrix& m1,
+                                          const ComplexMatrix& m2);
+    friend ComplexMatrix operator*(const std::complex<double>& s, const ComplexMatrix& m);
+    friend ComplexMatrix operator*(const ComplexMatrix& m, const std::complex<double>& s);
+    friend ComplexMatrix operator*(const double& s, const ComplexMatrix& m);
+    friend ComplexMatrix operator*(const ComplexMatrix& m, const double& s);
+    inline ComplexMatrix& operator=(const ComplexMatrix& m);
+    inline ComplexMatrix& operator*=(const std::complex<double>& s);
+    inline ComplexMatrix& operator+=(const ComplexMatrix& m);
+    inline ComplexMatrix& operator-=(const ComplexMatrix& m);
+
+    inline ComplexMatrix& add(const ComplexMatrix& m1, const ComplexMatrix& m2);
+    inline ComplexMatrix& minus(const ComplexMatrix& m1, const ComplexMatrix& m2);
+    inline ComplexMatrix& multiply(const ComplexMatrix& m1, const ComplexMatrix& m2);
+    inline ComplexMatrix& multiply_special(const char Symmetry,
+                                           const char Side,
+                                           const char Uplo,
+                                           const ComplexMatrix& m1,
+                                           const ComplexMatrix& m2);
+
+    //==================
+    // member function:
+    //==================
+    inline void create(const int nrow, const int ncol);
+    inline void zero_out(void);
+    static int& getMCount(void)
+    {
+        return mCount;
+    }
+    inline void set_as_identity_matrix(void);
+
+  private:
+    static int mCount;
+    inline void freemem(void); // mohan add 2007-11-20
+    inline void init(const int nrows, const int ncols);
 };
 
-inline std::complex<double> trace(const ComplexMatrix &m);
+inline std::complex<double> trace(const ComplexMatrix& m);
 // mohan add 2008-7-1
-inline double abs2_row(const ComplexMatrix &m,const int ir);
+inline double abs2_row(const ComplexMatrix& m, const int ir);
 // mohan add 2008-7-1
-inline double abs2_column(const ComplexMatrix &m,const int ic);
-inline double abs2(const ComplexMatrix &m);
-inline double abs2(const int nmat, ComplexMatrix **m);
-inline double abs2_triangle_matrix( const char Uplo, ComplexMatrix &m);
-
-inline ComplexMatrix transpose(const ComplexMatrix &m, const bool &conjugate);
-
-inline void scale_accumulate(
-		const std::complex<double> &s, 
-		const ComplexMatrix &min, 
-		ComplexMatrix &mout);
-
-inline void scale_accumulate(
-		const int &nmat, 
-		const std::complex<double> &s, 
-		ComplexMatrix **min, 
-		ComplexMatrix **mout);
-
-inline void scaled_sum(
-		const std::complex<double> &s1, 
-		const ComplexMatrix &m1, 
-		const std::complex<double> &s2, 
-		const ComplexMatrix &m2, 
-		ComplexMatrix &mout);
-
-inline void scaled_sum(
-		const int &nmat,
-		const std::complex<double> &s1, 
-		ComplexMatrix **m1, 
-		const std::complex<double> &s2,
-		ComplexMatrix **m2, 
-		ComplexMatrix **mout);
-		
-inline ComplexMatrix multiply_special(
-		const char Symmetry, 
-		const char Side, 
-		const char Uplo, 
-		const ComplexMatrix &m1, 
-		const ComplexMatrix &m2);
-		
+inline double abs2_column(const ComplexMatrix& m, const int ic);
+inline double abs2(const ComplexMatrix& m);
+inline double abs2(const int nmat, ComplexMatrix** m);
+inline double abs2_triangle_matrix(const char Uplo, ComplexMatrix& m);
+
+inline ComplexMatrix transpose(const ComplexMatrix& m, const bool& conjugate);
+
+inline void scale_accumulate(const std::complex<double>& s, const ComplexMatrix& min, ComplexMatrix& mout);
+
+inline void scale_accumulate(const int& nmat, const std::complex<double>& s, ComplexMatrix** min, ComplexMatrix** mout);
+
+inline void scaled_sum(const std::complex<double>& s1,
+                       const ComplexMatrix& m1,
+                       const std::complex<double>& s2,
+                       const ComplexMatrix& m2,
+                       ComplexMatrix& mout);
+
+inline void scaled_sum(const int& nmat,
+                       const std::complex<double>& s1,
+                       ComplexMatrix** m1,
+                       const std::complex<double>& s2,
+                       ComplexMatrix** m2,
+                       ComplexMatrix** mout);
+
+inline ComplexMatrix multiply_special(const char Symmetry,
+                                      const char Side,
+                                      const char Uplo,
+                                      const ComplexMatrix& m1,
+                                      const ComplexMatrix& m2);
+
 #endif
diff --git a/tools/01_NAO_generation/SIAB/src_tools/complexmatrix_inline.h b/tools/01_NAO_generation/SIAB/src_tools/complexmatrix_inline.h
index e4698810a89..68d3e3527e9 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/complexmatrix_inline.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/complexmatrix_inline.h
@@ -6,472 +6,469 @@
 #ifndef COMPLEXMATRIX_INLINE_H
 #define COMPLEXMATRIX_INLINE_H
 
+#include "complexarray.h"
+#include "complexmatrix.h"
+#include "lapack_connector.h"
+
 #include <cassert>
-#include <new>
 #include <cstdlib>
-#include <iostream>
-#include <exception>
 #include <cstring>
-#include "complexmatrix.h"
-#include "lapack_connector.h"
-#include "complexarray.h"
+#include <exception>
+#include <iostream>
+#include <new>
 
 // constructor with sizes
 inline ComplexMatrix::ComplexMatrix(int nrows, int ncols)
 {
-	this->init(nrows, ncols);
-	++mCount;
+    this->init(nrows, ncols);
+    ++mCount;
 }
 
 // zero out the ComplexMatrix
 inline void ComplexMatrix::zero_out(void)
 {
-//	for (int i=0; i<size; i++) c[i] = std::complex<double>(0.0,0.0);
-	memset ( c, 0, size*sizeof(complex<double>) );		// Peize Lin update 2015-12-08
+    //	for (int i=0; i<size; i++) c[i] = std::complex<double>(0.0,0.0);
+    memset(c, 0, size * sizeof(complex<double>)); // Peize Lin update 2015-12-08
 }
 
 /*
 void need_more_memory()
 {
-	cout << "\n Sorry to crash... but the running need more momory! Exit." << endl;
-	exit(0);
+    cout << "\n Sorry to crash... but the running need more momory! Exit." << endl;
+    exit(0);
 }
 */
 
-inline void ComplexMatrix::init(const int nrows,const int ncols)
+inline void ComplexMatrix::init(const int nrows, const int ncols)
 {
-	this->nr = nrows;
-	this->nc = ncols;
-	this->size = nr * nc;
-//	std::set_new_handler( need_more_memory );
-	this->c = new std::complex<double>[size];
-	assert(c != 0);
-	this->zero_out();
-}//mohan the principle:as simple as possible:modify inline 2007-10-13
+    this->nr = nrows;
+    this->nc = ncols;
+    this->size = nr * nc;
+    //	std::set_new_handler( need_more_memory );
+    this->c = new std::complex<double>[size];
+    assert(c != 0);
+    this->zero_out();
+} // mohan the principle:as simple as possible:modify inline 2007-10-13
 
 // Copy constructor
-inline ComplexMatrix::ComplexMatrix(const ComplexMatrix &m1)
+inline ComplexMatrix::ComplexMatrix(const ComplexMatrix& m1)
 {
-	this->init(m1.nr, m1.nc);
-	++mCount;
-	for(int i=0; i<this->size; i++) c[i] = m1.c[i];
+    this->init(m1.nr, m1.nc);
+    ++mCount;
+    for (int i = 0; i < this->size; i++)
+        c[i] = m1.c[i];
 }
 
 // deconstructor
 inline ComplexMatrix::~ComplexMatrix()
 {
-	this->freemem();
+    this->freemem();
 }
 
 // Free up memory for ComplexMatrix
 inline void ComplexMatrix::freemem(void)
 {
-	delete[] c;
-	c = NULL;
+    delete[] c;
+    c = NULL;
 }
 
 // reallocate memory for Complex Matrix
-inline void ComplexMatrix::create(const int nrow,const int ncol)
+inline void ComplexMatrix::create(const int nrow, const int ncol)
 {
-	// because c has been 'new' in  init function.
-	delete[] c;
-	this->init(nrow, ncol);
-	return;
+    // because c has been 'new' in  init function.
+    delete[] c;
+    this->init(nrow, ncol);
+    return;
 }
 
 inline void ComplexMatrix::set_as_identity_matrix()
 {
-/*	for(int i=0; i<nr; i++)
-	{
-		for(int j=0; j<nc; j++)
-		{
-			if(i==j) c[nc * i + j] = std::complex<double>(1.0, 0.0);  
-			else c[nc * i + j] = std::complex<double>(0.0, 0.0); 
-		}
-	}
-*/
-	// Peize Lin update 2015-12-08
-	zero_out();
-	const int nrc_min = min(nr,nc);
-	for( int i=0; i<nrc_min; ++i)
-	{
-		c[nc*i+i] = std::complex<double>(1.0,0.0);
-	}
-	return;
+    /*	for(int i=0; i<nr; i++)
+        {
+            for(int j=0; j<nc; j++)
+            {
+                if(i==j) c[nc * i + j] = std::complex<double>(1.0, 0.0);
+                else c[nc * i + j] = std::complex<double>(0.0, 0.0);
+            }
+        }
+    */
+    // Peize Lin update 2015-12-08
+    zero_out();
+    const int nrc_min = min(nr, nc);
+    for (int i = 0; i < nrc_min; ++i)
+    {
+        c[nc * i + i] = std::complex<double>(1.0, 0.0);
+    }
+    return;
 }
 
-// Adding matrices, as a friend 
-inline ComplexMatrix operator+(const ComplexMatrix &m1, const ComplexMatrix &m2)
+// Adding matrices, as a friend
+inline ComplexMatrix operator+(const ComplexMatrix& m1, const ComplexMatrix& m2)
 {
-	ComplexMatrix tm(m1);
-	tm+=m2;
-	return tm;
+    ComplexMatrix tm(m1);
+    tm += m2;
+    return tm;
 }
 
-// Subtracting matrices, as a friend 
-inline ComplexMatrix operator-(const ComplexMatrix &m1, const ComplexMatrix &m2)
+// Subtracting matrices, as a friend
+inline ComplexMatrix operator-(const ComplexMatrix& m1, const ComplexMatrix& m2)
 {
-	ComplexMatrix tm(m1);
-	tm-=m2;
-	return tm;
+    ComplexMatrix tm(m1);
+    tm -= m2;
+    return tm;
 }
 
 // Multiplying matrices, as a friend
 // mprod = m1 * m2
-inline ComplexMatrix operator*(const ComplexMatrix &m1, const ComplexMatrix &m2)
-{
-	assert(m1.nc == m2.nr);	
-	ComplexMatrix mprod(m1.nr, m2.nc);
-
-/*							
-	complex<double> z;
-	for (int i = 0;i < m1.nr;i++)
-	{
-		for (int j = 0;j < m2.nc;j++)
-		{
-			z = std::complex<double>(0,0);
-			for (int k = 0;k < m1.nc;k++)
-			{
-				z += m1(i, k) * m2(k, j);
-			}
-			mprod(i, j) = z;
-		}
-	}
-*/
-	// Peize Lin update 2015-12-07
-	LapackConnector::zgemm('N', 'N',
-		m1.nr, m2.nc, m1.nc,
-		1.0, m1.c, m1.nc, m2.c, m2.nc,
-		0.0, mprod.c, mprod.nc);
-	return mprod;
+inline ComplexMatrix operator*(const ComplexMatrix& m1, const ComplexMatrix& m2)
+{
+    assert(m1.nc == m2.nr);
+    ComplexMatrix mprod(m1.nr, m2.nc);
+
+    /*
+        complex<double> z;
+        for (int i = 0;i < m1.nr;i++)
+        {
+            for (int j = 0;j < m2.nc;j++)
+            {
+                z = std::complex<double>(0,0);
+                for (int k = 0;k < m1.nc;k++)
+                {
+                    z += m1(i, k) * m2(k, j);
+                }
+                mprod(i, j) = z;
+            }
+        }
+    */
+    // Peize Lin update 2015-12-07
+    LapackConnector::zgemm('N', 'N', m1.nr, m2.nc, m1.nc, 1.0, m1.c, m1.nc, m2.c, m2.nc, 0.0, mprod.c, mprod.nc);
+    return mprod;
 }
 
 // Peize Lin add 2015-12-08
 // m1^T=m1('T') or m1^H=m1('H')
 // mprod=m1*m2('L'eft) or mprod=m2*m1('R'ight)
 // m1 stored in 'U'pper or 'L'ower
-inline ComplexMatrix multiply_special(const char Symmetry, const char Side, const char Uplo, const ComplexMatrix &m1, const ComplexMatrix &m2)
-{
-	assert(m1.nr == m1.nc);
-	assert(m1.nc == m2.nr);
-	ComplexMatrix mprod(m1.nr, m2.nc);
-	if('T'==Symmetry)
-	{
-		LapackConnector::zsymm(Side, Uplo,
-			  m1.nr, m2.nc,
-			  1.0, m1.c, m1.nc, m2.c, m2.nc,
-			  0.0, mprod.c, mprod.nc);
-	}
-	else if('H'==Symmetry)
-	{
-		LapackConnector::zhemm(Side, Uplo,
-			  m1.nr, m2.nc,
-			  1.0, m1.c, m1.nc, m2.c, m2.nc,
-			  0.0, mprod.c, mprod.nc);
-	}
-	else
-	{
-		throw invalid_argument("Symmetry must be 'T' or 'H'");
-	}
-	return mprod;	
+inline ComplexMatrix multiply_special(const char Symmetry,
+                                      const char Side,
+                                      const char Uplo,
+                                      const ComplexMatrix& m1,
+                                      const ComplexMatrix& m2)
+{
+    assert(m1.nr == m1.nc);
+    assert(m1.nc == m2.nr);
+    ComplexMatrix mprod(m1.nr, m2.nc);
+    if ('T' == Symmetry)
+    {
+        LapackConnector::zsymm(Side, Uplo, m1.nr, m2.nc, 1.0, m1.c, m1.nc, m2.c, m2.nc, 0.0, mprod.c, mprod.nc);
+    }
+    else if ('H' == Symmetry)
+    {
+        LapackConnector::zhemm(Side, Uplo, m1.nr, m2.nc, 1.0, m1.c, m1.nc, m2.c, m2.nc, 0.0, mprod.c, mprod.nc);
+    }
+    else
+    {
+        throw invalid_argument("Symmetry must be 'T' or 'H'");
+    }
+    return mprod;
 }
 
 // Scale a ComplexMatrix
-inline ComplexMatrix operator*(const std::complex<double> &c,const ComplexMatrix &m)
+inline ComplexMatrix operator*(const std::complex<double>& c, const ComplexMatrix& m)
 {
-	ComplexMatrix sm(m);
-	for (int i=0 ;i<m.size; i++) sm.c[i] *= c;
-	return sm;
+    ComplexMatrix sm(m);
+    for (int i = 0; i < m.size; i++)
+        sm.c[i] *= c;
+    return sm;
 }
 
 // ComplexMatrix scalar
-inline ComplexMatrix operator*(const ComplexMatrix &m,const std::complex<double> &c)
+inline ComplexMatrix operator*(const ComplexMatrix& m, const std::complex<double>& c)
 {
-	ComplexMatrix sm(m);
-	for (int i = 0;i < m.size;i++) sm.c[i] *= c;
-	return sm;
+    ComplexMatrix sm(m);
+    for (int i = 0; i < m.size; i++)
+        sm.c[i] *= c;
+    return sm;
 }
 
-inline ComplexMatrix operator*(const double &r,const ComplexMatrix &m)
+inline ComplexMatrix operator*(const double& r, const ComplexMatrix& m)
 {
-	ComplexMatrix sm(m);
-	for(int i=0; i<m.size; i++) sm.c[i]*= r;
-	return sm;
+    ComplexMatrix sm(m);
+    for (int i = 0; i < m.size; i++)
+        sm.c[i] *= r;
+    return sm;
 }
 
-inline ComplexMatrix operator*(const ComplexMatrix &m,const double &r)
+inline ComplexMatrix operator*(const ComplexMatrix& m, const double& r)
 {
-	ComplexMatrix sm(m);
-	for (int i=0; i<m.size; i++) sm.c[i] *= r;
-	return sm;
+    ComplexMatrix sm(m);
+    for (int i = 0; i < m.size; i++)
+        sm.c[i] *= r;
+    return sm;
 }
 
-inline ComplexMatrix &ComplexMatrix::add(const ComplexMatrix &m1, const ComplexMatrix &m2)
+inline ComplexMatrix& ComplexMatrix::add(const ComplexMatrix& m1, const ComplexMatrix& m2)
 {
-	for(int i=0; i<size; ++i) this->c[i] = m1.c[i] + m2.c[i];
-	return *this;
+    for (int i = 0; i < size; ++i)
+        this->c[i] = m1.c[i] + m2.c[i];
+    return *this;
 }
 
-inline ComplexMatrix &ComplexMatrix::minus(const ComplexMatrix &m1, const ComplexMatrix &m2)
+inline ComplexMatrix& ComplexMatrix::minus(const ComplexMatrix& m1, const ComplexMatrix& m2)
 {
-	for(int i=0; i<size; ++i) this->c[i] = m1.c[i] - m2.c[i];
-	return *this;
+    for (int i = 0; i < size; ++i)
+        this->c[i] = m1.c[i] - m2.c[i];
+    return *this;
 }
-			   
-inline ComplexMatrix &ComplexMatrix::multiply(const ComplexMatrix &m1, const ComplexMatrix &m2)
+
+inline ComplexMatrix& ComplexMatrix::multiply(const ComplexMatrix& m1, const ComplexMatrix& m2)
 {
-	LapackConnector::zgemm('N', 'N',
-		m1.nr, m2.nc, m1.nc,
-		1.0, m1.c, m1.nc, m2.c, m2.nc,
-		0.0, this->c, this->nc);
-	return *this;
+    LapackConnector::zgemm('N', 'N', m1.nr, m2.nc, m1.nc, 1.0, m1.c, m1.nc, m2.c, m2.nc, 0.0, this->c, this->nc);
+    return *this;
 }
-		
+
 // Peize Lin add 2015-12-08
 // m1^T=m1('T') or m1^H=m1('H')
 // mprod=m1*m2('L'eft) or mprod=m2*m1('R'ight)
-// m1 stored in 'U'pper or 'L'ower		
-inline ComplexMatrix &ComplexMatrix::multiply_special(const char Symmetry, const char Side, const char Uplo, const ComplexMatrix &m1, const ComplexMatrix &m2)
-{
-	if('T'==Symmetry)
-	{
-		LapackConnector::zsymm(Side, Uplo,
-			  m1.nr, m2.nc,
-			  1.0, m1.c, m1.nc, m2.c, m2.nc,
-			  0.0, this->c, this->nc);
-	}
-	else if('H'==Symmetry)
-	{
-		LapackConnector::zhemm(Side, Uplo,
-			  m1.nr, m2.nc,
-			  1.0, m1.c, m1.nc, m2.c, m2.nc,
-			  0.0, this->c, this->nc);
-	}
-	else
-	{
-		throw invalid_argument("Symmetry must be 'T' or 'H'");
-	}
-	return *this;	
-}
-
-inline ComplexMatrix& ComplexMatrix::operator=(const ComplexMatrix &m)
-{
-	if(m.nr!=this->nr || m.nc!=this->nc)
-	{
-		cout << "\n row/col number can't match in ComplexMatrix '=' operator\n";
-		cout << " this nr = " << this->nr;
-		cout << " this nc = " << this->nc;
-		cout << " in nr = " << m.nr;
-		cout << " in nc = " << m.nc;
-		exit(0);
-	}
-/*	else {
-		this->create(m.nr, m.nc);
-	} */	// why ?
-
-//	for(int i=0;i<this->size;i++) c[i] = m.c[i];
-	memcpy ( c, m.c, this->size * sizeof(complex<double>) );	// Peize Lin update 2015-12-08
-	return *this;
-}
-
-inline ComplexMatrix& ComplexMatrix::operator*=(const std::complex<double> &s)
-{
-	for (int i = 0;i < this->size;i++) c[i] *= s;
-	return *this;
+// m1 stored in 'U'pper or 'L'ower
+inline ComplexMatrix& ComplexMatrix::multiply_special(const char Symmetry,
+                                                      const char Side,
+                                                      const char Uplo,
+                                                      const ComplexMatrix& m1,
+                                                      const ComplexMatrix& m2)
+{
+    if ('T' == Symmetry)
+    {
+        LapackConnector::zsymm(Side, Uplo, m1.nr, m2.nc, 1.0, m1.c, m1.nc, m2.c, m2.nc, 0.0, this->c, this->nc);
+    }
+    else if ('H' == Symmetry)
+    {
+        LapackConnector::zhemm(Side, Uplo, m1.nr, m2.nc, 1.0, m1.c, m1.nc, m2.c, m2.nc, 0.0, this->c, this->nc);
+    }
+    else
+    {
+        throw invalid_argument("Symmetry must be 'T' or 'H'");
+    }
+    return *this;
+}
+
+inline ComplexMatrix& ComplexMatrix::operator=(const ComplexMatrix& m)
+{
+    if (m.nr != this->nr || m.nc != this->nc)
+    {
+        cout << "\n row/col number can't match in ComplexMatrix '=' operator\n";
+        cout << " this nr = " << this->nr;
+        cout << " this nc = " << this->nc;
+        cout << " in nr = " << m.nr;
+        cout << " in nc = " << m.nc;
+        exit(0);
+    }
+    /*	else {
+            this->create(m.nr, m.nc);
+        } */	// why ?
+
+    //	for(int i=0;i<this->size;i++) c[i] = m.c[i];
+    memcpy(c, m.c, this->size * sizeof(complex<double>)); // Peize Lin update 2015-12-08
+    return *this;
+}
+
+inline ComplexMatrix& ComplexMatrix::operator*=(const std::complex<double>& s)
+{
+    for (int i = 0; i < this->size; i++)
+        c[i] *= s;
+    return *this;
 }
 
 // Accumulate to a ComplexMatrix in place
-inline ComplexMatrix& ComplexMatrix::operator+=(const ComplexMatrix &m)
+inline ComplexMatrix& ComplexMatrix::operator+=(const ComplexMatrix& m)
 {
-	for(int i=0; i<size; i++) this->c[i] += m.c[i];
-	return *this;
+    for (int i = 0; i < size; i++)
+        this->c[i] += m.c[i];
+    return *this;
 }
 
 // decumulate to a ComplexMatrix in place
-inline ComplexMatrix& ComplexMatrix::operator-=(const ComplexMatrix &m)
+inline ComplexMatrix& ComplexMatrix::operator-=(const ComplexMatrix& m)
 {
-	for(int i=0; i<size; i++) this->c[i] -= m.c[i];
-	return *this;
+    for (int i = 0; i < size; i++)
+        this->c[i] -= m.c[i];
+    return *this;
 }
 
 // Returns trace of ComplexMatrix
-inline std::complex<double> trace(const ComplexMatrix &m)
+inline std::complex<double> trace(const ComplexMatrix& m)
 {
-	complex<double> tr=complex<double>(0,0);
-	assert(m.nr == m.nc);
-	for (int i=0; i<m.nr; i++) tr += m(i, i);
-	return tr;
+    complex<double> tr = complex<double>(0, 0);
+    assert(m.nr == m.nc);
+    for (int i = 0; i < m.nr; i++)
+        tr += m(i, i);
+    return tr;
 }
 
 // Do mout += s*min
-inline void scale_accumulate(const std::complex<double> &s,
-                      const ComplexMatrix &min,
-                      ComplexMatrix &mout)
+inline void scale_accumulate(const std::complex<double>& s, const ComplexMatrix& min, ComplexMatrix& mout)
 {
-	assert(min.nr == mout.nr);
-	assert(min.nc == mout.nc);
-	for (int j=0; j<min.size; j++)
-	{
-		mout.c[j] += s * min.c[j];
-	}
-	return;
+    assert(min.nr == mout.nr);
+    assert(min.nc == mout.nc);
+    for (int j = 0; j < min.size; j++)
+    {
+        mout.c[j] += s * min.c[j];
+    }
+    return;
 }
 
 // Do mout[i] += s*min[i]
-inline void scale_accumulate(const int &nmat,
-                      const std::complex<double> &s,
-                      ComplexMatrix **min,
-                      ComplexMatrix **mout)
+inline void scale_accumulate(const int& nmat, const std::complex<double>& s, ComplexMatrix** min, ComplexMatrix** mout)
 {
-	assert(nmat>=0);
-	for (int i=0; i<nmat; i++)
-	{
-		scale_accumulate(s, *min[i], *mout[i]);
-	}
-	return;
+    assert(nmat >= 0);
+    for (int i = 0; i < nmat; i++)
+    {
+        scale_accumulate(s, *min[i], *mout[i]);
+    }
+    return;
 }
 
 // Do mout = s1*m1 + s2*m2
-inline void scaled_sum(const std::complex<double> &s1,
-                const ComplexMatrix &m1,
-                const std::complex<double> &s2,
-                const ComplexMatrix &m2,
-                ComplexMatrix &mout)
+inline void scaled_sum(const std::complex<double>& s1,
+                       const ComplexMatrix& m1,
+                       const std::complex<double>& s2,
+                       const ComplexMatrix& m2,
+                       ComplexMatrix& mout)
 {
-	assert(m1.nr == m2.nr);
-	assert(m1.nr == mout.nr);
-	assert(m1.nc == m2.nc);
-	assert(m1.nc == mout.nc);
+    assert(m1.nr == m2.nr);
+    assert(m1.nr == mout.nr);
+    assert(m1.nc == m2.nc);
+    assert(m1.nc == mout.nc);
 
-	for(int i=0; i<m1.size; i++)
-	{
-		mout.c[i] = s1 * m1.c[i] + s2 * m2.c[i];
-	}
-	return;
+    for (int i = 0; i < m1.size; i++)
+    {
+        mout.c[i] = s1 * m1.c[i] + s2 * m2.c[i];
+    }
+    return;
 }
 
 // Does mout[i] = s1*m1[i] + s2*m2[i]
-inline void scaled_sum(const int &nmat,
-                const std::complex<double> &s1,
-                ComplexMatrix **m1,
-                const std::complex<double> &s2,
-                ComplexMatrix **m2,
-                ComplexMatrix **mout)
+inline void scaled_sum(const int& nmat,
+                       const std::complex<double>& s1,
+                       ComplexMatrix** m1,
+                       const std::complex<double>& s2,
+                       ComplexMatrix** m2,
+                       ComplexMatrix** mout)
 {
-	assert(nmat>0);
-	for(int i=0; i<nmat; i++)
-	{
-		scaled_sum(s1, *m1[i], s2, *m2[i], *mout[i]);
-	}
-	return;
+    assert(nmat > 0);
+    for (int i = 0; i < nmat; i++)
+    {
+        scaled_sum(s1, *m1[i], s2, *m2[i], *mout[i]);
+    }
+    return;
 }
 
-inline double abs2_row(const ComplexMatrix &m,const int ir)
+inline double abs2_row(const ComplexMatrix& m, const int ir)
 {
-	double r=0.0;
-	complex<double> z;
-	for(int ic=0;ic<m.nc;ic++)
-	{
-		z = m.c[ m.nc*ir + ic];
-		r += z.real()*z.real() + z.imag()*z.imag();
-	}
-	return r;
+    double r = 0.0;
+    complex<double> z;
+    for (int ic = 0; ic < m.nc; ic++)
+    {
+        z = m.c[m.nc * ir + ic];
+        r += z.real() * z.real() + z.imag() * z.imag();
+    }
+    return r;
 }
 
-inline double abs2_column(const ComplexMatrix &m,const int ic)
+inline double abs2_column(const ComplexMatrix& m, const int ic)
 {
-	double r=0.0;
-	complex<double> z;
-	for(int ir=0;ir<m.nr;ir++)
-	{
-		z = m.c[ m.nc*ir + ic ];
-		r += z.real()*z.real() + z.imag()*z.imag();
-	}
-	return r;
+    double r = 0.0;
+    complex<double> z;
+    for (int ir = 0; ir < m.nr; ir++)
+    {
+        z = m.c[m.nc * ir + ic];
+        r += z.real() * z.real() + z.imag() * z.imag();
+    }
+    return r;
 }
 
 // returns absolute square magnitude of sum of all ComplexMatrix elements
-inline double abs2(const ComplexMatrix &m)
+inline double abs2(const ComplexMatrix& m)
 {
-	double r=0.0;
-	complex<double> z;
-	
-	for (int i = 0;i < m.size;i++)
-	{
-		z = m.c[i];
-		r += z.real() * z.real() + z.imag() * z.imag();
-	}
-	return r;
+    double r = 0.0;
+    complex<double> z;
+
+    for (int i = 0; i < m.size; i++)
+    {
+        z = m.c[i];
+        r += z.real() * z.real() + z.imag() * z.imag();
+    }
+    return r;
 }
 
 // Same for an array of matrices
-inline double abs2(const int nmat, ComplexMatrix **m)
-{
-	double r = 0.0;
-	for (int i = 0;i < nmat;i++)
-	{
-		r += abs2(*m[i]);
-	}
-	return r;
-}
-
-inline double abs2_triangle_matrix( const char Uplo, ComplexMatrix &m)
-{
-	double r(0.0);
-	complex<double> z;
-	
-	if('U'==Uplo)
-	{
-		for (int i=0; i<m.nr; ++i)
-		{
-			for (int j=i+1; j<m.nc; ++j)
-			{
-				z = m(i,j);
-				r += z.real() * z.real() + z.imag() * z.imag();
-			}
-		}		
-	}
-	else if('L'==Uplo)
-	{
-		for (int i=0; i<m.nr; ++i)
-		{
-			for (int j=0; j<i; ++j)
-			{
-				z = m(i,j);
-				r += z.real() * z.real() + z.imag() * z.imag();				
-			}
-		}
-	}
-	else
-	{
-		throw invalid_argument("Uplo must be 'U' or 'L'");
-	}
-}
-
-inline ComplexMatrix transpose(const ComplexMatrix &m, const bool &conjugate)
-{
-	ComplexMatrix tm(m.nc, m.nr);
-	if(conjugate)
-	{
-		for (int i = 0;i < m.nr;i++)
-		{
-			for (int j = 0;j < m.nc;j++)
-			{
-				tm(j, i) = conj ( m(i, j) );
-			}
-		}
-	}
-	else
-	{
-		for (int i = 0;i < m.nr;i++)
-		{
-			for (int j = 0;j < m.nc;j++)
-			{
-				tm(j, i) = m(i, j);
-			}
-		}
-	}
-	return tm;
+inline double abs2(const int nmat, ComplexMatrix** m)
+{
+    double r = 0.0;
+    for (int i = 0; i < nmat; i++)
+    {
+        r += abs2(*m[i]);
+    }
+    return r;
+}
+
+inline double abs2_triangle_matrix(const char Uplo, ComplexMatrix& m)
+{
+    double r(0.0);
+    complex<double> z;
+
+    if ('U' == Uplo)
+    {
+        for (int i = 0; i < m.nr; ++i)
+        {
+            for (int j = i + 1; j < m.nc; ++j)
+            {
+                z = m(i, j);
+                r += z.real() * z.real() + z.imag() * z.imag();
+            }
+        }
+    }
+    else if ('L' == Uplo)
+    {
+        for (int i = 0; i < m.nr; ++i)
+        {
+            for (int j = 0; j < i; ++j)
+            {
+                z = m(i, j);
+                r += z.real() * z.real() + z.imag() * z.imag();
+            }
+        }
+    }
+    else
+    {
+        throw invalid_argument("Uplo must be 'U' or 'L'");
+    }
+}
+
+inline ComplexMatrix transpose(const ComplexMatrix& m, const bool& conjugate)
+{
+    ComplexMatrix tm(m.nc, m.nr);
+    if (conjugate)
+    {
+        for (int i = 0; i < m.nr; i++)
+        {
+            for (int j = 0; j < m.nc; j++)
+            {
+                tm(j, i) = conj(m(i, j));
+            }
+        }
+    }
+    else
+    {
+        for (int i = 0; i < m.nr; i++)
+        {
+            for (int j = 0; j < m.nc; j++)
+            {
+                tm(j, i) = m(i, j);
+            }
+        }
+    }
+    return tm;
 }
 
 #endif
\ No newline at end of file
diff --git a/tools/01_NAO_generation/SIAB/src_tools/intarray.cpp b/tools/01_NAO_generation/SIAB/src_tools/intarray.cpp
index 0f8cbf486cf..ed987d1438e 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/intarray.cpp
+++ b/tools/01_NAO_generation/SIAB/src_tools/intarray.cpp
@@ -1,212 +1,274 @@
 /*******************************************
  * ESCP:Electro-Structure Calculate Package.
  ********************************************/
-#include <cstdlib>
 #include "intarray.h"
 
+#include <cstdlib>
+
 int IntArray::arrayCount = 0;
 
 void IntArrayAlloc()
 {
-	cout << "\n Allocation error for IntArray " << endl;
-	exit(0);
+    cout << "\n Allocation error for IntArray " << endl;
+    exit(0);
 }
 
-IntArray::IntArray(const int d1,const int d2)
+IntArray::IntArray(const int d1, const int d2)
 {
-	dim = 2;
-	bound1 = (d1 <= 0) ? 1 : d1;
-	bound2 = (d2 <= 0) ? 1 : d2;
-	bound3 = bound4 = 0;
-	size = bound1 * bound2;
-	ptr = new int[size];
-	assert( ptr != 0);
-	++arrayCount;
+    dim = 2;
+    bound1 = (d1 <= 0) ? 1 : d1;
+    bound2 = (d2 <= 0) ? 1 : d2;
+    bound3 = bound4 = 0;
+    size = bound1 * bound2;
+    ptr = new int[size];
+    assert(ptr != 0);
+    ++arrayCount;
 }
 
-IntArray::IntArray(const int d1,const int d2,const int d3)
+IntArray::IntArray(const int d1, const int d2, const int d3)
 {
-	dim = 3;
-	bound1 = (d1 <= 0) ? 1 : d1;
-	bound2 = (d2 <= 0) ? 1 : d2;
-	bound3 = (d3 <= 0) ? 1 : d3;
-	bound4 = 0;
-	set_new_handler(IntArrayAlloc);
-	size = bound1 * bound2 * bound3 ;	//* sizeof(float);
-	ptr = new int[size]();assert(ptr != 0);
-	zero_out();
-	++arrayCount;
+    dim = 3;
+    bound1 = (d1 <= 0) ? 1 : d1;
+    bound2 = (d2 <= 0) ? 1 : d2;
+    bound3 = (d3 <= 0) ? 1 : d3;
+    bound4 = 0;
+    set_new_handler(IntArrayAlloc);
+    size = bound1 * bound2 * bound3; //* sizeof(float);
+    ptr = new int[size]();
+    assert(ptr != 0);
+    zero_out();
+    ++arrayCount;
 }
 
-IntArray::IntArray(const int d1,const int d2,const int d3,const int d4)
+IntArray::IntArray(const int d1, const int d2, const int d3, const int d4)
 {
-	dim = 4;
-	bound1 = (d1 <= 0) ? 1 : d1;
-	bound2 = (d2 <= 0) ? 1 : d2;
-	bound3 = (d3 <= 0) ? 1 : d3;
-	bound4 = (d4 <= 0) ? 1 : d4;
-	set_new_handler(IntArrayAlloc);
-	size = bound1 * bound2 * bound3 * bound4 ;	//* sizeof(float);
-	ptr = new int[size]();assert(ptr != 0);
-	zero_out();
-	++arrayCount;
+    dim = 4;
+    bound1 = (d1 <= 0) ? 1 : d1;
+    bound2 = (d2 <= 0) ? 1 : d2;
+    bound3 = (d3 <= 0) ? 1 : d3;
+    bound4 = (d4 <= 0) ? 1 : d4;
+    set_new_handler(IntArrayAlloc);
+    size = bound1 * bound2 * bound3 * bound4; //* sizeof(float);
+    ptr = new int[size]();
+    assert(ptr != 0);
+    zero_out();
+    ++arrayCount;
 }
 
-IntArray::IntArray(const int d1,const int d2,const int d3,
-		const int d4,const int d5)
+IntArray::IntArray(const int d1, const int d2, const int d3, const int d4, const int d5)
 {
-	dim = 5;
-	bound1 = (d1 <= 0) ? 1 : d1;
-	bound2 = (d2 <= 0) ? 1 : d2;
-	bound3 = (d3 <= 0) ? 1 : d3;
-	bound4 = (d4 <= 0) ? 1 : d4;
-	bound5 = (d5 <= 0) ? 1 : d5;
-	set_new_handler(IntArrayAlloc);
-	size = bound1 * bound2 * bound3 * bound4 * bound5;
-	ptr = new int[size]();assert(ptr != 0);
-	zero_out();
-	++arrayCount;
+    dim = 5;
+    bound1 = (d1 <= 0) ? 1 : d1;
+    bound2 = (d2 <= 0) ? 1 : d2;
+    bound3 = (d3 <= 0) ? 1 : d3;
+    bound4 = (d4 <= 0) ? 1 : d4;
+    bound5 = (d5 <= 0) ? 1 : d5;
+    set_new_handler(IntArrayAlloc);
+    size = bound1 * bound2 * bound3 * bound4 * bound5;
+    ptr = new int[size]();
+    assert(ptr != 0);
+    zero_out();
+    ++arrayCount;
 }
 
-IntArray::IntArray(const int d1,const int d2,const int d3,
-		const int d4,const int d5,const int d6)
+IntArray::IntArray(const int d1, const int d2, const int d3, const int d4, const int d5, const int d6)
 {
-	dim = 6;
-	bound1 = (d1 <= 0) ? 1 : d1;
+    dim = 6;
+    bound1 = (d1 <= 0) ? 1 : d1;
     bound2 = (d2 <= 0) ? 1 : d2;
     bound3 = (d3 <= 0) ? 1 : d3;
     bound4 = (d4 <= 0) ? 1 : d4;
     bound5 = (d5 <= 0) ? 1 : d5;
-	bound6 = (d6 <= 0) ? 1 : d6;
+    bound6 = (d6 <= 0) ? 1 : d6;
     set_new_handler(IntArrayAlloc);
     size = bound1 * bound2 * bound3 * bound4 * bound5 * bound6;
-	ptr = new int[size]();assert(ptr != 0);
-	zero_out();
-	++arrayCount;
+    ptr = new int[size]();
+    assert(ptr != 0);
+    zero_out();
+    ++arrayCount;
 }
 
 //********************************
 // Destructor for class IntArray
 //********************************
-IntArray ::~IntArray()
+IntArray::~IntArray()
 {
     freemem();
 }
 
 void IntArray::freemem()
 {
-	delete [] ptr;
-	ptr = NULL;
+    delete[] ptr;
+    ptr = NULL;
 }
 
-void IntArray::create(const int d1,const int d2,const int d3,const int d4,const int d5,const int d6)
+void IntArray::create(const int d1, const int d2, const int d3, const int d4, const int d5, const int d6)
 {
-	size = d1 * d2 * d3 * d4 * d5 * d6;assert(size>0);
-	dim = 6;
-	bound1 = d1;bound2 = d2;bound3 = d3;bound4 = d4;bound5 = d5;bound6 = d6;
-	delete[] ptr; ptr = new int[size];
-	assert(ptr != 0);zero_out();
+    size = d1 * d2 * d3 * d4 * d5 * d6;
+    assert(size > 0);
+    dim = 6;
+    bound1 = d1;
+    bound2 = d2;
+    bound3 = d3;
+    bound4 = d4;
+    bound5 = d5;
+    bound6 = d6;
+    delete[] ptr;
+    ptr = new int[size];
+    assert(ptr != 0);
+    zero_out();
 }
 
-void IntArray::create(const int d1,const int d2,const int d3,const int d4,const int d5)
+void IntArray::create(const int d1, const int d2, const int d3, const int d4, const int d5)
 {
-	size = d1 * d2 * d3 * d4 * d5;assert(size>0);
-	dim = 5;
-	bound1 = d1;bound2 = d2;bound3 = d3;bound4 = d4;bound5 = d5;
-	delete[] ptr; ptr = new int[size];
-	assert(ptr != 0);zero_out();
+    size = d1 * d2 * d3 * d4 * d5;
+    assert(size > 0);
+    dim = 5;
+    bound1 = d1;
+    bound2 = d2;
+    bound3 = d3;
+    bound4 = d4;
+    bound5 = d5;
+    delete[] ptr;
+    ptr = new int[size];
+    assert(ptr != 0);
+    zero_out();
 }
 
-void IntArray::create(const int d1,const int d2,const int d3,const int d4)
+void IntArray::create(const int d1, const int d2, const int d3, const int d4)
 {
-	size = d1 * d2 * d3 * d4;assert(size>0);
-	dim = 4;
-	bound1 = d1;bound2 = d2;bound3 = d3;bound4 = d4;
-	delete[] ptr; ptr = new int[size];
-	assert(ptr != 0);zero_out();
+    size = d1 * d2 * d3 * d4;
+    assert(size > 0);
+    dim = 4;
+    bound1 = d1;
+    bound2 = d2;
+    bound3 = d3;
+    bound4 = d4;
+    delete[] ptr;
+    ptr = new int[size];
+    assert(ptr != 0);
+    zero_out();
 }
 
-void IntArray::create(const int d1,const int d2,const int d3)
+void IntArray::create(const int d1, const int d2, const int d3)
 {
-	size = d1 * d2 * d3;assert(size>0);
-	dim = 3;
-	bound1 = d1;bound2 = d2;bound3 = d3;bound4 = 1;
-	delete [] ptr;ptr = new int[size];
-	assert(ptr != 0);zero_out();
+    size = d1 * d2 * d3;
+    assert(size > 0);
+    dim = 3;
+    bound1 = d1;
+    bound2 = d2;
+    bound3 = d3;
+    bound4 = 1;
+    delete[] ptr;
+    ptr = new int[size];
+    assert(ptr != 0);
+    zero_out();
 }
 
 void IntArray::create(const int d1, const int d2)
 {
-	size = d1 * d2;assert(size>0);
-	dim = 2;
-	bound1 = d1;bound2 = d2;bound3 = bound4 = 1;
-	delete[] ptr;ptr = new int[size];
-	assert(ptr !=0 );zero_out();
+    size = d1 * d2;
+    assert(size > 0);
+    dim = 2;
+    bound1 = d1;
+    bound2 = d2;
+    bound3 = bound4 = 1;
+    delete[] ptr;
+    ptr = new int[size];
+    assert(ptr != 0);
+    zero_out();
 }
 
-const IntArray &IntArray::operator=(const IntArray &right)
+const IntArray& IntArray::operator=(const IntArray& right)
 {
-	for (int i = 0;i < size;i++) ptr[i] = right.ptr[i];
-	return *this;// enables x = y = z;
+    for (int i = 0; i < size; i++)
+        ptr[i] = right.ptr[i];
+    return *this; // enables x = y = z;
 }
 
-const IntArray &IntArray::operator=(const int &value)
+const IntArray& IntArray::operator=(const int& value)
 {
-	for (int i = 0;i < size;i++) ptr[i] = value;
-	return *this;// enables x = y = z;
+    for (int i = 0; i < size; i++)
+        ptr[i] = value;
+    return *this; // enables x = y = z;
 }
 
 //********************************************************
 // overloaded subscript operator for const Int Array
 // const reference return create an cvakue
 //********************************************************
-const int &IntArray::operator()
-(const int ind1,const int ind2)const
-{return ptr[ ind1 * bound2 + ind2 ];}
+const int& IntArray::operator()(const int ind1, const int ind2) const
+{
+    return ptr[ind1 * bound2 + ind2];
+}
 
-const int &IntArray::operator()
-(const int ind1,const int ind2,const int ind3)const
-{return ptr[ (ind1 * bound2 + ind2) * bound3 + ind3 ];}
+const int& IntArray::operator()(const int ind1, const int ind2, const int ind3) const
+{
+    return ptr[(ind1 * bound2 + ind2) * bound3 + ind3];
+}
 
-const int &IntArray::operator()
-(const int ind1,const int ind2,const int ind3,const int ind4)const
-{return ptr[ ((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4 ];}
+const int& IntArray::operator()(const int ind1, const int ind2, const int ind3, const int ind4) const
+{
+    return ptr[((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4];
+}
 
-const int &IntArray::operator()
-(const int ind1,const int ind2,const int ind3,const int ind4,const int ind5)const
-{return ptr[ (((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4) * bound5 + ind5 ];}
+const int& IntArray::operator()(const int ind1, const int ind2, const int ind3, const int ind4, const int ind5) const
+{
+    return ptr[(((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4) * bound5 + ind5];
+}
 
-const int &IntArray::operator()
-(const int ind1,const int ind2,const int ind3,const int ind4,const int ind5,const int ind6)const
-{return ptr[ ((((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4) * bound5 + ind5) * bound6 + ind6 ];}
+const int& IntArray::operator()(const int ind1,
+                                const int ind2,
+                                const int ind3,
+                                const int ind4,
+                                const int ind5,
+                                const int ind6) const
+{
+    return ptr[((((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4) * bound5 + ind5) * bound6 + ind6];
+}
 
 //********************************************************
 // overloaded subscript operator for non-const Int Array
 // const reference return creates an lvakue
 //********************************************************
-int &IntArray::operator()(const int ind1,const int ind2)
-{return ptr[ind1 * bound2 + ind2];}
+int& IntArray::operator()(const int ind1, const int ind2)
+{
+    return ptr[ind1 * bound2 + ind2];
+}
 
-int &IntArray::operator()(const int ind1,const int ind2,const int ind3)
-{return ptr[ (ind1 * bound2 + ind2) * bound3 + ind3 ];}
+int& IntArray::operator()(const int ind1, const int ind2, const int ind3)
+{
+    return ptr[(ind1 * bound2 + ind2) * bound3 + ind3];
+}
 
-int &IntArray::operator()(const int ind1,const int ind2,const int ind3,const int ind4)
-{return ptr[ ((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4 ];}
+int& IntArray::operator()(const int ind1, const int ind2, const int ind3, const int ind4)
+{
+    return ptr[((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4];
+}
 
-int &IntArray::operator()
-(const int ind1,const int ind2,const int ind3,const int ind4,const int ind5)
-{return ptr[ (((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4) * bound5 + ind5 ];}
+int& IntArray::operator()(const int ind1, const int ind2, const int ind3, const int ind4, const int ind5)
+{
+    return ptr[(((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4) * bound5 + ind5];
+}
 
-int &IntArray::operator()
-(const int ind1,const int ind2,const int ind3,const int ind4,const int ind5,const int ind6)
-{return ptr[ ((((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4) * bound5 + ind5) * bound6 + ind6 ];}
+int& IntArray::operator()(const int ind1,
+                          const int ind2,
+                          const int ind3,
+                          const int ind4,
+                          const int ind5,
+                          const int ind6)
+{
+    return ptr[((((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4) * bound5 + ind5) * bound6 + ind6];
+}
 
 //****************************
 // zeroes out the whole array
 //****************************
 void IntArray::zero_out(void)
 {
-	if (size <= 0) return;
-	for (int i = 0;i < size; i++) ptr[i] = 0;
-	return;
+    if (size <= 0)
+        return;
+    for (int i = 0; i < size; i++)
+        ptr[i] = 0;
+    return;
 }
diff --git a/tools/01_NAO_generation/SIAB/src_tools/intarray.h b/tools/01_NAO_generation/SIAB/src_tools/intarray.h
index 3318578da1a..de929211858 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/intarray.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/intarray.h
@@ -5,68 +5,94 @@
 #ifndef INTARRAY_H
 #define INTARRAY_H
 
-#include <iostream>
+#include <cassert>
 #include <fstream>
 #include <iomanip>
-#include <cassert>
+#include <iostream>
 
 using namespace std;
 
 class IntArray
 {
-public:
-	int *ptr;
+  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);
+    // 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();
+    ~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);
+    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);
+    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);
+    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;
+    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);
+    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;}
+    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;
+    }
 
-	int getArrayCount(void)
-	{ return arrayCount;}
+    int getArrayCount(void)
+    {
+        return arrayCount;
+    }
 
-private:
-	int size;
-	int dim;
-	int bound1, bound2, bound3, bound4, bound5, bound6;
-	static int arrayCount;
-	void freemem();
+  private:
+    int size;
+    int dim;
+    int bound1, bound2, bound3, bound4, bound5, bound6;
+    static int arrayCount;
+    void freemem();
 };
 
-#endif	// IntArray class
+#endif // IntArray class
diff --git a/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix.cpp b/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix.cpp
index ab9c73f635f..1d1c92bdc09 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix.cpp
+++ b/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix.cpp
@@ -1,62 +1,66 @@
 #include "inverse_matrix.h"
+
 #include "../src_spillage/Inverse_Matrix_S.h"
 #include "complexmatrix_inline.h"
 
-Inverse_Matrix::Inverse_Matrix(){}
-Inverse_Matrix::~Inverse_Matrix(){}
-void Inverse_Matrix::init(const int &dim_in)
+Inverse_Matrix::Inverse_Matrix()
+{
+}
+Inverse_Matrix::~Inverse_Matrix()
+{
+}
+void Inverse_Matrix::init(const int& dim_in)
 {
     this->dim = dim_in;
     this->e = new double[dim];
-    this->lwork = 2*dim;
+    this->lwork = 2 * dim;
     this->work2 = new std::complex<double>[lwork];
-    this->rwork = new double[3*dim-2];
+    this->rwork = new double[3 * dim - 2];
     this->info = 0;
     this->A.create(dim, dim);
     this->EA.create(dim, dim);
     return;
 }
 
-void Inverse_Matrix::using_zheev( const ComplexMatrix &Sin, ComplexMatrix &Sout)
+void Inverse_Matrix::using_zheev(const ComplexMatrix& Sin, ComplexMatrix& Sout)
 {
-    timer::tick("Inverse","using_zheev");
+    timer::tick("Inverse", "using_zheev");
     this->A = Sin;
 
-    //LapackConnector::zhegv( 1, 'V', 'U', nwan , A ,  nwan , B , nwan , e, work2 , 80 , rwork , &info );
+    // LapackConnector::zhegv( 1, 'V', 'U', nwan , A ,  nwan , B , nwan , e, work2 , 80 , rwork , &info );
     LapackConnector::zheev('V', 'U', dim, this->A, dim, e, work2, lwork, rwork, &info);
 
-    for(int i=0; i<dim; i++)
+    for (int i = 0; i < dim; i++)
     {
-        for(int j=0; j<dim; j++)
+        for (int j = 0; j < dim; j++)
         {
-            EA(i,j)= conj( this->A(j,i) ) / e[i] ;
+            EA(i, j) = conj(this->A(j, i)) / e[i];
         }
     }
 
     Sout = this->A * this->EA;
-    timer::tick("Inverse","using_zheev");
+    timer::tick("Inverse", "using_zheev");
     return;
 }
 
-void Inverse_Matrix::using_zpotrf( const ComplexMatrix &Sin)
+void Inverse_Matrix::using_zpotrf(const ComplexMatrix& Sin)
 {
-//  timer::tick("Inverse","using_zpotrf");
+    //  timer::tick("Inverse","using_zpotrf");
 
-    for(int i=0; i<dim; i++)
+    for (int i = 0; i < dim; i++)
     {
-        for(int j=i; j<dim; j++)
+        for (int j = i; j < dim; j++)
         {
-            A(i,j) = Sin(i,j);
+            A(i, j) = Sin(i, j);
         }
     }
 
-    LapackConnector::zpotrf('U',dim,A,dim,info);	
-//  cout << "\n info=" << info << endl;
-    LapackConnector::zpotri('U',dim,A,dim,info);
-//  cout << "\n info=" << info << endl;
-    //PRINTCM("A_zpotrf",A);
+    LapackConnector::zpotrf('U', dim, A, dim, info);
+    //  cout << "\n info=" << info << endl;
+    LapackConnector::zpotri('U', dim, A, dim, info);
+    //  cout << "\n info=" << info << endl;
+    // PRINTCM("A_zpotrf",A);
 
-//  timer::tick("Inverse","using_zpotrf");
+    //  timer::tick("Inverse","using_zpotrf");
     return;
 }
-
diff --git a/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix.h b/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix.h
index 17007e0cbdc..461aab9e6ce 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix.h
@@ -5,21 +5,21 @@
 
 class Inverse_Matrix
 {
-    public:
+  public:
     ComplexMatrix A;
-    void using_zpotrf( const ComplexMatrix &Sin);
+    void using_zpotrf(const ComplexMatrix& Sin);
 
-    void using_zheev(const ComplexMatrix &in, ComplexMatrix &out);
-    void init( const int &dim_in);
+    void using_zheev(const ComplexMatrix& in, ComplexMatrix& out);
+    void init(const int& dim_in);
 
     Inverse_Matrix();
     ~Inverse_Matrix();
 
-    private:
+  private:
     int dim;
-    double *e;
+    double* e;
     int lwork;
-    std::complex<double> *work2;
+    std::complex<double>* work2;
     double* rwork;
     int info;
 
diff --git a/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_iterate.cpp b/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_iterate.cpp
index 448b818233e..8a897dac7aa 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_iterate.cpp
+++ b/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_iterate.cpp
@@ -1,6 +1,3 @@
 #include "inverse_matrix_iterate.h"
-#include "inverse_matrix_iterate_inline.h"
-
-
-
 
+#include "inverse_matrix_iterate_inline.h"
diff --git a/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_iterate.h b/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_iterate.h
index 19ef02ee38d..18cfefa95eb 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_iterate.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_iterate.h
@@ -3,53 +3,88 @@
 
 #include "complexmatrix.h"
 
-#include<cstddef>
-#include<exception>
-#include<limits>
+#include <cstddef>
+#include <exception>
+#include <limits>
 
 // input A, X_0
 // cal A_I
 // iterate: X_{n+1} = X_n (2I - A*X_n )
- 
+
 class Inverse_Matrix_Iterate
 {
-	public:
-	inline void init(size_t dim_in);
-	inline Inverse_Matrix_Iterate &set_X0(const ComplexMatrix &X0_in);
-	Inverse_Matrix_Iterate &set_A(ComplexMatrix &A_in) { this->A_ptr=&A_in; return *this;};
-	inline const ComplexMatrix &cal_inverse();
-	const ComplexMatrix &get_inverse() { return *this->inverse_ptr; }
-
-	private:
-	inline void iterate(const ComplexMatrix &old_X, ComplexMatrix &new_X);	
-	
-	private:
-	ComplexMatrix *A_ptr = NULL;
-	ComplexMatrix X_store_1;
-	ComplexMatrix X_store_2;
-	const ComplexMatrix *X0_ptr = NULL;
-	const ComplexMatrix *inverse_ptr = NULL;
-	ComplexMatrix identity;
-	ComplexMatrix matrix_tmp_1;
-	ComplexMatrix matrix_tmp_2;
-	ComplexMatrix X0_tmp;
-	size_t iterate_num{20};
-	double tolerate_error2{1E-20};
-	double fluctrate_error2{0.0};
-	double last_error2{std::numeric_limits<double>::max()};
-	struct Matrix_Info{ size_t dim; char symmetry{NULL}; char uplo{NULL};} matrix_info;
-
-	public:
-	Inverse_Matrix_Iterate &set_iterate_num(size_t iterate_num_in) { this->iterate_num=iterate_num_in; return *this;}
-	Inverse_Matrix_Iterate &set_tolerate_error(double tolerate_error_in) { this->tolerate_error2=tolerate_error_in*tolerate_error_in; return *this;}
-	Inverse_Matrix_Iterate &set_fluctrate_error(double fluctrate_error_in) { this->fluctrate_error2=fluctrate_error_in*fluctrate_error_in; return *this;}
-	Inverse_Matrix_Iterate &set_symmetric_info(char symmetry_in, char uplo_in) { this->matrix_info.symmetry = symmetry_in; this->matrix_info.uplo = uplo_in; return *this;}
-	
-	private:
-	class Little_Error: public exception {};
-	public:
-	class Unstable: public exception {};
-	private:
+  public:
+    inline void init(size_t dim_in);
+    inline Inverse_Matrix_Iterate& set_X0(const ComplexMatrix& X0_in);
+    Inverse_Matrix_Iterate& set_A(ComplexMatrix& A_in)
+    {
+        this->A_ptr = &A_in;
+        return *this;
+    };
+    inline const ComplexMatrix& cal_inverse();
+    const ComplexMatrix& get_inverse()
+    {
+        return *this->inverse_ptr;
+    }
+
+  private:
+    inline void iterate(const ComplexMatrix& old_X, ComplexMatrix& new_X);
+
+  private:
+    ComplexMatrix* A_ptr = NULL;
+    ComplexMatrix X_store_1;
+    ComplexMatrix X_store_2;
+    const ComplexMatrix* X0_ptr = NULL;
+    const ComplexMatrix* inverse_ptr = NULL;
+    ComplexMatrix identity;
+    ComplexMatrix matrix_tmp_1;
+    ComplexMatrix matrix_tmp_2;
+    ComplexMatrix X0_tmp;
+    size_t iterate_num{20};
+    double tolerate_error2{1E-20};
+    double fluctrate_error2{0.0};
+    double last_error2{std::numeric_limits<double>::max()};
+    struct Matrix_Info
+    {
+        size_t dim;
+        char symmetry{NULL};
+        char uplo{NULL};
+    } matrix_info;
+
+  public:
+    Inverse_Matrix_Iterate& set_iterate_num(size_t iterate_num_in)
+    {
+        this->iterate_num = iterate_num_in;
+        return *this;
+    }
+    Inverse_Matrix_Iterate& set_tolerate_error(double tolerate_error_in)
+    {
+        this->tolerate_error2 = tolerate_error_in * tolerate_error_in;
+        return *this;
+    }
+    Inverse_Matrix_Iterate& set_fluctrate_error(double fluctrate_error_in)
+    {
+        this->fluctrate_error2 = fluctrate_error_in * fluctrate_error_in;
+        return *this;
+    }
+    Inverse_Matrix_Iterate& set_symmetric_info(char symmetry_in, char uplo_in)
+    {
+        this->matrix_info.symmetry = symmetry_in;
+        this->matrix_info.uplo = uplo_in;
+        return *this;
+    }
+
+  private:
+    class Little_Error : public exception
+    {
+    };
+
+  public:
+    class Unstable : public exception
+    {
+    };
+
+  private:
 };
 
 #endif
diff --git a/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_iterate_bak.cpp b/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_iterate_bak.cpp
index 6ba253892b7..4409006abd9 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_iterate_bak.cpp
+++ b/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_iterate_bak.cpp
@@ -2,79 +2,95 @@
 
 void Inverse_Matrix_Iterate::init(size_t dim_in)
 {
-	this->matrix_info.dim = dim_in;
-	this->X_store_1.create(dim_in,dim_in);
-	this->X_store_2.create(dim_in,dim_in);
-	
-	this->X0_tmp.create(dim_in,dim_in);
-	this->identity.create(dim_in,dim_in);
-	this->identity.set_as_identity_matrix();
-	this->error_matrix.create(dim_in,dim_in);
+    this->matrix_info.dim = dim_in;
+    this->X_store_1.create(dim_in, dim_in);
+    this->X_store_2.create(dim_in, dim_in);
+
+    this->X0_tmp.create(dim_in, dim_in);
+    this->identity.create(dim_in, dim_in);
+    this->identity.set_as_identity_matrix();
+    this->error_matrix.create(dim_in, dim_in);
 }
 
-Inverse_Matrix_Iterate &Inverse_Matrix_Iterate::set_X0(const ComplexMatrix &X0_in)
+Inverse_Matrix_Iterate& Inverse_Matrix_Iterate::set_X0(const ComplexMatrix& X0_in)
 {
-	if( &X0_in == &this->X_store_1 )		// without this will cause a bug when istep==1
-	{
-		this->X0_ptr = &this->X0_tmp;
-		this->X0_tmp = X0_in;
-	}
-	else
-	{
-		this->X0_ptr = &X0_in;
-	}
-	return *this;
+    if (&X0_in == &this->X_store_1) // without this will cause a bug when istep==1
+    {
+        this->X0_ptr = &this->X0_tmp;
+        this->X0_tmp = X0_in;
+    }
+    else
+    {
+        this->X0_ptr = &X0_in;
+    }
+    return *this;
 }
 
-void Inverse_Matrix_Iterate::iterate(const ComplexMatrix &old_X, ComplexMatrix &new_X)
+void Inverse_Matrix_Iterate::iterate(const ComplexMatrix& old_X, ComplexMatrix& new_X)
 {
-	// I - A * X_n
-	if(this->matrix_info.symmetry) 
-	{ this->error_matrix = this->identity - multiply_special(this->matrix_info.symmetry, 'L', this->matrix_info.uplo, *(this->A_ptr), old_X); }
-	else 
-	{ this->error_matrix = this->identity - *(this->A_ptr) * old_X; }
-	
-	const double error_norm(sqrt(abs2(this->error_matrix)));
-	if ( error_norm < this->tolerate_error)
-	{
-		this->inverse_ptr = &old_X;
-		throw Little_Error();
-	}
-	else if ( error_norm > this->last_error + this->fluctrate_error)
-	{	throw Unstable();	}
-	else { this->last_error = error_norm; }
-	
-	// X_{n+1} = X_n (I + (I-A*X_n) ) = X_n ( 2I - A*X_n )
-	for( size_t i(0); i<this->matrix_info.dim; ++i) { this->error_matrix(i,i) += 1.0; }
-	new_X = old_X * this->error_matrix;
-//	new_X = old_X * (this->identity + this->error_matrix);
+    // I - A * X_n
+    if (this->matrix_info.symmetry)
+    {
+        this->error_matrix
+            = this->identity
+              - multiply_special(this->matrix_info.symmetry, 'L', this->matrix_info.uplo, *(this->A_ptr), old_X);
+    }
+    else
+    {
+        this->error_matrix = this->identity - *(this->A_ptr) * old_X;
+    }
+
+    const double error_norm(sqrt(abs2(this->error_matrix)));
+    if (error_norm < this->tolerate_error)
+    {
+        this->inverse_ptr = &old_X;
+        throw Little_Error();
+    }
+    else if (error_norm > this->last_error + this->fluctrate_error)
+    {
+        throw Unstable();
+    }
+    else
+    {
+        this->last_error = error_norm;
+    }
+
+    // X_{n+1} = X_n (I + (I-A*X_n) ) = X_n ( 2I - A*X_n )
+    for (size_t i(0); i < this->matrix_info.dim; ++i)
+    {
+        this->error_matrix(i, i) += 1.0;
+    }
+    new_X = old_X * this->error_matrix;
+    //	new_X = old_X * (this->identity + this->error_matrix);
 
-	this->inverse_ptr = &new_X;
+    this->inverse_ptr = &new_X;
 }
 
-const ComplexMatrix &Inverse_Matrix_Iterate::cal_inverse()
+const ComplexMatrix& Inverse_Matrix_Iterate::cal_inverse()
 {
-	try
-	{
-		this->inverse_ptr = this->X0_ptr;
-		this->last_error = std::numeric_limits<double>::max();
-		if(this->iterate_num >= 1)
-		{
-			iterate( *this->X0_ptr, this->X_store_1 );
-		}
-		for( size_t istep(2); istep<=this->iterate_num; ++istep)
-		{
-			if( istep&1 )		// istep odd, X_store_2 -> X_store_1
-			{
-				iterate( this->X_store_2, this->X_store_1 );
-			}
-			else				// istep even, X_store_1 -> X_store_2
-			{
-				iterate( this->X_store_1, this->X_store_2);
-			}	
-		}		
-	}
-	catch(const Little_Error &e){}
+    try
+    {
+        this->inverse_ptr = this->X0_ptr;
+        this->last_error = std::numeric_limits<double>::max();
+        if (this->iterate_num >= 1)
+        {
+            iterate(*this->X0_ptr, this->X_store_1);
+        }
+        for (size_t istep(2); istep <= this->iterate_num; ++istep)
+        {
+            if (istep & 1) // istep odd, X_store_2 -> X_store_1
+            {
+                iterate(this->X_store_2, this->X_store_1);
+            }
+            else // istep even, X_store_1 -> X_store_2
+            {
+                iterate(this->X_store_1, this->X_store_2);
+            }
+        }
+    }
+    catch (const Little_Error& e)
+    {
+    }
 
-	return get_inverse();
+    return get_inverse();
 }
\ No newline at end of file
diff --git a/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_iterate_inline.h b/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_iterate_inline.h
index 6f7457d2cfc..9aae1c7791c 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_iterate_inline.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_iterate_inline.h
@@ -1,87 +1,106 @@
 #ifndef INVERSE_MATRIX_ITERATE_INLINE_H
 #define INVERSE_MATRIX_ITERATE_INLINE_H
 
-#include "inverse_matrix_iterate.h"
 #include "complexmatrix.h"
 #include "complexmatrix_inline.h"
+#include "inverse_matrix_iterate.h"
 
 inline void Inverse_Matrix_Iterate::init(size_t dim_in)
 {
-	this->matrix_info.dim = dim_in;
-	this->X_store_1.create(dim_in,dim_in);
-	this->X_store_2.create(dim_in,dim_in);
-	
-	this->X0_tmp.create(dim_in,dim_in);
-	this->identity.create(dim_in,dim_in);
-	this->identity.set_as_identity_matrix();
-	this->matrix_tmp_1.create(dim_in,dim_in);
-	this->matrix_tmp_2.create(dim_in,dim_in);
+    this->matrix_info.dim = dim_in;
+    this->X_store_1.create(dim_in, dim_in);
+    this->X_store_2.create(dim_in, dim_in);
+
+    this->X0_tmp.create(dim_in, dim_in);
+    this->identity.create(dim_in, dim_in);
+    this->identity.set_as_identity_matrix();
+    this->matrix_tmp_1.create(dim_in, dim_in);
+    this->matrix_tmp_2.create(dim_in, dim_in);
 }
 
-inline Inverse_Matrix_Iterate &Inverse_Matrix_Iterate::set_X0(const ComplexMatrix &X0_in)
+inline Inverse_Matrix_Iterate& Inverse_Matrix_Iterate::set_X0(const ComplexMatrix& X0_in)
 {
-	if( &X0_in == &this->X_store_1 )		// without this will cause a bug when istep==1
-	{
-		this->X0_ptr = &this->X0_tmp;
-		this->X0_tmp = X0_in;
-	}
-	else
-	{
-		this->X0_ptr = &X0_in;
-	}
-	return *this;
+    if (&X0_in == &this->X_store_1) // without this will cause a bug when istep==1
+    {
+        this->X0_ptr = &this->X0_tmp;
+        this->X0_tmp = X0_in;
+    }
+    else
+    {
+        this->X0_ptr = &X0_in;
+    }
+    return *this;
 }
 
-inline void Inverse_Matrix_Iterate::iterate(const ComplexMatrix &old_X, ComplexMatrix &new_X)
+inline void Inverse_Matrix_Iterate::iterate(const ComplexMatrix& old_X, ComplexMatrix& new_X)
 {
-	// I - A * X_n
-	if(this->matrix_info.symmetry) 
-	{ this->matrix_tmp_1.minus( this->identity, this->matrix_tmp_2.multiply_special(this->matrix_info.symmetry, 'L', this->matrix_info.uplo, *(this->A_ptr), old_X) ); }
-	else 
-	{ this->matrix_tmp_1.minus( this->identity, this->matrix_tmp_2.multiply(*(this->A_ptr), old_X ) ); }
-	
-	const double error_norm2(abs2(this->matrix_tmp_1));
-	if ( error_norm2 < this->tolerate_error2)
-	{
-		this->inverse_ptr = &old_X;
-		throw Little_Error();
-	}
-	else if ( error_norm2 > this->last_error2 + this->fluctrate_error2)
-	{	throw Unstable();	}
-	else { this->last_error2 = error_norm2; }
-	
-	// X_{n+1} = X_n (I + (I-A*X_n) ) = X_n ( 2I - A*X_n )
-	for( size_t i(0); i<this->matrix_info.dim; ++i) { this->matrix_tmp_1(i,i) += 1.0; }
-	new_X.multiply(old_X, this->matrix_tmp_1);
-//	new_X = old_X * (this->identity + this->matrix_tmp_1);
+    // I - A * X_n
+    if (this->matrix_info.symmetry)
+    {
+        this->matrix_tmp_1.minus(this->identity,
+                                 this->matrix_tmp_2.multiply_special(this->matrix_info.symmetry,
+                                                                     'L',
+                                                                     this->matrix_info.uplo,
+                                                                     *(this->A_ptr),
+                                                                     old_X));
+    }
+    else
+    {
+        this->matrix_tmp_1.minus(this->identity, this->matrix_tmp_2.multiply(*(this->A_ptr), old_X));
+    }
+
+    const double error_norm2(abs2(this->matrix_tmp_1));
+    if (error_norm2 < this->tolerate_error2)
+    {
+        this->inverse_ptr = &old_X;
+        throw Little_Error();
+    }
+    else if (error_norm2 > this->last_error2 + this->fluctrate_error2)
+    {
+        throw Unstable();
+    }
+    else
+    {
+        this->last_error2 = error_norm2;
+    }
+
+    // X_{n+1} = X_n (I + (I-A*X_n) ) = X_n ( 2I - A*X_n )
+    for (size_t i(0); i < this->matrix_info.dim; ++i)
+    {
+        this->matrix_tmp_1(i, i) += 1.0;
+    }
+    new_X.multiply(old_X, this->matrix_tmp_1);
+    //	new_X = old_X * (this->identity + this->matrix_tmp_1);
 
-	this->inverse_ptr = &new_X;
+    this->inverse_ptr = &new_X;
 }
 
-inline const ComplexMatrix &Inverse_Matrix_Iterate::cal_inverse()
+inline const ComplexMatrix& Inverse_Matrix_Iterate::cal_inverse()
 {
-	try
-	{
-		this->inverse_ptr = this->X0_ptr;
-		this->last_error2 = std::numeric_limits<double>::max();
-		if(this->iterate_num >= 1)
-		{
-			iterate( *this->X0_ptr, this->X_store_1 );
-		}
-		for( size_t istep(2); istep<=this->iterate_num; ++istep)
-		{
-			if( istep&1 )		// istep odd, X_store_2 -> X_store_1
-			{
-				iterate( this->X_store_2, this->X_store_1 );
-			}
-			else				// istep even, X_store_1 -> X_store_2
-			{
-				iterate( this->X_store_1, this->X_store_2);
-			}	
-		}		
-	}
-	catch(const Little_Error &e){}
+    try
+    {
+        this->inverse_ptr = this->X0_ptr;
+        this->last_error2 = std::numeric_limits<double>::max();
+        if (this->iterate_num >= 1)
+        {
+            iterate(*this->X0_ptr, this->X_store_1);
+        }
+        for (size_t istep(2); istep <= this->iterate_num; ++istep)
+        {
+            if (istep & 1) // istep odd, X_store_2 -> X_store_1
+            {
+                iterate(this->X_store_2, this->X_store_1);
+            }
+            else // istep even, X_store_1 -> X_store_2
+            {
+                iterate(this->X_store_1, this->X_store_2);
+            }
+        }
+    }
+    catch (const Little_Error& e)
+    {
+    }
 
-	return get_inverse();
+    return get_inverse();
 }
 #endif
\ No newline at end of file
diff --git a/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_perturbe.cpp b/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_perturbe.cpp
index 8726ccd7aae..043d134eb64 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_perturbe.cpp
+++ b/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_perturbe.cpp
@@ -1,26 +1,39 @@
 #include "inverse_matrix_perturbe.h"
-#include<iostream>
-ComplexMatrix &Inverse_Matrix_Perturbe::cal_inverse()
+
+#include <iostream>
+ComplexMatrix& Inverse_Matrix_Perturbe::cal_inverse()
 {
-//std::cout<<"@";
-	this->I = *this->A_inverse;
-	const ComplexMatrix b_a_inverse(*this->B * *this->A_inverse);
-	ComplexMatrix add_item(*this->A_inverse);
-	for(size_t i(1); i<=iterate_num; ++i)
-	{
-		add_item = add_item * b_a_inverse;
-		if(plus_minus_state==2) { this->I += add_item; }
-		else if(plus_minus_state==1)
-		{
-			if((i&1) == 0) { this->I += add_item; }
-			else { this->I -= add_item; }
-		}	
-		else { throw invalid_argument("set plus_minus_state"); }
-	}
-	return this->I;
+    // std::cout<<"@";
+    this->I = *this->A_inverse;
+    const ComplexMatrix b_a_inverse(*this->B * *this->A_inverse);
+    ComplexMatrix add_item(*this->A_inverse);
+    for (size_t i(1); i <= iterate_num; ++i)
+    {
+        add_item = add_item * b_a_inverse;
+        if (plus_minus_state == 2)
+        {
+            this->I += add_item;
+        }
+        else if (plus_minus_state == 1)
+        {
+            if ((i & 1) == 0)
+            {
+                this->I += add_item;
+            }
+            else
+            {
+                this->I -= add_item;
+            }
+        }
+        else
+        {
+            throw invalid_argument("set plus_minus_state");
+        }
+    }
+    return this->I;
 }
 
-void Inverse_Matrix_Perturbe::init(const int &dim)
+void Inverse_Matrix_Perturbe::init(const int& dim)
 {
     this->I.create(dim, dim);
     return;
diff --git a/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_perturbe.h b/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_perturbe.h
index 5a461edbdf3..217a5b9c9a7 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_perturbe.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/inverse_matrix_perturbe.h
@@ -9,27 +9,58 @@
 
 class Inverse_Matrix_Perturbe
 {
-    public:
-	void init(const int &dim);
-	ComplexMatrix &cal_inverse();
-	ComplexMatrix &get_I() { return I; }
-	
-	private:
-    const ComplexMatrix *A_inverse{nullptr};
-	const ComplexMatrix *B{nullptr};
-	ComplexMatrix I;
-	size_t iterate_num;
-	int plus_minus_state{0};		// 1 plus, 2 minus
-	
-	public:
-	const ComplexMatrix *get_A_inverse() { return A_inverse; }
-	Inverse_Matrix_Perturbe &set_A_inverse(ComplexMatrix *A_inverse_in) { A_inverse = A_inverse_in; return *this;}
-	const ComplexMatrix *get_B() { return B; }
-	Inverse_Matrix_Perturbe &set_B(ComplexMatrix *B_in) { B = B_in; return *this;}
-	const size_t &get_iterate_num() { return iterate_num; }
-	Inverse_Matrix_Perturbe &set_iterate_num(size_t iterate_num_in) { iterate_num = iterate_num_in; return *this;}
-	const int &get_plus_minus_state() { return plus_minus_state; }
-	Inverse_Matrix_Perturbe &set_plus_minus_state(int plus_minus_state_in) { plus_minus_state = plus_minus_state_in; return *this;}
+  public:
+    void init(const int& dim);
+    ComplexMatrix& cal_inverse();
+    ComplexMatrix& get_I()
+    {
+        return I;
+    }
+
+  private:
+    const ComplexMatrix* A_inverse{nullptr};
+    const ComplexMatrix* B{nullptr};
+    ComplexMatrix I;
+    size_t iterate_num;
+    int plus_minus_state{0}; // 1 plus, 2 minus
+
+  public:
+    const ComplexMatrix* get_A_inverse()
+    {
+        return A_inverse;
+    }
+    Inverse_Matrix_Perturbe& set_A_inverse(ComplexMatrix* A_inverse_in)
+    {
+        A_inverse = A_inverse_in;
+        return *this;
+    }
+    const ComplexMatrix* get_B()
+    {
+        return B;
+    }
+    Inverse_Matrix_Perturbe& set_B(ComplexMatrix* B_in)
+    {
+        B = B_in;
+        return *this;
+    }
+    const size_t& get_iterate_num()
+    {
+        return iterate_num;
+    }
+    Inverse_Matrix_Perturbe& set_iterate_num(size_t iterate_num_in)
+    {
+        iterate_num = iterate_num_in;
+        return *this;
+    }
+    const int& get_plus_minus_state()
+    {
+        return plus_minus_state;
+    }
+    Inverse_Matrix_Perturbe& set_plus_minus_state(int plus_minus_state_in)
+    {
+        plus_minus_state = plus_minus_state_in;
+        return *this;
+    }
 };
 
 #endif
\ No newline at end of file
diff --git a/tools/01_NAO_generation/SIAB/src_tools/lapack_connector.h b/tools/01_NAO_generation/SIAB/src_tools/lapack_connector.h
index 91eafad567f..412ba6f5d06 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/lapack_connector.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/lapack_connector.h
@@ -4,7 +4,7 @@
 // File:            LapackConnector.hpp
 // Author:          sltk
 // Comment:         LapackConnector provide the connector to the fortran Lapack routine.
-// Warning:         
+// Warning:
 // Start time:      2007-03-08
 // Last modified:   2008-08-12 ywcui : add zhegvx
 // 					2008-08-13 mohan : find bug,test.
@@ -16,71 +16,60 @@
 
 #include <new>
 
-//#include <exception>
-#include <stdexcept>
-
+// #include <exception>
 #include <complex>
+#include <stdexcept>
 using namespace std;
 
 #include "complexmatrix.h"
 
-
 extern "C"
 {
 
-	void zpotrf_( const char* uplo,
-				const int* n,
-				const std::complex<double> *A,
-				const int* lda,
-				const int* infp);
-
-	void zpotri_( const char* uplo,
-				const int* n,
-				const std::complex<double> *A,
-				const int* lda,
-				const int* infp);
-
-	//=======================================================================
-	//ZGESV computes the solution to a complex system of linear equations
- 	//     A * X = B,
-	//  where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
-	//
-	//  The LU decomposition with partial pivoting and row interchanges is
-	//  used to factor A as
-	//     A = P * L * U,
-	//  where P is a permutation matrix, L is unit lower triangular, and U is
-	//  upper triangular.  The factored form of A is then used to solve the
-	//  system of equations A * X = B.
-	//
-	//  A = LDA * N
-	//  B = LDB * NRHS (exit X = B)
-	//========================================================================
-	void zgesv_(const int* n,
-				const int* nrhs,
-				const std::complex<double> *A,
-				const int *lda,
-				const int *ipiv,
-				const std::complex<double> *B,
-				const int* ldb,
-				int* info
-				);
-
-	/*zhegv computes all the eigenvalues, and optionally, the eigenvectors
-	  of a complex generalized Hermitian-definite eigenproblem, of the form
-	  A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
-	  Here A and B are assumed to be Hermitian and B is also
-	  positive definite.
-	*/
-
-	/*zhegvx computes selected eigenvalues, and optionally, eigenvectors
-	  of a complex generalized Hermitian-definite eigenproblem, of the form
-	  A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
-	  B are assumed to be Hermitian and B is also positive definite.
-	  Eigenvalues and eigenvectors can be selected by specifying either a
-	  range of values or a range of indices for the desired eigenvalues.
-	*/
-
-	void zhegv_(const int* itype,
+    void zpotrf_(const char* uplo, const int* n, const std::complex<double>* A, const int* lda, const int* infp);
+
+    void zpotri_(const char* uplo, const int* n, const std::complex<double>* A, const int* lda, const int* infp);
+
+    //=======================================================================
+    // ZGESV computes the solution to a complex system of linear equations
+    //     A * X = B,
+    //  where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
+    //
+    //  The LU decomposition with partial pivoting and row interchanges is
+    //  used to factor A as
+    //     A = P * L * U,
+    //  where P is a permutation matrix, L is unit lower triangular, and U is
+    //  upper triangular.  The factored form of A is then used to solve the
+    //  system of equations A * X = B.
+    //
+    //  A = LDA * N
+    //  B = LDB * NRHS (exit X = B)
+    //========================================================================
+    void zgesv_(const int* n,
+                const int* nrhs,
+                const std::complex<double>* A,
+                const int* lda,
+                const int* ipiv,
+                const std::complex<double>* B,
+                const int* ldb,
+                int* info);
+
+    /*zhegv computes all the eigenvalues, and optionally, the eigenvectors
+      of a complex generalized Hermitian-definite eigenproblem, of the form
+      A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
+      Here A and B are assumed to be Hermitian and B is also
+      positive definite.
+    */
+
+    /*zhegvx computes selected eigenvalues, and optionally, eigenvectors
+      of a complex generalized Hermitian-definite eigenproblem, of the form
+      A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
+      B are assumed to be Hermitian and B is also positive definite.
+      Eigenvalues and eigenvectors can be selected by specifying either a
+      range of values or a range of indices for the desired eigenvalues.
+    */
+
+    void zhegv_(const int* itype,
                 const char* jobz,
                 const char* uplo,
                 const int* n,
@@ -92,126 +81,149 @@ extern "C"
                 std::complex<double>* work,
                 int* lwork,
                 double* rwork,
-                int* info	);
-
-	void zheev_(const char* jobz,
-				const char* uplo,
-				const int* n,
-				complex<double> *a,
-				const int* lda,
-				double* w,
-				complex<double >* work,
-				const int* lwork,
-				double* rwork,
-				int* info	);
-
-	//===================================================================================
-	// zhegvx	:
-	// itype	:	1:	A*x = (lambda)*B*x
-	// 				2:	A*B*x = (lambda)*x
-	// 				3:	B*A*x = (lambda)*x
-	// jobz		:	'N':  Compute eigenvalues only
-	// 				'V':  Compute eigenvalues and eigenvectors.
-	// range	:	'A': all eigenvalues will be found.
-	//				'V': all eigenvalues in the half-open interval (VL,VU] will be found.
-	//				'I': the IL-th through IU-th eigenvalues will be found.
-	// uplo		:	'U':  Upper triangles of A and B are stored
-	// 				'L':  Lower triangles of A and B are stored
-	// N		: 	The order of the matrices A and B.  N >= 0.
-	// A		:	(input/output) COMPLEX*16 array, dimension (LDA, N)
+                int* info);
+
+    void zheev_(const char* jobz,
+                const char* uplo,
+                const int* n,
+                complex<double>* a,
+                const int* lda,
+                double* w,
+                complex<double>* work,
+                const int* lwork,
+                double* rwork,
+                int* info);
+
+    //===================================================================================
+    // zhegvx	:
+    // itype	:	1:	A*x = (lambda)*B*x
+    // 				2:	A*B*x = (lambda)*x
+    // 				3:	B*A*x = (lambda)*x
+    // jobz		:	'N':  Compute eigenvalues only
+    // 				'V':  Compute eigenvalues and eigenvectors.
+    // range	:	'A': all eigenvalues will be found.
+    //				'V': all eigenvalues in the half-open interval (VL,VU] will be found.
+    //				'I': the IL-th through IU-th eigenvalues will be found.
+    // uplo		:	'U':  Upper triangles of A and B are stored
+    // 				'L':  Lower triangles of A and B are stored
+    // N		: 	The order of the matrices A and B.  N >= 0.
+    // A		:	(input/output) COMPLEX*16 array, dimension (LDA, N)
     //            	On entry, the Hermitian matrix A.  If UPLO = 'U',  the  leading
     //           	N-by-N upper triangular part of A contains the upper triangular
     //           	part of the matrix A.  If UPLO = 'L', the leading N-by-N  lower
     //           	triangular  part of A contains the lower triangular part of the
     //          	matrix A.
-	// LDA		:	The leading dimension of the array A.  LDA >= max(1,N).
-	// LDB		:	On  entry,  the Hermitian matrix B.  If UPLO = 'U', the leading
+    // LDA		:	The leading dimension of the array A.  LDA >= max(1,N).
+    // LDB		:	On  entry,  the Hermitian matrix B.  If UPLO = 'U', the leading
     //            	N-by-N upper triangular part of B contains the upper triangular
     //           	part  of the matrix B.  If UPLO = 'L', the leading N-by-N lower
     //           	triangular part of B contains the lower triangular part of  the
     //           	matrix B.
-	// vl,vu	:	the lower and upper bounds of the interval to be searched for 
-	//				eigenvalues. VL < VU.  Not referenced if RANGE = 'A' or 'I'.
-	// il,iu	:	the indices (in ascending order) of the smallest and largest 
-	//				eigenvalues to be  returned. 1  <= IL <= IU <= N, if N > 0; 
-	//				IL = 1 and IU = 0 if N = 0.  Not referenced if RANGE = 'A' or 'V'.
-	// abstol	:	The absolute error tolerance for the eigenvalues
-	// m		:	The  total number of eigenvalues found ,0 <= M <= N.  If RANGE
-	//              = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
-	// w		:	The first  M  elements  contain  the  selected  eigenvalues  in
-	//              ascending order
-	// z		:   If  JOBZ  = 'N', then Z is not referenced.  If JOBZ = 'V', then
+    // vl,vu	:	the lower and upper bounds of the interval to be searched for
+    //				eigenvalues. VL < VU.  Not referenced if RANGE = 'A' or 'I'.
+    // il,iu	:	the indices (in ascending order) of the smallest and largest
+    //				eigenvalues to be  returned. 1  <= IL <= IU <= N, if N > 0;
+    //				IL = 1 and IU = 0 if N = 0.  Not referenced if RANGE = 'A' or 'V'.
+    // abstol	:	The absolute error tolerance for the eigenvalues
+    // m		:	The  total number of eigenvalues found ,0 <= M <= N.  If RANGE
+    //              = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
+    // w		:	The first  M  elements  contain  the  selected  eigenvalues  in
+    //              ascending order
+    // z		:   If  JOBZ  = 'N', then Z is not referenced.  If JOBZ = 'V', then
     //           	if INFO = 0, the first M columns of Z contain  the  orthonormal
     //          	eigenvectors  of  the  matrix  A  corresponding to the selected
     //          	eigenvalues, with the i-th column of Z holding the  eigenvector
     //           	associated  with W(i).  The eigenvectors are normalized as fol-
     //           	lows: if ITYPE  =  1  or  2,  Z**T*B*Z  =  I;  if  ITYPE  =  3,
     //           	Z**T*inv(B)*Z = I.
-	//				If  an  eigenvector  fails  to  converge, then that column of Z
+    //				If  an  eigenvector  fails  to  converge, then that column of Z
     //		        contains the latest approximation to the eigenvector,  and  the
     //	            index  of the eigenvector is returned in IFAIL.  Note: the user
     //              must ensure that at least max(1,M) columns are supplied in  the
     //              array  Z;  if RANGE = 'V', the exact value of M is not known in
-    //              advance and an upper bound must be used 
-	// LDZ		:	The leading dimension of the array Z.  LDZ >= 1, and if JOBZ  =
+    //              advance and an upper bound must be used
+    // LDZ		:	The leading dimension of the array Z.  LDZ >= 1, and if JOBZ  =
     //              'V', LDZ >= max(1,N).
-	// work		:	On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
-	// lwork	:	LWORK >= max(1,2*N). For optimal efficiency, LWORK >= (NB+1)*N, 
-	//				where NB is the blocksize for ZHETRD returned by ILAENV.
-	// rwork	:   dimension (7*N)
-	// iwork	:	dimension (5*N)
-	// info		:	= 0:  successful exit
-	// 				< 0:  if INFO = -i, the i-th argument had an illegal value
-	// 				> 0:  ZPOTRF or ZHEEVX returned an error code:
-	//===================================================================================
-
-	/*
-	void zhegvx_(const int* itype,
-				 const char* jobz,
-				 const char* range,	
-				 const char* uplo,
-				 const int* N,
-				 std::complex<double> *A,
-				 const int* LDA,
-				 std::complex<double> *B,
-				 const int* LDB,
-				 const double* vl,
-				 const double* vu,
-				 const int* il,
-				 const int* iu,
-				 const double* abstol,
-				 const int* m,
-				 double* w,
-				 std::complex<double> *z,
-				 const int *LDZ,
-				 std::complex<double> *work,
-				 const int* lwork,
-				 double* rwork,
-				 int* iwork,
-				 int* ifail,
-				 int* info
-				 );
-	*/
-
-	
-	// Peize Lin add 2015-11-20	
-	void daxpy_(int *N, double *alpha, double *X, int *incX, double *Y, int *incY);
-	void dscal_(int *N, double *alpha, double *X, int *incX);
-
-	// Peize Lin add 2015-12-07
-	void zgemm_(const char *TransA, const char *TransB,
-                       const int *M, const int *N, const int *K,
-					   const double *alpha, const void *A, const int *lda, const void *B, const int *ldb,
-                       const double *beta, void *C, const int *ldc);
-	void zsymm_(const char *Side, const char *Uplo,
-                       const int *M, const int *N,
-					   const double *alpha, const void *A, const int *lda, const void *B, const int *ldb,
-                       const double *beta, void *C, const int *ldc);					   
-	void zhemm_(const char *Side, const char *Uplo,
-                       const int *M, const int *N,
-					   const double *alpha, const void *A, const int *lda, const void *B, const int *ldb,
-                       const double *beta, void *C, const int *ldc);					   
-	
+    // work		:	On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+    // lwork	:	LWORK >= max(1,2*N). For optimal efficiency, LWORK >= (NB+1)*N,
+    //				where NB is the blocksize for ZHETRD returned by ILAENV.
+    // rwork	:   dimension (7*N)
+    // iwork	:	dimension (5*N)
+    // info		:	= 0:  successful exit
+    // 				< 0:  if INFO = -i, the i-th argument had an illegal value
+    // 				> 0:  ZPOTRF or ZHEEVX returned an error code:
+    //===================================================================================
+
+    /*
+    void zhegvx_(const int* itype,
+                 const char* jobz,
+                 const char* range,
+                 const char* uplo,
+                 const int* N,
+                 std::complex<double> *A,
+                 const int* LDA,
+                 std::complex<double> *B,
+                 const int* LDB,
+                 const double* vl,
+                 const double* vu,
+                 const int* il,
+                 const int* iu,
+                 const double* abstol,
+                 const int* m,
+                 double* w,
+                 std::complex<double> *z,
+                 const int *LDZ,
+                 std::complex<double> *work,
+                 const int* lwork,
+                 double* rwork,
+                 int* iwork,
+                 int* ifail,
+                 int* info
+                 );
+    */
+
+    // Peize Lin add 2015-11-20
+    void daxpy_(int* N, double* alpha, double* X, int* incX, double* Y, int* incY);
+    void dscal_(int* N, double* alpha, double* X, int* incX);
+
+    // Peize Lin add 2015-12-07
+    void zgemm_(const char* TransA,
+                const char* TransB,
+                const int* M,
+                const int* N,
+                const int* K,
+                const double* alpha,
+                const void* A,
+                const int* lda,
+                const void* B,
+                const int* ldb,
+                const double* beta,
+                void* C,
+                const int* ldc);
+    void zsymm_(const char* Side,
+                const char* Uplo,
+                const int* M,
+                const int* N,
+                const double* alpha,
+                const void* A,
+                const int* lda,
+                const void* B,
+                const int* ldb,
+                const double* beta,
+                void* C,
+                const int* ldc);
+    void zhemm_(const char* Side,
+                const char* Uplo,
+                const int* M,
+                const int* N,
+                const double* alpha,
+                const void* A,
+                const int* lda,
+                const void* B,
+                const int* ldb,
+                const double* beta,
+                void* C,
+                const int* ldc);
 };
 
 // Class LapackConnector provide the connector to fortran lapack routine.
@@ -219,268 +231,263 @@ extern "C"
 // Usage example:	LapackConnector::functionname(parameter list).
 class LapackConnector
 {
-private:
-	// Transpose the complex matrix to the fortran-form real-complex array.
-	static inline
-	complex<double>* transpose(const ComplexMatrix& a, const int n, const int lda)
-	{	
-		complex<double>* aux = new std::complex<double>[lda*n];
-		for(int i = 0; i < n; ++i)
-		{	for(int j = 0; j < lda; ++j)
-			{	
-				aux[i*lda+j] = a(j,i);		// aux[i*lda+j] means aux[i][j] in semantic, not in syntax!
-			}
-		}
-		return aux;
-	}
-
-	// Transpose the fortran-form real-complex array to the complex matrix.
-	static inline
-	void transpose(const std::complex<double>* aux, ComplexMatrix& a, const int n, const int lda)
-	{	
-		for(int i = 0; i < n; ++i)
-		{	for(int j = 0; j < lda; ++j)
-			{	
-				a(j, i) = aux[i*lda+j];		// aux[i*lda+j] means aux[i][j] in semantic, not in syntax!
-			}
-		}
-	}
-
-	static inline
-	char change_side(const char &side)
-	{
-		switch(side)
-		{
-			case 'L': return 'R';
-			case 'R': return 'L';
-			default: throw invalid_argument("Side must be 'L' or 'R'");
-		}		
-	}
-	
-	static inline
-	char change_uplo(const char &uplo)
-	{
-		switch(uplo)
-		{
-			case 'U': return 'L';
-			case 'L': return 'U';
-			default: throw invalid_argument("Uplo must be 'U' or 'L'");
-		}		
-	}
-
-
-	
-public:
-	static inline
-	void zgesv (const int n,
-				const int nrhs,
-				ComplexMatrix &A,
-				const int lda,
-				const int *ipiv,
-				complex<double> *B,
-				const int ldb,
-				int* info
-				)
-	{
-		complex<double> *aux = LapackConnector::transpose(A, n, lda);
-
-		zgesv_(&n, &nrhs, aux, &lda, ipiv, B, &ldb, info);
-
-		LapackConnector::transpose(aux, A, n, lda);
-	}
-
-	static inline
-	void zpotrf(char uplo,
-				int n,
-				ComplexMatrix &a,
-				const int lda,
-				int& info)
-	{
-/*		complex<double> *aux = LapackConnector::transpose(a, n, lda);
-		zpotrf_( &uplo, &n, aux, &lda, &info);
-		LapackConnector::transpose(aux, a, n, lda);
-		delete[] aux;
-*/		// Peize Lin update 2015-12-27
-		const char uplo_changed(change_uplo(uplo));
-		zpotrf_( &uplo_changed, &n, a.c, &lda, &info);
-	}
-
-	static inline
-	void zpotri(char uplo,
-				int n,
-				ComplexMatrix &a,
-				const int lda,
-				int& info)
-	{
-/*		complex<double> *aux = LapackConnector::transpose(a, n, lda);
-		zpotri_( &uplo, &n, aux, &lda, &info);
-		LapackConnector::transpose(aux, a, n, lda);
-		delete[] aux;
-*/		// Peize Lin update 2015-12-27
-		const char uplo_changed(change_uplo(uplo));
-		zpotri_( &uplo_changed, &n, a.c, &lda, &info);		
-	}
-
-
-	// wrap function of fortran lapack routine zhegv.
-    static inline
-	void zhegv(	const int itype,
-                const char jobz,
-                const char uplo,
-                const int n,
-                ComplexMatrix& a,
-                const int lda,
-                ComplexMatrix& b,
-                const int ldb,
-                double* w,
-                std::complex<double>* work,
-                int lwork,
-                double* rwork,
-                int& info	)
-	{	// Transpose the complex matrix to the fortran-form real-complex array.
-		complex<double>* aux = LapackConnector::transpose(a, n, lda);
-		complex<double>* bux = LapackConnector::transpose(b, n, ldb);
-		
-		// call the fortran routine
-		zhegv_(&itype, &jobz, &uplo, &n, aux, &lda, bux, &ldb, w, work, &lwork, rwork, &info);
-
-		// Transpose the fortran-form real-complex array to the complex matrix.
-		LapackConnector::transpose(aux, a, n, lda);
-		LapackConnector::transpose(bux, b, n, ldb);
-
-		// free the memory.
-		delete[] aux;
-		delete[] bux;
-	}
-
-	// wrap function of fortran lapack routine zheev.
-	static inline
-	void zheev( const char jobz,
-				const char uplo,
-				const int n,
-				ComplexMatrix& a,
-				const int lda,
-				double* w,
-				complex< double >* work,
-				const int lwork,
-				double* rwork,
-				int* info	)
-	{	// Transpose the complex matrix to the fortran-form real-complex array.
-		complex<double> *aux = LapackConnector::transpose(a, n, lda);
-
-		// call the fortran routine
-		zheev_(&jobz, &uplo, &n, aux, &lda, w, work, &lwork, rwork, info);
-
-		// Transpose the fortran-form real-complex array to the complex matrix.
-		LapackConnector::transpose(aux, a, n, lda);
-
-		// free the memory.
-		delete[] aux;
-	}
-	
-	// wrap function of fortran lapack routine zheev.	
-	static inline
-	void zhegvx( const int itype,
-				 const char jobz,
-				 const char range,
-				 const char uplo,
-				 const int n,
-				 ComplexMatrix& a,
-				 const int lda,
-				 ComplexMatrix& b,
-				 const int ldb,
-				 const double vl,
-				 const double vu,
-				 const int il,
-				 const int iu,
-				 const double abstol,
-				 const int m,
-				 double* w,
-				 ComplexMatrix& z,
-				 const int ldz,
-				 std::complex<double>* work,
-				 const int lwork,
-				 double* rwork,
-				 int* iwork,
-				 int* ifail,
-				 int& info)
-	{
-		 // Transpose the complex matrix to the fortran-form real-complex array.
-		complex<double>* aux = LapackConnector::transpose(a, n, lda);
-		complex<double>* bux = LapackConnector::transpose(b, n, ldb);
-		complex<double>* zux = LapackConnector::transpose(z, n, ldz);
-		
-		// call the fortran routine
-//	    zhegvx_(&itype, &jobz, &range, &uplo, &n, aux, &lda, bux, &ldb, &vl,
-//	    &vu, &il,&iu, &abstol, &m, w, zux, &ldz, work, &lwork, rwork, iwork, ifail, &info);
-	    
-	    // Transpose the fortran-form real-complex array to the complex matrix
-		LapackConnector::transpose(zux, z, n, ldz);	 
-           
-		// free the memory.
-		delete[] aux;
-		delete[] bux;   
- 		delete[] zux;
-		
-     }             
-
-	
-	// Peize Lin add 2015-11-20
-	static inline 
-	void daxpy( int N,  double alpha,  double *X,  int incX, double *Y,  int incY)
-	{
-		daxpy_(&N, &alpha, X, &incX, Y, &incY);
-	}	
-	
-	static inline 
-	void dscal( int N,  double alpha, double *X,  int incX)
-	{
-		dscal_(&N, &alpha, X, &incX);
-	}
-
-	// Peize Lin add 2015-12-07
-	// multiply two matrix. C=A*B
-	static inline 
-	void zgemm(const char TransA, const char TransB,
-			   const int M, const int N, const int K,
-			   const double alpha, const void *A, const int lda, const void *B, const int ldb,
-			   const double beta, void *C, const int ldc)				   
-	{
-		zgemm_(&TransA, &TransB,
-			   &N, &M, &K,
-			   &alpha, B, &ldb, A, &lda,
-			   &beta, C, &ldc);
-	}
-
-	// multiply two matrix. C=A*B('L'eft) or C=B*A('R'ight). A^T=A, stored in 'U'pper or 'L'ower
-	static inline 
-	void zsymm(const char side, const char uplo,
-			   const int M, const int N,
-			   const double alpha, const void *A, const int lda, const void *B, const int ldb,
-			   const double beta, void *C, const int ldc)				   
-	{
-		const char &side_changed(change_side(side));
-		const char &uplo_changed(change_uplo(uplo));
-		zsymm_(&side_changed, &uplo_changed,
-			   &N, &M,
-			   &alpha, A, &lda, B, &ldb,
-			   &beta, C, &ldc);
-	}
-	
-	// multiply two matrix. C=A*B('L'eft) or C=B*A('R'ight). A^H=A, stored in 'U'pper or 'L'ower
-	static inline 
-	void zhemm(const char side, const char uplo,
-			   const int M, const int N,
-			   const double alpha, const void *A, const int lda, const void *B, const int ldb,
-			   const double beta, void *C, const int ldc)				   
-	{
-		const char &side_changed(change_side(side));
-		const char &uplo_changed(change_uplo(uplo));
-		zhemm_(&side_changed, &uplo_changed,
-			   &N, &M,
-			   &alpha, A, &lda, B, &ldb,
-			   &beta, C, &ldc);
-	}
+  private:
+    // Transpose the complex matrix to the fortran-form real-complex array.
+    static inline complex<double>* transpose(const ComplexMatrix& a, const int n, const int lda)
+    {
+        complex<double>* aux = new std::complex<double>[lda * n];
+        for (int i = 0; i < n; ++i)
+        {
+            for (int j = 0; j < lda; ++j)
+            {
+                aux[i * lda + j] = a(j, i); // aux[i*lda+j] means aux[i][j] in semantic, not in syntax!
+            }
+        }
+        return aux;
+    }
+
+    // Transpose the fortran-form real-complex array to the complex matrix.
+    static inline void transpose(const std::complex<double>* aux, ComplexMatrix& a, const int n, const int lda)
+    {
+        for (int i = 0; i < n; ++i)
+        {
+            for (int j = 0; j < lda; ++j)
+            {
+                a(j, i) = aux[i * lda + j]; // aux[i*lda+j] means aux[i][j] in semantic, not in syntax!
+            }
+        }
+    }
+
+    static inline char change_side(const char& side)
+    {
+        switch (side)
+        {
+        case 'L':
+            return 'R';
+        case 'R':
+            return 'L';
+        default:
+            throw invalid_argument("Side must be 'L' or 'R'");
+        }
+    }
+
+    static inline char change_uplo(const char& uplo)
+    {
+        switch (uplo)
+        {
+        case 'U':
+            return 'L';
+        case 'L':
+            return 'U';
+        default:
+            throw invalid_argument("Uplo must be 'U' or 'L'");
+        }
+    }
+
+  public:
+    static inline void zgesv(const int n,
+                             const int nrhs,
+                             ComplexMatrix& A,
+                             const int lda,
+                             const int* ipiv,
+                             complex<double>* B,
+                             const int ldb,
+                             int* info)
+    {
+        complex<double>* aux = LapackConnector::transpose(A, n, lda);
+
+        zgesv_(&n, &nrhs, aux, &lda, ipiv, B, &ldb, info);
+
+        LapackConnector::transpose(aux, A, n, lda);
+    }
+
+    static inline void zpotrf(char uplo, int n, ComplexMatrix& a, const int lda, int& info)
+    {
+        /*		complex<double> *aux = LapackConnector::transpose(a, n, lda);
+                zpotrf_( &uplo, &n, aux, &lda, &info);
+                LapackConnector::transpose(aux, a, n, lda);
+                delete[] aux;
+        */		// Peize Lin update 2015-12-27
+        const char uplo_changed(change_uplo(uplo));
+        zpotrf_(&uplo_changed, &n, a.c, &lda, &info);
+    }
+
+    static inline void zpotri(char uplo, int n, ComplexMatrix& a, const int lda, int& info)
+    {
+        /*		complex<double> *aux = LapackConnector::transpose(a, n, lda);
+                zpotri_( &uplo, &n, aux, &lda, &info);
+                LapackConnector::transpose(aux, a, n, lda);
+                delete[] aux;
+        */		// Peize Lin update 2015-12-27
+        const char uplo_changed(change_uplo(uplo));
+        zpotri_(&uplo_changed, &n, a.c, &lda, &info);
+    }
+
+    // wrap function of fortran lapack routine zhegv.
+    static inline void zhegv(const int itype,
+                             const char jobz,
+                             const char uplo,
+                             const int n,
+                             ComplexMatrix& a,
+                             const int lda,
+                             ComplexMatrix& b,
+                             const int ldb,
+                             double* w,
+                             std::complex<double>* work,
+                             int lwork,
+                             double* rwork,
+                             int& info)
+    { // Transpose the complex matrix to the fortran-form real-complex array.
+        complex<double>* aux = LapackConnector::transpose(a, n, lda);
+        complex<double>* bux = LapackConnector::transpose(b, n, ldb);
+
+        // call the fortran routine
+        zhegv_(&itype, &jobz, &uplo, &n, aux, &lda, bux, &ldb, w, work, &lwork, rwork, &info);
+
+        // Transpose the fortran-form real-complex array to the complex matrix.
+        LapackConnector::transpose(aux, a, n, lda);
+        LapackConnector::transpose(bux, b, n, ldb);
+
+        // free the memory.
+        delete[] aux;
+        delete[] bux;
+    }
+
+    // wrap function of fortran lapack routine zheev.
+    static inline void zheev(const char jobz,
+                             const char uplo,
+                             const int n,
+                             ComplexMatrix& a,
+                             const int lda,
+                             double* w,
+                             complex<double>* work,
+                             const int lwork,
+                             double* rwork,
+                             int* info)
+    { // Transpose the complex matrix to the fortran-form real-complex array.
+        complex<double>* aux = LapackConnector::transpose(a, n, lda);
+
+        // call the fortran routine
+        zheev_(&jobz, &uplo, &n, aux, &lda, w, work, &lwork, rwork, info);
+
+        // Transpose the fortran-form real-complex array to the complex matrix.
+        LapackConnector::transpose(aux, a, n, lda);
+
+        // free the memory.
+        delete[] aux;
+    }
+
+    // wrap function of fortran lapack routine zheev.
+    static inline void zhegvx(const int itype,
+                              const char jobz,
+                              const char range,
+                              const char uplo,
+                              const int n,
+                              ComplexMatrix& a,
+                              const int lda,
+                              ComplexMatrix& b,
+                              const int ldb,
+                              const double vl,
+                              const double vu,
+                              const int il,
+                              const int iu,
+                              const double abstol,
+                              const int m,
+                              double* w,
+                              ComplexMatrix& z,
+                              const int ldz,
+                              std::complex<double>* work,
+                              const int lwork,
+                              double* rwork,
+                              int* iwork,
+                              int* ifail,
+                              int& info)
+    {
+        // Transpose the complex matrix to the fortran-form real-complex array.
+        complex<double>* aux = LapackConnector::transpose(a, n, lda);
+        complex<double>* bux = LapackConnector::transpose(b, n, ldb);
+        complex<double>* zux = LapackConnector::transpose(z, n, ldz);
+
+        // call the fortran routine
+        //	    zhegvx_(&itype, &jobz, &range, &uplo, &n, aux, &lda, bux, &ldb, &vl,
+        //	    &vu, &il,&iu, &abstol, &m, w, zux, &ldz, work, &lwork, rwork, iwork, ifail, &info);
+
+        // Transpose the fortran-form real-complex array to the complex matrix
+        LapackConnector::transpose(zux, z, n, ldz);
+
+        // free the memory.
+        delete[] aux;
+        delete[] bux;
+        delete[] zux;
+    }
+
+    // Peize Lin add 2015-11-20
+    static inline void daxpy(int N, double alpha, double* X, int incX, double* Y, int incY)
+    {
+        daxpy_(&N, &alpha, X, &incX, Y, &incY);
+    }
+
+    static inline void dscal(int N, double alpha, double* X, int incX)
+    {
+        dscal_(&N, &alpha, X, &incX);
+    }
+
+    // Peize Lin add 2015-12-07
+    // multiply two matrix. C=A*B
+    static inline void zgemm(const char TransA,
+                             const char TransB,
+                             const int M,
+                             const int N,
+                             const int K,
+                             const double alpha,
+                             const void* A,
+                             const int lda,
+                             const void* B,
+                             const int ldb,
+                             const double beta,
+                             void* C,
+                             const int ldc)
+    {
+        zgemm_(&TransA, &TransB, &N, &M, &K, &alpha, B, &ldb, A, &lda, &beta, C, &ldc);
+    }
+
+    // multiply two matrix. C=A*B('L'eft) or C=B*A('R'ight). A^T=A, stored in 'U'pper or 'L'ower
+    static inline void zsymm(const char side,
+                             const char uplo,
+                             const int M,
+                             const int N,
+                             const double alpha,
+                             const void* A,
+                             const int lda,
+                             const void* B,
+                             const int ldb,
+                             const double beta,
+                             void* C,
+                             const int ldc)
+    {
+        const char& side_changed(change_side(side));
+        const char& uplo_changed(change_uplo(uplo));
+        zsymm_(&side_changed, &uplo_changed, &N, &M, &alpha, A, &lda, B, &ldb, &beta, C, &ldc);
+    }
+
+    // multiply two matrix. C=A*B('L'eft) or C=B*A('R'ight). A^H=A, stored in 'U'pper or 'L'ower
+    static inline void zhemm(const char side,
+                             const char uplo,
+                             const int M,
+                             const int N,
+                             const double alpha,
+                             const void* A,
+                             const int lda,
+                             const void* B,
+                             const int ldb,
+                             const double beta,
+                             void* C,
+                             const int ldc)
+    {
+        const char& side_changed(change_side(side));
+        const char& uplo_changed(change_uplo(uplo));
+        zhemm_(&side_changed, &uplo_changed, &N, &M, &alpha, A, &lda, B, &ldb, &beta, C, &ldc);
+    }
 };
 
-#endif  // LAPACKCONNECTOR_HPP
+#endif // LAPACKCONNECTOR_HPP
diff --git a/tools/01_NAO_generation/SIAB/src_tools/mathzone.cpp b/tools/01_NAO_generation/SIAB/src_tools/mathzone.cpp
index a42bbbae093..2c4b101ef79 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/mathzone.cpp
+++ b/tools/01_NAO_generation/SIAB/src_tools/mathzone.cpp
@@ -1,6 +1,5 @@
 #include "mathzone.h"
 
-
 double Mathzone::PI = 3.14159265;
 
 Mathzone::Mathzone()
@@ -11,25 +10,23 @@ Mathzone::~Mathzone()
 {
 }
 
-double Mathzone::BESSJY(double x, double xnu, double *rj, double *ry, double *rjp, double *ryp)
+double Mathzone::BESSJY(double x, double xnu, double* rj, double* ry, double* rjp, double* ryp)
 {
-	const int XMIN = 2;
-	const double FPMIN = 1.0e-30;
-	const double EPS = 1.0e-10;
-	const int MAXIT = 10000;
+    const int XMIN = 2;
+    const double FPMIN = 1.0e-30;
+    const double EPS = 1.0e-10;
+    const int MAXIT = 10000;
 
-	int i, isign, l, nl;
-    double a, b, br, bi, c, cr, ci, d, del, del1, den, di, dlr, dli, dr, e, f, fact, fact2,
-    fact3, ff, gam, gam1, gam2, gammi, gampl, h, p, pimu, pimu2, q, r, rjl,
-    rjl1, rjmu, rjp1, rjpl, rjtemp, ry1, rymu, rymup, rytemp, sum, sum1,
-    temp, w, x2, xi, xi2;
+    int i, isign, l, nl;
+    double a, b, br, bi, c, cr, ci, d, del, del1, den, di, dlr, dli, dr, e, f, fact, fact2, fact3, ff, gam, gam1, gam2,
+        gammi, gampl, h, p, pimu, pimu2, q, r, rjl, rjl1, rjmu, rjp1, rjpl, rjtemp, ry1, rymu, rymup, rytemp, sum, sum1,
+        temp, w, x2, xi, xi2;
 
     if (x <= 0.0 || xnu < 0.0)
-	{
-//		WARNING_QUIT("Mathzone::BESSJY","bad arguments");
-		exit(0);
-	}
-
+    {
+        //		WARNING_QUIT("Mathzone::BESSJY","bad arguments");
+        exit(0);
+    }
 
     nl = (x < XMIN ? (int)(xnu + 0.5) : IMAX(0, (int)(xnu - x + 1.5)));
     const double xmu = xnu - nl;
@@ -39,24 +36,27 @@ double Mathzone::BESSJY(double x, double xnu, double *rj, double *ry, double *rj
     w = xi2 / PI;
     isign = 1;
     h = xnu * xi;
-    if (h < FPMIN) h = FPMIN;
+    if (h < FPMIN)
+        h = FPMIN;
 
     b = xi2 * xnu;
 
     d = 0.0;
 
     c = h;
-	
-	for (i = 1;i <= MAXIT;i++)
+
+    for (i = 1; i <= MAXIT; i++)
     {
         b += xi2;
         d = b - d;
 
-        if (fabs(d) < FPMIN) d = FPMIN;
+        if (fabs(d) < FPMIN)
+            d = FPMIN;
 
         c = b - 1.0 / c;
 
-        if (fabs(c) < FPMIN) c = FPMIN;
+        if (fabs(c) < FPMIN)
+            c = FPMIN;
 
         d = 1.0 / d;
 
@@ -64,12 +64,15 @@ double Mathzone::BESSJY(double x, double xnu, double *rj, double *ry, double *rj
 
         h = del * h;
 
-        if (d < 0.0) isign = -isign;
+        if (d < 0.0)
+            isign = -isign;
 
-        if (fabs(del - 1.0) < EPS) break;
+        if (fabs(del - 1.0) < EPS)
+            break;
     }
 
-    if (i > MAXIT) cout << "x too large in bessjy; try asymptotic expansion" << endl;
+    if (i > MAXIT)
+        cout << "x too large in bessjy; try asymptotic expansion" << endl;
 
     rjl = isign * FPMIN;
 
@@ -81,7 +84,7 @@ double Mathzone::BESSJY(double x, double xnu, double *rj, double *ry, double *rj
 
     fact = xnu * xi;
 
-	for (l = nl;l >= 1;l--)
+    for (l = nl; l >= 1; l--)
     {
         rjtemp = fact * rjl + rjpl;
         fact -= xi;
@@ -89,7 +92,8 @@ double Mathzone::BESSJY(double x, double xnu, double *rj, double *ry, double *rj
         rjl = rjtemp;
     }
 
-    if (rjl == 0.0) rjl = EPS;
+    if (rjl == 0.0)
+        rjl = EPS;
 
     f = rjpl / rjl;
 
@@ -101,7 +105,7 @@ double Mathzone::BESSJY(double x, double xnu, double *rj, double *ry, double *rj
         d = -log(x2);
         e = xmu * d;
         fact2 = (fabs(e) < EPS ? 1.0 : sinh(e) / e);
-		// call BESCHB
+        // call BESCHB
         BESCHB(xmu, &gam1, &gam2, &gampl, &gammi);
         ff = 2.0 / PI * fact * (gam1 * cosh(e) + gam2 * fact2 * d);
         e = exp(e);
@@ -115,7 +119,7 @@ double Mathzone::BESSJY(double x, double xnu, double *rj, double *ry, double *rj
         sum = ff + r * q;
         sum1 = p;
 
-		for (i = 1;i <= MAXIT;i++)
+        for (i = 1; i <= MAXIT; i++)
         {
             ff = (i * ff + p + q) / (i * i - xmu2);
             c *= (d / i);
@@ -126,10 +130,12 @@ double Mathzone::BESSJY(double x, double xnu, double *rj, double *ry, double *rj
             del1 = c * p - i * del;
             sum1 += del1;
 
-            if (fabs(del) < (1.0 + fabs(sum))*EPS) break;
+            if (fabs(del) < (1.0 + fabs(sum)) * EPS)
+                break;
         }
 
-        if (i > MAXIT) cout << "bessy series failed to converge";
+        if (i > MAXIT)
+            cout << "bessy series failed to converge";
 
         rymu = -sum;
 
@@ -140,7 +146,7 @@ double Mathzone::BESSJY(double x, double xnu, double *rj, double *ry, double *rj
         rjmu = w / (rymup - f * rymu);
     }
 
-	else
+    else
     {
         a = 0.25 - xmu2;
         p = -0.5 * xi;
@@ -159,14 +165,15 @@ double Mathzone::BESSJY(double x, double xnu, double *rj, double *ry, double *rj
         q = p * dli + q * dlr;
         p = temp;
 
-		for (i = 2;i <= MAXIT;i++)
+        for (i = 2; i <= MAXIT; i++)
         {
             a += 2 * (i - 1);
             bi += 2.0;
             dr = a * dr + br;
             di = a * di + bi;
 
-            if (fabs(dr) + fabs(di) < FPMIN) dr = FPMIN;
+            if (fabs(dr) + fabs(di) < FPMIN)
+                dr = FPMIN;
 
             fact = a / (cr * cr + ci * ci);
 
@@ -174,7 +181,8 @@ double Mathzone::BESSJY(double x, double xnu, double *rj, double *ry, double *rj
 
             ci = bi - ci * fact;
 
-            if (fabs(cr) + fabs(ci) < FPMIN) cr = FPMIN;
+            if (fabs(cr) + fabs(ci) < FPMIN)
+                cr = FPMIN;
 
             den = dr * dr + di * di;
 
@@ -192,17 +200,21 @@ double Mathzone::BESSJY(double x, double xnu, double *rj, double *ry, double *rj
 
             p = temp;
 
-            if (fabs(dlr - 1.0) + fabs(dli) < EPS) break;
+            if (fabs(dlr - 1.0) + fabs(dli) < EPS)
+                break;
         }
 
-		if (i > MAXIT) cout << "cf2 failed in bessjy";
+        if (i > MAXIT)
+            cout << "cf2 failed in bessjy";
 
         gam = (p - f) / q;
 
         rjmu = sqrt(w / ((p - f) * gam + q));
 
-		if(rjl >=0 ) rjmu = fabs(rjmu);
-		else rjmu = -fabs(rjmu);
+        if (rjl >= 0)
+            rjmu = fabs(rjmu);
+        else
+            rjmu = -fabs(rjmu);
 
         rymu = rjmu * gam;
 
@@ -216,7 +228,7 @@ double Mathzone::BESSJY(double x, double xnu, double *rj, double *ry, double *rj
     *rj = rjl1 * fact;
     *rjp = rjp1 * fact;
 
-    for (i = 1;i <= nl;i++)
+    for (i = 1; i <= nl; i++)
     {
         rytemp = (xmu + i) * xi2 * ry1 - rymu;
         rymu = ry1;
@@ -230,24 +242,34 @@ double Mathzone::BESSJY(double x, double xnu, double *rj, double *ry, double *rj
 
 int Mathzone::IMAX(int a, int b)
 {
-    if (a > b) return a;
-    else return b;
+    if (a > b)
+        return a;
+    else
+        return b;
 }
 
-void Mathzone::BESCHB(double x, double *gam1, double *gam2, double *gampl, double *gammi)
+void Mathzone::BESCHB(double x, double* gam1, double* gam2, double* gampl, double* gammi)
 {
-	const int NUSE1 = 7;
-	const int NUSE2 = 8;
+    const int NUSE1 = 7;
+    const int NUSE2 = 8;
     double xx;
-    static double c1[] = {   -1.142022680371168e0, 6.5165112670737e-3,
-                             3.087090173086e-4, -3.4706269649e-6,
-                             6.9437664e-9, 3.67795e-11, -1.356e-13
-                         };
-    static double c2[] = {   1.843740587300905e0, -7.68528408447867e-2,
-                             1.2719271366546e-3, -4.9717367042e-6, -3.31261198e-8,
-                             2.423096e-10, -1.702e-13, -1.49e-15
-                         };
-    xx = 8.0 * x * x - 1.0; //Multiply x by 2 to make range be .1 to 1,and then apply transformation for evaluating even     Chebyshev series.
+    static double c1[] = {-1.142022680371168e0,
+                          6.5165112670737e-3,
+                          3.087090173086e-4,
+                          -3.4706269649e-6,
+                          6.9437664e-9,
+                          3.67795e-11,
+                          -1.356e-13};
+    static double c2[] = {1.843740587300905e0,
+                          -7.68528408447867e-2,
+                          1.2719271366546e-3,
+                          -4.9717367042e-6,
+                          -3.31261198e-8,
+                          2.423096e-10,
+                          -1.702e-13,
+                          -1.49e-15};
+    xx = 8.0 * x * x - 1.0; // Multiply x by 2 to make range be .1 to 1,and then apply transformation for evaluating
+                            // even     Chebyshev series.
     *gam1 = CHEBEV(-1.0, 1.0, c1, NUSE1, xx);
     *gam2 = CHEBEV(-1.0, 1.0, c2, NUSE2, xx);
     *gampl = *gam2 - x * (*gam1);
@@ -259,384 +281,361 @@ double Mathzone::CHEBEV(double a, double b, double c[], int m, double x)
     double d = 0.0, dd = 0.0, sv, y, y2;
     int j;
 
-    if ((x - a)*(x - b) > 0.0) cout << "x not in range in routine chebev" << endl;
+    if ((x - a) * (x - b) > 0.0)
+        cout << "x not in range in routine chebev" << endl;
 
     y2 = 2.0 * (y = (2.0 * x - a - b) / (b - a));
 
-    for (j = m - 1;j >= 1;j--)
+    for (j = m - 1; j >= 1; j--)
     {
         sv = d;
         d = y2 * d - dd + c[j];
         dd = sv;
     }
 
-    return y*d - dd + 0.5*c[0];
+    return y * d - dd + 0.5 * c[0];
 }
 
-
-double Mathzone::Spherical_Bessel_7(const int n, const double &x)
+double Mathzone::Spherical_Bessel_7(const int n, const double& x)
 {
-	if(x==0)
-	{
-		if(n!=0) return 0;
-		if(n==0) return 1;
-	}
-	double order, rj, rjp, ry, ryp;
+    if (x == 0)
+    {
+        if (n != 0)
+            return 0;
+        if (n == 0)
+            return 1;
+    }
+    double order, rj, rjp, ry, ryp;
 
-	if (n < 0 || x <= 0.0) 
-	{
-		exit(0);
-//		WARNING_QUIT("Mathzone::Spherical_Bessel_7","bad arguments in sphbes");
-	}
+    if (n < 0 || x <= 0.0)
+    {
+        exit(0);
+        //		WARNING_QUIT("Mathzone::Spherical_Bessel_7","bad arguments in sphbes");
+    }
 
-	order = n + 0.5;
+    order = n + 0.5;
 
-	// call BESSSJY
-	BESSJY(x, order, &rj, &ry, &rjp, &ryp);
+    // call BESSSJY
+    BESSJY(x, order, &rj, &ry, &rjp, &ryp);
 
-	const double RTPIO2=1.2533141; // mohan remian question: here is INT????????????
+    const double RTPIO2 = 1.2533141; // mohan remian question: here is INT????????????
 
-	const double factor = RTPIO2 / sqrt(x);
+    const double factor = RTPIO2 / sqrt(x);
 
-	return factor*rj;
+    return factor * rj;
 }
 
-
-void Mathzone::Spherical_Bessel_Roots
-(
-	const int &num,
-	const int &l,	
-	const double &epsilon,
-	double* eigenvalue,
-	const double &rcut
-)
+void Mathzone::Spherical_Bessel_Roots(const int& num,
+                                      const int& l,
+                                      const double& epsilon,
+                                      double* eigenvalue,
+                                      const double& rcut)
 {
-//	TITLE("Mathzone","Spherical_Bessel_Roots");
-//	if(num<=0) WARNING_QUIT("Mathzone::Spherical_Bessel_Roots","num<=0");
-
-	const double max = 2*PI + (num + (l+0.5)/2 + 0.75)*PI/2 + 
-		sqrt((num + (l+0.5)/2+0.75)*(num + (l+0.5)/2+0.75)*PI*PI/4-(l+0.5)*(l+0.5)/2);
-
-	// magic number !!
-	// only need to > 1
-	const int msh = 10 * num;
-//	cout<<"\n msh = "<<msh;
-
-	// delta don't need to be small,
-	// it only needs to make sure can find the eigenstates
-	const double delta = max / static_cast<double>(msh);
-
-//	cout<<"\n delta = "<<delta;
-
-	double *r = new double[msh];
-	for(int i=0; i<msh; i++) 
-	{
-		r[i] = i*delta;
-	}
-//	cout<<"\n rmax(r) = "<<r[msh-1];
-//	cout<<"\n Enumber = "<<num;
-
-	double *jl = new double[msh];	
-	Mathzone::Spherical_Bessel(msh, r, 1, l, jl);
-
-	int n=0;
-	for(int i=0; i<msh && n<num; i++)
-	{
-		if(jl[i]*jl[i+1] < 0.0)
-		{
-			double y_1 = jl[i];
-			double y_2 = jl[i+1];
-			double x_1 = r[i];
-			double x_2 = r[i+1];
-			double acc = abs( y_2 - y_1 );
-			while(acc > epsilon)
-			{
-				double *rad = new double[100];
-				double *jl_new = new double[100];
-
-				// if not enough accurate, divide again.
-				const double new_delta = (x_2 - x_1)/99.0;
-				for(int j=0;j<100;j++)
-				{
-					rad[j] = x_1 + j*new_delta;
-				}
-				Mathzone::Spherical_Bessel(100,rad,1,l,jl_new);
-				
-				int j=0;
-				for(;j<100;j++)
-				{	
-					if(jl_new[j]*jl_new[j+1]<0)break;
-				}
-
-				x_1 = rad[j];
-				x_2 = rad[j+1];
-				y_1 = jl_new[j];
-				y_2 = jl_new[j+1];
-				acc = abs( y_2 - y_1 );
-				delete[] rad;
-				delete[] jl_new;
-			}	
-			eigenvalue[n]=(x_2 + x_1)*0.5/rcut;
-//			cout << "\n e="<< eigenvalue[n] << " n="<<n;
-			n++;
-		}
-	}
-
-//	cout << "\n emax="<< eigenvalue[num-1] << " num="<<num;
-//	cout << "\n rmax( rcut * r ) = " << r[msh-1]/rcut << "(a.u.)" << endl;
-
-	delete[] r;
-	delete[] jl;
-	return;
+    //	TITLE("Mathzone","Spherical_Bessel_Roots");
+    //	if(num<=0) WARNING_QUIT("Mathzone::Spherical_Bessel_Roots","num<=0");
+
+    const double max
+        = 2 * PI + (num + (l + 0.5) / 2 + 0.75) * PI / 2
+          + sqrt((num + (l + 0.5) / 2 + 0.75) * (num + (l + 0.5) / 2 + 0.75) * PI * PI / 4 - (l + 0.5) * (l + 0.5) / 2);
+
+    // magic number !!
+    // only need to > 1
+    const int msh = 10 * num;
+    //	cout<<"\n msh = "<<msh;
+
+    // delta don't need to be small,
+    // it only needs to make sure can find the eigenstates
+    const double delta = max / static_cast<double>(msh);
+
+    //	cout<<"\n delta = "<<delta;
+
+    double* r = new double[msh];
+    for (int i = 0; i < msh; i++)
+    {
+        r[i] = i * delta;
+    }
+    //	cout<<"\n rmax(r) = "<<r[msh-1];
+    //	cout<<"\n Enumber = "<<num;
+
+    double* jl = new double[msh];
+    Mathzone::Spherical_Bessel(msh, r, 1, l, jl);
+
+    int n = 0;
+    for (int i = 0; i < msh && n < num; i++)
+    {
+        if (jl[i] * jl[i + 1] < 0.0)
+        {
+            double y_1 = jl[i];
+            double y_2 = jl[i + 1];
+            double x_1 = r[i];
+            double x_2 = r[i + 1];
+            double acc = abs(y_2 - y_1);
+            while (acc > epsilon)
+            {
+                double* rad = new double[100];
+                double* jl_new = new double[100];
+
+                // if not enough accurate, divide again.
+                const double new_delta = (x_2 - x_1) / 99.0;
+                for (int j = 0; j < 100; j++)
+                {
+                    rad[j] = x_1 + j * new_delta;
+                }
+                Mathzone::Spherical_Bessel(100, rad, 1, l, jl_new);
+
+                int j = 0;
+                for (; j < 100; j++)
+                {
+                    if (jl_new[j] * jl_new[j + 1] < 0)
+                        break;
+                }
+
+                x_1 = rad[j];
+                x_2 = rad[j + 1];
+                y_1 = jl_new[j];
+                y_2 = jl_new[j + 1];
+                acc = abs(y_2 - y_1);
+                delete[] rad;
+                delete[] jl_new;
+            }
+            eigenvalue[n] = (x_2 + x_1) * 0.5 / rcut;
+            //			cout << "\n e="<< eigenvalue[n] << " n="<<n;
+            n++;
+        }
+    }
+
+    //	cout << "\n emax="<< eigenvalue[num-1] << " num="<<num;
+    //	cout << "\n rmax( rcut * r ) = " << r[msh-1]/rcut << "(a.u.)" << endl;
+
+    delete[] r;
+    delete[] jl;
+    return;
 }
 
 // from sph_bes.f90
-void Mathzone::Spherical_Bessel
-(
-	const int &msh,	//number of grid points
-	const double *r,//radial grid
-	const double &q,	//
-	const int &l,	//angular momentum
-	double *jl		//jl(1:msh) = j_l(q*r(i)),spherical bessel function
+void Mathzone::Spherical_Bessel(const int& msh,  // number of grid points
+                                const double* r, // radial grid
+                                const double& q, //
+                                const int& l,    // angular momentum
+                                double* jl       // jl(1:msh) = j_l(q*r(i)),spherical bessel function
 )
 {
-//	timer::tick("Mathzone","Spherical_Bessel");
-	double x1;
-	int i, ir, ir0;
-
-	if(l>=7)
-	{
-		for(int ir=0; ir<msh; ir++)
-		{
-			x1 = q * r[ir];
-			jl[ir] = Spherical_Bessel_7(l, x1);
-		}
-		return;
-	}
-	
-	if (abs(q) < 1.0e-8)
-	{
-		if (l == -1)
-		{
-			cout << "\n sph_bes, j_{-1}(0) ????";
-		}
-		else if (l == 0)
-		{
-			for (i = 0;i < msh;i++)
-			{
-				jl[i] = 1.0;
-			}
-		}
-		else
-		{
-			for (i = 0;i < msh;i++)
-			{
-				jl[i] = 0.0;
-			}
-		}
-	}
-	else
-	{
-		if (abs(q * r [0]) > 1.0e-8)
-		{
-			ir0 = 0;//mohan modify 2007-10-13
-		}
-		else
-		{
-			if (l == -1)
-			{
-				cout << "\n sph_bes, j_{-1}(0) ?//?";
-			}
-			else if (l == 0)
-			{
-				jl [0] = 1.0;//mohan modify 2007-10-13
-			}
-			else
-			{
-				jl [0] = 0.0;//mohan modify 2007-10-13
-			}
-			ir0 = 1;//mohan modify 2007-10-13
-		}
-		if (l == - 1)
-		{
-			for (ir = ir0;ir < msh; ir++)
-			{
-				x1 = q * r[ir];
-				jl [ir] = cos(x1) / x1;
-			}
-		}
-		else if (l == 0)
-		{
-			for (ir = ir0;ir < msh;ir++)
-			{
-				x1 = q * r[ir];
-				jl [ir] = sin(x1) / x1;
-			}
-		}
-		else if (l == 1)
-		{
-			for (ir = ir0;ir < msh;ir++)
-			{
-				x1 = q * r[ir];
-				const double sinx = sin(x1);
-				const double cosx = cos(x1);
-				jl [ir] = (sinx / x1 - cosx) / x1;
-			}
-		}
-		else if (l == 2)
-		{
-			for (ir = ir0;ir < msh;ir++)
-			{
-				const double x1 = q * r[ir];
-				const double sinx = sin(x1);
-				const double cosx = cos(x1);
-				jl [ir] = ((3.0 / x1  - x1) * sinx
-				           - 3.0 * cosx) / (x1 * x1);
-			}
-		}
-		else if (l == 3)
-		{
-			for (ir = ir0;ir < msh;ir++)
-			{
-				x1 = q * r[ir];
-				jl [ir] = (sin(x1) * (15.0 / x1 - 6.0 * x1) +
-				           cos(x1) * (x1 * x1 - 15.0)) / pow(x1, 3);//mohan modify 2007-10-13
-			}
-		}
-		else if (l == 4)
-		{
-			for (ir = ir0;ir < msh;ir++)
-			{
-				const double x1 = q * r[ir];
-				const double x2 = x1 * x1;
-				const double x3 = x1 * x2;
-				const double x4 = x1 * x3;
-				const double x5 = x1 * x4;
-				jl [ir] = (sin(x1) * (105.0 - 45.0 * x2 + x4) +
-				           cos(x1)  * (10.0 * x3 - 105.0 * x1)) / x5;   // mohan modify 2007-10-13
-			}
-		}
-		else if (l == 5)
-		{
-			for (ir = ir0;ir < msh;ir++)
-			{
-				x1 = q * r[ir];
-
-				if (x1 < 0.14)
-				{
-					jl[ir] = 0;//mohan add 2007-10-15
-				}
-				else
-				{
-					double cx1 = cos(x1);
-					double sx1 = sin(x1);
-					jl [ir] = (-cx1 -
-					           (945.0 * cx1) / pow(x1, 4) +
-					           (105.0 * cx1) / (x1 * x1)  +
-					           (945.0 * sx1) / pow(x1, 5) -
-					           (420.0 * sx1) / pow(x1, 3) +
-					           (15.0 * sx1) / x1) / x1;
-
-				}
-			}
-		}
-		else if (l == 6)
-		{
-			for (ir = ir0;ir < msh;ir++)
-			{
-				x1 = q * r[ir];
-
-				if (x1 < 0.29)
-				{
-					jl[ir] = 0;//mohan add 2007-10-15
-				}
-				else
-				{
-					double cx1 = cos(x1);
-					double sx1 = sin(x1);
-					jl [ir] = ((-10395.0 * cx1) / pow(x1, 5) +
-					           (1260.0 * cx1) / pow(x1, 3) -
-					           (21.0 * cx1) / x1 - sx1 +
-					           (10395.0 * sx1) / pow(x1, 6) -
-					           (4725.0 * sx1) / pow(x1, 4) +
-					           (210.0 * sx1) / (x1 * x1)) / x1;
-				}
-			}
-		}//mohan modify 2007-11-20 reduce cos , sin , q*r[ir] times;
-		else
-		{
-			cout << "\n error in sph_bes, l out of {-1 ... 6},l = " << l ;
-		}
-	}
-//	printr1_d(ofs,"jl[] :",jl,msh);
-//	timer::tick("Mathzone","Spherical_Bessel");
-	return;
-} 
-
-void Mathzone::Pointwise_Product
-( 
-	const size_t mesh,
-	const double * func1,
-	const double * func2,
-	double * mul_func
-)
+    //	timer::tick("Mathzone","Spherical_Bessel");
+    double x1;
+    int i, ir, ir0;
+
+    if (l >= 7)
+    {
+        for (int ir = 0; ir < msh; ir++)
+        {
+            x1 = q * r[ir];
+            jl[ir] = Spherical_Bessel_7(l, x1);
+        }
+        return;
+    }
+
+    if (abs(q) < 1.0e-8)
+    {
+        if (l == -1)
+        {
+            cout << "\n sph_bes, j_{-1}(0) ????";
+        }
+        else if (l == 0)
+        {
+            for (i = 0; i < msh; i++)
+            {
+                jl[i] = 1.0;
+            }
+        }
+        else
+        {
+            for (i = 0; i < msh; i++)
+            {
+                jl[i] = 0.0;
+            }
+        }
+    }
+    else
+    {
+        if (abs(q * r[0]) > 1.0e-8)
+        {
+            ir0 = 0; // mohan modify 2007-10-13
+        }
+        else
+        {
+            if (l == -1)
+            {
+                cout << "\n sph_bes, j_{-1}(0) ?//?";
+            }
+            else if (l == 0)
+            {
+                jl[0] = 1.0; // mohan modify 2007-10-13
+            }
+            else
+            {
+                jl[0] = 0.0; // mohan modify 2007-10-13
+            }
+            ir0 = 1; // mohan modify 2007-10-13
+        }
+        if (l == -1)
+        {
+            for (ir = ir0; ir < msh; ir++)
+            {
+                x1 = q * r[ir];
+                jl[ir] = cos(x1) / x1;
+            }
+        }
+        else if (l == 0)
+        {
+            for (ir = ir0; ir < msh; ir++)
+            {
+                x1 = q * r[ir];
+                jl[ir] = sin(x1) / x1;
+            }
+        }
+        else if (l == 1)
+        {
+            for (ir = ir0; ir < msh; ir++)
+            {
+                x1 = q * r[ir];
+                const double sinx = sin(x1);
+                const double cosx = cos(x1);
+                jl[ir] = (sinx / x1 - cosx) / x1;
+            }
+        }
+        else if (l == 2)
+        {
+            for (ir = ir0; ir < msh; ir++)
+            {
+                const double x1 = q * r[ir];
+                const double sinx = sin(x1);
+                const double cosx = cos(x1);
+                jl[ir] = ((3.0 / x1 - x1) * sinx - 3.0 * cosx) / (x1 * x1);
+            }
+        }
+        else if (l == 3)
+        {
+            for (ir = ir0; ir < msh; ir++)
+            {
+                x1 = q * r[ir];
+                jl[ir] = (sin(x1) * (15.0 / x1 - 6.0 * x1) + cos(x1) * (x1 * x1 - 15.0))
+                         / pow(x1, 3); // mohan modify 2007-10-13
+            }
+        }
+        else if (l == 4)
+        {
+            for (ir = ir0; ir < msh; ir++)
+            {
+                const double x1 = q * r[ir];
+                const double x2 = x1 * x1;
+                const double x3 = x1 * x2;
+                const double x4 = x1 * x3;
+                const double x5 = x1 * x4;
+                jl[ir] = (sin(x1) * (105.0 - 45.0 * x2 + x4) + cos(x1) * (10.0 * x3 - 105.0 * x1))
+                         / x5; // mohan modify 2007-10-13
+            }
+        }
+        else if (l == 5)
+        {
+            for (ir = ir0; ir < msh; ir++)
+            {
+                x1 = q * r[ir];
+
+                if (x1 < 0.14)
+                {
+                    jl[ir] = 0; // mohan add 2007-10-15
+                }
+                else
+                {
+                    double cx1 = cos(x1);
+                    double sx1 = sin(x1);
+                    jl[ir] = (-cx1 - (945.0 * cx1) / pow(x1, 4) + (105.0 * cx1) / (x1 * x1) + (945.0 * sx1) / pow(x1, 5)
+                              - (420.0 * sx1) / pow(x1, 3) + (15.0 * sx1) / x1)
+                             / x1;
+                }
+            }
+        }
+        else if (l == 6)
+        {
+            for (ir = ir0; ir < msh; ir++)
+            {
+                x1 = q * r[ir];
+
+                if (x1 < 0.29)
+                {
+                    jl[ir] = 0; // mohan add 2007-10-15
+                }
+                else
+                {
+                    double cx1 = cos(x1);
+                    double sx1 = sin(x1);
+                    jl[ir] = ((-10395.0 * cx1) / pow(x1, 5) + (1260.0 * cx1) / pow(x1, 3) - (21.0 * cx1) / x1 - sx1
+                              + (10395.0 * sx1) / pow(x1, 6) - (4725.0 * sx1) / pow(x1, 4) + (210.0 * sx1) / (x1 * x1))
+                             / x1;
+                }
+            }
+        } // mohan modify 2007-11-20 reduce cos , sin , q*r[ir] times;
+        else
+        {
+            cout << "\n error in sph_bes, l out of {-1 ... 6},l = " << l;
+        }
+    }
+    //	printr1_d(ofs,"jl[] :",jl,msh);
+    //	timer::tick("Mathzone","Spherical_Bessel");
+    return;
+}
+
+void Mathzone::Pointwise_Product(const size_t mesh, const double* func1, const double* func2, double* mul_func)
 {
-	for( size_t i(0); i<mesh; ++i)
-	{
-		mul_func[i] = func1[i] * func2[i];
-	}
-	return;
+    for (size_t i(0); i < mesh; ++i)
+    {
+        mul_func[i] = func1[i] * func2[i];
+    }
+    return;
 }
 
-void Mathzone::Simpson_Integral
-(
-	const int mesh,
-	const double *func,
-	const double *rab,
-	double &asum
-)
+void Mathzone::Simpson_Integral(const int mesh, const double* func, const double* rab, double& asum)
 {
-  /*     simpson's rule integration. On input:
-  !      mesh = mhe number of grid points (should be odd)
-  !      func(i)= function to be integrated
-  !      rab(i) = r(i) * dr(i)/di * di
-  !      For the logarithmic grid not including r=0 :
-  !      r(i) = r_0*exp((i-1)*dx) ==> rab(i)=r(i)*dx
-  !      For the logarithmic grid including r=0 :
-  !      r(i) = a(exp((i-1)*dx)-1) ==> rab(i)=(r(i)+a)*dx
-  !      Output in asum = \sum_i c_i f(i)*rab(i) = \int_0^\infty f(r) dr
-  !      where c_i are alternativaly 2/3, 4/3 except c_1 = c_mesh = 1/3
-  */
-	//  simpson's rule integrator for function stored on the
-	//  radial logarithmic mesh
-	//	routine assumes that mesh is an odd number so run check
-	if (mesh % 2 == 0)
-	{
-		cout << "\n error in subroutine simpson ";
-		cout << "\n routine assumes mesh is odd but mesh = "
-		<< mesh << endl;
-		return;
-	}
-
-	asum = 0.00;
-	const double r12 = 1.00 / 12.00;
-	double f3 = func [0] * rab [0] * r12;
-	for (int i = 1;i < mesh;i += 2)
-	{
-		const double f1 = f3;
-		const double f2 = func [i] * rab [i] * r12;
-		f3 = func [i + 1] * rab [i + 1] * r12;
-		asum += 4.00 * f1 + 16.00 * f2 + 4.00 * f3;
-	}
-	return;
-}// end subroutine simpson
-
-double Mathzone::Boltzmann_dist(const double &value, const double &t)
+    /*     simpson's rule integration. On input:
+    !      mesh = mhe number of grid points (should be odd)
+    !      func(i)= function to be integrated
+    !      rab(i) = r(i) * dr(i)/di * di
+    !      For the logarithmic grid not including r=0 :
+    !      r(i) = r_0*exp((i-1)*dx) ==> rab(i)=r(i)*dx
+    !      For the logarithmic grid including r=0 :
+    !      r(i) = a(exp((i-1)*dx)-1) ==> rab(i)=(r(i)+a)*dx
+    !      Output in asum = \sum_i c_i f(i)*rab(i) = \int_0^\infty f(r) dr
+    !      where c_i are alternativaly 2/3, 4/3 except c_1 = c_mesh = 1/3
+    */
+    //  simpson's rule integrator for function stored on the
+    //  radial logarithmic mesh
+    //	routine assumes that mesh is an odd number so run check
+    if (mesh % 2 == 0)
+    {
+        cout << "\n error in subroutine simpson ";
+        cout << "\n routine assumes mesh is odd but mesh = " << mesh << endl;
+        return;
+    }
+
+    asum = 0.00;
+    const double r12 = 1.00 / 12.00;
+    double f3 = func[0] * rab[0] * r12;
+    for (int i = 1; i < mesh; i += 2)
+    {
+        const double f1 = f3;
+        const double f2 = func[i] * rab[i] * r12;
+        f3 = func[i + 1] * rab[i + 1] * r12;
+        asum += 4.00 * f1 + 16.00 * f2 + 4.00 * f3;
+    }
+    return;
+} // end subroutine simpson
+
+double Mathzone::Boltzmann_dist(const double& value, const double& t)
 {
-	static const double small = 1.0e-20;
-	assert(t > small);
-	assert(value > 0.0);
-//	cout << "\n value = " << value;
-//	cout << "\n Boltzmann_dist = " << exp( -value/t );
-	return exp(-value / t);
+    static const double small = 1.0e-20;
+    assert(t > small);
+    assert(value > 0.0);
+    //	cout << "\n value = " << value;
+    //	cout << "\n Boltzmann_dist = " << exp( -value/t );
+    return exp(-value / t);
 }
\ No newline at end of file
diff --git a/tools/01_NAO_generation/SIAB/src_tools/mathzone.h b/tools/01_NAO_generation/SIAB/src_tools/mathzone.h
index a1e9b93f426..d130525c6be 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/mathzone.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/mathzone.h
@@ -1,69 +1,47 @@
 #ifndef MATHZONE_H
 #define MATHZONE_H
 
-#include <cstdlib>
+#include <cassert>
 #include <cmath>
-
+#include <cstdlib>
 #include <fstream>
+#include <iomanip>
 #include <iostream>
 #include <string>
-#include <iomanip>
-#include <cassert>
 
 using namespace std;
 
-class Mathzone 
+class Mathzone
 {
-public:
-
-	Mathzone();
-	~Mathzone();
-
-	//========================================================
-	// Spherical Bessel
-	//========================================================
-	static void Spherical_Bessel
-	(
-		const int &msh,	//number of grid points
-		const double *r,//radial grid
-		const double &q,	//
-		const int &l,	//angular momentum
-		double *jl		//jl(1:msh) = j_l(q*r(i)),spherical bessel function
-	);
-	static void Spherical_Bessel_Roots
-	(
-		const int &num,
-		const int &l,
-		const double &epsilon,
-		double* eigenvalue,
-		const double &rcut
-	);
-	static void Pointwise_Product( 
-		const size_t mesh,
-		const double * func1,
-		const double * func2,
-		double * mul_func
-	);
-	static void Simpson_Integral(
-		const int mesh,
-		const double *func,
-		const double *rab,
-		double &asum
-	);
-	static double Boltzmann_dist( 
-		const double &value, 
-		const double &t 
-	);
-	static double PI;
-
-private:
-	static double Spherical_Bessel_7(const int n, const double &x);
-	static double BESSJY(double x, double xnu, double *rj, double *ry, double *rjp, double *ryp);
-	static void BESCHB(double x, double *gam1, double *gam2, double *gampl, double *gammi);
-	static double CHEBEV(double a, double b, double c[], int m, double x);
-	static int IMAX(int a, int b);
-
-	
+  public:
+    Mathzone();
+    ~Mathzone();
+
+    //========================================================
+    // Spherical Bessel
+    //========================================================
+    static void Spherical_Bessel(const int& msh,  // number of grid points
+                                 const double* r, // radial grid
+                                 const double& q, //
+                                 const int& l,    // angular momentum
+                                 double* jl       // jl(1:msh) = j_l(q*r(i)),spherical bessel function
+    );
+    static void Spherical_Bessel_Roots(const int& num,
+                                       const int& l,
+                                       const double& epsilon,
+                                       double* eigenvalue,
+                                       const double& rcut);
+    static void Pointwise_Product(const size_t mesh, const double* func1, const double* func2, double* mul_func);
+    static void Simpson_Integral(const int mesh, const double* func, const double* rab, double& asum);
+    static double Boltzmann_dist(const double& value, const double& t);
+    static double PI;
+
+  private:
+    static double Spherical_Bessel_7(const int n, const double& x);
+    static double BESSJY(double x, double xnu, double* rj, double* ry, double* rjp, double* ryp);
+    static void BESCHB(double x, double* gam1, double* gam2, double* gampl, double* gammi);
+    static double CHEBEV(double a, double b, double c[], int m, double x);
+    static int IMAX(int a, int b);
 };
 
 #endif
diff --git a/tools/01_NAO_generation/SIAB/src_tools/matrix.cpp b/tools/01_NAO_generation/SIAB/src_tools/matrix.cpp
index 5b6df2600d5..f2aa507eef0 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/matrix.cpp
+++ b/tools/01_NAO_generation/SIAB/src_tools/matrix.cpp
@@ -1,12 +1,11 @@
 /* matrix.cpp file */
-#include <iostream>
-#include <fstream>
-#include <iomanip>
-#include <cstdio>
 #include <cassert>
-
 #include <cmath>
+#include <cstdio>
 #include <cstdlib>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
 
 using namespace std;
 #include "matrix.h"
@@ -21,40 +20,41 @@ int matrix::mCount = 0;
 
 void matrixAlloc()
 {
-	cout << "\n Allocation error for Matrix " << endl;
-	abort();
+    cout << "\n Allocation error for Matrix " << endl;
+    abort();
 }
 
-void matrix::init(const int nrows,const int ncols)
+void matrix::init(const int nrows, const int ncols)
 {
-	nr = nrows;
-	nc = ncols;
-	c = NULL;
-//	hermetian = 0;
-	set_new_handler(matrixAlloc);
+    nr = nrows;
+    nc = ncols;
+    c = NULL;
+    //	hermetian = 0;
+    set_new_handler(matrixAlloc);
 
-	if(nr*nc == 0) c = NULL;
-	else
-	{
-		c = new double[nrows*ncols]();
-		this->zero_out();
-		assert(c!=0);// terminate if memory not allocated
-	}
+    if (nr * nc == 0)
+        c = NULL;
+    else
+    {
+        c = new double[nrows * ncols]();
+        this->zero_out();
+        assert(c != 0); // terminate if memory not allocated
+    }
 
-	mCount++;
+    mCount++;
 }
 
 // Free up memory for matrix
 void matrix::freemem(void)
 {
-	delete [] c;
-	c = NULL;
+    delete[] c;
+    c = NULL;
 }
 
 // constructor with sizes
-matrix::matrix(const int nrows,const int ncols)
+matrix::matrix(const int nrows, const int ncols)
 {
-	this->init(nrows, ncols);
+    this->init(nrows, ncols);
 }
 
 //******************
@@ -62,12 +62,13 @@ matrix::matrix(const int nrows,const int ncols)
 // Copy constructor
 //
 //******************
-matrix::matrix(const matrix &m1)
+matrix::matrix(const matrix& m1)
 {
-	init(m1.nr, m1.nc);
-	//hermetian = m1.hermetian;
-	// Copy over m1 contents
-	for (int i = 0; i < nr*nc; i++) c[i] = m1.c[i];
+    init(m1.nr, m1.nc);
+    // hermetian = m1.hermetian;
+    //  Copy over m1 contents
+    for (int i = 0; i < nr * nc; i++)
+        c[i] = m1.c[i];
 }
 
 //*************
@@ -77,147 +78,153 @@ matrix::matrix(const matrix &m1)
 //*************
 matrix::~matrix()
 {
-	this->freemem();
+    this->freemem();
 }
 
 //******************************
 // reallocate memory for matrix
 //******************************
-void matrix::create(const int nrow,const int ncol)
+void matrix::create(const int nrow, const int ncol)
 {
-	delete [] c;
-	nr = nrow;
-	nc = ncol;
-	c = new double[nrow * ncol]();
-	return;
+    delete[] c;
+    nr = nrow;
+    nc = ncol;
+    c = new double[nrow * ncol]();
+    return;
 }
 
 /* Assignment:  nonstandard in that it returns void.  To make it standard,
  * replace void -> matrix and uncomment the return *this; */
 
-void matrix::operator=(const matrix & m1)
+void matrix::operator=(const matrix& m1)
 {
-	for (int i = 0; i < nr*nc; i++) c[i] = m1.c[i];
+    for (int i = 0; i < nr * nc; i++)
+        c[i] = m1.c[i];
 }
 
 /* Adding matrices, as a friend */
-matrix operator+(const matrix &m1, const matrix &m2)
+matrix operator+(const matrix& m1, const matrix& m2)
 {
-	assert(m1.nr == m2.nr);
-	assert(m2.nc == m2.nc);
+    assert(m1.nr == m2.nr);
+    assert(m2.nc == m2.nc);
 
-	matrix tm(m1);
-	for (int i = 0; i < m1.nr*m1.nc; i++) tm.c[i] += m2.c[i];
-	return tm;
+    matrix tm(m1);
+    for (int i = 0; i < m1.nr * m1.nc; i++)
+        tm.c[i] += m2.c[i];
+    return tm;
 }
 
 /* Subtracting matrices, as a friend */
-matrix operator-(const matrix &m1, const matrix &m2)
+matrix operator-(const matrix& m1, const matrix& m2)
 {
-	assert(m1.nr == m2.nr);
-	assert(m2.nc == m2.nc);
+    assert(m1.nr == m2.nr);
+    assert(m2.nc == m2.nc);
 
-	matrix tm(m1);
-	for(int i = 0; i < m1.nr*m1.nc; i++) tm.c[i] -= m2.c[i];
-	return tm;
+    matrix tm(m1);
+    for (int i = 0; i < m1.nr * m1.nc; i++)
+        tm.c[i] -= m2.c[i];
+    return tm;
 }
 
 //***************************************
 //
-// Multiplying matrices, as a friend 
+// Multiplying matrices, as a friend
 //
 // *************************************
-matrix operator*(const matrix &m1, const matrix &m2)
+matrix operator*(const matrix& m1, const matrix& m2)
 {
-	assert(m1.nr == m2.nr);
-	assert(m2.nc == m2.nc);
-	
-	// allocate the result and zero it out
-	matrix mprod(m2.nr, m1.nc);
+    assert(m1.nr == m2.nr);
+    assert(m2.nc == m2.nc);
 
-	mprod.zero_out();
+    // allocate the result and zero it out
+    matrix mprod(m2.nr, m1.nc);
 
-	// do the multiply and return
-	for (int i = 0;i < m2.nr;i++)
-	{
-		for (int j = 0;j < m1.nc;j++)
-		{
-			for (int k = 0;k < m2.nc;k++)
-			{
-				mprod(i, j) += m2(i, k) * m1(k, j);
-			}
-		}
-	}
+    mprod.zero_out();
 
-	return mprod;
+    // do the multiply and return
+    for (int i = 0; i < m2.nr; i++)
+    {
+        for (int j = 0; j < m1.nc; j++)
+        {
+            for (int k = 0; k < m2.nc; k++)
+            {
+                mprod(i, j) += m2(i, k) * m1(k, j);
+            }
+        }
+    }
+
+    return mprod;
 }
 
 /* Scale a matrix */
-matrix operator*(const double &s, const matrix &m)
+matrix operator*(const double& s, const matrix& m)
 {
-	matrix sm(m);
-	for (int i = 0; i < m.nr*m.nc; i++) sm.c[i] *= s;
-	return sm;
+    matrix sm(m);
+    for (int i = 0; i < m.nr * m.nc; i++)
+        sm.c[i] *= s;
+    return sm;
 }
 
 /* matrix * double */
-matrix operator*(const matrix &m,const double &s)
+matrix operator*(const matrix& m, const double& s)
 {
-	matrix sm(m);
-	for (int i = 0; i < m.nr*m.nc; i++)sm.c[i] *= s;
-	return sm;
+    matrix sm(m);
+    for (int i = 0; i < m.nr * m.nc; i++)
+        sm.c[i] *= s;
+    return sm;
 }
 
 /* Scale a matrix in place */
-void matrix::operator*=(const double &s)
+void matrix::operator*=(const double& s)
 {
-	for (int i = 0; i < nc*nr; i++) c[i] *= s;
+    for (int i = 0; i < nc * nr; i++)
+        c[i] *= s;
 }
 
 /* Accumulate to a matrix in place */
-void matrix::operator+=(const matrix & m)
+void matrix::operator+=(const matrix& m)
 {
-	if (nr != m.nr || nc != m.nc)
-	{
-		cout << " void matrix::operator+=(const matrix &m) has size mismatch\n";
-		abort();
-	}
-	for (int i = 0; i < nc*nr; i++) c[i] += m.c[i];
+    if (nr != m.nr || nc != m.nc)
+    {
+        cout << " void matrix::operator+=(const matrix &m) has size mismatch\n";
+        abort();
+    }
+    for (int i = 0; i < nc * nr; i++)
+        c[i] += m.c[i];
 }
 
-
 /* decumulate to a matrix in place */
-void matrix::operator-=(const matrix & m)
+void matrix::operator-=(const matrix& m)
 {
-	if (nr != m.nr || nc != m.nc)
-	{
-		cout << "void matrix::operator-=(const matrix &m) has size mismatch\n";
-		abort();
-	}
-	for (int i = 0; i < nc*nr; i++) c[i] -= m.c[i];
+    if (nr != m.nr || nc != m.nc)
+    {
+        cout << "void matrix::operator-=(const matrix &m) has size mismatch\n";
+        abort();
+    }
+    for (int i = 0; i < nc * nr; i++)
+        c[i] -= m.c[i];
 }
 
 /* zero out the matrix */
 void matrix::zero_out(void)
 {
-	for(int i = 0; i < nr*nc; i++)
-	{
-		c[i] = 0.0;
-	}
-	return;
+    for (int i = 0; i < nr * nc; i++)
+    {
+        c[i] = 0.0;
+    }
+    return;
 }
 
-
 matrix transpose(matrix m)
 {
-	matrix tm(m.nc, m.nr);
+    matrix tm(m.nc, m.nr);
 
-	for (int i = 0;i < m.nr;i++)
-	{
-		for (int j = 0;j < m.nc;j++)
-		{
-			tm(j, i) = m(i, j);
-		}
-	}
-	return tm;
+    for (int i = 0; i < m.nr; i++)
+    {
+        for (int j = 0; j < m.nc; j++)
+        {
+            tm(j, i) = m(i, j);
+        }
+    }
+    return tm;
 }
diff --git a/tools/01_NAO_generation/SIAB/src_tools/matrix.h b/tools/01_NAO_generation/SIAB/src_tools/matrix.h
index 5c3e896f555..2232e497c12 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/matrix.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/matrix.h
@@ -2,53 +2,57 @@
 #define MATRIX_H
 class matrix
 {
-	/* data */
-public:
-	double *c;    /* Holds the data */
-	int nr;
-	int nc;   /* Number of rows and columns */
-	int hermetian;
-	static int mCount;
-
-	/* Constructors and destructor */
-	matrix(const int nrows = 1,const int ncols = 1);
-	matrix(const matrix &m1); /* copy constructor */
-	~matrix();
-
-	void create(const int nrow,const int ncol);
-	void operator=(const matrix &m1); /* Nonstandard: returns void */
-
-	double &operator()(const int ir,const int ic)
-	{ return c[nc*ir + ic]; };
-
-	const double &operator()(const int ir,const int ic)const
-	{ return c[nc*ir + ic]; };
-
-//	inline double &operator()(int i,int j) const
-//	    { return c[nc*i+j]; }
-
-	friend matrix operator+(const matrix &m1, const matrix &m2);
-	friend matrix operator-(const matrix &m1, const matrix &m2);
-	friend matrix operator*(const matrix &m1, const matrix &m2);
-	friend matrix operator*(const double &s, const matrix &m);
-	friend matrix operator*(const matrix &m, const double &s);
-	void operator*=(const double &s);
-	void operator+=(const matrix &m);
-	void operator-=(const matrix &m);
-
-	/* member function: */
-	matrix Inverse(void);
-
-	double det(void);
-
-	/* Does the memory allocations of the constructor */
-	void init(int nrows, int ncols);
-
-	/* Free up memory */
-	void freemem(void);
-
-	/* zero out all the entries */
-	void zero_out(void);
+    /* data */
+  public:
+    double* c; /* Holds the data */
+    int nr;
+    int nc; /* Number of rows and columns */
+    int hermetian;
+    static int mCount;
+
+    /* Constructors and destructor */
+    matrix(const int nrows = 1, const int ncols = 1);
+    matrix(const matrix& m1); /* copy constructor */
+    ~matrix();
+
+    void create(const int nrow, const int ncol);
+    void operator=(const matrix& m1); /* Nonstandard: returns void */
+
+    double& operator()(const int ir, const int ic)
+    {
+        return c[nc * ir + ic];
+    };
+
+    const double& operator()(const int ir, const int ic) const
+    {
+        return c[nc * ir + ic];
+    };
+
+    //	inline double &operator()(int i,int j) const
+    //	    { return c[nc*i+j]; }
+
+    friend matrix operator+(const matrix& m1, const matrix& m2);
+    friend matrix operator-(const matrix& m1, const matrix& m2);
+    friend matrix operator*(const matrix& m1, const matrix& m2);
+    friend matrix operator*(const double& s, const matrix& m);
+    friend matrix operator*(const matrix& m, const double& s);
+    void operator*=(const double& s);
+    void operator+=(const matrix& m);
+    void operator-=(const matrix& m);
+
+    /* member function: */
+    matrix Inverse(void);
+
+    double det(void);
+
+    /* Does the memory allocations of the constructor */
+    void init(int nrows, int ncols);
+
+    /* Free up memory */
+    void freemem(void);
+
+    /* zero out all the entries */
+    void zero_out(void);
 };
 
 matrix transpose(matrix m);
diff --git a/tools/01_NAO_generation/SIAB/src_tools/matrix3.cpp b/tools/01_NAO_generation/SIAB/src_tools/matrix3.cpp
index 67ccc3c8f45..b4e2782b239 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/matrix3.cpp
+++ b/tools/01_NAO_generation/SIAB/src_tools/matrix3.cpp
@@ -1,212 +1,247 @@
 #include "matrix3.h"
 
-Matrix3::Matrix3(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)
-{
-	e11 = r11;e12 = r12;e13 = r13;
-	e21 = r21;e22 = r22;e23 = r23;
-	e31 = r31;e32 = r32;e33 = r33;
+Matrix3::Matrix3(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)
+{
+    e11 = r11;
+    e12 = r12;
+    e13 = r13;
+    e21 = r21;
+    e22 = r22;
+    e23 = r23;
+    e31 = r31;
+    e32 = r32;
+    e33 = r33;
 }
 
 void Matrix3::Reset(void)
 {
-	e11 = 1;e12 = 0;e13 = 0;
-	e21 = 0;e22 = 1;e23 = 0;
-	e31 = 0;e32 = 0;e33 = 1;
+    e11 = 1;
+    e12 = 0;
+    e13 = 0;
+    e21 = 0;
+    e22 = 1;
+    e23 = 0;
+    e31 = 0;
+    e32 = 0;
+    e33 = 1;
 }
 
 void Matrix3::Identity(void)
 {
-	e11 = 1;e12 = 0;e13 = 0;
-	e21 = 0;e22 = 1;e23 = 0;
-	e31 = 0;e32 = 0;e33 = 1;
+    e11 = 1;
+    e12 = 0;
+    e13 = 0;
+    e21 = 0;
+    e22 = 1;
+    e23 = 0;
+    e31 = 0;
+    e32 = 0;
+    e33 = 1;
 }
 
-double Matrix3::Det(void) const 
+double Matrix3::Det(void) const
 {
-	return	e11*e22*e33 -
-	        e11*e32*e23 +
-	        e21*e32*e13 -
-	        e21*e12*e33 +
-	        e31*e12*e23 -
-	        e31*e22*e13;
+    return e11 * e22 * e33 - e11 * e32 * e23 + e21 * e32 * e13 - e21 * e12 * e33 + e31 * e12 * e23 - e31 * e22 * e13;
 }
 
 Matrix3 Matrix3::Transpose(void) const
 {
-	return Matrix3(e11, e21, e31, e12, e22, e32, e13, e23, e33);
+    return Matrix3(e11, e21, e31, e12, e22, e32, e13, e23, e33);
 }
 
 Matrix3 Matrix3::Inverse(void) const
 {
-	double d = this->Det();
-
-	if(d == 0)d = 1;
-
-	return Matrix3((e22*e33 - e23*e32) / d,
-	               -(e12*e33 - e13*e32) / d,
-	               (e12*e23 - e13*e22) / d,
-	               -(e21*e33 - e23*e31) / d,
-	               (e11*e33 - e13*e31) / d,
-	               -(e11*e23 - e13*e21) / d,
-	               (e21*e32 - e22*e31) / d,
-	               -(e11*e32 - e12*e31) / d,
-	               (e11*e22 - e12*e21) / d);
-}
-
-Matrix3& Matrix3::operator = (const Matrix3 &m)
-{
-	e11 = m.e11;e12 = m.e12;e13 = m.e13;
-	e21 = m.e21;e22 = m.e22;e23 = m.e23;
-	e31 = m.e31;e32 = m.e32;e33 = m.e33;
-	return *this;
+    double d = this->Det();
+
+    if (d == 0)
+        d = 1;
+
+    return Matrix3((e22 * e33 - e23 * e32) / d,
+                   -(e12 * e33 - e13 * e32) / d,
+                   (e12 * e23 - e13 * e22) / d,
+                   -(e21 * e33 - e23 * e31) / d,
+                   (e11 * e33 - e13 * e31) / d,
+                   -(e11 * e23 - e13 * e21) / d,
+                   (e21 * e32 - e22 * e31) / d,
+                   -(e11 * e32 - e12 * e31) / d,
+                   (e11 * e22 - e12 * e21) / d);
+}
+
+Matrix3& Matrix3::operator=(const Matrix3& m)
+{
+    e11 = m.e11;
+    e12 = m.e12;
+    e13 = m.e13;
+    e21 = m.e21;
+    e22 = m.e22;
+    e23 = m.e23;
+    e31 = m.e31;
+    e32 = m.e32;
+    e33 = m.e33;
+    return *this;
+}
+
+Matrix3& Matrix3::operator+=(const Matrix3& m)
+{
+    e11 += m.e11;
+    e12 += m.e12;
+    e13 += m.e13;
+    e21 += m.e21;
+    e22 += m.e22;
+    e23 += m.e23;
+    e31 += m.e31;
+    e32 += m.e32;
+    e33 += m.e33;
+    return *this;
 }
 
-Matrix3& Matrix3::operator +=(const Matrix3 &m)
+Matrix3& Matrix3::operator-=(const Matrix3& m)
 {
-	e11 += m.e11;e12 += m.e12;e13 += m.e13;
-	e21 += m.e21;e22 += m.e22;e23 += m.e23;
-	e31 += m.e31;e32 += m.e32;e33 += m.e33;
-	return *this;
+    e11 -= m.e11;
+    e12 -= m.e12;
+    e13 -= m.e13;
+    e21 -= m.e21;
+    e22 -= m.e22;
+    e23 -= m.e23;
+    e31 -= m.e31;
+    e32 -= m.e32;
+    e33 -= m.e33;
+    return *this;
 }
 
-Matrix3& Matrix3::operator -=(const Matrix3 &m)
+Matrix3& Matrix3::operator*=(const double& s)
 {
-	e11 -= m.e11;e12 -= m.e12;e13 -= m.e13;
-	e21 -= m.e21;e22 -= m.e22;e23 -= m.e23;
-	e31 -= m.e31;e32 -= m.e32;e33 -= m.e33;
-	return *this;
+    e11 *= s;
+    e12 *= s;
+    e13 *= s;
+    e21 *= s;
+    e22 *= s;
+    e23 *= s;
+    e31 *= s;
+    e32 *= s;
+    e33 *= s;
+    return *this;
 }
 
-Matrix3& Matrix3::operator *=(const double &s)
+Matrix3& Matrix3::operator/=(const double& s)
 {
-	e11 *= s;e12 *= s;e13 *= s;
-	e21 *= s;e22 *= s;e23 *= s;
-	e31 *= s;e32 *= s;e33 *= s;
-	return *this;
+    e11 /= s;
+    e12 /= s;
+    e13 /= s;
+    e21 /= s;
+    e22 /= s;
+    e23 /= s;
+    e31 /= s;
+    e32 /= s;
+    e33 /= s;
+    return *this;
 }
 
-Matrix3& Matrix3::operator /=(const double &s)
+// m1+m2
+Matrix3 operator+(const Matrix3& m1, const Matrix3& m2)
 {
-	e11 /= s;e12 /= s;e13 /= s;
-	e21 /= s;e22 /= s;e23 /= s;
-	e31 /= s;e32 /= s;e33 /= s;
-	return *this;
+    return Matrix3(m1.e11 + m2.e11,
+                   m1.e12 + m2.e12,
+                   m1.e13 + m2.e13,
+                   m1.e21 + m2.e21,
+                   m1.e22 + m2.e22,
+                   m1.e23 + m2.e23,
+                   m1.e31 + m2.e31,
+                   m1.e32 + m2.e32,
+                   m1.e33 + m2.e33);
 }
 
-//m1+m2
-Matrix3 operator +(const Matrix3 &m1, const Matrix3 &m2)
+// m1-m2
+Matrix3 operator-(const Matrix3& m1, const Matrix3& m2)
 {
-	return Matrix3(m1.e11 + m2.e11,
-	               m1.e12 + m2.e12,
-	               m1.e13 + m2.e13,
-	               m1.e21 + m2.e21,
-	               m1.e22 + m2.e22,
-	               m1.e23 + m2.e23,
-	               m1.e31 + m2.e31,
-	               m1.e32 + m2.e32,
-	               m1.e33 + m2.e33);
+    return Matrix3(m1.e11 + m2.e11,
+                   m1.e12 - m2.e12,
+                   m1.e13 - m2.e13,
+                   m1.e21 - m2.e21,
+                   m1.e22 - m2.e22,
+                   m1.e23 - m2.e23,
+                   m1.e31 - m2.e31,
+                   m1.e32 - m2.e32,
+                   m1.e33 - m2.e33);
 }
 
-//m1-m2
-Matrix3 operator -(const Matrix3 &m1, const Matrix3 &m2)
+// m/s
+Matrix3 operator/(const Matrix3& m, const double& s)
 {
-	return Matrix3(m1.e11 + m2.e11,
-	               m1.e12 - m2.e12,
-	               m1.e13 - m2.e13,
-	               m1.e21 - m2.e21,
-	               m1.e22 - m2.e22,
-	               m1.e23 - m2.e23,
-	               m1.e31 - m2.e31,
-	               m1.e32 - m2.e32,
-	               m1.e33 - m2.e33);
+    return Matrix3(m.e11 / s, m.e12 / s, m.e13 / s, m.e21 / s, m.e22 / s, m.e23 / s, m.e31 / s, m.e32 / s, m.e33 / s);
 }
 
-//m/s
-Matrix3 operator /(const Matrix3 &m, const double &s)
+// m1*m2
+Matrix3 operator*(const Matrix3& m1, const Matrix3& m2)
 {
-	return Matrix3(m.e11 / s, m.e12 / s, m.e13 / s,
-	               m.e21 / s, m.e22 / s, m.e23 / s,
-	               m.e31 / s, m.e32 / s, m.e33 / s);
+    return Matrix3(m1.e11 * m2.e11 + m1.e12 * m2.e21 + m1.e13 * m2.e31,
+                   m1.e11 * m2.e12 + m1.e12 * m2.e22 + m1.e13 * m2.e32,
+                   m1.e11 * m2.e13 + m1.e12 * m2.e23 + m1.e13 * m2.e33,
+                   m1.e21 * m2.e11 + m1.e22 * m2.e21 + m1.e23 * m2.e31,
+                   m1.e21 * m2.e12 + m1.e22 * m2.e22 + m1.e23 * m2.e32,
+                   m1.e21 * m2.e13 + m1.e22 * m2.e23 + m1.e23 * m2.e33,
+                   m1.e31 * m2.e11 + m1.e32 * m2.e21 + m1.e33 * m2.e31,
+                   m1.e31 * m2.e12 + m1.e32 * m2.e22 + m1.e33 * m2.e32,
+                   m1.e31 * m2.e13 + m1.e32 * m2.e23 + m1.e33 * m2.e33);
 }
 
-//m1*m2
-Matrix3 operator *(const Matrix3 &m1, const Matrix3 &m2)
+// m*s
+Matrix3 operator*(const Matrix3& m, const double& s)
 {
-	return Matrix3(m1.e11*m2.e11 + m1.e12*m2.e21 + m1.e13*m2.e31,
-	               m1.e11*m2.e12 + m1.e12*m2.e22 + m1.e13*m2.e32,
-	               m1.e11*m2.e13 + m1.e12*m2.e23 + m1.e13*m2.e33,
-	               m1.e21*m2.e11 + m1.e22*m2.e21 + m1.e23*m2.e31,
-	               m1.e21*m2.e12 + m1.e22*m2.e22 + m1.e23*m2.e32,
-	               m1.e21*m2.e13 + m1.e22*m2.e23 + m1.e23*m2.e33,
-	               m1.e31*m2.e11 + m1.e32*m2.e21 + m1.e33*m2.e31,
-	               m1.e31*m2.e12 + m1.e32*m2.e22 + m1.e33*m2.e32,
-	               m1.e31*m2.e13 + m1.e32*m2.e23 + m1.e33*m2.e33);
+    return Matrix3(m.e11 * s, m.e12 * s, m.e13 * s, m.e21 * s, m.e22 * s, m.e23 * s, m.e31 * s, m.e32 * s, m.e33 * s);
 }
 
-//m*s
-Matrix3 operator *(const Matrix3 &m,const double &s)
+// s*m
+Matrix3 operator*(const double& s, const Matrix3& m)
 {
-	return Matrix3(m.e11*s, m.e12*s, m.e13*s,
-	               m.e21*s, m.e22*s, m.e23*s,
-	               m.e31*s, m.e32*s, m.e33*s);
+    return Matrix3(m.e11 * s, m.e12 * s, m.e13 * s, m.e21 * s, m.e22 * s, m.e23 * s, m.e31 * s, m.e32 * s, m.e33 * s);
 }
 
-//s*m
-Matrix3 operator *(const double &s, const Matrix3 &m)
+// m*u
+Vector3<double> operator*(const Matrix3& m, const Vector3<double>& u)
 {
-	return Matrix3(m.e11*s, m.e12*s, m.e13*s,
-	               m.e21*s, m.e22*s, m.e23*s,
-	               m.e31*s, m.e32*s, m.e33*s);
+    return Vector3<double>(m.e11 * u.x + m.e12 * u.y + m.e13 * u.z,
+                           m.e21 * u.x + m.e22 * u.y + m.e23 * u.z,
+                           m.e31 * u.x + m.e32 * u.y + m.e33 * u.z);
 }
 
-//m*u
-Vector3<double> operator *(const Matrix3 &m, const Vector3<double> &u)
+// u*m
+Vector3<double> operator*(const Vector3<double>& u, const Matrix3& m)
 {
-	return Vector3<double>(m.e11*u.x + m.e12*u.y + m.e13*u.z,
-	                       m.e21*u.x + m.e22*u.y + m.e23*u.z,
-	                       m.e31*u.x + m.e32*u.y + m.e33*u.z);
+    return Vector3<double>(u.x * m.e11 + u.y * m.e21 + u.z * m.e31,
+                           u.x * m.e12 + u.y * m.e22 + u.z * m.e32,
+                           u.x * m.e13 + u.y * m.e23 + u.z * m.e33);
 }
 
-//u*m
-Vector3<double> operator *(const Vector3<double> &u, const Matrix3 &m)
+// whether m1==m2
+bool operator==(const Matrix3& m1, const Matrix3& m2)
 {
-	return Vector3<double>(u.x*m.e11 + u.y*m.e21 + u.z*m.e31,
-	                       u.x*m.e12 + u.y*m.e22 + u.z*m.e32,
-	                       u.x*m.e13 + u.y*m.e23 + u.z*m.e33);
+    if (m1.e11 == m2.e11 && m1.e12 == m2.e12 && m1.e13 == m2.e13 && m1.e21 == m2.e21 && m1.e22 == m2.e22
+        && m1.e23 == m2.e23 && m1.e31 == m2.e31 && m1.e32 == m2.e32 && m1.e33 == m2.e33)
+    {
+        return true;
+    }
+    return false;
 }
 
-
-// whether m1==m2
-bool operator==(const Matrix3 &m1, const Matrix3 &m2)
-{
-	if(m1.e11 == m2.e11 &&
-	   m1.e12 == m2.e12 &&
-	   m1.e13 == m2.e13 &&
-	   m1.e21 == m2.e21 &&
-	   m1.e22 == m2.e22 &&
-	   m1.e23 == m2.e23 &&
-	   m1.e31 == m2.e31 &&
-	   m1.e32 == m2.e32 &&
-	   m1.e33 == m2.e33)
-	{
-		return true;
-	}
-	return false;
-}
-
-//whether m1 != m2
-bool operator!=(const Matrix3 &m1, const Matrix3 &m2)
+// whether m1 != m2
+bool operator!=(const Matrix3& m1, const Matrix3& m2)
 {
     return !(m1 == m2); //!= defined in terms of operator ==
 }
 
-
-void Matrix3::print(void)const
+void Matrix3::print(void) const
 {
-	cout << e11 << setw(15) << e12 << setw(15) << e13 << endl ;
-	cout << e21 << setw(15) << e22 << setw(15) << e23 << endl ;
-	cout << e31 << setw(15) << e32 << setw(15) << e33 << endl ;
-	return;
+    cout << e11 << setw(15) << e12 << setw(15) << e13 << endl;
+    cout << e21 << setw(15) << e22 << setw(15) << e23 << endl;
+    cout << e31 << setw(15) << e32 << setw(15) << e33 << endl;
+    return;
 }
diff --git a/tools/01_NAO_generation/SIAB/src_tools/matrix3.h b/tools/01_NAO_generation/SIAB/src_tools/matrix3.h
index 024f19870d7..9fb09ca20ce 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/matrix3.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/matrix3.h
@@ -13,42 +13,47 @@
 
 class Matrix3
 {
-	/* data */
-public:
-	// element eij:i_column,j_row
-	double e11, e12, e13, e21, e22, e23, e31, e32, e33;
-
-	/* Constructors and destructor */
-	Matrix3(const double &r11 = 1,const double &r12 = 0,const double &r13 = 0,
-	        const double &r21 = 0,const double &r22 = 1,const double &r23 = 0,
-	        const double &r31 = 0,const double &r32 = 0,const double &r33 = 1);
-
-	void Reset(void);
-	void Identity(void);
-	double Det(void) const ;
-	Matrix3	Transpose(void) const ;
-	Matrix3	Inverse(void) const ;
-
-	Matrix3& operator=(const Matrix3 &m);
-	Matrix3& operator+=(const Matrix3 &m);
-	Matrix3& operator-=(const Matrix3 &m);
-	Matrix3& operator*=(const double &s);
-	Matrix3& operator/=(const double &s);
-
-	void print(void)const;
+    /* data */
+  public:
+    // element eij:i_column,j_row
+    double e11, e12, e13, e21, e22, e23, e31, e32, e33;
+
+    /* Constructors and destructor */
+    Matrix3(const double& r11 = 1,
+            const double& r12 = 0,
+            const double& r13 = 0,
+            const double& r21 = 0,
+            const double& r22 = 1,
+            const double& r23 = 0,
+            const double& r31 = 0,
+            const double& r32 = 0,
+            const double& r33 = 1);
+
+    void Reset(void);
+    void Identity(void);
+    double Det(void) const;
+    Matrix3 Transpose(void) const;
+    Matrix3 Inverse(void) const;
+
+    Matrix3& operator=(const Matrix3& m);
+    Matrix3& operator+=(const Matrix3& m);
+    Matrix3& operator-=(const Matrix3& m);
+    Matrix3& operator*=(const double& s);
+    Matrix3& operator/=(const double& s);
+
+    void print(void) const;
 };
 
-Matrix3 operator +(const Matrix3 &m1, const Matrix3 &m2);	//m1+m2
-Matrix3 operator -(const Matrix3 &m1, const Matrix3 &m2);	//m1-m2
-Matrix3 operator /(const Matrix3 &m,const double &s);		//m/s
-Matrix3 operator *(const Matrix3 &m1,const  Matrix3 &m2);	//m1*m2
-Matrix3 operator *(const Matrix3 &m,const double &s);		//m*s
-Matrix3 operator *(double &s, const Matrix3 &m);		//s*m
-Vector3<double> operator *(const Matrix3 &m, const Vector3<double> &u);	//m*u
-Vector3<double> operator *(const Vector3<double> &u, const Matrix3 &m);	//u*m
+Matrix3 operator+(const Matrix3& m1, const Matrix3& m2);               // m1+m2
+Matrix3 operator-(const Matrix3& m1, const Matrix3& m2);               // m1-m2
+Matrix3 operator/(const Matrix3& m, const double& s);                  // m/s
+Matrix3 operator*(const Matrix3& m1, const Matrix3& m2);               // m1*m2
+Matrix3 operator*(const Matrix3& m, const double& s);                  // m*s
+Matrix3 operator*(double& s, const Matrix3& m);                        // s*m
+Vector3<double> operator*(const Matrix3& m, const Vector3<double>& u); // m*u
+Vector3<double> operator*(const Vector3<double>& u, const Matrix3& m); // u*m
 
-bool operator ==(const Matrix3 &m1, const Matrix3 &m2); //whether m1 == m2
-bool operator !=(const Matrix3 &m1, const Matrix3 &m2); //whethor m1 != m2
+bool operator==(const Matrix3& m1, const Matrix3& m2); // whether m1 == m2
+bool operator!=(const Matrix3& m1, const Matrix3& m2); // whethor m1 != m2
 
 #endif // MATRIX3_H
-
diff --git a/tools/01_NAO_generation/SIAB/src_tools/realarray.cpp b/tools/01_NAO_generation/SIAB/src_tools/realarray.cpp
index 86c17471dbe..cd5040ff8fd 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/realarray.cpp
+++ b/tools/01_NAO_generation/SIAB/src_tools/realarray.cpp
@@ -2,53 +2,56 @@
  * ESCP:Electro-Structure Calculate Package.
  ********************************************/
 
-#include <cstdlib>
 #include "realarray.h"
 
+#include <cstdlib>
+
 int realArray::arrayCount = 0;
 
 void realArrayAlloc()
 {
-	cout << "\n Allocation error for realArray " << endl;
-	exit(0);
+    cout << "\n Allocation error for realArray " << endl;
+    exit(0);
 }
 
-realArray::realArray(const int d1,const int d2,const int d3)
+realArray::realArray(const int d1, const int d2, const int d3)
 {
-	dim = 3;
-	bound1 = (d1 <= 0) ? 1 : d1;
-	bound2 = (d2 <= 0) ? 1 : d2;
-	bound3 = (d3 <= 0) ? 1 : d3;
-	bound4 = 0;
+    dim = 3;
+    bound1 = (d1 <= 0) ? 1 : d1;
+    bound2 = (d2 <= 0) ? 1 : d2;
+    bound3 = (d3 <= 0) ? 1 : d3;
+    bound4 = 0;
 
-	set_new_handler(realArrayAlloc);
+    set_new_handler(realArrayAlloc);
 
-	size = bound1 * bound2 * bound3 ;	//* sizeof(float);
+    size = bound1 * bound2 * bound3; //* sizeof(float);
 
-	ptr = new double[size]();
-	assert(ptr != 0);
-	for(int i=0; i<size; i++) ptr[i] = 0.0;
+    ptr = new double[size]();
+    assert(ptr != 0);
+    for (int i = 0; i < size; i++)
+        ptr[i] = 0.0;
 
-	++arrayCount;
+    ++arrayCount;
 }
 
-realArray::realArray(const int d1,const int d2,const int d3,const int d4)
+realArray::realArray(const int d1, const int d2, const int d3, const int d4)
 {
-	dim = 4;
-	bound1 = (d1 <= 0) ? 1 : d1;
-	bound2 = (d2 <= 0) ? 1 : d2;
-	bound3 = (d3 <= 0) ? 1 : d3;
-	bound4 = (d4 <= 0) ? 1 : d4;
+    dim = 4;
+    bound1 = (d1 <= 0) ? 1 : d1;
+    bound2 = (d2 <= 0) ? 1 : d2;
+    bound3 = (d3 <= 0) ? 1 : d3;
+    bound4 = (d4 <= 0) ? 1 : d4;
 
-	set_new_handler(realArrayAlloc);
+    set_new_handler(realArrayAlloc);
 
-	size = bound1 * bound2 * bound3 * bound4 ;	//* sizeof(float);
+    size = bound1 * bound2 * bound3 * bound4; //* sizeof(float);
 
-	ptr = new double[size]();
-	assert(ptr != 0);
-	for(int i=0; i<size; i++) ptr[i] = 0.0;
+    ptr = new double[size]();
+    assert(ptr != 0);
+    for (int i = 0; i < size; i++)
+        ptr[i] = 0.0;
 
-	++arrayCount;
+    ++arrayCount;
 }
 
 //********************************
@@ -56,63 +59,67 @@ realArray::realArray(const int d1,const int d2,const int d3,const int d4)
 // Destructor for class realArray
 //
 //********************************
-realArray ::~realArray()
+realArray::~realArray()
 {
     freemem();
 }
 
 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)
+void realArray::create(const int d1, const int d2, const int d3, const int d4)
 {
-	size = d1 * d2 * d3 * d4;
-	assert(size>0);
+    size = d1 * d2 * d3 * d4;
+    assert(size > 0);
 
-	dim = 4;
+    dim = 4;
 
-	bound1 = d1;
-	bound2 = d2;
-	bound3 = d3;
-	bound4 = d4;
+    bound1 = d1;
+    bound2 = d2;
+    bound3 = d3;
+    bound4 = d4;
 
-	delete [] ptr;
-	ptr = new double[size]();
-	assert(ptr != 0);
-	for(int i=0; i<size; i++) ptr[i] = 0.0;
+    delete[] ptr;
+    ptr = new double[size]();
+    assert(ptr != 0);
+    for (int i = 0; i < size; i++)
+        ptr[i] = 0.0;
 }
 
-void realArray::create(const int d1,const int d2,const int d3)
+void realArray::create(const int d1, const int d2, const int d3)
 {
-	size = d1 * d2 * d3;
-	assert(size>0);
+    size = d1 * d2 * d3;
+    assert(size > 0);
 
-	dim = 3;
+    dim = 3;
 
-	bound1 = d1;
-	bound2 = d2;
-	bound3 = d3;
-	bound4 = 1;
+    bound1 = d1;
+    bound2 = d2;
+    bound3 = d3;
+    bound4 = 1;
 
-	delete [] ptr;
-	ptr = new double[size]();
-	assert(ptr != 0);
-	for(int i=0; i<size; i++) ptr[i] = 0.0;
+    delete[] ptr;
+    ptr = new double[size]();
+    assert(ptr != 0);
+    for (int i = 0; i < size; i++)
+        ptr[i] = 0.0;
 }
 
-const realArray &realArray::operator=(const realArray &right)
+const realArray& realArray::operator=(const realArray& right)
 {
-	for (int i = 0;i < size;i++) ptr[i] = right.ptr[i];
-	return *this;// enables x = y = z;
+    for (int i = 0; i < size; i++)
+        ptr[i] = right.ptr[i];
+    return *this; // enables x = y = z;
 }
 
-const realArray &realArray::operator=(const double &right)
+const realArray& realArray::operator=(const double& right)
 {
-	for (int i = 0;i < size;i++) ptr[i] = right;
-	return *this;// enables x = y = z;
+    for (int i = 0; i < size; i++)
+        ptr[i] = right;
+    return *this; // enables x = y = z;
 }
 
 //********************************************************
@@ -121,25 +128,22 @@ const realArray &realArray::operator=(const double &right)
 // const reference return create an cvakue
 //
 //********************************************************
-const double &realArray::operator()
-(const int ind1,const int ind2,const int ind3)const
+const double& realArray::operator()(const int ind1, const int ind2, const int ind3) const
 {
-	const int ind = (ind1 * bound2 + ind2) * bound3 + ind3 ;
-	return ptr[ind];
+    const int ind = (ind1 * bound2 + ind2) * bound3 + ind3;
+    return ptr[ind];
 }
 
-
 //********************************************************
 //
 // overloaded subscript operator for const real Array
 // const reference return creates an cvakue
 //
 //********************************************************
-const double &realArray::operator()
-(const int ind1,const int ind2,const int ind3,const int ind4)const
+const double& realArray::operator()(const int ind1, const int ind2, const int ind3, const int ind4) const
 {
-	const int ind = ((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4;
-	return ptr[ind];
+    const int ind = ((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4;
+    return ptr[ind];
 }
 
 //********************************************************
@@ -148,10 +152,10 @@ const double &realArray::operator()
 // const reference return creates an lvakue
 //
 //********************************************************
-double &realArray::operator()(const int ind1,const int ind2,const int ind3)
+double& realArray::operator()(const int ind1, const int ind2, const int ind3)
 {
-	const int ind = (ind1 * bound2 + ind2) * bound3 + ind3;
-	return ptr[ind]; // reference return
+    const int ind = (ind1 * bound2 + ind2) * bound3 + ind3;
+    return ptr[ind]; // reference return
 }
 
 //********************************************************
@@ -160,10 +164,10 @@ double &realArray::operator()(const int ind1,const int ind2,const int ind3)
 // const reference return creates an lvakue
 //
 //********************************************************
-double &realArray::operator()(const int ind1,const int ind2,const int ind3,const int ind4)
+double& realArray::operator()(const int ind1, const int ind2, const int ind3, const int ind4)
 {
-	const int ind = ((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4;
-	return ptr[ind];// reference return
+    const int ind = ((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4;
+    return ptr[ind]; // reference return
 }
 
 //****************************
@@ -173,7 +177,9 @@ double &realArray::operator()(const int ind1,const int ind2,const int ind3,const
 //****************************
 void realArray::zero_out(void)
 {
-	if (size <= 0) return;
-	for (int i = 0;i < size; i++) ptr[i] = 0;
-	return;
+    if (size <= 0)
+        return;
+    for (int i = 0; i < size; i++)
+        ptr[i] = 0;
+    return;
 }
diff --git a/tools/01_NAO_generation/SIAB/src_tools/realarray.h b/tools/01_NAO_generation/SIAB/src_tools/realarray.h
index 8c9e917f982..1ff521e6400 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/realarray.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/realarray.h
@@ -5,74 +5,89 @@
 #ifndef REALARRAY_H
 #define REALARRAY_H
 
-#include <iostream>
+#include <cassert>
 #include <fstream>
 #include <iomanip>
-#include <cassert>
+#include <iostream>
 
 using namespace std;
 
 class realArray
 {
-public:
-	double *ptr;
+  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();
+    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);
+    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);
+    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);
+    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;
+    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);
+    void zero_out(void);
 
-	int getSize() const
-	{ return size;}
+    int getSize() const
+    {
+        return size;
+    }
 
-	int getDim() const
-	{ return dim;}
+    int getDim() const
+    {
+        return dim;
+    }
 
-	int getBound1() const
-	{ return bound1;}
+    int getBound1() const
+    {
+        return bound1;
+    }
 
-	int getBound2() const
-	{ return bound2;}
+    int getBound2() const
+    {
+        return bound2;
+    }
 
-	int getBound3() const
-	{ return bound3;}
+    int getBound3() const
+    {
+        return bound3;
+    }
 
-	int getBound4() const
-	{ return bound4;}
+    int getBound4() const
+    {
+        return bound4;
+    }
 
-	int getArrayCount(void)
-	{ return arrayCount;}
+    int getArrayCount(void)
+    {
+        return arrayCount;
+    }
 
-private:
-	int size;
-	int dim;
-	int bound1, bound2, bound3, bound4;
-	static int arrayCount;
+  private:
+    int size;
+    int dim;
+    int bound1, bound2, bound3, bound4;
+    static int arrayCount;
 
-	void freemem();
+    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;
 }
 
-#endif	// realArray class
+#endif // realArray class
diff --git a/tools/01_NAO_generation/SIAB/src_tools/timer.cpp b/tools/01_NAO_generation/SIAB/src_tools/timer.cpp
index ec1660d9ce3..62688968e7f 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/timer.cpp
+++ b/tools/01_NAO_generation/SIAB/src_tools/timer.cpp
@@ -7,7 +7,7 @@
 using namespace std;
 
 //----------------------------------------------------------
-// EXPLAIN :   
+// EXPLAIN :
 //----------------------------------------------------------
 bool timer::disabled = false;
 int timer::n_clock = 100;
@@ -18,28 +18,24 @@ double* timer::cpu_start;
 string* timer::name;
 string* timer::class_name;
 int* timer::calls;
-bool timer::delete_flag=false;
+bool timer::delete_flag = false;
 
-timer::timer()
-{
-};
+timer::timer() {};
 
-timer::~timer()
-{
-};
+timer::~timer() {};
 
 void timer::finish()
 {
-	timer::tick("","total");
-	print_all();
-	if(delete_flag)
-	{
-		delete[] cpu_second;
-		delete[] cpu_start;
-		delete[] name;
-		delete[] class_name;
-		delete[] calls;
-	}
+    timer::tick("", "total");
+    print_all();
+    if (delete_flag)
+    {
+        delete[] cpu_second;
+        delete[] cpu_start;
+        delete[] name;
+        delete[] class_name;
+        delete[] calls;
+    }
 }
 
 //----------------------------------------------------------
@@ -47,198 +43,194 @@ void timer::finish()
 //----------------------------------------------------------
 void timer::start(void)
 {
-	cpu_second = new double[n_clock]();
-	cpu_start = new double[n_clock]();
-	name = new string[n_clock]();
-	class_name = new string[n_clock]();
-	calls = new int[n_clock]();
-	delete_flag = true;
-
-	for (int i = 0; i < n_clock; i++)
-	{
-		cpu_start[i] = (double)start_flag;
-		name[i]= "\0";
-		class_name[i]= "\0";
-		calls[i] = 0;
-	}
-
-	// first init ,then we can use tick
-	timer::tick("","total");
-	return;
+    cpu_second = new double[n_clock]();
+    cpu_start = new double[n_clock]();
+    name = new string[n_clock]();
+    class_name = new string[n_clock]();
+    calls = new int[n_clock]();
+    delete_flag = true;
+
+    for (int i = 0; i < n_clock; i++)
+    {
+        cpu_start[i] = (double)start_flag;
+        name[i] = "\0";
+        class_name[i] = "\0";
+        calls[i] = 0;
+    }
+
+    // first init ,then we can use tick
+    timer::tick("", "total");
+    return;
 }
 
 double timer::cpu_time(void)
 {
-	clock_t t1 = 0;
-//----------------------------------------------------------
-// EXPLAIN : here static is important !!
-// only first call can let t0 = 0,clock begin
-// when enter this function second time , t0 > 0
-//----------------------------------------------------------
-	static clock_t t0 = 0;
-	if (t0 == 0) 
-	{
-		t0 = clock();
-	}
-	t1 = clock() - t0;
-
-	if(t1 < 0) return 0; // mohan add, abandon the cross point time 2^32 ~ -2^32 .
-	else return (double)t1/CLOCKS_PER_SEC;
+    clock_t t1 = 0;
+    //----------------------------------------------------------
+    // EXPLAIN : here static is important !!
+    // only first call can let t0 = 0,clock begin
+    // when enter this function second time , t0 > 0
+    //----------------------------------------------------------
+    static clock_t t0 = 0;
+    if (t0 == 0)
+    {
+        t0 = clock();
+    }
+    t1 = clock() - t0;
+
+    if (t1 < 0)
+        return 0; // mohan add, abandon the cross point time 2^32 ~ -2^32 .
+    else
+        return (double)t1 / CLOCKS_PER_SEC;
 }
 
-void timer::tick(const string &tagc,const string &tag)
+void timer::tick(const string& tagc, const string& tag)
 {
-//----------------------------------------------------------
-// EXPLAIN : if timer is disabled , return
-//----------------------------------------------------------
-	if (disabled) 
-	{
-		return;
-	}
-
-	int find_clock=0;//counter
-//----------------------------------------------------------
-// EXPLAIN :  find if the tag has been used 
-//----------------------------------------------------------
-	for(find_clock=0; find_clock<n_now; find_clock++)
-	{
-		if (tag == name[find_clock] && tagc == class_name[find_clock]) 
-		{
-			break;
-		}
-	}
-
-//----------------------------------------------------------
-// EXPLAIN : if it's a new tag; add a new tag to list;
-// add this new tag to name list;
-//----------------------------------------------------------
-	if (find_clock == n_now)
-	{
-		n_now++;
-		name[find_clock] = tag;
-		class_name[find_clock] = tagc;
-	}
-
-//----------------------------------------------------------
-// EXPLAIN : if exceed the the uplimits of list 
-//----------------------------------------------------------
-	if (n_now >= n_clock) 
-	{
-		cout << "\nError! Too many timer!";
-		return;
-	}
-
-//----------------------------------------------------------
-// CALL MEMBER FUNCTION :
-// NAME : cpu_time
-//
-// EXPLAIN : start_flag is minus 1, 
-// so if cpu_start == start_flag,means a new clock counting
-// begin, hence we record the start time of this clock 
-// counting , if cpu_start != start_flag, means it's
-// the end of this counting, so we add the time during
-// this two 'time point'  to the clock time storage.
-//----------------------------------------------------------
-	if (cpu_start[find_clock] == start_flag )
-	{
-		cpu_start[find_clock] = cpu_time();
-		calls[find_clock]++;
-	}
-	else
-	{
-		cpu_second[find_clock] += cpu_time() - cpu_start[find_clock];
-		cpu_start[find_clock] = start_flag;
-	}
-	return;
+    //----------------------------------------------------------
+    // EXPLAIN : if timer is disabled , return
+    //----------------------------------------------------------
+    if (disabled)
+    {
+        return;
+    }
+
+    int find_clock = 0; // counter
+    //----------------------------------------------------------
+    // EXPLAIN :  find if the tag has been used
+    //----------------------------------------------------------
+    for (find_clock = 0; find_clock < n_now; find_clock++)
+    {
+        if (tag == name[find_clock] && tagc == class_name[find_clock])
+        {
+            break;
+        }
+    }
+
+    //----------------------------------------------------------
+    // EXPLAIN : if it's a new tag; add a new tag to list;
+    // add this new tag to name list;
+    //----------------------------------------------------------
+    if (find_clock == n_now)
+    {
+        n_now++;
+        name[find_clock] = tag;
+        class_name[find_clock] = tagc;
+    }
+
+    //----------------------------------------------------------
+    // EXPLAIN : if exceed the the uplimits of list
+    //----------------------------------------------------------
+    if (n_now >= n_clock)
+    {
+        cout << "\nError! Too many timer!";
+        return;
+    }
+
+    //----------------------------------------------------------
+    // CALL MEMBER FUNCTION :
+    // NAME : cpu_time
+    //
+    // EXPLAIN : start_flag is minus 1,
+    // so if cpu_start == start_flag,means a new clock counting
+    // begin, hence we record the start time of this clock
+    // counting , if cpu_start != start_flag, means it's
+    // the end of this counting, so we add the time during
+    // this two 'time point'  to the clock time storage.
+    //----------------------------------------------------------
+    if (cpu_start[find_clock] == start_flag)
+    {
+        cpu_start[find_clock] = cpu_time();
+        calls[find_clock]++;
+    }
+    else
+    {
+        cpu_second[find_clock] += cpu_time() - cpu_start[find_clock];
+        cpu_start[find_clock] = start_flag;
+    }
+    return;
 }
 
-
-
 void timer::enable(void)
 {
-	disabled = false;
-	return;
+    disabled = false;
+    return;
 }
 
 void timer::disable(void)
 {
-	disabled = true;
-	return;
+    disabled = true;
+    return;
 }
 
-double timer::print(const string &tag)
+double timer::print(const string& tag)
 {
-	int index = 0;
-	for(int i=0; i<n_now; i++)
-	{
-		if (tag == name[i]) 
-		{
-			index = i;
-			break;
-		}
-	}
-//	cout << "\n " << name[i] 
-//		 << " time : " 
-//		 << cpu_second[i] << " (sec)" << endl;
-	return cpu_second[index];
+    int index = 0;
+    for (int i = 0; i < n_now; i++)
+    {
+        if (tag == name[i])
+        {
+            index = i;
+            break;
+        }
+    }
+    //	cout << "\n " << name[i]
+    //		 << " time : "
+    //		 << cpu_second[i] << " (sec)" << endl;
+    return cpu_second[index];
 }
 
 long double timer::print_until_now(void)
 {
-	// stop the clock
-	timer::tick("","total");
-	// start again
-	timer::tick("","total");
-	return print("total");
+    // stop the clock
+    timer::tick("", "total");
+    // start again
+    timer::tick("", "total");
+    return print("total");
 }
 
 void timer::print_all()
 {
-//	cout<<"\n timer::print_all()"<<endl;
-	const double small = 0.1; // cpu = 10^6
-	// if want to print > 1s , set small = 10^6
-
-	cout<<"\n--------- CLASS NAME--------------- NAME ----- TIME(sec) -- CALLS ----- AVG ----- PER%";
-	bool *print_flag = new bool[n_clock];
-	for(int i=0; i<n_clock; i++) print_flag[i] = false;
-	
-	for (int i=0; i<n_clock; i++)
-	{
-		int k = 0;
-		double tmp = -1.0;
-		for(int j=0; j<n_clock; j++)
-		{
-			if(print_flag[j]) 
-			{
-				continue;
-			}
-			if(tmp < cpu_second[j])
-			{
-				k = j;
-		//		cout << "\n k = " << k;
-				tmp = cpu_second[j];
-			}
-		}
-		print_flag[k]=true;
-	
-		if ((cpu_second[k] >= 0 && cpu_second[k] < small) ||
-		        (cpu_second[k] <= 0 && cpu_second[k] > -small))
-		{
-			continue;
-		}
-		const long double spend_time = cpu_second[k];
-		const double average_spend_time = spend_time/calls[k];
-
-		cout  << "\n" 
-			 << setw(20) << class_name[k]
-			 << setw(20) << name[k]
-			 << setw(15) << spend_time
-			 << setw(10) << calls[k]
-			 << setw(10) << setprecision(2) << average_spend_time
-			 << setw(10) << spend_time / cpu_second[0] * 100 << "%";
-			
-	}
-	cout<<"\n--------------------------------------------------------------------------------------"<<endl;
-	delete[] print_flag;
-	return;
+    //	cout<<"\n timer::print_all()"<<endl;
+    const double small = 0.1; // cpu = 10^6
+    // if want to print > 1s , set small = 10^6
+
+    cout << "\n--------- CLASS NAME--------------- NAME ----- TIME(sec) -- CALLS ----- AVG ----- PER%";
+    bool* print_flag = new bool[n_clock];
+    for (int i = 0; i < n_clock; i++)
+        print_flag[i] = false;
+
+    for (int i = 0; i < n_clock; i++)
+    {
+        int k = 0;
+        double tmp = -1.0;
+        for (int j = 0; j < n_clock; j++)
+        {
+            if (print_flag[j])
+            {
+                continue;
+            }
+            if (tmp < cpu_second[j])
+            {
+                k = j;
+                //		cout << "\n k = " << k;
+                tmp = cpu_second[j];
+            }
+        }
+        print_flag[k] = true;
+
+        if ((cpu_second[k] >= 0 && cpu_second[k] < small) || (cpu_second[k] <= 0 && cpu_second[k] > -small))
+        {
+            continue;
+        }
+        const long double spend_time = cpu_second[k];
+        const double average_spend_time = spend_time / calls[k];
+
+        cout << "\n"
+             << setw(20) << class_name[k] << setw(20) << name[k] << setw(15) << spend_time << setw(10) << calls[k]
+             << setw(10) << setprecision(2) << average_spend_time << setw(10) << spend_time / cpu_second[0] * 100
+             << "%";
+    }
+    cout << "\n--------------------------------------------------------------------------------------" << endl;
+    delete[] print_flag;
+    return;
 }
diff --git a/tools/01_NAO_generation/SIAB/src_tools/timer.h b/tools/01_NAO_generation/SIAB/src_tools/timer.h
index 826e56ca28c..43070b581a4 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/timer.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/timer.h
@@ -6,10 +6,10 @@
 #define TIMER_H
 
 #include <ctime>
-#include <string>
 #include <fstream>
-#include <iostream>
 #include <iomanip>
+#include <iostream>
+#include <string>
 using namespace std;
 
 //==========================================================
@@ -18,62 +18,59 @@ using namespace std;
 //==========================================================
 class timer
 {
-	timer();
-	~timer();
+    timer();
+    ~timer();
 
-	public:
-//==========================================================
-// MEMBER FUNCTIONS :
-// NAME : tick(use twice at a time)
-// NAME : start
-// NAME : finish
-// NAME : enable
-// NAME : disable
-// NAME : print
-// NAME : print_all
-//==========================================================
-	static void tick(const string &class_name_in,const string &name_in);
+  public:
+    //==========================================================
+    // MEMBER FUNCTIONS :
+    // NAME : tick(use twice at a time)
+    // NAME : start
+    // NAME : finish
+    // NAME : enable
+    // NAME : disable
+    // NAME : print
+    // NAME : print_all
+    //==========================================================
+    static void tick(const string& class_name_in, const string& name_in);
 
-	static void start();
-	static void finish();
+    static void start();
+    static void finish();
 
-	static void enable();
-	static void disable();
+    static void enable();
+    static void disable();
 
-	static void print_all();
-	static long double print_until_now();
-	static double print(const string &name_in);
+    static void print_all();
+    static long double print_until_now();
+    static double print(const string& name_in);
 
-	private:
-
-//==========================================================
-// MEMBER VARIABLES : 
-// NAME : disabled(if disabled , timer can't work)
-// NAME : n_clock(uplimit number of clock numbers)
-// NAME : n_now(the index of clocks) 
-// NAME : start_flag
-// NAME : cpu_start
-// NAME : cput_second(record cpu time)
-// NAME : name(record clock name)
-// NAME : class_name(record another name for clock)
-//==========================================================
-	static bool disabled;
-	static int n_clock;
-	static int n_now;
-	static int start_flag;
-	static double *cpu_start;
-	static double *cpu_second;
-	static string *name;
-	static string *class_name;
-	static int *calls;
-
-//==========================================================
-// MEMBER FUNCTIONS :
-// NAME : cpu_time(calculate time)
-//==========================================================
-	static double cpu_time();
-	static bool delete_flag;
+  private:
+    //==========================================================
+    // MEMBER VARIABLES :
+    // NAME : disabled(if disabled , timer can't work)
+    // NAME : n_clock(uplimit number of clock numbers)
+    // NAME : n_now(the index of clocks)
+    // NAME : start_flag
+    // NAME : cpu_start
+    // NAME : cput_second(record cpu time)
+    // NAME : name(record clock name)
+    // NAME : class_name(record another name for clock)
+    //==========================================================
+    static bool disabled;
+    static int n_clock;
+    static int n_now;
+    static int start_flag;
+    static double* cpu_start;
+    static double* cpu_second;
+    static string* name;
+    static string* class_name;
+    static int* calls;
 
+    //==========================================================
+    // MEMBER FUNCTIONS :
+    // NAME : cpu_time(calculate time)
+    //==========================================================
+    static double cpu_time();
+    static bool delete_flag;
 };
 #endif
-
diff --git a/tools/01_NAO_generation/SIAB/src_tools/vector3.h b/tools/01_NAO_generation/SIAB/src_tools/vector3.h
index 945f117855f..f4d76ba9071 100644
--- a/tools/01_NAO_generation/SIAB/src_tools/vector3.h
+++ b/tools/01_NAO_generation/SIAB/src_tools/vector3.h
@@ -9,229 +9,242 @@
 #endif
 
 #include <cmath>
-#include <iostream>
 #include <iomanip>
+#include <iostream>
 using namespace std;
 
 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);
-
-	Vector3<T>& operator=(const Vector3<T> &u);
-	Vector3<T>& operator+=(const Vector3<T> &u);
-	Vector3<T>& operator-=(const Vector3<T> &u);
-	Vector3<T>& operator*=(const Vector3<T> &u);
-	Vector3<T>& operator*=(const T &s);
-	Vector3<T>& operator/=(const Vector3<T> &u);
-	Vector3<T>& operator/=(const T &s);
-
-	void set(const T &x1, const T &y1,const T &z1);
-
-	inline T norm(void) const
-	{ return sqrt(x*x + y*y + z*z); }
-
-	T norm2(void)const;
-	void normalize(void);
-	void reverse(void);
-
-	// mohan add 2009-11-29
-	void print(void)const ;
+  public:
+    T x;
+    T y;
+    T z;
+    Vector3(const T& x1 = 0, const T& y1 = 0, const T& z1 = 0);
+
+    Vector3<T>& operator=(const Vector3<T>& u);
+    Vector3<T>& operator+=(const Vector3<T>& u);
+    Vector3<T>& operator-=(const Vector3<T>& u);
+    Vector3<T>& operator*=(const Vector3<T>& u);
+    Vector3<T>& operator*=(const T& s);
+    Vector3<T>& operator/=(const Vector3<T>& u);
+    Vector3<T>& operator/=(const T& s);
+
+    void set(const T& x1, const T& y1, const T& z1);
+
+    inline T norm(void) const
+    {
+        return sqrt(x * x + y * y + z * z);
+    }
+
+    T norm2(void) const;
+    void normalize(void);
+    void reverse(void);
+
+    // mohan add 2009-11-29
+    void print(void) const;
 };
 
 template <class T>
-Vector3<T> operator+(const Vector3<T> &u,const Vector3<T> &v);
+Vector3<T> operator+(const Vector3<T>& u, const Vector3<T>& v);
 template <class T>
-Vector3<T> operator-(const Vector3<T> &u,const Vector3<T> &v);
+Vector3<T> operator-(const Vector3<T>& u, const Vector3<T>& v);
 // | 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
+// u.v=(uy*vz-uz*vy)i+(-ux*vz+uz*vx)j+(ux*vy-uy*vx)k
 template <class T>
-Vector3<T> operator^(const Vector3<T> &u,const Vector3<T> &v);
-//u.v=(ux*vx)+(uy*vy)+(uz*vz)
+Vector3<T> operator^(const Vector3<T>& u, const Vector3<T>& v);
+// u.v=(ux*vx)+(uy*vy)+(uz*vz)
 template <class T>
-T operator*(const Vector3<T> &u,const Vector3<T> &v);
+T operator*(const Vector3<T>& u, const Vector3<T>& v);
 template <class T>
-Vector3<T> operator*(const T &s,const Vector3<T> &u);
+Vector3<T> operator*(const T& s, const Vector3<T>& u);
 template <class T>
-Vector3<T> operator/(const Vector3<T> &u,const T &s);
+Vector3<T> operator/(const Vector3<T>& u, const T& s);
 template <class T>
-Vector3<T> cross(const Vector3<T> &u,const Vector3<T> &v);
-//u.v=(ux*vx)+(uy*vy)+(uz*vz)
+Vector3<T> cross(const Vector3<T>& u, const Vector3<T>& v);
+// u.v=(ux*vx)+(uy*vy)+(uz*vz)
 template <class T>
-T dot(const Vector3<T> &u,const Vector3<T> &v);
-//template <class T>
-//T TripleSalarProduct(Vector3<T> u, Vector3<T> v, Vector3<T> w);
+T dot(const Vector3<T>& u, const Vector3<T>& v);
+// template <class T>
+// T TripleSalarProduct(Vector3<T> u, Vector3<T> v, Vector3<T> w);
 
 //==========================================================
 // Define Constructor
 //==========================================================
 template <class T>
-Vector3<T>::Vector3(const T &x1,const T &y1, const T &z1)
-:x(x1),y(y1),z(z1)
-{}
+Vector3<T>::Vector3(const T& x1, const T& y1, const T& z1) : x(x1), y(y1), z(z1)
+{
+}
 
 template <class T>
-void Vector3<T>::set(const T &x1,const T &y1,const T &z1)
+void Vector3<T>::set(const T& x1, const T& y1, const T& z1)
 {
-	x = x1; y = y1; z = z1;
+    x = x1;
+    y = y1;
+    z = z1;
 }
 
 template <class T>
-T Vector3<T>::norm2(void)const
+T Vector3<T>::norm2(void) const
 {
-	return x*x + y*y + z*z;
+    return x * x + y * y + z * z;
 }
 
 //|v| = sqrt(x*x+y*y+z*z)=1
 template <class T>
 void Vector3<T>::normalize(void)
 {
-	T m = sqrt(x*x + y*y + z*z);
-	x /= m;y /= m;z /= m;
+    T m = sqrt(x * x + y * y + z * z);
+    x /= m;
+    y /= m;
+    z /= m;
 }
 
 template <class T>
 void Vector3<T>::reverse(void)
 {
-	x = -x;y = -y;z = -z;
+    x = -x;
+    y = -y;
+    z = -z;
 }
 
 template <class T>
-Vector3<T>& Vector3<T>::operator=(const Vector3<T> &u)
+Vector3<T>& Vector3<T>::operator=(const Vector3<T>& u)
 {
-	x = u.x;y = u.y;z = u.z;
-	return *this;
+    x = u.x;
+    y = u.y;
+    z = u.z;
+    return *this;
 }
 
 template <class T>
-Vector3<T>& Vector3<T>::operator+=(const Vector3<T> &u)
+Vector3<T>& Vector3<T>::operator+=(const Vector3<T>& u)
 {
-	x += u.x;y += u.y;z += u.z;
-	return *this;
+    x += u.x;
+    y += u.y;
+    z += u.z;
+    return *this;
 }
 
 template <class T>
-Vector3<T>& Vector3<T>::operator-=(const Vector3<T> &u)
+Vector3<T>& Vector3<T>::operator-=(const Vector3<T>& u)
 {
-	x -= u.x;y -= u.y;z -= u.z;
-	return *this;
+    x -= u.x;
+    y -= u.y;
+    z -= u.z;
+    return *this;
 }
 
 template <class T>
-Vector3<T>& Vector3<T>::operator*=(const T &s)
+Vector3<T>& Vector3<T>::operator*=(const T& s)
 {
-	x *= s;y *= s;z *= s;
-	return *this;
+    x *= s;
+    y *= s;
+    z *= s;
+    return *this;
 }
 
 template <class T>
-Vector3<T>& Vector3<T>::operator/=(const T &s)
+Vector3<T>& Vector3<T>::operator/=(const T& s)
 {
-	x /= s;y /= s;z /= s;
-	return *this;
+    x /= s;
+    y /= s;
+    z /= s;
+    return *this;
 }
 
 template <class T>
-Vector3<T> operator+(const Vector3<T> &u,const Vector3<T> &v)
+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);
+    return Vector3<T>(u.x + v.x, u.y + v.y, u.z + v.z);
 }
 
 template <class T>
-Vector3<T> operator-(const Vector3<T> &u,const Vector3<T> &v)
+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);
+    return Vector3<T>(u.x - v.x, u.y - v.y, u.z - v.z);
 }
 
 //	u.v=(uy*vz-uz*vy)i+(-ux*vz+uz*vx)j+(ux*vy-uy*vzx)k
 template <class T>
-Vector3<T> operator^(const Vector3<T> &u,const Vector3<T> &v)
+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);
+    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);
 }
 
-//u.v=(ux*vx)+(uy*vy)+(uz*vz)
+// u.v=(ux*vx)+(uy*vy)+(uz*vz)
 template <class T>
-T operator*(const Vector3<T> &u,const Vector3<T> &v)
+T operator*(const Vector3<T>& u, const Vector3<T>& v)
 {
-	return (u.x * v.x + u.y * v.y + u.z * v.z);
+    return (u.x * v.x + u.y * v.y + u.z * v.z);
 }
 
-template <class T>// mohan add 2009-5-10
-Vector3<T> operator*(const Vector3<T> &u, const T &s)
+template <class T> // mohan add 2009-5-10
+Vector3<T> operator*(const Vector3<T>& u, const T& s)
 {
-	return Vector3<T>(u.x * s, u.y * s, u.z * s);
+    return Vector3<T>(u.x * s, u.y * s, u.z * s);
 }
 
 template <class T>
-Vector3<T> operator*(const T &s,const Vector3<T> &u)
+Vector3<T> operator*(const T& s, const Vector3<T>& u)
 {
-	return Vector3<T>(u.x * s, u.y * s, u.z * s);
+    return Vector3<T>(u.x * s, u.y * s, u.z * s);
 }
 
 template <class T>
-Vector3<T> operator /(const Vector3<T> &u,const T &s)
+Vector3<T> operator/(const Vector3<T>& u, const T& s)
 {
-	return Vector3<T>(u.x / s, u.y / s, u.z / s);
+    return Vector3<T>(u.x / s, u.y / s, u.z / s);
 }
 
 //	u.v=(uy*vz-uz*vy)i+(-ux*vz+uz*vx)j+(ux*vy-uy*vzx)k
 template <class T>
-Vector3<T> cross(const Vector3<T> &u, const Vector3<T> &v)
+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);
+    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);
 }
 
 template <class T>
-T dot(const Vector3<T> &u,const Vector3<T> &v)
+T dot(const Vector3<T>& u, const Vector3<T>& v)
 {
-	return (u.x * v.x + u.y * v.y + u.z * v.z);
+    return (u.x * v.x + u.y * v.y + u.z * v.z);
 }
 
-//whether u == v
+// whether u == v
 template <class T>
-bool operator ==(const Vector3<T> &u, const Vector3<T> &v)
+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;
+    if (u.x == v.x && u.y == v.y && u.z == v.z)
+    {
+        return true;
+    }
+    return false;
 }
 
-//whether m1 != m2
+// whether m1 != m2
 template <class T>
-bool operator !=(const Vector3<T> &u, const Vector3<T> &v)
+bool operator!=(const Vector3<T>& u, const Vector3<T>& v)
 {
     return !(u == v);
 }
 
 template <class T>
-void Vector3<T>::print(void)const
+void Vector3<T>::print(void) const
 {
-	cout.precision(5) ;
-	cout << "(" << setw(10) << x << "," << setw(10) << y << ","
-	<< setw(10) << z  << ")"  << endl ;
-	return ;
+    cout.precision(5);
+    cout << "(" << setw(10) << x << "," << setw(10) << y << "," << setw(10) << z << ")" << endl;
+    return;
 }
 
-//s = u.(v x w)
-//template <class T>
-//T TripleScalarProduct(Vector3<T> u, Vector3<T> v, Vector3<T> w)
+// 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)));
-//}
+// }
 
 #endif