From 927c9f03e0822b7ab06f44743fa17dee66af8b6f Mon Sep 17 00:00:00 2001 From: SY Wang Date: Mon, 19 Oct 2026 12:38:41 +0800 Subject: [PATCH] Precommit: Prettify tools directory (Part I) --- .../SIAB/src_parallel/parallel_common.cpp | 135 ++-- .../SIAB/src_parallel/parallel_common.h | 28 +- .../SIAB/src_parallel/parallel_global.cpp | 81 +-- .../SIAB/src_parallel/parallel_global.h | 9 +- .../SIAB/src_parallel/parallel_kpoints.cpp | 486 +++++++------- .../SIAB/src_parallel/parallel_kpoints.h | 42 +- .../SIAB/src_parallel/parallel_reduce.cpp | 222 +++---- .../SIAB/src_parallel/parallel_reduce.h | 50 +- .../SIAB/src_pw/bessel_basis.cpp | 265 ++++---- .../SIAB/src_pw/bessel_basis.h | 45 +- .../SIAB/src_pw/heapsort.cpp | 76 +-- .../01_NAO_generation/SIAB/src_pw/heapsort.h | 6 +- .../SIAB/src_pw/memory_calculation.cpp | 97 +-- .../SIAB/src_pw/memory_calculation.h | 9 +- .../SIAB/src_pw/numerical_basis.cpp | 22 +- .../SIAB/src_pw/numerical_basis.h | 16 +- .../SIAB/src_pw/pw_basis.cpp | 601 +++++++++--------- .../01_NAO_generation/SIAB/src_pw/pw_basis.h | 203 +++--- .../SIAB/src_pw/pw_complement.cpp | 138 ++-- .../SIAB/src_pw/pw_complement.h | 32 +- .../SIAB/src_pw/ylm_real.cpp | 243 ++++--- .../01_NAO_generation/SIAB/src_pw/ylm_real.h | 14 +- .../SIAB/src_unittest/src_tools/common_test.h | 72 +-- .../inverse_matrix_iterate_unittest.h | 45 +- .../inverse_matrix_perturbe_unittest.h | 134 ++-- .../src_tools/inverse_matrix_unittest.h | 28 +- .../src_unittest/src_tools/matrix_multiply.h | 77 ++- tools/deepmd_input_nlist_probe.cpp | 29 +- 28 files changed, 1581 insertions(+), 1624 deletions(-) diff --git a/tools/01_NAO_generation/SIAB/src_parallel/parallel_common.cpp b/tools/01_NAO_generation/SIAB/src_parallel/parallel_common.cpp index bc86f2c2fda..061ebe775ff 100644 --- a/tools/01_NAO_generation/SIAB/src_parallel/parallel_common.cpp +++ b/tools/01_NAO_generation/SIAB/src_parallel/parallel_common.cpp @@ -1,101 +1,106 @@ #include "parallel_common.h" + #include //this is needed because strcpy #ifdef __MPI -void Parallel_Common::bcast_string(string &object) +void Parallel_Common::bcast_string(string& object) { - char swap[100]; - if(MY_RANK == 0) strcpy(swap, object.c_str() ); - MPI_Bcast(swap, 100, MPI_CHAR, 0, MPI_COMM_WORLD); - if(MY_RANK != 0) object = static_cast( swap ); -// cout<<"\n"<(swap); + // cout<<"\n"<( swap ); - } - return; + for (int i = 0; i < n; i++) + { + char swap[100]; + if (MY_RANK == 0) + strcpy(swap, object[i].c_str()); + MPI_Bcast(swap, 100, MPI_CHAR, 0, MPI_COMM_WORLD); + if (MY_RANK != 0) + object[i] = static_cast(swap); + } + return; } -void Parallel_Common::bcast_complex_double(complex &object) +void Parallel_Common::bcast_complex_double(complex& object) { - double a = object.real(); - double b = object.imag(); - Parallel_Common::bcast_double(a); - Parallel_Common::bcast_double(b); - object = std::complex( a, b); - return; + double a = object.real(); + double b = object.imag(); + Parallel_Common::bcast_double(a); + Parallel_Common::bcast_double(b); + object = std::complex(a, b); + return; } -void Parallel_Common::bcast_complex_double(complex *object, const int n) +void Parallel_Common::bcast_complex_double(complex* object, const int n) { - double *a = new double[n]; - double *b = new double[n]; - - for(int i=0; i( a[i], b[i]); - } - delete[] a; - delete[] b; + double* a = new double[n]; + double* b = new double[n]; + + for (int i = 0; i < n; i++) + { + a[i] = object[i].real(); + b[i] = object[i].imag(); + } + + Parallel_Common::bcast_double(a, n); + Parallel_Common::bcast_double(b, n); + + for (int i = 0; i < n; i++) + { + object[i] = std::complex(a[i], b[i]); + } + delete[] a; + delete[] b; } -void Parallel_Common::bcast_double(double &object) +void Parallel_Common::bcast_double(double& object) { - MPI_Bcast(&object, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); -// cout<<"\n"<( swap ); -// cout<<"\n"<(swap); + // cout<<"\n"< *object, const int n); - void bcast_string(string *object,const int n); - void bcast_double(double *object,const int n); - void bcast_int(int *object,const int n); - void bcast_char(char *object,const int n); - - //(2) bcast single - void bcast_complex_double( std::complex &object); - void bcast_string(string &object); - void bcast_double(double &object); - void bcast_int(int &object); - void bcast_bool(bool &object); +//(1) bcast array +void bcast_complex_double(std::complex* object, const int n); +void bcast_string(string* object, const int n); +void bcast_double(double* object, const int n); +void bcast_int(int* object, const int n); +void bcast_char(char* object, const int n); -} +//(2) bcast single +void bcast_complex_double(std::complex& object); +void bcast_string(string& object); +void bcast_double(double& object); +void bcast_int(int& object); +void bcast_bool(bool& object); + +} // namespace Parallel_Common #endif diff --git a/tools/01_NAO_generation/SIAB/src_parallel/parallel_global.cpp b/tools/01_NAO_generation/SIAB/src_parallel/parallel_global.cpp index a6f2323a27d..858df1d0b08 100644 --- a/tools/01_NAO_generation/SIAB/src_parallel/parallel_global.cpp +++ b/tools/01_NAO_generation/SIAB/src_parallel/parallel_global.cpp @@ -3,62 +3,63 @@ // DATE : 2009-11-08 //========================================================== #include "parallel_global.h" + +#include "../src_spillage/tools.h" #include "parallel_common.h" #include "parallel_reduce.h" -#include "../src_spillage/tools.h" using namespace std; #if defined __MPI MPI_Comm POOL_WORLD; -void Parallel_Global::myProd(complex *in, std::complex *inout,int *len,MPI_Datatype *dptr) +void Parallel_Global::myProd(complex* in, std::complex* inout, int* len, MPI_Datatype* dptr) { - for(int i=0;i<*len;i++) - { - (*inout).real()=(*inout).real()+(*in).real(); - (*inout).imag()=(*inout).imag()+(*in).imag(); - in++; - inout++; - } - return; + for (int i = 0; i < *len; i++) + { + (*inout).real() = (*inout).real() + (*in).real(); + (*inout).imag() = (*inout).imag() + (*in).imag(); + in++; + inout++; + } + return; } #endif -void Parallel_Global::read_pal_param(int argc,char **argv) +void Parallel_Global::read_pal_param(int argc, char** argv) { #if defined __MPI -//for test -/* - cout << "\n Hello! Test MPI NOW : argc = "< *in, std::complex *inout,int *len,MPI_Datatype *dptr); +// void myProd(complex *in, std::complex *inout,int *len,MPI_Datatype *dptr); namespace Parallel_Global { - void read_pal_param(int argc, char **argv); +void read_pal_param(int argc, char** argv); #ifdef __MPI - void myProd(complex *in, std::complex *inout,int *len,MPI_Datatype *dptr); +void myProd(complex* in, std::complex* inout, int* len, MPI_Datatype* dptr); #endif -} - +} // namespace Parallel_Global #endif // GMPI diff --git a/tools/01_NAO_generation/SIAB/src_parallel/parallel_kpoints.cpp b/tools/01_NAO_generation/SIAB/src_parallel/parallel_kpoints.cpp index ef0b281f82b..9581cb5453a 100644 --- a/tools/01_NAO_generation/SIAB/src_parallel/parallel_kpoints.cpp +++ b/tools/01_NAO_generation/SIAB/src_parallel/parallel_kpoints.cpp @@ -1,76 +1,61 @@ #include "parallel_kpoints.h" -#include "parallel_global.h" + #include "parallel_common.h" +#include "parallel_global.h" Parallel_Kpoints::Parallel_Kpoints() { - nks_pool = new int[1]; - startk_pool = new int[1]; - whichpool = new int[1]; + nks_pool = new int[1]; + startk_pool = new int[1]; + whichpool = new int[1]; } - + Parallel_Kpoints::~Parallel_Kpoints() { - delete[] nks_pool; - delete[] startk_pool; - delete[] whichpool; + delete[] nks_pool; + delete[] startk_pool; + delete[] whichpool; } void Parallel_Kpoints::init(void) { - //TITLE("Parallel_Kpoints","init"); + // TITLE("Parallel_Kpoints","init"); #ifdef __MPI -//---------------------------------------------------------- -// CALL Function : divide_pools -//---------------------------------------------------------- - this->divide_pools(); - -// for test - - if(MY_RANK==0) - { - cout << "\n " << setw(8) << "MY_RANK" - << setw(8) << "MY_POOL" - << setw(13) << "RANK_IN_POOL" - << setw(6) << "NPROC" - << setw(6) << "KPAR" - << setw(14) << "NPROC_IN_POOL" << endl; - } - - ofs_running << "\n " << setw(8) << "MY_RANK" - << setw(8) << "MY_POOL" - << setw(13) << "RANK_IN_POOL" - << setw(6) << "NPROC" - << setw(6) << "KPAR" - << setw(14) << "NPROC_IN_POOL" << endl; - - for(int i=0;idivide_pools(); + + // for test + + if (MY_RANK == 0) + { + cout << "\n " << setw(8) << "MY_RANK" << setw(8) << "MY_POOL" << setw(13) << "RANK_IN_POOL" << setw(6) + << "NPROC" << setw(6) << "KPAR" << setw(14) << "NPROC_IN_POOL" << endl; + } + + ofs_running << "\n " << setw(8) << "MY_RANK" << setw(8) << "MY_POOL" << setw(13) << "RANK_IN_POOL" << setw(6) + << "NPROC" << setw(6) << "KPAR" << setw(14) << "NPROC_IN_POOL" << endl; + + for (int i = 0; i < NPROC; i++) + { + if (MY_RANK == i) + { + cout << " I'm" << setw(8) << MY_RANK << setw(8) << MY_POOL << setw(13) << RANK_IN_POOL << setw(6) << NPROC + << setw(6) << KPAR << setw(14) << NPROC_IN_POOL << endl; + + ofs_running << " I'm" << setw(8) << MY_RANK << setw(8) << MY_POOL << setw(13) << RANK_IN_POOL << setw(6) + << NPROC << setw(6) << KPAR << setw(14) << NPROC_IN_POOL << endl; + } + MPI_Barrier(MPI_COMM_WORLD); + } + + if (MY_RANK != 0) + { + cout.rdbuf(NULL); + } + + return; #endif } @@ -78,226 +63,223 @@ void Parallel_Kpoints::divide_pools(void) { #ifdef __MPI - //cout<<"\n ==> mpi_split()"<= startpro_pool[i]) - { - MY_POOL=i; - } - } - - int key = 1; - RANK_IN_POOL = MY_RANK-startpro_pool[MY_POOL]; - - //======================================================== - // MPI_Comm_Split: Creates new communicators based on - // colors(2nd parameter) and keys(3rd parameter) - // Note: The color must be non-negative or MPI_UNDEFINED. - //======================================================== - MPI_Comm_split(MPI_COMM_WORLD,MY_POOL,key,&POOL_WORLD); + // cout<<"\n ==> mpi_split()"<= startpro_pool[i]) + { + MY_POOL = i; + } + } + + int key = 1; + RANK_IN_POOL = MY_RANK - startpro_pool[MY_POOL]; + + //======================================================== + // MPI_Comm_Split: Creates new communicators based on + // colors(2nd parameter) and keys(3rd parameter) + // Note: The color must be non-negative or MPI_UNDEFINED. + //======================================================== + MPI_Comm_split(MPI_COMM_WORLD, MY_POOL, key, &POOL_WORLD); #endif - - return; -} + return; +} // the kpoints here are reduced after symmetry applied. -void Parallel_Kpoints::kinfo(int &nkstot) +void Parallel_Kpoints::kinfo(int& nkstot) { #ifdef __MPI Parallel_Common::bcast_int(nkstot); - this->get_nks_pool(nkstot); - this->get_startk_pool(nkstot); - this->get_whichpool(nkstot); + this->get_nks_pool(nkstot); + this->get_startk_pool(nkstot); + this->get_whichpool(nkstot); #endif - return; + return; } -void Parallel_Kpoints::get_whichpool(const int &nkstot) +void Parallel_Kpoints::get_whichpool(const int& nkstot) { - delete[] whichpool; - this->whichpool = new int[nkstot]; - ZEROS(whichpool, nkstot); - - for(int i=0; inks_pool[i]; ik++) - { - const int k_now = ik + startk_pool[i]; - this->whichpool[k_now] = i; - //ofs_running << "\n whichpool[" << k_now <<"] = " << whichpool[k_now]; - } - } - - return; + delete[] whichpool; + this->whichpool = new int[nkstot]; + ZEROS(whichpool, nkstot); + + for (int i = 0; i < KPAR; i++) + { + for (int ik = 0; ik < this->nks_pool[i]; ik++) + { + const int k_now = ik + startk_pool[i]; + this->whichpool[k_now] = i; + // ofs_running << "\n whichpool[" << k_now <<"] = " << whichpool[k_now]; + } + } + + return; } -void Parallel_Kpoints::get_nks_pool(const int &nkstot) +void Parallel_Kpoints::get_nks_pool(const int& nkstot) { - delete[] nks_pool; - this->nks_pool = new int[KPAR]; - ZEROS(nks_pool, KPAR); - - const int nks_ave = nkstot/KPAR; - const int remain = nkstot%KPAR; - - //ofs_running << "\n nkstot = " << nkstot; - //ofs_running << "\n KPAR = " << KPAR; - //ofs_running << "\n nks_ave = " << nks_ave; - - for(int i=0; inks_pool[i] = nks_ave; - if(i &object) +void Parallel_Reduce::reduce_complex_double_all(complex& object) { #ifdef __MPI MPI_Allreduce(MPI_IN_PLACE, &object, 1, MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD); @@ -73,126 +71,128 @@ void Parallel_Reduce::reduce_complex_double_all(complex &object) return; } -void Parallel_Reduce::reduce_complex_double_pool(complex &object) +void Parallel_Reduce::reduce_complex_double_pool(complex& object) { #ifdef __MPI - MPI_Allreduce(MPI_IN_PLACE, &object, 1, MPI_DOUBLE_COMPLEX, MPI_SUM, POOL_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &object, 1, MPI_DOUBLE_COMPLEX, MPI_SUM, POOL_WORLD); #endif - return; + return; } -void Parallel_Reduce::reduce_complex_double_pool(complex *object, const int n) +void Parallel_Reduce::reduce_complex_double_pool(complex* object, const int n) { #ifdef __MPI - MPI_Allreduce(MPI_IN_PLACE, object, n, MPI_DOUBLE_COMPLEX, MPI_SUM, POOL_WORLD); + MPI_Allreduce(MPI_IN_PLACE, object, n, MPI_DOUBLE_COMPLEX, MPI_SUM, POOL_WORLD); #endif - return; + return; } -void Parallel_Reduce::gather_min_int_all(int &v) +void Parallel_Reduce::gather_min_int_all(int& v) { #ifdef __MPI - int *all = new int[NPROC]; - MPI_Allgather(&v, 1, MPI_INT, all, 1, MPI_INT, MPI_COMM_WORLD); - for(int i=0; iall[i]) - { - v = all[i]; - } - } - delete[] all; + int* all = new int[NPROC]; + MPI_Allgather(&v, 1, MPI_INT, all, 1, MPI_INT, MPI_COMM_WORLD); + for (int i = 0; i < NPROC; i++) + { + if (v > all[i]) + { + v = all[i]; + } + } + delete[] all; #endif } -void Parallel_Reduce::gather_max_double_all(double &v) +void Parallel_Reduce::gather_max_double_all(double& v) { #ifdef __MPI - double *value=new double[NPROC]; - ZEROS(value, NPROC); - MPI_Allgather(&v, 1, MPI_DOUBLE, value, 1, MPI_DOUBLE, MPI_COMM_WORLD); - for(int i=0; ivalue[i]) - { - v = value[i]; - } - } - delete[] value; + if (NPROC_IN_POOL == 1) + return; + double* value = new double[NPROC_IN_POOL]; + ZEROS(value, NPROC_IN_POOL); + MPI_Allgather(&v, 1, MPI_DOUBLE, value, 1, MPI_DOUBLE, POOL_WORLD); + for (int i = 0; i < NPROC_IN_POOL; i++) + { + if (v > value[i]) + { + v = value[i]; + } + } + delete[] value; #endif } -void Parallel_Reduce::gather_min_double_all(double &v) +void Parallel_Reduce::gather_min_double_all(double& v) { #ifdef __MPI - double *value=new double[NPROC]; - ZEROS(value, NPROC); - MPI_Allgather(&v, 1, MPI_DOUBLE, value, 1, MPI_DOUBLE, MPI_COMM_WORLD); - for(int i=0; ivalue[i]) - { - v = value[i]; - } - } - delete[] value; + double* value = new double[NPROC]; + ZEROS(value, NPROC); + MPI_Allgather(&v, 1, MPI_DOUBLE, value, 1, MPI_DOUBLE, MPI_COMM_WORLD); + for (int i = 0; i < NPROC; i++) + { + if (v > value[i]) + { + v = value[i]; + } + } + delete[] value; #endif } -bool Parallel_Reduce::check_if_equal(double &v) +bool Parallel_Reduce::check_if_equal(double& v) { #ifdef __MPI - double *all=new double[NPROC]; - MPI_Allgather(&v, 1, MPI_DOUBLE, all, 1, MPI_DOUBLE, MPI_COMM_WORLD); - for(int i=0; i 1.0e-9 ) - { - for(int j=0; jDk = dk; + TITLE("Bessel_Basis", "init"); + this->Dk = dk; - // to make a table - this->init_TableOne( SMOOTH, SIGMA, ECUT, RCUT, dr, Dk, LMAXALL, NE, TOLERENCE); + // to make a table + this->init_TableOne(SMOOTH, SIGMA, ECUT, RCUT, dr, Dk, LMAXALL, NE, TOLERENCE); - return; + return; } -double Bessel_Basis::Polynomial_Interpolation2 - (const int &l, const int &ie, const double &gnorm)const +double Bessel_Basis::Polynomial_Interpolation2(const int& l, const int& ie, const double& gnorm) const { - const double position = gnorm / this->Dk; - const int iq = static_cast(position); -// if(iq 0.0); - assert(dr > 0.0); - assert(dk > 0.0); - - // init kmesh - this->kmesh = static_cast(sqrt(ecutwfc) / dk) +1 + 4; - if (kmesh % 2 == 0)++kmesh; - cout << "\n kmesh = " << kmesh; - - // init Table One - this->TableOne.create(lmax+1, ecut_number, kmesh); - - // init rmesh - int rmesh = static_cast( rcut / dr ) + 4; - if (rmesh % 2 == 0) ++rmesh; - cout << "\n rmesh = " << rmesh; - - // allocate rmesh and Jlk and eigenvalue of Jlq - double *r = new double[rmesh]; - double *rab = new double[rmesh]; - double *jle = new double[rmesh]; - double *jlk = new double[rmesh]; - double *g = new double[rmesh]; // smooth function - double *function = new double[rmesh]; - double *en = new double[ecut_number]; - - for(int ir=0; ir(ir) * dr; - rab[ir] = dr; - if(smooth_in) - { - g[ir] = 1.0 - std::exp(-( (r[ir]-rcut)*(r[ir]-rcut)/2.0/sigma_in/sigma_in ) ); - } - } - - // init eigenvalue of Jl - for(int l=0; lTableOne(l, ie, ik) ); - } - - }// end ie - }// end ; - - delete[] en; - delete[] jle; - delete[] jlk; - delete[] rab; - delete[] g; - delete[] r; - delete[] function; - timer::tick("Spillage","TableOne"); - return; + TITLE("Bessel_Basis", "init_TableOne"); + timer::tick("Spillage", "TableOne"); + // check + assert(ecutwfc > 0.0); + assert(dr > 0.0); + assert(dk > 0.0); + + // init kmesh + this->kmesh = static_cast(sqrt(ecutwfc) / dk) + 1 + 4; + if (kmesh % 2 == 0) + ++kmesh; + cout << "\n kmesh = " << kmesh; + + // init Table One + this->TableOne.create(lmax + 1, ecut_number, kmesh); + + // init rmesh + int rmesh = static_cast(rcut / dr) + 4; + if (rmesh % 2 == 0) + ++rmesh; + cout << "\n rmesh = " << rmesh; + + // allocate rmesh and Jlk and eigenvalue of Jlq + double* r = new double[rmesh]; + double* rab = new double[rmesh]; + double* jle = new double[rmesh]; + double* jlk = new double[rmesh]; + double* g = new double[rmesh]; // smooth function + double* function = new double[rmesh]; + double* en = new double[ecut_number]; + + for (int ir = 0; ir < rmesh; ir++) + { + r[ir] = static_cast(ir) * dr; + rab[ir] = dr; + if (smooth_in) + { + g[ir] = 1.0 - std::exp(-((r[ir] - rcut) * (r[ir] - rcut) / 2.0 / sigma_in / sigma_in)); + } + } + + // init eigenvalue of Jl + for (int l = 0; l < lmax + 1; l++) + { + ZEROS(en, ecut_number); + ZEROS(jle, rmesh); + ZEROS(jlk, rmesh); + + // calculate eigenvalue for l + Mathzone::Spherical_Bessel_Roots(ecut_number, l, tolerence, en, rcut); + + // for each eigenvalue + for (int ie = 0; ie < ecut_number; ie++) + { + // calculate J_{l}( en[ir]*r) + Mathzone::Spherical_Bessel(rmesh, r, en[ie], l, jle); + + for (int ir = 0; ir < rmesh; ir++) + { + jle[ir] = jle[ir] * r[ir] * r[ir]; + } + + // mohan add 2009-08-28 + if (smooth_in) + { + for (int ir = 0; ir < rmesh; ir++) + { + jle[ir] *= g[ir]; + } + } + + for (int ik = 0; ik < kmesh; ik++) + { + // calculate J_{l}( ik*dk*r ) + Mathzone::Spherical_Bessel(rmesh, r, ik * dk, l, jlk); + + // calculate the function will be integrated + for (int ir = 0; ir < rmesh; ir++) + { + function[ir] = jle[ir] * jlk[ir]; + } + + // make table value + Mathzone::Simpson_Integral(rmesh, function, rab, this->TableOne(l, ie, ik)); + } + + } // end ie + } // end ; + + delete[] en; + delete[] jle; + delete[] jlk; + delete[] rab; + delete[] g; + delete[] r; + delete[] function; + timer::tick("Spillage", "TableOne"); + return; } diff --git a/tools/01_NAO_generation/SIAB/src_pw/bessel_basis.h b/tools/01_NAO_generation/SIAB/src_pw/bessel_basis.h index 0857edffa53..f2094873e53 100644 --- a/tools/01_NAO_generation/SIAB/src_pw/bessel_basis.h +++ b/tools/01_NAO_generation/SIAB/src_pw/bessel_basis.h @@ -13,35 +13,32 @@ //========================================================== class Bessel_Basis { -public: - Bessel_Basis(); - ~Bessel_Basis(); + public: + Bessel_Basis(); + ~Bessel_Basis(); - // used for a specific group of C4 coefficients. - void init( - const double &dk = 0.01, - const double &dr = 0.01); + // used for a specific group of C4 coefficients. + void init(const double& dk = 0.01, const double& dr = 0.01); - // get value from interpolation - double Polynomial_Interpolation2(const int &l, const int &ie, const double &gnorm)const; + // get value from interpolation + double Polynomial_Interpolation2(const int& l, const int& ie, const double& gnorm) const; - realArray TableOne; + realArray TableOne; -private: - int kmesh; - double Dk; + private: + int kmesh; + double Dk; - // init table, used for output overlap Q. - void init_TableOne( - const bool smooth_in, - const double &sigma_in, - const double &ecut, - const double &rcut, - const double &dr, - const double &dk, - const int &lmax, - const int &ecut_number, - const double &tolerence); + // init table, used for output overlap Q. + void init_TableOne(const bool smooth_in, + const double& sigma_in, + const double& ecut, + const double& rcut, + const double& dr, + const double& dk, + const int& lmax, + const int& ecut_number, + const double& tolerence); }; #endif diff --git a/tools/01_NAO_generation/SIAB/src_pw/heapsort.cpp b/tools/01_NAO_generation/SIAB/src_pw/heapsort.cpp index 67ecac6489b..6506e31b207 100644 --- a/tools/01_NAO_generation/SIAB/src_pw/heapsort.cpp +++ b/tools/01_NAO_generation/SIAB/src_pw/heapsort.cpp @@ -1,18 +1,21 @@ -#include "../src_spillage/common.h" #include "heapsort.h" -void heapAjust(double *r, int *ind, int s, int m) +#include "../src_spillage/common.h" + +void heapAjust(double* r, int* ind, int s, int m) { int j, ic; double rc; rc = r[s]; ic = ind[s]; - for (j = 2 * s;j <= m;j *= 2) + for (j = 2 * s; j <= m; j *= 2) { - if (j < m && (r[j] < r[j+1])) j++; + if (j < m && (r[j] < r[j + 1])) + j++; - if (!(rc < r[j])) break; + if (!(rc < r[j])) + break; r[s] = r[j]; @@ -27,26 +30,26 @@ void heapAjust(double *r, int *ind, int s, int m) return; } -void heapsort(const int n, double *r, int *ind) +void heapsort(const int n, double* r, int* ind) { - timer::tick("mymath","heapsort"); + timer::tick("mymath", "heapsort"); int i, ic; double rc; if (ind[0] == 0) { - for (i = 0;i < n;i++) + for (i = 0; i < n; i++) { ind[i] = i; } } - for (i = n / 2;i >= 0;i--) + for (i = n / 2; i >= 0; i--) { heapAjust(r, ind, i, n - 1); } - for (i= n - 1;i > 0;i--) + for (i = n - 1; i > 0; i--) { rc = r[0]; r[0] = r[i]; @@ -56,7 +59,7 @@ void heapsort(const int n, double *r, int *ind) ind[i] = ic; heapAjust(r, ind, 0, i - 1); } - timer::tick("mymath","heapsort"); + timer::tick("mymath", "heapsort"); return; } @@ -81,67 +84,68 @@ c adapted from Numerical Recipes pg. 329 (new edition) *********************************************************************/ // from hpsort.f90 -void hpsort(int n, double *ra, int *ind) +void hpsort(int n, double* ra, int* ind) { int i, ir, j, k, iind; double rra; if (ind[1] == 0) { - for (i = 1;i <= n;i++) + for (i = 1; i <= n; i++) ind[i] = i; } - if (n < 2) return; // nothing to order + if (n < 2) + return; // nothing to order - k = n / 2 + 1; + k = n / 2 + 1; ir = n; while (true) { - if (k > 1) // still in hiring phase + if (k > 1) // still in hiring phase { - k = k - 1; + k = k - 1; rra = ra[k]; iind = ind[k]; } - else // in retirement-promotion phase. + else // in retirement-promotion phase. { - rra = ra[ir]; // clear a space at the end of the array - iind = ind[ir]; // - ra[ir] = ra[1]; // retire the top of the heap into it - ind[ir] = ind[1]; // - ir = ir - 1; // decrease the size of the corporation + rra = ra[ir]; // clear a space at the end of the array + iind = ind[ir]; // + ra[ir] = ra[1]; // retire the top of the heap into it + ind[ir] = ind[1]; // + ir = ir - 1; // decrease the size of the corporation - if (ir == 1) // done with the last promotion + if (ir == 1) // done with the last promotion { - ra[1] = rra; // the least competent worker at all // - ind[1] = iind; // + ra[1] = rra; // the least competent worker at all // + ind[1] = iind; // return; } } - i = k; // wheter in hiring or promotion phase, we + i = k; // wheter in hiring or promotion phase, we - j = k + k; // set up to place rra in its proper level + j = k + k; // set up to place rra in its proper level while (j <= ir) { if (j < ir) { - if (ra[j] < ra[j+1]) // compare to better underling + if (ra[j] < ra[j + 1]) // compare to better underling { j = j + 1; } - else if (ra[j] == ra[j+1]) + else if (ra[j] == ra[j + 1]) { - if (ind[j] < ind[j+1]) + if (ind[j] < ind[j + 1]) j = j + 1; } } - if (rra < ra[j]) // demote rra + if (rra < ra[j]) // demote rra { ra[i] = ra[j]; ind[i] = ind[j]; @@ -150,7 +154,7 @@ void hpsort(int n, double *ra, int *ind) } else if (rra == ra[j]) { - if (iind < ind[j]) // demote rra + if (iind < ind[j]) // demote rra { ra[i] = ra[j]; ind[i] = ind[j]; @@ -158,10 +162,10 @@ void hpsort(int n, double *ra, int *ind) j = j + j; } else - j = ir + 1; // set j to terminate do-while loop + j = ir + 1; // set j to terminate do-while loop } - else // this is the right place for rra - j = ir + 1; // set j to terminate do-while loop + else // this is the right place for rra + j = ir + 1; // set j to terminate do-while loop } ra[i] = rra; diff --git a/tools/01_NAO_generation/SIAB/src_pw/heapsort.h b/tools/01_NAO_generation/SIAB/src_pw/heapsort.h index ae40098c5a2..96697fa6e4b 100644 --- a/tools/01_NAO_generation/SIAB/src_pw/heapsort.h +++ b/tools/01_NAO_generation/SIAB/src_pw/heapsort.h @@ -2,10 +2,8 @@ #define HEAPSROT_H #include "../src_spillage/common.h" -void heapsort(int n, double *r, int *ind); +void heapsort(int n, double* r, int* ind); void heapAjust(double r[], int ind[], int s, int m); -void hpsort(int n, double *ra, int *ind); +void hpsort(int n, double* ra, int* ind); #endif // HEAPSORT_H - - diff --git a/tools/01_NAO_generation/SIAB/src_pw/memory_calculation.cpp b/tools/01_NAO_generation/SIAB/src_pw/memory_calculation.cpp index 70ce802cd20..4ec3a584eaa 100644 --- a/tools/01_NAO_generation/SIAB/src_pw/memory_calculation.cpp +++ b/tools/01_NAO_generation/SIAB/src_pw/memory_calculation.cpp @@ -3,64 +3,65 @@ #include "../src_tools/complexmatrix_inline.h" Memory::Memory() -{} +{ +} Memory::~Memory() -{} +{ +} void Memory::calculation(void) { - TITLE("Memory","calculation"); - - const double dk = 0.01; - int kmesh = static_cast(sqrt(ECUT) / dk) +1 + 4; - if (kmesh % 2 == 0)++kmesh; + TITLE("Memory", "calculation"); + + const double dk = 0.01; + int kmesh = static_cast(sqrt(ECUT) / dk) + 1 + 4; + if (kmesh % 2 == 0) + ++kmesh; + + const int n1 = NTYPE; + const int n2 = LMAXUSED + 1; + const int n3 = NMAXUSED; + const int n4 = kmesh; + const int n1234 = n1 * n2 * n3 * n4; - const int n1 = NTYPE; - const int n2 = LMAXUSED + 1; - const int n3 = NMAXUSED; - const int n4 = kmesh; - const int n1234 = n1 * n2 * n3 * n4; + cout << "\n psi1d dimension(NTYPE, LMAXUSED+1, NMAXUSED, kmesh) : " << n1234 << " = " << n1 << " * " << n2 << " * " + << n3 << " * " << n4 << " = " << n1234 * 8 / (double)1024 / (double)1024 << " MB"; - cout << "\n psi1d dimension(NTYPE, LMAXUSED+1, NMAXUSED, kmesh) : " << n1234 - << " = " << n1 << " * " << n2 << " * " - << n3 << " * " << n4 << " = " << n1234 * 8 / (double)1024 / (double)1024 << " MB"; + for (int i = 0; i < mz.nlevel; i++) + { + int nwfc2 = 0; + for (int it = 0; it < NTYPE; it++) + { + for (int l = 0; l < mz.lmax_type[it] + 1; l++) + { + nwfc2 += NA[it] * (2 * l + 1); + cout << "\n na=" << NA[it] << " 2*l+1=" << 2 * l + 1 << " n=1"; + } + } + cout << "\n newfc2=" << nwfc2; + cout << "\n psi3d dimension(NKSTOT=" << NKSTOT << ", nwfc2=" << nwfc2 << ", npwx=" << PW.npwx + << ")=" << NKSTOT * nwfc2 * PW.npwx * 16 / (double)1024 / (double)1024 << " MB"; - for(int i=0; i sum = std::complex(0, 0); + for (int t = 0; t < nwfc2 * nwfc2; t++) + { + timer::tick("Memory", "test"); + for (int j = 0; j < PW.npwx; j++) + { + // sum += conj(haha[j]) * haha[j]; + sum += conj(haha(2, j)) * haha(2, j); + } + timer::tick("Memory", "test"); + } + } - return; + return; } diff --git a/tools/01_NAO_generation/SIAB/src_pw/memory_calculation.h b/tools/01_NAO_generation/SIAB/src_pw/memory_calculation.h index 6a669b10d46..b43929b07fe 100644 --- a/tools/01_NAO_generation/SIAB/src_pw/memory_calculation.h +++ b/tools/01_NAO_generation/SIAB/src_pw/memory_calculation.h @@ -4,12 +4,11 @@ class Memory { - public: - Memory(); - ~Memory(); - - void calculation(void); + public: + Memory(); + ~Memory(); + void calculation(void); }; #endif diff --git a/tools/01_NAO_generation/SIAB/src_pw/numerical_basis.cpp b/tools/01_NAO_generation/SIAB/src_pw/numerical_basis.cpp index af2eabf7b18..734b686c879 100644 --- a/tools/01_NAO_generation/SIAB/src_pw/numerical_basis.cpp +++ b/tools/01_NAO_generation/SIAB/src_pw/numerical_basis.cpp @@ -1,14 +1,19 @@ -#include "ylm_real.h" #include "numerical_basis.h" +#include "ylm_real.h" + Bessel_Basis Numerical_Basis::bessel_basis; -Numerical_Basis::Numerical_Basis(){} +Numerical_Basis::Numerical_Basis() +{ +} -Numerical_Basis::~Numerical_Basis(){} +Numerical_Basis::~Numerical_Basis() +{ +} //============================================================ -// MEMBER FUNCTION : +// MEMBER FUNCTION : // NAME : init // Firstly, use check(0) to call bessel_basis,init // to generate TableOne @@ -18,10 +23,9 @@ Numerical_Basis::~Numerical_Basis(){} // The function is called in run_fp.cpp. void Numerical_Basis::init_table(void) { - TITLE("Numerical_Basis","output_overlap"); + TITLE("Numerical_Basis", "output_overlap"); - Numerical_Basis::bessel_basis.init(); - - return; -} + Numerical_Basis::bessel_basis.init(); + return; +} diff --git a/tools/01_NAO_generation/SIAB/src_pw/numerical_basis.h b/tools/01_NAO_generation/SIAB/src_pw/numerical_basis.h index 12aab7f6745..5719b11551f 100644 --- a/tools/01_NAO_generation/SIAB/src_pw/numerical_basis.h +++ b/tools/01_NAO_generation/SIAB/src_pw/numerical_basis.h @@ -9,21 +9,19 @@ #include "bessel_basis.h" //========================================================== // CLASS : -// NAME : Numerical_Basis +// NAME : Numerical_Basis //========================================================== class Numerical_Basis { - public: - Numerical_Basis(); - ~Numerical_Basis(); + public: + Numerical_Basis(); + ~Numerical_Basis(); - void init_table(void); - - static Bessel_Basis bessel_basis; - - private: + void init_table(void); + static Bessel_Basis bessel_basis; + private: }; #endif diff --git a/tools/01_NAO_generation/SIAB/src_pw/pw_basis.cpp b/tools/01_NAO_generation/SIAB/src_pw/pw_basis.cpp index 20b81df804d..18dcb3f4fd4 100644 --- a/tools/01_NAO_generation/SIAB/src_pw/pw_basis.cpp +++ b/tools/01_NAO_generation/SIAB/src_pw/pw_basis.cpp @@ -1,9 +1,10 @@ #include "pw_basis.h" -#include "pw_complement.h" -#include "../src_spillage/tools.h" -#include "ylm_real.h" + #include "../src_parallel/parallel_reduce.h" +#include "../src_spillage/tools.h" #include "../src_tools/complexmatrix_inline.h" +#include "pw_complement.h" +#include "ylm_real.h" PW_Basis::PW_Basis() { @@ -27,40 +28,40 @@ PW_Basis::PW_Basis() ylm = new matrix[1]; Jlq = new std::complex[1]; - save = new std::complex[1]; - posmu = new int[1]; - posnu = new int[1]; + save = new std::complex[1]; + posmu = new int[1]; + posnu = new int[1]; } PW_Basis::~PW_Basis() { - delete [] gcar_global; - delete [] gdirect_global; - delete [] gg_global; + delete[] gcar_global; + delete[] gdirect_global; + delete[] gg_global; #ifdef __MPI - delete [] gcar; - delete [] gdirect; - delete [] gg; + delete[] gcar; + delete[] gdirect; + delete[] gg; #endif - delete [] ig1; - delete [] ig2; - delete [] ig3; + delete[] ig1; + delete[] ig2; + delete[] ig3; delete[] psi3d; delete[] ylm; delete[] Jlq; - delete[] save; - delete[] posmu; - delete[] posnu; + delete[] save; + delete[] posmu; + delete[] posnu; } void PW_Basis::init(void) { - TITLE("PW_Basis","init"); - timer::tick("PW_Basis","init"); + TITLE("PW_Basis", "init"); + timer::tick("PW_Basis", "init"); // (1) setup the energy cutoff, which is different in each pool, // because the number of k points is differnt. @@ -69,17 +70,17 @@ void PW_Basis::init(void) // (2) setup the FFT dimension according to energy cutoff. this->setup_FFT_dimension(); - OUT(ofs_running,"ggwfc2(Ry)",ggwfc2); + OUT(ofs_running, "ggwfc2(Ry)", ggwfc2); // (3) get the number of plane waves for wave functions. this->ngmc_g = PW_complement::get_total_pw_number(0.0, ggchg, nx, ny, nz, GGT); this->ngmw_g = PW_complement::get_total_pw_number(0.0, ggwfc2, nx, ny, nz, GGT); - OUT(ofs_running,"ngmc_g",ngmc_g); - OUT(ofs_running,"ngmw_g",ngmw_g); + OUT(ofs_running, "ngmc_g", ngmc_g); + OUT(ofs_running, "ngmw_g", ngmw_g); // (4) get the local number of plane waves. - assert(NPROC_IN_POOL>0); + assert(NPROC_IN_POOL > 0); const int remain = ngmw_g % NPROC_IN_POOL; this->ngmw = this->ngmw_g / NPROC_IN_POOL; if (RANK_IN_POOL < remain) @@ -93,23 +94,23 @@ void PW_Basis::init(void) } else { - ngmw_start = (ngmw+1) * remain + ngmw * (RANK_IN_POOL-remain); + ngmw_start = (ngmw + 1) * remain + ngmw * (RANK_IN_POOL - remain); } - - OUT(ofs_running,"ngmw",ngmw); - //ofs_running << "\n ngmw_start = " << ngmw_start; + + OUT(ofs_running, "ngmw", ngmw); + // ofs_running << "\n ngmw_start = " << ngmw_start; // (5) init |g|, g_direct, g_cartesian delete[] gdirect_global; delete[] gcar_global; delete[] gg_global; - gdirect_global = new Vector3[ngmc_g];// indices of G vectors + gdirect_global = new Vector3[ngmc_g]; // indices of G vectors gcar_global = new Vector3[ngmc_g]; - gg_global = new double[ngmc_g];// store the |G|^2 of the 1d array + gg_global = new double[ngmc_g]; // store the |G|^2 of the 1d array PW_complement::get_total_pw(gg_global, gdirect_global, 0.0, ggchg, nx, ny, nz, this->GGT, ngmc_g); - //PW_complement::get_total_pw(gg_global, gdirect_global, 0.0, ggwfc2, nx, ny, nz, this->GGT, ngmw_g); + // PW_complement::get_total_pw(gg_global, gdirect_global, 0.0, ggwfc2, nx, ny, nz, this->GGT, ngmw_g); PW_complement::setup_GVectors(this->G, ngmc_g, gg_global, gdirect_global, gcar_global); // (6) get ig1, ig2, ig3 @@ -118,28 +119,28 @@ void PW_Basis::init(void) #else this->get_GVectors(); #endif - // mohan add on 2011-07-23 - stringstream ss; - ss << "ORBITAL_PW_GK" << MY_RANK+1 << ".dat"; - ofstream ofs(ss.str().c_str()); - ofs << ggwfc2 << " (ggwfc2, Ry)" << endl; - ofs << ngmw << " (Number of plane wave basis for wave functions)" << endl; - for(int ig=0; igset_igk(); this->setup_structure_factor(); - timer::tick("PW_Basis","init"); + timer::tick("PW_Basis", "init"); return; } void PW_Basis::setup_gg(void) { - TITLE("PW_Basis","setup_gg"); + TITLE("PW_Basis", "setup_gg"); //(1.1) lattice parameters this->GT = LATVEC.Inverse(); @@ -148,38 +149,39 @@ void PW_Basis::setup_gg(void) // (1.2) tpiba and ecut assert(LAT0 > 0.0); - this->tpiba = TWO_PI/LAT0; + this->tpiba = TWO_PI / LAT0; this->tpiba2 = tpiba * tpiba; this->ggpsi = ECUT / tpiba2; - double wfac = 4.0; - this->ggchg = wfac * this->ggpsi; + double wfac = 4.0; + this->ggchg = wfac * this->ggpsi; // (1.3) get ggwfc2 - assert( NKS > 0); - assert( CARKX != NULL ); - assert( CARKY != NULL ); - assert( CARKZ != NULL ); + assert(NKS > 0); + assert(CARKX != NULL); + assert(CARKY != NULL); + assert(CARKZ != NULL); - for (int ik=0; ikggpsi) + k_mod; - const double tmp2 = tmp * tmp ; - if (this->ggwfc2 < tmp2) this->ggwfc2 = tmp2; + const double tmp2 = tmp * tmp; + if (this->ggwfc2 < tmp2) + this->ggwfc2 = tmp2; } - OUT(ofs_running,"NKS",NKS); - OUT(ofs_running,"ggpsi(Ry)",ggpsi*tpiba2); - OUT(ofs_running,"ggwfc2(Ry)",ggwfc2*tpiba2); + OUT(ofs_running, "NKS", NKS); + OUT(ofs_running, "ggpsi(Ry)", ggpsi * tpiba2); + OUT(ofs_running, "ggwfc2(Ry)", ggwfc2 * tpiba2); // (1.4) initialize itia2iat this->itia2iat = new int*[NTYPE]; - int iat=0; - for (int it=0; itnpwx = 0; this->ngk = new int[NKS]; ZEROS(this->ngk, NKS); - assert(this->ngmw>0); + assert(this->ngmw > 0); -// ofs_running << "\n\n ngmw = " << ngmw; -// ofs_running << "\n gg[ngmw-1]=" << gg[ngmw-1]; - for (int ik=0; ikigk = new int*[NKS]; - for (int ik=0; ikngk[ik]]; ZEROS(igk[ik], ngk[ik]); } - for (int ik=0; ik kcar(CARKX[ik], CARKY[ik], CARKZ[ik]); - int ng=0; + int ng = 0; const double k2 = kcar * kcar; - for (int ig=0; ig f = this->gcar[ig] + kcar; const double gk2 = f * f; @@ -267,34 +268,32 @@ void PW_Basis::set_igk(void) // Set up crystal structure parameters. void PW_Basis::setup_FFT_dimension() { - TITLE("PW_Basis","setup_FFT_dimension"); + TITLE("PW_Basis", "setup_FFT_dimension"); PW_complement::get_FFT_dimension(LATVEC, this->ggchg, nx, ny, nz); this->nxyz = nx * ny * nz; - OUT(ofs_running,"FFT wavefunctions",nx,ny,nz); + OUT(ofs_running, "FFT wavefunctions", nx, ny, nz); return; } - - #ifdef __MPI void PW_Basis::get_MPI_GVectors(void) { - TITLE("PW_Basis","get_MPI_GVectors"); + TITLE("PW_Basis", "get_MPI_GVectors"); - assert(ngmw>0); + assert(ngmw > 0); // (1) first part delete[] gdirect; delete[] gcar; delete[] gg; this->gdirect = new Vector3[ngmw]; - this->gcar = new Vector3[ngmw]; + this->gcar = new Vector3[ngmw]; this->gg = new double[ngmw]; - for (int ig=0; igig2 = new int[ngmw]; this->ig3 = new int[ngmw]; -// cout << "\n dim of ig1,2,3 = " << ngmw; + // cout << "\n dim of ig1,2,3 = " << ngmw; for (int i = 0; i < ngmw; i++) { @@ -323,12 +322,12 @@ void PW_Basis::get_MPI_GVectors(void) } return; -}//end setup_mpi_GVectors +} // end setup_mpi_GVectors #else void PW_Basis::get_GVectors(void) { - TITLE("PW_Basis","get_GVectors"); - timer::tick("PW_Basis","get_GVectors"); + TITLE("PW_Basis", "get_GVectors"); + timer::tick("PW_Basis", "get_GVectors"); //************************************************************ // g : Store the G vectors in 1d array (Cartian coordinate) @@ -345,7 +344,7 @@ void PW_Basis::get_GVectors(void) this->gg = this->gg_global; // (2) calculate ig1, ig2, ig3 - assert(ngmw>0); + assert(ngmw > 0); delete[] ig1; delete[] ig2; delete[] ig3; @@ -360,34 +359,33 @@ void PW_Basis::get_GVectors(void) this->ig3[i] = int(this->gdirect[i].z) + nz; } - timer::tick("PW_Basis","get_GVectors"); + timer::tick("PW_Basis", "get_GVectors"); return; -}//end get_GVectors; +} // end get_GVectors; #endif - // Calculate structure factor void PW_Basis::setup_structure_factor(void) { - TITLE("PW_Basis","setup_structure_factor"); - timer::tick("PW_Basis","setup_struc_factor"); + TITLE("PW_Basis", "setup_structure_factor"); + timer::tick("PW_Basis", "setup_struc_factor"); std::complex ci_tpi = NEG_IMAG_UNIT * TWO_PI; std::complex x; this->strucFac = new std::complex*[NTYPE]; - for (int it=0; itstrucFac[it] = new std::complex[ngmw]; - ZEROS( strucFac[it], ngmw); + ZEROS(strucFac[it], ngmw); } - for (int it=0; it< NTYPE; it++) + for (int it = 0; it < NTYPE; it++) { - for (int ig=0; igstrucFac[it][ig] = std::complex( sum_cos, -sum_sin ); + this->strucFac[it][ig] = std::complex(sum_cos, -sum_sin); - double tmpx = strucFac[it][ig].real() ; - double tmpy = strucFac[it][ig].imag() ; + double tmpx = strucFac[it][ig].real(); + double tmpy = strucFac[it][ig].imag(); } } int nat = 0; - for (int it=0; iteigts1.create(nat, 2*this->nx + 1); - this->eigts2.create(nat, 2*this->ny + 1); - this->eigts3.create(nat, 2*this->nz + 1); + int i, j, ng; + this->eigts1.create(nat, 2 * this->nx + 1); + this->eigts2.create(nat, 2 * this->ny + 1); + this->eigts3.create(nat, 2 * this->nz + 1); Vector3 gtau; Vector3 tau; @@ -431,78 +428,75 @@ void PW_Basis::setup_structure_factor(void) tau.y = CARPOSY[it][ia]; tau.z = CARPOSZ[it][ia]; - gtau = this->G * tau; //HLX: fixed on 10/13/2006 + gtau = this->G * tau; // HLX: fixed on 10/13/2006 for (int n1 = -nx; n1 <= nx; n1++) { double arg = n1 * gtau.x; - this->eigts1(iat, n1 + nx) = exp( ci_tpi*arg ); + this->eigts1(iat, n1 + nx) = exp(ci_tpi * arg); } for (int n2 = -ny; n2 <= ny; n2++) { double arg = n2 * gtau.y; - this->eigts2(iat, n2 + ny) = exp( ci_tpi*arg ); + this->eigts2(iat, n2 + ny) = exp(ci_tpi * arg); } for (int n3 = -nz; n3 <= nz; n3++) { double arg = n3 * gtau.z; - this->eigts3(iat, n3 + nz) = exp( ci_tpi*arg ); + this->eigts3(iat, n3 + nz) = exp(ci_tpi * arg); } iat++; } } - timer::tick("PW_Basis","setup_struc_factor"); + timer::tick("PW_Basis", "setup_struc_factor"); return; } -complex* PW_Basis::get_sk(const int ik, const int it, const int ia)const +complex* PW_Basis::get_sk(const int ik, const int it, const int ia) const { - timer::tick("PW_Basis","get_sk"); - const double arg = (CARKX[ik] * CARPOSX[it][ia] - + CARKY[ik] * CARPOSY[it][ia] - + CARKZ[ik] * CARPOSZ[it][ia] ) * TWO_PI; - const std::complex kphase = complex ( cos(arg), -sin(arg) ); - std::complex *sk = new std::complex[ this->ngk[ik] ]; + timer::tick("PW_Basis", "get_sk"); + const double arg + = (CARKX[ik] * CARPOSX[it][ia] + CARKY[ik] * CARPOSY[it][ia] + CARKZ[ik] * CARPOSZ[it][ia]) * TWO_PI; + const std::complex kphase = complex(cos(arg), -sin(arg)); + std::complex* sk = new std::complex[this->ngk[ik]]; const int iat = this->itia2iat[it][ia]; - for (int ig=0; ig< this->ngk[ik]; ig++) + for (int ig = 0; ig < this->ngk[ik]; ig++) { -// cout << "\n ig=" << ig << endl; + // cout << "\n ig=" << ig << endl; const int iig = this->igk[ik][ig]; - assert( iig < ngmw ); + assert(iig < ngmw); - if (!(ig1[iig] >= 0 && ig1[iig] < 2*nx + 1) ) + if (!(ig1[iig] >= 0 && ig1[iig] < 2 * nx + 1)) { cout << "\n ig1[iig] = " << ig1[iig] << endl; - QUIT(); + QUIT(); } - if (!(ig2[iig] >= 0 && ig2[iig] < 2*ny + 1) ) + if (!(ig2[iig] >= 0 && ig2[iig] < 2 * ny + 1)) { cout << "\n ig2[iig] = " << ig2[iig] << endl; - QUIT(); + QUIT(); } - if (!(ig3[iig] >= 0 && ig3[iig] < 2*nz + 1) ) + if (!(ig3[iig] >= 0 && ig3[iig] < 2 * nz + 1)) { cout << "\n iig = " << iig << endl; cout << "\n ig3[iig] = " << ig3[iig] << endl; - QUIT(); + QUIT(); } - sk[ig] = kphase - * this->eigts1(iat, this->ig1[iig]) - * this->eigts2(iat, this->ig2[iig]) + sk[ig] = kphase * this->eigts1(iat, this->ig1[iig]) * this->eigts2(iat, this->ig2[iig]) * this->eigts3(iat, this->ig3[iig]); } - timer::tick("PW_Basis","get_sk"); + timer::tick("PW_Basis", "get_sk"); return sk; } -Vector3 PW_Basis::get_1qvec_cartesian(const int ik,const int ig)const +Vector3 PW_Basis::get_1qvec_cartesian(const int ik, const int ig) const { Vector3 kvec = Vector3(CARKX[ik], CARKY[ik], CARKZ[ik]); - Vector3 qvec = kvec + this->gcar[ this->igk[ik][ig] ]; + Vector3 qvec = kvec + this->gcar[this->igk[ik][ig]]; return qvec; } @@ -511,19 +505,19 @@ void PW_Basis::table(void) TITLE("PW_Basis", "table"); NBasis.init_table(); - const int total_lm = ( LMAXUSED + 1) * ( LMAXUSED + 1); + const int total_lm = (LMAXUSED + 1) * (LMAXUSED + 1); delete[] ylm; this->ylm = new matrix[NKS]; - for (int ik=0; ikylm[ik].create(total_lm, npw); // (1) get the vector k+G - Vector3 *gk = new Vector3 [npw]; - for (int ig=0; ig* gk = new Vector3[npw]; + for (int ig = 0; ig < npw; ig++) { gk[ig] = this->get_1qvec_cartesian(ik, ig); } @@ -535,37 +529,35 @@ void PW_Basis::table(void) delete[] gk; } + // init index of iw00 + assert(NWFCALL > 0); + this->iwindex = new Way2iw[NWFCALL]; - // init index of iw00 - assert(NWFCALL>0); - this->iwindex = new Way2iw[NWFCALL]; - - int iw=0; - for(int it=0; itDk = 0.01; this->kmesh = static_cast(sqrt(ECUT) / Dk) + 1 + 4; - if (kmesh % 2 == 0) ++kmesh; - - this->psi1d.create(NTYPE, LMAXUSED+1, NMAXUSED, kmesh); - - delete[] save; - delete[] posmu; - delete[] posnu; - - const int dim = nwfc2 * nwfc2; - save = new std::complex[dim]; - posmu = new int[dim]; - posnu = new int[dim]; - ZEROS(save, dim); - ZEROS(posmu, dim); - ZEROS(posnu, dim); + if (kmesh % 2 == 0) + ++kmesh; + + this->psi1d.create(NTYPE, LMAXUSED + 1, NMAXUSED, kmesh); + + delete[] save; + delete[] posmu; + delete[] posnu; + + const int dim = nwfc2 * nwfc2; + save = new std::complex[dim]; + posmu = new int[dim]; + posnu = new int[dim]; + ZEROS(save, dim); + ZEROS(posmu, dim); + ZEROS(posnu, dim); return; } void PW_Basis::calculate_psi1d(void) { - timer::tick("PW_Basis","calculate_psi1d"); + timer::tick("PW_Basis", "calculate_psi1d"); psi1d.zero_out(); - for (int it=0; itpsi1d(it,l,n,ikm) += input.Coef.C4(it, l, n ,ie) - * Numerical_Basis::bessel_basis.TableOne(l, ie, ikm); + this->psi1d(it, l, n, ikm) + += input.Coef.C4(it, l, n, ie) * Numerical_Basis::bessel_basis.TableOne(l, ie, ikm); } } } } } - timer::tick("PW_Basis","calculate_psi1d"); + timer::tick("PW_Basis", "calculate_psi1d"); return; } -void PW_Basis::allocate_psi3d(const int &level) +void PW_Basis::allocate_psi3d(const int& level) { // allocate the number of wave functions. @@ -624,17 +617,16 @@ void PW_Basis::allocate_psi3d(const int &level) delete[] psi3d; psi3d = new ComplexMatrix[NKS]; - for (int ik=0; ikpsi3d[ik].create(nwfc2, npwx); } return; } - -void PW_Basis::calculate_psi3d(const int &ilevel, const int &ik) +void PW_Basis::calculate_psi3d(const int& ilevel, const int& ik) { - //TITLE("PW_Basis", "calculate_psi3d"); + // TITLE("PW_Basis", "calculate_psi3d"); timer::tick("PW_Basis", "calculate_psi3d"); // (1) get the vector k+G @@ -643,36 +635,36 @@ void PW_Basis::calculate_psi3d(const int &ilevel, const int &ik) const int npw = ngk[ik]; - Vector3 *gk = new Vector3 [npw]; - for (int ig=0; ig* gk = new Vector3[npw]; + for (int ig = 0; ig < npw; ig++) { gk[ig] = this->get_1qvec_cartesian(ik, ig); } - double *flq = new double[npw]; - for (int iw=0; iwnwfc2; iw++) + double* flq = new double[npw]; + for (int iw = 0; iw < this->nwfc2; iw++) { // for each wave function in nwfc2(it, i, l, n) const int it = mz.Level[ilevel].wayd[iw].type; const int ia = mz.Level[ilevel].wayd[iw].i; - const int l = mz.Level[ilevel].wayd[iw].L; - const int n = mz.Level[ilevel].wayd[iw].N; - const int m = mz.Level[ilevel].wayd[iw].m; + const int l = mz.Level[ilevel].wayd[iw].L; + const int n = mz.Level[ilevel].wayd[iw].N; + const int m = mz.Level[ilevel].wayd[iw].m; std::complex lphase = pow(IMAG_UNIT, l); // get flq from psi1d for each k point - for (int ig=0; igPolynomial_Interpolation(it, l, n, gk[ig].norm() * this->tpiba ); + flq[ig] = this->Polynomial_Interpolation(it, l, n, gk[ig].norm() * this->tpiba); } // get the structure wave functions for each k point. - std::complex *sk = this->get_sk(ik, it, ia); + std::complex* sk = this->get_sk(ik, it, ia); - const int lm = l*l+m; - for (int ig=0; ig PW_Basis::calculateS(const int &iw, const int &iw2, const int &ik) +complex PW_Basis::calculateS(const int& iw, const int& iw2, const int& ik) { std::complex overlap = 0.0; - for (int ig=0; igpsi3d[ik](iw, ig)) * this->psi3d[ik](iw2, ig); } @@ -698,11 +690,9 @@ complex PW_Basis::calculateS(const int &iw, const int &iw2, const int &i return overlap; } - -void PW_Basis::update_psi1d(const int &il, const int &ic, const int &ie, - const double &c4_now, const double &c4_old) +void PW_Basis::update_psi1d(const int& il, const int& ic, const int& ie, const double& c4_now, const double& c4_old) { - timer::tick("PW_Basis","update_psi1d"); + timer::tick("PW_Basis", "update_psi1d"); // a small trick // ALGORITHM : when we change psi1d, we don't know @@ -718,34 +708,34 @@ void PW_Basis::update_psi1d(const int &il, const int &ic, const int &ie, const int l_now = mz.Level[il].wayc[ic].L; const int n_now = mz.Level[il].wayc[ic].N; - for (int it=0; itpsi1d(it,l,n,ikm) = sum; + this->psi1d(it, l, n, ikm) = sum; + } + } + else if (it == it_now || l == l_now || n == n_now) + { + for (int ikm = 0; ikm < kmesh; ikm++) + { + // a small trick + this->psi1d(it, l, n, ikm) + += (c4_now - c4_old) * Numerical_Basis::bessel_basis.TableOne(l, ie, ikm); } } - else if(it==it_now || l==l_now || n==n_now) - { - for (int ikm=0; ikmpsi1d(it,l,n,ikm) += (c4_now-c4_old) * Numerical_Basis::bessel_basis.TableOne(l, ie, ikm); - } - } } } } @@ -754,27 +744,28 @@ void PW_Basis::update_psi1d(const int &il, const int &ic, const int &ie, l_last = l_now; n_last = n_now; - timer::tick("PW_Basis","update_psi1d"); + timer::tick("PW_Basis", "update_psi1d"); return; } -void PW_Basis::update_psi3d( const int &il, const int &ic, const int &ik) +void PW_Basis::update_psi3d(const int& il, const int& ic, const int& ik) { -// TITLE("PW_Basis","update_psi1d3d"); - timer::tick("PW_Basis","update_psi3d"); + // TITLE("PW_Basis","update_psi1d3d"); + timer::tick("PW_Basis", "update_psi3d"); const int npw = ngk[ik]; static int it_old; static int ia_old; static int l_old; static int n_old; - static bool allocate=false; - static int *ic_last; + static bool allocate = false; + static int* ic_last; if (!allocate) { ic_last = new int[NKS]; - for (int ik=0; ik *gk = new Vector3 [npw]; - for (int ig=0; ig* gk = new Vector3[npw]; + for (int ig = 0; ig < npw; ig++) { gk[ig] = this->get_1qvec_cartesian(ik, ig); } - double *flq = new double[npw]; - std::complex *sk = new std::complex[1]; - std::complex *samepart = new std::complex[npw]; - for (int iw=0; iwnwfc2; iw++) + double* flq = new double[npw]; + std::complex* sk = new std::complex[1]; + std::complex* samepart = new std::complex[npw]; + for (int iw = 0; iw < this->nwfc2; iw++) { const int ic1 = mz.Level[il].wayd[iw].ic; if (ic1 == ic || ic1 == ic_last[ik]) @@ -800,13 +790,13 @@ void PW_Basis::update_psi3d( const int &il, const int &ic, const int &ik) // for each wave function in nwfc2(it, i, l, n) const int it = mz.Level[il].wayd[iw].type; const int ia = mz.Level[il].wayd[iw].i; - const int l = mz.Level[il].wayd[iw].L; - const int n = mz.Level[il].wayd[iw].N; - const int m = mz.Level[il].wayd[iw].m; + const int l = mz.Level[il].wayd[iw].L; + const int n = mz.Level[il].wayd[iw].N; + const int m = mz.Level[il].wayd[iw].m; std::complex lphase = pow(IMAG_UNIT, l); // a small trick: get the structure wave functions for each k point. - if (it==it_old && ia==ia_old) + if (it == it_old && ia == ia_old) { // do nothing } @@ -817,47 +807,46 @@ void PW_Basis::update_psi3d( const int &il, const int &ic, const int &ik) } // small trick (2) : get flq from psi1d for each k point - if (it==it_old && l==l_old && n==n_old) + if (it == it_old && l == l_old && n == n_old) { // do nothing } else { - for (int ig=0; igPolynomial_Interpolation(it, l, n, gk[ig].norm() * this->tpiba ); + flq[ig] = this->Polynomial_Interpolation(it, l, n, gk[ig].norm() * this->tpiba); } } // small trick (3) - if (it==it_old && ia==ia_old && l==l_old && n==n_old) + if (it == it_old && ia == ia_old && l == l_old && n == n_old) { - //samepart is the same as the previous one. + // samepart is the same as the previous one. } else { - for (int ig=0; igDk; + const double position = gnorm / this->Dk; const int iq = static_cast(position); - assert(iq < kmesh-4); + assert(iq < kmesh - 4); const double x0 = position - static_cast(iq); const double x1 = 1.0 - x0; const double x2 = 2.0 - x0; const double x3 = 3.0 - x0; - const double y= - this->psi1d(it, l, n, iq) * x1 * x2 * x3 / 6.0 + - this->psi1d(it, l, n, iq) * x0 * x2 * x3 / 2.0 - - this->psi1d(it, l, n, iq) * x1 * x0 * x3 / 2.0 + - this->psi1d(it, l, n, iq) * x1 * x2 * x0 / 6.0 ; + const double y = this->psi1d(it, l, n, iq) * x1 * x2 * x3 / 6.0 + this->psi1d(it, l, n, iq) * x0 * x2 * x3 / 2.0 + - this->psi1d(it, l, n, iq) * x1 * x0 * x3 / 2.0 + this->psi1d(it, l, n, iq) * x1 * x2 * x0 / 6.0; return y; } - -void PW_Basis::calculate_Jlq(const int &ik, const int &iw, const int &ie) +void PW_Basis::calculate_Jlq(const int& ik, const int& iw, const int& ie) { - //TITLE("PW_Basis", "calculate_Jlq"); - timer::tick("PW_Basis","calculate_Jlq"); + // TITLE("PW_Basis", "calculate_Jlq"); + timer::tick("PW_Basis", "calculate_Jlq"); const int npw = this->ngk[ik]; - const int it= iwindex[iw].type; - const int ia= iwindex[iw].i; - const int l = iwindex[iw].L; - const int m = iwindex[iw].m; + const int it = iwindex[iw].type; + const int ia = iwindex[iw].i; + const int l = iwindex[iw].L; + const int m = iwindex[iw].m; // generate 1d Jlq. - Vector3 *gk = new Vector3 [npw]; - double* flq=new double[npw]; - for (int ig=0; ig* gk = new Vector3[npw]; + double* flq = new double[npw]; + for (int ig = 0; ig < npw; ig++) { gk[ig] = this->get_1qvec_cartesian(ik, ig); } - for (int ig=0; igtpiba ); + flq[ig] = Numerical_Basis::bessel_basis.Polynomial_Interpolation2(l, ie, gk[ig].norm() * this->tpiba); } // generate 3d Jlq. delete[] Jlq; Jlq = new std::complex[npw]; - std::complex *sk = this->get_sk(ik, it, ia); + std::complex* sk = this->get_sk(ik, it, ia); - const int lm = l*l+m; + const int lm = l * l + m; std::complex lphase = pow(IMAG_UNIT, l); - for (int ig=0; ig PW_Basis::calculate_Jlq_Phi(const int &ik, const int &mu) +complex PW_Basis::calculate_Jlq_Phi(const int& ik, const int& mu) { - const int npw = this->ngk[ik]; - complex overlap = std::complex(0,0); - for(int ig=0; igngk[ik]; + complex overlap = std::complex(0, 0); + for (int ig = 0; ig < npw; ig++) + { + overlap += conj(Jlq[ig]) * psi3d[ik](mu, ig); + } #ifdef __MPI Parallel_Reduce::reduce_complex_double_pool(overlap); #endif - return overlap; + return overlap; } diff --git a/tools/01_NAO_generation/SIAB/src_pw/pw_basis.h b/tools/01_NAO_generation/SIAB/src_pw/pw_basis.h index 21e7b57f8ce..6bf7b1af3bc 100644 --- a/tools/01_NAO_generation/SIAB/src_pw/pw_basis.h +++ b/tools/01_NAO_generation/SIAB/src_pw/pw_basis.h @@ -11,125 +11,122 @@ struct Way2iw // dimension : WFCALL { int type; - int i; // atom + int i; // atom int L; - int m; + int m; }; class PW_Basis { -public: - - // about lattice. - Matrix3 G; - Matrix3 GT; - Matrix3 GGT; - double tpiba; // 2pi/lat0 - double tpiba2; // tpiba^2 - - // about k point - int* ngk; - int** igk; - int** itia2iat; - int npwx; - - // about fft grid. - int nx; - int ny; - int nz; - int nxyz; - - double ggpsi; // planewave cut off for the wavefunctions, unit (2*PI/lat0)^2 - double ggwfc2; // ggwav=wfact*ggpsi, default value: wfact=4.0 - double ggchg; - - std::complex** strucFac; // StrucFac [ntype,ngmax] - - int ngmw; //(= ngmax) / num. of G vectors within ggfft - int ngmw_g; - int ngmw_start; - - int ngmc; - int ngmc_g; - int ngmc_start; - - Vector3 *gdirect; //(= *G1d) ; // ig = new Vector igc[ngmc], - Vector3 *gdirect_global; //(= *G1d) ; // ig = new Vector igc[ngmc], - Vector3 *gcar; //G vectors in cartesian corrdinate - Vector3 *gcar_global; //G vectors in cartesian corrdinate - double *gg; // modulus (G^2) of G vectors - double *gg_global; // modulus (G^2) of G vectors - - ComplexMatrix eigts1; // - ComplexMatrix eigts2; //the phases e^{-iG*tau_s} - ComplexMatrix eigts3; // - int *ig1; // - int *ig2; // the indices of G components - int *ig3; // + public: + // about lattice. + Matrix3 G; + Matrix3 GT; + Matrix3 GGT; + double tpiba; // 2pi/lat0 + double tpiba2; // tpiba^2 + + // about k point + int* ngk; + int** igk; + int** itia2iat; + int npwx; + + // about fft grid. + int nx; + int ny; + int nz; + int nxyz; + + double ggpsi; // planewave cut off for the wavefunctions, unit (2*PI/lat0)^2 + double ggwfc2; // ggwav=wfact*ggpsi, default value: wfact=4.0 + double ggchg; + + std::complex** strucFac; // StrucFac [ntype,ngmax] + + int ngmw; //(= ngmax) / num. of G vectors within ggfft + int ngmw_g; + int ngmw_start; + + int ngmc; + int ngmc_g; + int ngmc_start; + + Vector3* gdirect; //(= *G1d) ; // ig = new Vector igc[ngmc], + Vector3* gdirect_global; //(= *G1d) ; // ig = new Vector igc[ngmc], + Vector3* gcar; // G vectors in cartesian corrdinate + Vector3* gcar_global; // G vectors in cartesian corrdinate + double* gg; // modulus (G^2) of G vectors + double* gg_global; // modulus (G^2) of G vectors + + ComplexMatrix eigts1; // + ComplexMatrix eigts2; // the phases e^{-iG*tau_s} + ComplexMatrix eigts3; // + int* ig1; // + int* ig2; // the indices of G components + int* ig3; // PW_Basis(); ~PW_Basis(); void init(void); - void setup_structure_factor(void); // Calculate structur factors - void get_sk(void); - complex* get_sk(const int ik, const int it, const int ia)const; - Vector3 get_1qvec_cartesian(const int ik,const int ig)const; - - void table(void); - Numerical_Basis NBasis; - - void allocate_psi1d(const int &il); - void allocate_psi3d(const int &level); - realArray psi1d; - ComplexMatrix *psi3d; - matrix *ylm; - - void calculate_psi1d(void); - // FUNCTION: - // use psi1d to generate 3D wave functions. - // psi3d(it, ia, l, n, m) for each k point. - void calculate_psi3d(const int &ilevel, const int &ik); - - void update_psi1d(const int &il, const int &ic, const int &ie, - const double &c4_now, const double &c4_old); - - int nwfc2; - - void update_psi3d( const int &il, const int &ic, const int &ik); - - int kmesh; - double Dk; - - // FNCTION: - // calculate the overlap between psi(it, ia, l, n, m) - // be called in SpillageStep::init_QS_matrix - complex calculateS(const int &iw, const int &iw2, const int &ik); - - void calculate_Jlq(const int &ik, const int &iw, const int &ie); - complex *Jlq; - complex calculate_Jlq_Phi(const int &ik, const int &mu); - - complex *save; - int* posmu; - int* posnu; - -private: - void set_igk(void); + void setup_structure_factor(void); // Calculate structur factors + void get_sk(void); + complex* get_sk(const int ik, const int it, const int ia) const; + Vector3 get_1qvec_cartesian(const int ik, const int ig) const; + + void table(void); + Numerical_Basis NBasis; + + void allocate_psi1d(const int& il); + void allocate_psi3d(const int& level); + realArray psi1d; + ComplexMatrix* psi3d; + matrix* ylm; + + void calculate_psi1d(void); + // FUNCTION: + // use psi1d to generate 3D wave functions. + // psi3d(it, ia, l, n, m) for each k point. + void calculate_psi3d(const int& ilevel, const int& ik); + + void update_psi1d(const int& il, const int& ic, const int& ie, const double& c4_now, const double& c4_old); + + int nwfc2; + + void update_psi3d(const int& il, const int& ic, const int& ik); + + int kmesh; + double Dk; + + // FNCTION: + // calculate the overlap between psi(it, ia, l, n, m) + // be called in SpillageStep::init_QS_matrix + complex calculateS(const int& iw, const int& iw2, const int& ik); + + void calculate_Jlq(const int& ik, const int& iw, const int& ie); + complex* Jlq; + complex calculate_Jlq_Phi(const int& ik, const int& mu); + + complex* save; + int* posmu; + int* posnu; + + private: + void set_igk(void); void setup_gg(void); - void setup_FFT_dimension(void); // set up FFT dimensions + void setup_FFT_dimension(void); // set up FFT dimensions #ifdef __MPI void get_MPI_GVectors(void); #else void get_GVectors(void); #endif - // FUNCTION: - // use interpolation scheme to get the one dimension wave function for - // each k point - double Polynomial_Interpolation(const int &it, const int &l, const int &n, const double &gnorm)const; + // FUNCTION: + // use interpolation scheme to get the one dimension wave function for + // each k point + double Polynomial_Interpolation(const int& it, const int& l, const int& n, const double& gnorm) const; - Way2iw* iwindex; - + Way2iw* iwindex; }; -#endif //PlaneWave class +#endif // PlaneWave class diff --git a/tools/01_NAO_generation/SIAB/src_pw/pw_complement.cpp b/tools/01_NAO_generation/SIAB/src_pw/pw_complement.cpp index a2ab75cadad..8e12b6a15a9 100644 --- a/tools/01_NAO_generation/SIAB/src_pw/pw_complement.cpp +++ b/tools/01_NAO_generation/SIAB/src_pw/pw_complement.cpp @@ -1,21 +1,21 @@ #include "pw_complement.h" + #include "heapsort.h" -int PW_complement::get_total_pw_number( - const double& ggcut_start, - const double& ggcut_end, - const int& nx, - const int& ny, - const int& nz, - const Matrix3& GGT) +int PW_complement::get_total_pw_number(const double& ggcut_start, + const double& ggcut_end, + const int& nx, + const int& ny, + const int& nz, + const Matrix3& GGT) { - TITLE("PW_complement","get_total_pw_number"); - if (ggcut_end<=0.0) + TITLE("PW_complement", "get_total_pw_number"); + if (ggcut_end <= 0.0) { - WARNING_QUIT("PW_complement::get_total_pw_number","ggcut <= 0.0"); + WARNING_QUIT("PW_complement::get_total_pw_number", "ggcut <= 0.0"); } // First, figure out the number of G-vectors within the cutoff G2max. - int ibox[3]={0,0,0}; + int ibox[3] = {0, 0, 0}; // mohan modify 2008-3-25 nx -> ncx ,ny-> ncy, nz->ncz ibox[0] = int(nx / 2) + 1; @@ -24,7 +24,7 @@ int PW_complement::get_total_pw_number( // first create the auxiliary arrays for the 1D G vectors Vector3 f; - int ngm = 0 ; + int ngm = 0; for (int i = -ibox[0]; i <= ibox[0]; i++) { for (int j = -ibox[1]; j <= ibox[1]; j++) @@ -34,9 +34,9 @@ int PW_complement::get_total_pw_number( f.x = i; f.y = j; f.z = k; - double g2 = f * (GGT * f); //G2= |f|^2 in the unit of (2Pi/lat0)^2 - - if (g2 < ggcut_end && g2 >= ggcut_start) + double g2 = f * (GGT * f); // G2= |f|^2 in the unit of (2Pi/lat0)^2 + + if (g2 < ggcut_end && g2 >= ggcut_start) { ngm++; } @@ -46,26 +46,25 @@ int PW_complement::get_total_pw_number( return ngm; } -void PW_complement::get_total_pw( - double* gg, - Vector3 *ig, - const double& ggcut_start, - const double& ggcut_end, - const int& nx, - const int& ny, - const int& nz, - const Matrix3& GGT, // GGT = G*GT. - int& ngm// number of total plane waves. +void PW_complement::get_total_pw(double* gg, + Vector3* ig, + const double& ggcut_start, + const double& ggcut_end, + const int& nx, + const int& ny, + const int& nz, + const Matrix3& GGT, // GGT = G*GT. + int& ngm // number of total plane waves. ) { - TITLE("PW_complement","get_total_pw"); - timer::tick("PW_complement","get_total_pw"); - if (ggcut_end<=0.0) + TITLE("PW_complement", "get_total_pw"); + timer::tick("PW_complement", "get_total_pw"); + if (ggcut_end <= 0.0) { - WARNING_QUIT("PW_complement::get_total_pw","ggcut <= 0.0"); + WARNING_QUIT("PW_complement::get_total_pw", "ggcut <= 0.0"); } // First, figure out the number of G-vectors within the cutoff G2max. - int ibox[3]={0,0,0}; + int ibox[3] = {0, 0, 0}; // mohan modify 2008-3-25 nx -> ncx ,ny-> ncy, nz->ncz ibox[0] = int(nx / 2) + 1; @@ -84,39 +83,34 @@ void PW_complement::get_total_pw( f.x = i; f.y = j; f.z = k; - double g2 = f * (GGT * f); //G2= |f|^2 in the unit of (2Pi/lat0)^2 + double g2 = f * (GGT * f); // G2= |f|^2 in the unit of (2Pi/lat0)^2 if (g2 < ggcut_end && g2 >= ggcut_start) { /*** g vector indices f=(i,j,k) ***/ - ig[ng] = f ; + ig[ng] = f; gg[ng] = g2; - //cout<0,use the input // FFT grid, otherwise generate the FFT grid in the code. int i = 0; Vector3 lat; - int ibox[3]={0,0,0}; + int ibox[3] = {0, 0, 0}; // ibox[i] are the minimal FFT dimensions, lat.x = latvec.e11; @@ -133,12 +127,12 @@ void PW_complement::get_FFT_dimension( lat.y = latvec.e32; lat.z = latvec.e33; ibox[2] = 2 * int(sqrt(ggcut) * sqrt(lat * lat)) + 1; - //lat*lat=lat.x*lat.x+lat.y*lat.y+lat.z+lat.z + // lat*lat=lat.x*lat.x+lat.y*lat.y+lat.z+lat.z /* Find the minimal FFT box size the factors into the primes (2,3,5,7). */ -// log << "\n ibox[0]= "<< ibox[0] << " ibox[1] = "< *igsort = new Vector3[ngmc_g_in]; - for (int i=0;i* igsort = new Vector3[ngmc_g_in]; + for (int i = 0; i < ngmc_g_in; i++) { igsort[i] = ig[ind[i]]; -// cout << i << " " << ind[i] << " " << ig[ind[i]].x << " " << ig[ind[i]].y << " " << ig[ind[i]].z << endl; + // cout << i << " " << ind[i] << " " << ig[ind[i]].x << " " << ig[ind[i]].y << " " << ig[ind[i]].z << endl; } - for (int i=0; i *igsort, - const double& ggcut_start, const double& ggcut_end, - const int& nx, const int& ny, const int& nz, const Matrix3& GGT, // GGT = G*GT - int& ngm); // number of total plane waves. +int get_total_pw_number(const double& ggcut_start, + const double& ggcut_end, + const int& nx, + const int& ny, + const int& nz, + const Matrix3& GGT); + +void get_total_pw(double* ig2sort, + Vector3* igsort, + const double& ggcut_start, + const double& ggcut_end, + const int& nx, + const int& ny, + const int& nz, + const Matrix3& GGT, // GGT = G*GT + int& ngm); // number of total plane waves. // output nx, ny, nz according to input: latve and ggcut. -void get_FFT_dimension(const Matrix3& latvec, const double &ggcut, int &nx, int &ny, int &nz); +void get_FFT_dimension(const Matrix3& latvec, const double& ggcut, int& nx, int& ny, int& nz); //========================================================== // MEMBER FUNCTION : @@ -30,10 +39,7 @@ void get_FFT_dimension(const Matrix3& latvec, const double &ggcut, int &nx, int // the 3D G-grid and the FFT grid. // generate : gg_global, g_global, ig_global //========================================================== -void setup_GVectors(const Matrix3& G, const int &ngmc_g, double* gg, - Vector3* ig, Vector3* g); -} - +void setup_GVectors(const Matrix3& G, const int& ngmc_g, double* gg, Vector3* ig, Vector3* g); +} // namespace PW_complement #endif - diff --git a/tools/01_NAO_generation/SIAB/src_pw/ylm_real.cpp b/tools/01_NAO_generation/SIAB/src_pw/ylm_real.cpp index 5fcc7724f4a..2c57ba77bf5 100644 --- a/tools/01_NAO_generation/SIAB/src_pw/ylm_real.cpp +++ b/tools/01_NAO_generation/SIAB/src_pw/ylm_real.cpp @@ -1,4 +1,5 @@ #include "ylm_real.h" + #include "../src_parallel/parallel_reduce.h" #include "../src_spillage/tools.h" @@ -7,33 +8,31 @@ // NAME : YLM_REAL(Real spherical harmonics ylm(G) up to l=lmax // Use Numerical recursive algorithm as given in Numerical Recipes //========================================================== -void Ylm_Real -( - const int lmax2, // lmax2 = (lmax+1)^2 - const int ng, // - const Vector3 *g, // g_cartesian_vec(x,y,z) - matrix &ylm // output +void Ylm_Real(const int lmax2, // lmax2 = (lmax+1)^2 + const int ng, // + const Vector3* g, // g_cartesian_vec(x,y,z) + matrix& ylm // output ) { -// TITLE("Mathzone","Ylm_Real"); -// timer::tick("Mathzone","Ylm_Real"); + // TITLE("Mathzone","Ylm_Real"); + // timer::tick("Mathzone","Ylm_Real"); - if (ng<1 || lmax2<1) + if (ng < 1 || lmax2 < 1) { - cout << "\n YLM_REAL: ng < 1 or lmax2 < 1"; - QUIT(); - timer::tick("Mathzone","Ylm_Real"); + cout << "\n YLM_REAL: ng < 1 or lmax2 < 1"; + QUIT(); + timer::tick("Mathzone", "Ylm_Real"); return; } -//---------------------------------------------------------- -// EXPLAIN : find out lmax -//---------------------------------------------------------- + //---------------------------------------------------------- + // EXPLAIN : find out lmax + //---------------------------------------------------------- bool out_of_range = true; int lmax = 0; - for (int l= 0; l< 30; l++) + for (int l = 0; l < 30; l++) { - if ((l+1)*(l+1) == lmax2) + if ((l + 1) * (l + 1) == lmax2) { lmax = l; out_of_range = false; @@ -42,15 +41,15 @@ void Ylm_Real } if (out_of_range) { - WARNING_QUIT("YLM_REAL","l>30 or l<0"); + WARNING_QUIT("YLM_REAL", "l>30 or l<0"); } -//---------------------------------------------------------- -// EXPLAIN : if lmax = 1,only use Y00 , output result. -//---------------------------------------------------------- + //---------------------------------------------------------- + // EXPLAIN : if lmax = 1,only use Y00 , output result. + //---------------------------------------------------------- if (lmax == 0) { - for (int i=0;i 1.0e-9) @@ -89,59 +88,58 @@ void Ylm_Real } else { - phi[ig] = PI_HALF * ((g[ig].y >= 0.0) ? 1.0 : -1.0); //HLX: modified on 10/13/2006 + phi[ig] = PI_HALF * ((g[ig].y >= 0.0) ? 1.0 : -1.0); // HLX: modified on 10/13/2006 } // end if } // enddo -//========================================================== -// NAME : p(Legendre Polynomials) (0 <= m <= l) -//========================================================== - realArray p(lmax+1,lmax+1,ng); + //========================================================== + // NAME : p(Legendre Polynomials) (0 <= m <= l) + //========================================================== + realArray p(lmax + 1, lmax + 1, ng); int m; int i; double x1, x2; int lm = -1; - for (int l=0; l<=lmax; l++) + for (int l = 0; l <= lmax; l++) { - const double c = sqrt((2*l+1) / FOUR_PI); + const double c = sqrt((2 * l + 1) / FOUR_PI); if (l == 0) { - for (i=0;i(l) / 2.0) ;//mohan modify 2007-10-13 - if (l%2 == 1) + p(l, l, i) = Semi_Fact(l3) * pow(x2, static_cast(l) / 2.0); // mohan modify 2007-10-13 + if (l % 2 == 1) { p(l, l, i) = -p(l, l, i); } @@ -150,90 +148,81 @@ void Ylm_Real // Y_lm, m = 0 ++lm; - for (i=0;i 0 - const double same = c * sqrt - ( - static_cast(Fact(l - m)) / - static_cast(Fact(l + m)) - ) - *SQRT2; + const double same = c * sqrt(static_cast(Fact(l - m)) / static_cast(Fact(l + m))) * SQRT2; ++lm; - for (i=0;i1; i--) + for (int i = n; i > 1; i--) { f *= i; -// if(n>16) -// { -// cout<<"\n n="<16) + // { + // cout<<"\n n="<16) QUIT(); + // if(n>16) QUIT(); return f; } @@ -286,10 +274,9 @@ long double Fact(const int n) int Semi_Fact(const int n) { int semif = 1; - for (int i=n; i>2; i -= 2) + for (int i = n; i > 2; i -= 2) { semif *= i; } return semif; } - diff --git a/tools/01_NAO_generation/SIAB/src_pw/ylm_real.h b/tools/01_NAO_generation/SIAB/src_pw/ylm_real.h index 84a344392e7..7bbb98609f4 100644 --- a/tools/01_NAO_generation/SIAB/src_pw/ylm_real.h +++ b/tools/01_NAO_generation/SIAB/src_pw/ylm_real.h @@ -1,16 +1,14 @@ #ifndef YLM_REAL_H #define YLM_REAL_H -#include "../src_tools/vector3.h" -#include "../src_tools/realarray.h" #include "../src_tools/matrix.h" +#include "../src_tools/realarray.h" +#include "../src_tools/vector3.h" using namespace std; -void Ylm_Real -( - const int lmax2, // lmax2 = (lmax+1)^2 - const int ng, // - const Vector3 *g, // g_cartesian_vec(x,y,z) - matrix &ylm // output +void Ylm_Real(const int lmax2, // lmax2 = (lmax+1)^2 + const int ng, // + const Vector3* g, // g_cartesian_vec(x,y,z) + matrix& ylm // output ); long double Fact(const int n); diff --git a/tools/01_NAO_generation/SIAB/src_unittest/src_tools/common_test.h b/tools/01_NAO_generation/SIAB/src_unittest/src_tools/common_test.h index 67e3c09bdab..b2f02614742 100644 --- a/tools/01_NAO_generation/SIAB/src_unittest/src_tools/common_test.h +++ b/tools/01_NAO_generation/SIAB/src_unittest/src_tools/common_test.h @@ -1,56 +1,56 @@ #ifndef COMMON_TEST #define COMMON_TEST -#include -#include -#include +#include +#include +#include using std::cout; using std::endl; -#include +#include extern double sum_t; // gettimeofday( &time_now, NULL); -inline double get_time(const timeval &begin, const timeval &end) +inline double get_time(const timeval& begin, const timeval& end) { - return (double)(end.tv_sec-begin.tv_sec) + (double)(end.tv_usec-begin.tv_usec)/1000000.0; + return (double)(end.tv_sec - begin.tv_sec) + (double)(end.tv_usec - begin.tv_usec) / 1000000.0; } -inline void cout_matrix( const ComplexMatrix &matrix) +inline void cout_matrix(const ComplexMatrix& matrix) { - cout<< std::setprecision(10)<< std::setiosflags(ios::fixed); - const size_t nr(matrix.nr), nc(matrix.nc); - for(size_t i(0); i