diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Calculate_C4.cpp b/tools/01_NAO_generation/SIAB/src_spillage/Calculate_C4.cpp index 9a96e259c03..91f4aba870c 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Calculate_C4.cpp +++ b/tools/01_NAO_generation/SIAB/src_spillage/Calculate_C4.cpp @@ -5,39 +5,47 @@ // 2009-05-19 change file name //========================================================== #include "Calculate_C4.h" + #include "tools.h" int Calculate_C4::test = 0; -Calculate_C4::Calculate_C4(){} -Calculate_C4::~Calculate_C4(){} +Calculate_C4::Calculate_C4() +{ +} +Calculate_C4::~Calculate_C4() +{ +} //=================================================================== // FUNCTION : Find_Eigenvalues -// DESCRIPTION : Find 'Enumber' number of eigenvalues of +// DESCRIPTION : Find 'Enumber' number of eigenvalues of // spherical bessel function from '0~lmax', and put the -// result eigenvalues in 'en' matrix. +// result eigenvalues in 'en' matrix. // P.S: (1) The eigenvalues are calculated '< tolerence' // accuracy. // P.S: (2) The rcut is important for cacluating eigenvalues // because it makes jl( rcut * eigenvalue ) = 0.00 //==================================================================== -void Calculate_C4::Find_Eigenvalues(const int &Enumber, matrix &en, - const int &lmax, const double &rcut, const double &tolerence) +void Calculate_C4::Find_Eigenvalues(const int& Enumber, + matrix& en, + const int& lmax, + const double& rcut, + const double& tolerence) { - double *e = new double[Enumber]; - for(int l=0; l( rcut / dr)+1; - if(mesh % 2 == 0) mesh++; - cout << "\n\n Mesh for test = " << mesh << " L = " << L; - - double *r = new double[mesh]; - double *rab = new double[mesh]; - double *e = new double[Enumber]; - double *jle = new double[mesh]; - double *jtest = new double[mesh]; - double *function = new double[mesh]; - - for(int i=0; i(rcut / dr) + 1; + if (mesh % 2 == 0) + mesh++; + cout << "\n\n Mesh for test = " << mesh << " L = " << L; + + double* r = new double[mesh]; + double* rab = new double[mesh]; + double* e = new double[Enumber]; + double* jle = new double[mesh]; + double* jtest = new double[mesh]; + double* function = new double[mesh]; + + for (int i = 0; i < Enumber; i++) + { + e[i] = en(L, i); + } + + for (int ir = 0; ir < mesh; ir++) + { + r[ir] = ir * dr; + rab[ir] = dr; + } + + Mathzone::Spherical_Bessel(mesh, r, e[test_ei], L, jtest); + + cout << "\n r=" << r[mesh - 1] << "(a.u.) jtest=" << jtest[mesh - 1]; + + double C4 = 0.0; + for (int i = 0; i < Enumber; i++) + { + Mathzone::Spherical_Bessel(mesh, r, e[i], L, jle); + for (int ir = 0; ir < mesh; ir++) + { + function[ir] = jle[ir] * jtest[ir] * r[ir] * r[ir]; + } + Mathzone::Simpson_Integral(mesh, function, rab, C4); + // if( abs(C4) < 10e-5 ) + // { + // cout << "\n" << setw(5) << i << setw(15) << "0"; + // } + // else + // { + cout << "\n" << setw(5) << i << setw(15) << C4; + // } + } + + delete[] function; + delete[] jtest; + delete[] r; + delete[] rab; + delete[] e; + delete[] jle; + return; } //=================================================================== @@ -120,441 +133,450 @@ void Calculate_C4::Test_Orthogonal(const double &rcut, const double &dr, // (5) ues int_{0}^{Rc} J1(q1*r)*Jl(q1*r)*r*r dr ==> norm // (6) C4 = C4_firstly / norm //=================================================================== -void Calculate_C4::key_solver(const string &name, const int &L, const matrix &en, const int &Enumber, double *C4) +void Calculate_C4::key_solver(const string& name, const int& L, const matrix& en, const int& Enumber, double* C4) { - ifstream ifs(name.c_str()); - if(!ifs) - { - cout << "\n Can't find file : " << name; - } - else - { - cout << "\n Open file " << name; - } - - string title; - - int mesh; - ifs >> title >> mesh; - assert(mesh > 0); - cout << "\t mesh = "<> title >> title >> title; - - for(int i=0; i> r[i] >> psi[i] >> rab[i]; - } - - double norm = 0.0; - for(int i=0; i> title >> mesh; + assert(mesh > 0); + cout << "\t mesh = " << mesh; + + double* r = new double[mesh]; + double* psi = new double[mesh]; + double* rab = new double[mesh]; + double* e = new double[Enumber]; + double* jle = new double[mesh]; + double* function = new double[mesh]; + double* function2 = new double[mesh]; + + for (int i = 0; i < Enumber; i++) + { + e[i] = en(L, i); + } + + ifs >> title >> title >> title; + + for (int i = 0; i < mesh; i++) + { + ifs >> r[i] >> psi[i] >> rab[i]; + } + + double norm = 0.0; + for (int i = 0; i < Enumber; i++) + { + Mathzone::Spherical_Bessel(mesh, r, e[i], L, jle); + for (int ir = 0; ir < mesh; ir++) + { + function[ir] = jle[ir] * psi[ir] * r[ir]; + } + + for (int ir = 0; ir < mesh; ir++) + { + function2[ir] = jle[ir] * jle[ir] * r[ir] * r[ir]; + } + + Mathzone::Simpson_Integral(mesh, function, rab, C4[i]); + Mathzone::Simpson_Integral(mesh, function2, rab, norm); + assert(norm != 0.0); + C4[i] /= norm; + // cout << "\n norm = "<(rcut / dr)+1; //mohan update 2009-10-10 - if(mesh % 2 == 0) mesh++; - - if(test==1) - { - cout << "\n norm : mesh = " << mesh; - cout << "\n type = " << ntype; - cout << "\n lmax = " << lmax; - cout << "\n enumber = " << enumber; - } - - // init radial mesh - double *r = new double[mesh]; - for(int i=0; i(rcut / dr) + 1; // mohan update 2009-10-10 + if (mesh % 2 == 0) + mesh++; + + if (test == 1) + { + cout << "\n norm : mesh = " << mesh; + cout << "\n type = " << ntype; + cout << "\n lmax = " << lmax; + cout << "\n enumber = " << enumber; + } + + // init radial mesh + double* r = new double[mesh]; + for (int i = 0; i < mesh; i++) + { + r[i] = i * dr; + } + + // init relta r, used in simpson integral + double* rab = new double[mesh]; + for (int i = 0; i < mesh; i++) + { + rab[i] = dr; + } + + double* e = new double[enumber]; + double* psi = new double[mesh]; + double* function = new double[mesh]; + double* jle = new double[mesh]; + + for (int it = 0; it < ntype; it++) + { + for (int l = 0; l < lmax + 1; l++) + { + // init eigenvalue of order l + for (int i = 0; i < enumber; i++) + { + e[i] = en(l, i); + } + for (int n = 0; n < nmax; n++) + { + ZEROS(jle, mesh); + ZEROS(psi, mesh); + for (int ie = 0; ie < enumber; ie++) + { + Mathzone::Spherical_Bessel(mesh, r, e[ie], l, jle); + for (int ir = 0; ir < mesh; ir++) + { + psi[ir] += C4(it, l, n, ie) * jle[ir]; + } + } + + for (int ir = 0; ir < mesh; ir++) + { + function[ir] = psi[ir] * psi[ir] * r[ir] * r[ir]; + } + double norm = 0.0; + Mathzone::Simpson_Integral(mesh, function, rab, norm); + + if (test == 2) + cout << "it=" << it << " L=" << l << " N=" << n << " norm = " << norm << endl; + norm = sqrt(norm); + assert(norm > 1.0e-5); + for (int ie = 0; ie < enumber; ie++) + { + C4(it, l, n, ie) /= norm; + accept_rate(it, l, n, ie) /= norm; // Peize Lin add 2015-12-24 + } + } + } + } + + delete[] r; + delete[] psi; + delete[] e; + delete[] jle; + return; } -void Calculate_C4::plot(ifstream &ifs) +void Calculate_C4::plot(ifstream& ifs) { - TITLE("Calculate_C4","plot"); - int nfiles; - double dr; - READ_VALUE( ifs, nfiles); - READ_VALUE( ifs, dr); - cout << "\n Number of files to plot : " << nfiles; - cout << "\n dr = " << dr; - assert( nfiles > 0); - assert( dr > 0.0 ); - - int l; - double ecut; - double rcut; - int enumber = 0; - double tolerence; - double *c4; - double *r; - double *e; - double *psi; - double *jle; - for(int i=0; i") ) - { - string filedir; - READ_VALUE( ifs, filedir ); - READ_VALUE( ifs, l ); - READ_VALUE( ifs, ecut ); - READ_VALUE( ifs, rcut ); - READ_VALUE( ifs, enumber ); - READ_VALUE( ifs, tolerence ); - cout << "\n l = " << l; - cout << "\n ecut = " << ecut; - cout << "\n enumber = " << enumber; - SCAN_END(ifs, ""); - } - - // init c4 - c4 = new double[enumber]; - ZEROS(c4, enumber); - if( SCAN_BEGIN(ifs,"") ) - { - assert(enumber>0); - for(int i=0; i> id >> c4[i]; - // cout << "\n c4[" << i << "]=" << c4[i]; - } - SCAN_END(ifs, ""); - } - - // init mesh - int mesh = static_cast(rcut / dr)+1; - if(mesh % 2 == 0) mesh++; - cout << "\n mesh = " << mesh; - - // init r - r = new double[mesh]; - for(int i=0; i")) + { + string filedir; + READ_VALUE(ifs, filedir); + READ_VALUE(ifs, l); + READ_VALUE(ifs, ecut); + READ_VALUE(ifs, rcut); + READ_VALUE(ifs, enumber); + READ_VALUE(ifs, tolerence); + cout << "\n l = " << l; + cout << "\n ecut = " << ecut; + cout << "\n enumber = " << enumber; + SCAN_END(ifs, ""); + } + + // init c4 + c4 = new double[enumber]; + ZEROS(c4, enumber); + if (SCAN_BEGIN(ifs, "")) + { + assert(enumber > 0); + for (int i = 0; i < enumber; i++) + { + int id; + ifs >> id >> c4[i]; + // cout << "\n c4[" << i << "]=" << c4[i]; + } + SCAN_END(ifs, ""); + } + + // init mesh + int mesh = static_cast(rcut / dr) + 1; + if (mesh % 2 == 0) + mesh++; + cout << "\n mesh = " << mesh; + + // init r + r = new double[mesh]; + for (int i = 0; i < mesh; i++) + { + r[i] = dr * i; + } + + // init eigenvalue en + assert(l < 3); + matrix en(l + 1, enumber); + Find_Eigenvalues(enumber, en, l, rcut, tolerence); + + // init e + e = new double[enumber]; + for (int i = 0; i < enumber; i++) + { + e[i] = en(l, i); + } + + // init psi + psi = new double[mesh]; + ZEROS(psi, mesh); + + // init jle + jle = new double[mesh]; + + // sum all jle + for (int ie = 0; ie < enumber; ie++) + { + ZEROS(jle, mesh); + Mathzone::Spherical_Bessel(mesh, r, e[ie], l, jle); + for (int ir = 0; ir < mesh; ir++) + { + psi[ir] += c4[ie] * jle[ir]; + } + } + + // unit + double norm; + double* function = new double[mesh]; + double* rab = new double[mesh]; + for (int ir = 0; ir < mesh; ir++) + { + rab[ir] = dr; + function[ir] = psi[ir] * psi[ir] * r[ir] * r[ir]; + } + + Mathzone::Simpson_Integral(mesh, function, rab, norm); + cout << "\n norm = " << norm; + assert(norm > 1.0e-5); // small number, + norm = sqrt(norm); + + stringstream ss; + ss << filenames << ".plot.txt"; + + ofstream ofs(ss.str().c_str()); + for (int ir = 0; ir < mesh; ir++) + { + ofs << ir * dr << " " << psi[ir] / norm << endl; + } + + delete[] c4; + delete[] r; + delete[] e; + delete[] psi; + delete[] jle; + delete[] function; + delete[] rab; + } + + exit(0); + return; } -void Calculate_C4::Test(ifstream &ifs, const double &ecut, const double &tolerence ) +void Calculate_C4::Test(ifstream& ifs, const double& ecut, const double& tolerence) { - TITLE("Calculate_C4","Test"); - double rcut;// rcut for input orbitals - double dr;// dr for simpson integrals( only used in test case), - // in output C4 case : dr is read from file as 'rab' - int test_ei;// the eigenvalue we want to test - int lmax; // max value of L - - READ_VALUE( ifs, rcut ); - READ_VALUE( ifs, dr ); - READ_VALUE( ifs, test_ei); - READ_VALUE( ifs, lmax ); - - cout << "\n Energy cutoff = "<< ecut; - cout << "\n tolerence = " << tolerence; - cout << "\n rcut = "< get_C4 Program Start." << endl; - - cout << "\n rcut = " << rcut << " (a.u.)"; - cout << "\n ecut = " << ecut << " (Ry)"; - cout << "\n tolerence = " << tolerence; - cout << "\n lmax = " << lmax; - - int nw;//number of orbitals(files) to calculate C4. - string *name;// file name, containing psi(r)*r, rab, radial - int *L;//angular momentum for each file( each orbital); - - READ_VALUE( ifs, nw ); - assert(nw>0); - cout << "\n Files = " << nw; - - name = new string[nw]; - L = new int[nw]; - //================================================ - // read in the file name and the corresponding L. - //================================================ - for(int i=0; i( sqrt( ecut * 2) * rcut / Mathzone::PI ); - cout << "\n Eigenvalue number(Jl) = " << Enumber; - - //================================================ - // Calculate eigenvalue of Jl - // CALL FUNCTION : Find_Eigenvalues - //================================================ - matrix en( lmax+1, Enumber); - Find_Eigenvalues(Enumber, en, lmax, rcut, tolerence); - - for(int i=0; ikey_solver(name[i], L[i], en, Enumber, C4); - stringstream ss; - ss << name[i] << ".C4"; - ofstream ofs( ss.str().c_str() ); - - ofs << "Calculate_C4_Program"; - ofs << "\n"; - ofs << "\n AUTHOR : Mohan Chen"; - ofs << "\n Date : 2009-4-1"; - ofs << "\n LOCATION : LQCC, Hefei, China"; - ofs << "\n EMAIL : mohan@mail.ustc.edu.cn"; - ofs << "\n Description : Calculate the coefficients C4 of f(r) in Spherical Bessel Basis J(qr)."; - ofs << "\n Formula : C4 = integral(r)[ f(r)*jl(qr)*r^{2} ]dr "; - ofs << "\n P.S. : We default consider f(r) read from file is in the form : ( f(r) * r )."; - ofs << "\n"; - - ofs << "\n\n"; - ofs << "\n" << name[i] << " ( f(r) read from this file)"; - ofs << "\n" << setw(8) << L[i] << " angular momentum for this orbital."; - ofs << "\n" << setw(8) << ecut << " Energy cutoff(Ry.)."; - ofs << "\n" << setw(8) << rcut << " rcut (a.u.)"; - ofs << "\n" << setw(8) << Enumber << " eigenvalue number( sqrt(ecut*2)*rcut/PI )."; - ofs << "\n" << setw(8) << tolerence << " tolerence to calculate eigenvalue."; - ofs << "\n"; - - ofs << "\n\n"; - for(int j=0; j"; - ofs.close(); - delete[] C4; - } - - delete[] name; - delete[] L; + cout << "\n\n ==> get_C4 Program Start." << endl; + + cout << "\n rcut = " << rcut << " (a.u.)"; + cout << "\n ecut = " << ecut << " (Ry)"; + cout << "\n tolerence = " << tolerence; + cout << "\n lmax = " << lmax; + + int nw; // number of orbitals(files) to calculate C4. + string* name; // file name, containing psi(r)*r, rab, radial + int* L; // angular momentum for each file( each orbital); + + READ_VALUE(ifs, nw); + assert(nw > 0); + cout << "\n Files = " << nw; + + name = new string[nw]; + L = new int[nw]; + //================================================ + // read in the file name and the corresponding L. + //================================================ + for (int i = 0; i < nw; i++) + { + READ_VALUE(ifs, name[i]); + READ_VALUE(ifs, L[i]); + } + + //================================================ + // Calculate the eigenvalues number from ecut. + //================================================ + const int Enumber = static_cast(sqrt(ecut * 2) * rcut / Mathzone::PI); + cout << "\n Eigenvalue number(Jl) = " << Enumber; + + //================================================ + // Calculate eigenvalue of Jl + // CALL FUNCTION : Find_Eigenvalues + //================================================ + matrix en(lmax + 1, Enumber); + Find_Eigenvalues(Enumber, en, lmax, rcut, tolerence); + + for (int i = 0; i < nw; i++) + { + double* C4 = new double[Enumber]; + //=========================================== + // CALL FUNCTION : Calculate + // calculate C4. + //=========================================== + this->key_solver(name[i], L[i], en, Enumber, C4); + stringstream ss; + ss << name[i] << ".C4"; + ofstream ofs(ss.str().c_str()); + + ofs << "Calculate_C4_Program"; + ofs << "\n"; + ofs << "\n AUTHOR : Mohan Chen"; + ofs << "\n Date : 2009-4-1"; + ofs << "\n LOCATION : LQCC, Hefei, China"; + ofs << "\n EMAIL : mohan@mail.ustc.edu.cn"; + ofs << "\n Description : Calculate the coefficients C4 of f(r) in Spherical Bessel Basis J(qr)."; + ofs << "\n Formula : C4 = integral(r)[ f(r)*jl(qr)*r^{2} ]dr "; + ofs << "\n P.S. : We default consider f(r) read from file is in the form : ( f(r) * r )."; + ofs << "\n"; + + ofs << "\n\n"; + ofs << "\n" << name[i] << " ( f(r) read from this file)"; + ofs << "\n" << setw(8) << L[i] << " angular momentum for this orbital."; + ofs << "\n" << setw(8) << ecut << " Energy cutoff(Ry.)."; + ofs << "\n" << setw(8) << rcut << " rcut (a.u.)"; + ofs << "\n" << setw(8) << Enumber << " eigenvalue number( sqrt(ecut*2)*rcut/PI )."; + ofs << "\n" << setw(8) << tolerence << " tolerence to calculate eigenvalue."; + ofs << "\n"; + + ofs << "\n\n"; + for (int j = 0; j < Enumber; j++) + { + ofs << "\n" << setw(5) << j << setiosflags(ios::scientific) << setprecision(8) << setw(25) << C4[j]; + } + ofs << "\n"; + ofs.close(); + delete[] C4; + } + + delete[] name; + delete[] L; } diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Calculate_C4.h b/tools/01_NAO_generation/SIAB/src_spillage/Calculate_C4.h index d50a2fbda5c..7f8056966bc 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Calculate_C4.h +++ b/tools/01_NAO_generation/SIAB/src_spillage/Calculate_C4.h @@ -5,35 +5,42 @@ class Calculate_C4 { - public: - - Calculate_C4(); - ~Calculate_C4(); - - void init( ifstream &ifs, const double &ecut, const double &rcut, const double &tolerence, const int &lmax ); - void Test( ifstream &ifs, const double &ecut, const double &tolerence); - void plot( ifstream &ifs ); - - // Peize Lin add accept_rate 2015-12-24 - static void norm_ic( realArray &C4, realArray &accept_rate, const int &ntype, const int &lmax, - const int &nmax, const int &enumber, const double &tolerence, - const double &rcut, const double &dr); - - static void Find_Eigenvalues(const int &Enumber, matrix &en, - const int &lmax, const double &rcut, const double &tolerence); - - private: - - void Test_Orthogonal(const double &rcut, const double &dr, - const int &L, const matrix &en, const int &Enumber, const int &test_ei); - - void key_solver(const string &name, const int &L, - const matrix &en, const int &Enumber, double *C4); - - static int test; - - - + public: + Calculate_C4(); + ~Calculate_C4(); + + void init(ifstream& ifs, const double& ecut, const double& rcut, const double& tolerence, const int& lmax); + void Test(ifstream& ifs, const double& ecut, const double& tolerence); + void plot(ifstream& ifs); + + // Peize Lin add accept_rate 2015-12-24 + static void norm_ic(realArray& C4, + realArray& accept_rate, + const int& ntype, + const int& lmax, + const int& nmax, + const int& enumber, + const double& tolerence, + const double& rcut, + const double& dr); + + static void Find_Eigenvalues(const int& Enumber, + matrix& en, + const int& lmax, + const double& rcut, + const double& tolerence); + + private: + void Test_Orthogonal(const double& rcut, + const double& dr, + const int& L, + const matrix& en, + const int& Enumber, + const int& test_ei); + + void key_solver(const string& name, const int& L, const matrix& en, const int& Enumber, double* C4); + + static int test; }; #endif diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Coefficients.cpp b/tools/01_NAO_generation/SIAB/src_spillage/Coefficients.cpp index 5d7cf79c257..7b7dc135226 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Coefficients.cpp +++ b/tools/01_NAO_generation/SIAB/src_spillage/Coefficients.cpp @@ -1,4 +1,5 @@ #include "Coefficients.h" + #include "../src_parallel/parallel_common.h" bool Coefficients::write_c4_flag; @@ -7,25 +8,25 @@ Coefficients::Coefficients() { test = 0; output_c4_name = new string[1]; - accumulate_num = 0;//mohan add 2009-09-25 + accumulate_num = 0; // mohan add 2009-09-25 } Coefficients::~Coefficients() { - if (TEST1) cout << "\n ~Coefficients()" << endl; + if (TEST1) + cout << "\n ~Coefficients()" << endl; delete[] output_c4_name; } -void Coefficients::init( - const int &ntype_in, - const int &lmax_in, - const int &nmax_in, - const double &ecut_in, - const double &ecut_jlq_in, // mohan add 2009-07-23 - const double &rcut_in, - const double &tolerence_in) +void Coefficients::init(const int& ntype_in, + const int& lmax_in, + const int& nmax_in, + const double& ecut_in, + const double& ecut_jlq_in, // mohan add 2009-07-23 + const double& rcut_in, + const double& tolerence_in) { - TITLE(ofs_running,"Coefficients", "Constructor"); + TITLE(ofs_running, "Coefficients", "Constructor"); ntype = ntype_in; lmax = lmax_in; @@ -33,21 +34,21 @@ void Coefficients::init( ecut = ecut_in; ecut_jlq = ecut_jlq_in; // mohan add 2009-07-23 rcut = rcut_in; -// in Rydberg Unit - enumber = static_cast( sqrt( ecut )* rcut/PI ); - enumber_jlq = static_cast ( sqrt( ecut_jlq) * rcut/PI ); // mohan add 2009-07-23 - tolerence = tolerence_in; - - OUT(ofs_running,"ntype",ntype); - OUT(ofs_running,"lmax",lmax); - OUT(ofs_running,"nmax",nmax); - OUT(ofs_running,"ecut(Ry)",ecut); - OUT(ofs_running,"ecut_jlq(Ry)",ecut_jlq); - OUT(ofs_running,"rcut(Bohr)",rcut); - OUT(ofs_running,"enumber",enumber); - OUT(ofs_running,"enumber_jlq",enumber_jlq); - OUT(ofs_running,"tolerence",tolerence); - + // in Rydberg Unit + enumber = static_cast(sqrt(ecut) * rcut / PI); + enumber_jlq = static_cast(sqrt(ecut_jlq) * rcut / PI); // mohan add 2009-07-23 + tolerence = tolerence_in; + + OUT(ofs_running, "ntype", ntype); + OUT(ofs_running, "lmax", lmax); + OUT(ofs_running, "nmax", nmax); + OUT(ofs_running, "ecut(Ry)", ecut); + OUT(ofs_running, "ecut_jlq(Ry)", ecut_jlq); + OUT(ofs_running, "rcut(Bohr)", rcut); + OUT(ofs_running, "enumber", enumber); + OUT(ofs_running, "enumber_jlq", enumber_jlq); + OUT(ofs_running, "tolerence", tolerence); + assert(ntype > 0); assert(lmax >= 0); assert(nmax > 0); @@ -69,39 +70,43 @@ void Coefficients::init( void Coefficients::allocate_C4(void) { - if (test==1)TITLE("Coefficients", "allocate_C4"); + if (test == 1) + TITLE("Coefficients", "allocate_C4"); this->C4.create(ntype, lmax + 1, nmax, enumber); this->C4_old.create(ntype, lmax + 1, nmax, enumber); - this->C4_accumulate.create(ntype, lmax + 1, nmax, enumber); // mohan add 2009-09-25 -// this->acc_count.create(ntype, lmax+1, nmax, enumber); - this->accept_number.create(ntype, lmax+1, nmax, enumber); // mohan add 2009-10-16 - this->accept_rate.create(ntype, lmax+1, nmax, enumber);// mohan add 2009-10-31 + this->C4_accumulate.create(ntype, + lmax + 1, + nmax, + enumber); // mohan add 2009-09-25 + // this->acc_count.create(ntype, lmax+1, nmax, enumber); + this->accept_number.create(ntype, lmax + 1, nmax, enumber); // mohan add 2009-10-16 + this->accept_rate.create(ntype, lmax + 1, nmax, enumber); // mohan add 2009-10-31 // init accept_rate - for (int i=0; iaccept_rate.ptr[i] = 1.0; } for (int it = 0; it < ntype; it++) { - for (int il = 0; il < lmax+1; il++) + for (int il = 0; il < lmax + 1; il++) { for (int in = 0; in < nmax; in++) { - ofs_running << "\n Init C4: T=" << it+1 << " L=" << il << " N=" << in; + ofs_running << "\n Init C4: T=" << it + 1 << " L=" << il << " N=" << in; for (int ie = 0; ie < enumber; ie++) { // mohan modify 2009-08-26 int ne_cut; int ne_cut_min = -1; // must initalized - if (BLOCK_NE_MIN>0) + if (BLOCK_NE_MIN > 0) { ne_cut_min = BLOCK_NE_MIN; } - if (BLOCK_NE<0) + if (BLOCK_NE < 0) { ne_cut = enumber_jlq; } @@ -110,7 +115,7 @@ void Coefficients::allocate_C4(void) ne_cut = BLOCK_NE; } - if (iene_cut_min) + if (ie < ne_cut && ie > ne_cut_min) { this->C4(it, il, in, ie) = Random::between0and1(); this->C4_old(it, il, in, ie) = this->C4(it, il, in, ie); @@ -120,7 +125,8 @@ void Coefficients::allocate_C4(void) this->C4(it, il, in, ie) = 0.0; this->C4_old(it, il, in, ie) = 0.0; } - if (ie%4==0)ofs_running< 1. if (RESTART) { - if (MY_RANK==0) + if (MY_RANK == 0) { - ifstream ifs( "ORBITAL_RESULTS.txt" ); + ifstream ifs("ORBITAL_RESULTS.txt"); if (!ifs) { - WARNING_QUIT("Coefficients::allocate_C4()","can't find the restart file: ORBITAL_RESULTS.txt"); + WARNING_QUIT("Coefficients::allocate_C4()", "can't find the restart file: ORBITAL_RESULTS.txt"); } else { @@ -149,7 +155,7 @@ void Coefficients::allocate_C4(void) double ecut_in; double rcut_in; int ne_in; - if ( SCAN_BEGIN(ifs,"") ) + if (SCAN_BEGIN(ifs, "")) { READ_VALUE(ifs, ecut_in); READ_VALUE(ifs, rcut_in); @@ -158,11 +164,11 @@ void Coefficients::allocate_C4(void) { cout << "\n ne_in = " << ne_in; cout << "\n enumber = " << enumber; - WARNING_QUIT("Coefficients::Allocate_C4()","enumber != ne_in"); + WARNING_QUIT("Coefficients::Allocate_C4()", "enumber != ne_in"); } } - if ( SCAN_BEGIN(ifs,"") ) + if (SCAN_BEGIN(ifs, "")) { int nchi_in; READ_VALUE(ifs, nchi_in); @@ -174,7 +180,7 @@ void Coefficients::allocate_C4(void) int count_nchi = 0; for (int it = 0; it < ntype; it++) { - for (int il = 0; il < lmax+1; il++) + for (int il = 0; il < lmax + 1; il++) { for (int in = 0; in < nmax; in++) { @@ -184,11 +190,11 @@ void Coefficients::allocate_C4(void) assert(no_use1 == "Type"); assert(no_use2 == "L"); assert(no_use3 == "Zeta-Orbital"); - ifs >> it2>> il2>> in2; - it2 -= 1; - in2 -= 1; + ifs >> it2 >> il2 >> in2; + it2 -= 1; + in2 -= 1; } - if (it2!=it || il2!=il || in2!=in) + if (it2 != it || il2 != il || in2 != in) { hold_on = true; continue; @@ -196,68 +202,66 @@ void Coefficients::allocate_C4(void) else { hold_on = false; - cout << "\n\n Read in T=" << it+1 << " L=" << il << " Zeta-Orbital=" << in; - for (int ie=0; ieC4.getSize(); i++) + for (int i = 0; i < this->C4.getSize(); i++) { this->C4_old.ptr[i] = this->C4.ptr[i]; } return; } - -void Coefficients::trial_c4( const int &t, const int &l, const int &n, const int &ie ) +void Coefficients::trial_c4(const int& t, const int& l, const int& n, const int& ie) { // (1) this->C4_old(t, l, n, ie) = this->C4(t, l, n, ie); // update between -0.1~0.1 // (2) - this->C4(t, l, n, ie) += Random::betweenMinus1and1()/10.0* this->accept_rate(t, l, n, ie); + this->C4(t, l, n, ie) += Random::betweenMinus1and1() / 10.0 * this->accept_rate(t, l, n, ie); #ifdef __MPI - // don't need to bcast ! it's a random seed. + // don't need to bcast ! it's a random seed. // ofs_running << "\n C4=" << this->C4(t, l, n, ie); // Parallel_Common::bcast_double( this->C4(t, l, n, ie) ); #endif -// if( abs(this->C4(t, l, n, ie)) < 1.0e-3) //mohan modify 2009-10-10 -// { -// this->C4(t, l, n, ie) = 0.0; -// } -// this->C4(t, l, n, ie) += Random::betweenMinus1and1(); + // if( abs(this->C4(t, l, n, ie)) < 1.0e-3) //mohan modify 2009-10-10 + // { + // this->C4(t, l, n, ie) = 0.0; + // } + // this->C4(t, l, n, ie) += Random::betweenMinus1and1(); // this->C4(t, l, n, ie) *= Random::betweenMinus2and2(); // cout << "\n ic = " << ic << " ie = " << ie << " C4 = " << C4(t, l, n, ie); @@ -267,26 +271,26 @@ void Coefficients::trial_c4( const int &t, const int &l, const int &n, const int /* void Coefficients::accumulate_num_zero(void) { - cout << "\n Accumulate_num = " << this->accumulate_num; - this->accumulate_num = 0; + cout << "\n Accumulate_num = " << this->accumulate_num; + this->accumulate_num = 0; } */ /* void Coefficients::accumulating_C4( const int &t, const int &l, const int &n, const int &ie ) { - // mohan add 2009-09-25 - double mix = this->C4_accumulate(t, l, n, ie) * this->acc_count(t,l,n,ie) + this->C4_old(t, l, n, ie); - ++acc_count(t,l,n,ie); - ++accumulate_num; + // mohan add 2009-09-25 + double mix = this->C4_accumulate(t, l, n, ie) * this->acc_count(t,l,n,ie) + this->C4_old(t, l, n, ie); + ++acc_count(t,l,n,ie); + ++accumulate_num; this->C4_accumulate(t, l, n, ie) = mix / acc_count(t,l,n,ie); // cout << "\n Accumulate_num = " << this->accumulate_num; - return; + return; } */ -void Coefficients::update_c4( const int &t, const int &l, const int &n, const int &ie ) +void Coefficients::update_c4(const int& t, const int& l, const int& n, const int& ie) { //=============================================== // store the new c4 parameters(random accepted!! @@ -294,12 +298,12 @@ void Coefficients::update_c4( const int &t, const int &l, const int &n, const in //=============================================== this->C4_old(t, l, n, ie) = this->C4(t, l, n, ie); - ++accept_number(t,l,n,ie); + ++accept_number(t, l, n, ie); return; } -void Coefficients::go_back_c4( const int &t, const int &l, const int &n, const int &ie ) +void Coefficients::go_back_c4(const int& t, const int& l, const int& n, const int& ie) { //================================== // give up the random added number; @@ -311,40 +315,40 @@ void Coefficients::go_back_c4( const int &t, const int &l, const int &n, const i /* void Coefficients::kinetic_energy( const int &t, const int &l, const int &n, const int &ie) { - if(!ke_flip_flag) - { - if(ke_up_flag) - { - // C4>0.0 and C4_old>0.0 - this->ke_total_up += ( C4(t, l, n, ie) - C4_old(t, l, n, ie) ) * ie * ie; - } - else - { - // C4<0.0 and C4_old>0.0 - this->ke_total_down += ( C4_old(t, l, n, ie) - C4(t, l, n, ie) ) * ie * ie; - } - } - else - { - // (1) C4>0.0 and C4_old<0.0, - // so ke_total_up energy increase - // ke_total_down energy decrease - // - // (2) C4<0.0 and C4_old>0.0, - // so ke_total_up energy decrease - // and ke_total_down energy increase - this->ke_total_up += C4(t, l, n, ie) * ie * ie; - this->ke_total_down += C4_old(t, l, n, ie) * ie * ie; - } - - this->ke_total_old = this->ke_total_new; - this->ke_total_new = abs(abs(ke_total_up) - abs(ke_total_down)); + if(!ke_flip_flag) + { + if(ke_up_flag) + { + // C4>0.0 and C4_old>0.0 + this->ke_total_up += ( C4(t, l, n, ie) - C4_old(t, l, n, ie) ) * ie * ie; + } + else + { + // C4<0.0 and C4_old>0.0 + this->ke_total_down += ( C4_old(t, l, n, ie) - C4(t, l, n, ie) ) * ie * ie; + } + } + else + { + // (1) C4>0.0 and C4_old<0.0, + // so ke_total_up energy increase + // ke_total_down energy decrease + // + // (2) C4<0.0 and C4_old>0.0, + // so ke_total_up energy decrease + // and ke_total_down energy increase + this->ke_total_up += C4(t, l, n, ie) * ie * ie; + this->ke_total_down += C4_old(t, l, n, ie) * ie * ie; + } + + this->ke_total_old = this->ke_total_new; + this->ke_total_new = abs(abs(ke_total_up) - abs(ke_total_down)); // cout << "\n ke_total_up = " << ke_total_up; // cout << "\n ke_total_down = " << ke_total_down; // cout << "\n ke_total_old = " << ke_total_old; // cout << "\n ke_total_new = " << ke_total_new << endl; - return; + return; } */ diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Coefficients.h b/tools/01_NAO_generation/SIAB/src_spillage/Coefficients.h index 32c692e3452..81e49efaa51 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Coefficients.h +++ b/tools/01_NAO_generation/SIAB/src_spillage/Coefficients.h @@ -4,51 +4,55 @@ #include "common.h" class Coefficients { - public: - //============================ - // C4(type, L, N, ie) - //============================ - realArray C4; - realArray C4_old; - realArray C4_accumulate; - //IntArray acc_count; - IntArray accept_number; - realArray accept_rate; - - // constructor, available for SZ, DZ, DZP... - Coefficients(); - ~Coefficients(); - - void init( const int &ntype, const int &lmax, const int &namx, - const double &ecut, const double &ecut_jlq, - const double &rcut, const double &tolerence); - void trial_c4( const int &t, const int &l, const int &n, const int &ie ); - - void update_c4( const int &t, const int &l, const int &n, const int &ie ); - void go_back_c4( const int &t, const int &l, const int &n, const int &ie ); - void copy_c4_to_old(void); // for normalize - void kinetic_energy( const int &t, const int &l, const int &n, const int &ie);// mohan add 2009-07-23 - void accumulate_num_zero(void); - void accumulating_C4( const int &t, const int &l, const int &n, const int &ie ); - - int ntype; - int lmax; - int nmax; - int enumber; - double rcut; - double tolerence; - int fix; - - string *output_c4_name; - static bool write_c4_flag; - - private: - void allocate_C4(); - double ecut; - double ecut_jlq; // mohan add 2009-07-23 - int enumber_jlq; // mohan add 2009-07-23 - int test; - int accumulate_num; + public: + //============================ + // C4(type, L, N, ie) + //============================ + realArray C4; + realArray C4_old; + realArray C4_accumulate; + // IntArray acc_count; + IntArray accept_number; + realArray accept_rate; + + // constructor, available for SZ, DZ, DZP... + Coefficients(); + ~Coefficients(); + + void init(const int& ntype, + const int& lmax, + const int& namx, + const double& ecut, + const double& ecut_jlq, + const double& rcut, + const double& tolerence); + void trial_c4(const int& t, const int& l, const int& n, const int& ie); + + void update_c4(const int& t, const int& l, const int& n, const int& ie); + void go_back_c4(const int& t, const int& l, const int& n, const int& ie); + void copy_c4_to_old(void); // for normalize + void kinetic_energy(const int& t, const int& l, const int& n, const int& ie); // mohan add 2009-07-23 + void accumulate_num_zero(void); + void accumulating_C4(const int& t, const int& l, const int& n, const int& ie); + + int ntype; + int lmax; + int nmax; + int enumber; + double rcut; + double tolerence; + int fix; + + string* output_c4_name; + static bool write_c4_flag; + + private: + void allocate_C4(); + double ecut; + double ecut_jlq; // mohan add 2009-07-23 + int enumber_jlq; // mohan add 2009-07-23 + int test; + int accumulate_num; }; #endif diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Inverse_Matrix_S.cpp b/tools/01_NAO_generation/SIAB/src_spillage/Inverse_Matrix_S.cpp index d4a38a0e5b9..9a7f28f82ce 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Inverse_Matrix_S.cpp +++ b/tools/01_NAO_generation/SIAB/src_spillage/Inverse_Matrix_S.cpp @@ -1,48 +1,51 @@ #include "Inverse_Matrix_S.h" -#include "Inverse_Matrix_S_inline.h" + #include "../src_tools/inverse_matrix_iterate_inline.h" +#include "Inverse_Matrix_S_inline.h" -void Inverse_Matrix_S::init( int dim_in) +void Inverse_Matrix_S::init(int dim_in) { - this->dim = dim_in; - this->inverse.init(this->dim); - this->inverse_iterate.init(this->dim); - this->inverse_iterate.set_iterate_num(20); // test, undetermined - this->inverse_iterate.set_tolerate_error(1e-10); // test, undetermined - this->inverse_iterate.set_symmetric_info('H', 'U'); - this->set_exact_loop_num(1); // test, undetermined - return; + this->dim = dim_in; + this->inverse.init(this->dim); + this->inverse_iterate.init(this->dim); + this->inverse_iterate.set_iterate_num(20); // test, undetermined + this->inverse_iterate.set_tolerate_error(1e-10); // test, undetermined + this->inverse_iterate.set_symmetric_info('H', 'U'); + this->set_exact_loop_num(1); // test, undetermined + return; } -const ComplexMatrix & Inverse_Matrix_S::cal_inverse(ComplexMatrix &S) +const ComplexMatrix& Inverse_Matrix_S::cal_inverse(ComplexMatrix& S) { - if(0==loop_count) - { - inverse.using_zpotrf(S); - } - else - { - inverse_iterate.set_X0(*this->inverse_old) - .set_A(S); - try - { - if (1==loop_count) { copy_matrix_U2L(const_cast(*this->inverse_old)); } -// copy_matrix_U2L(const_cast(*this->inverse_old)); -// copy_matrix_U2L(S); + if (0 == loop_count) + { + inverse.using_zpotrf(S); + } + else + { + inverse_iterate.set_X0(*this->inverse_old).set_A(S); + try + { + if (1 == loop_count) + { + copy_matrix_U2L(const_cast(*this->inverse_old)); + } + // copy_matrix_U2L(const_cast(*this->inverse_old)); + // copy_matrix_U2L(S); + + inverse_iterate.cal_inverse(); + + // zero_matrix_L(const_cast(*this->inverse_old)); + // zero_matrix_L(S); + } + catch (const Inverse_Matrix_Iterate::Unstable& e) + { + // zero_matrix_L(const_cast(*this->inverse_old)); + // zero_matrix_L(S); - inverse_iterate.cal_inverse(); - -// zero_matrix_L(const_cast(*this->inverse_old)); -// zero_matrix_L(S); - } - catch(const Inverse_Matrix_Iterate::Unstable &e) - { -// zero_matrix_L(const_cast(*this->inverse_old)); -// zero_matrix_L(S); - - loop_count = 0; - inverse.using_zpotrf(S); - } - } - return this->get_inverse(); + loop_count = 0; + inverse.using_zpotrf(S); + } + } + return this->get_inverse(); } \ No newline at end of file diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Inverse_Matrix_S.h b/tools/01_NAO_generation/SIAB/src_spillage/Inverse_Matrix_S.h index f177b534da8..9dbd168423f 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Inverse_Matrix_S.h +++ b/tools/01_NAO_generation/SIAB/src_spillage/Inverse_Matrix_S.h @@ -6,30 +6,49 @@ class Inverse_Matrix_S { -public: - Inverse_Matrix inverse; - Inverse_Matrix_Iterate inverse_iterate; + public: + Inverse_Matrix inverse; + Inverse_Matrix_Iterate inverse_iterate; - void init( int dim_in); - const ComplexMatrix & cal_inverse(ComplexMatrix &S); - const ComplexMatrix & get_inverse() { return (0==loop_count) ? inverse.A : inverse_iterate.get_inverse(); } - inline void update(); - void reset_loop() { loop_count = 0; } + void init(int dim_in); + const ComplexMatrix& cal_inverse(ComplexMatrix& S); + const ComplexMatrix& get_inverse() + { + return (0 == loop_count) ? inverse.A : inverse_iterate.get_inverse(); + } + inline void update(); + void reset_loop() + { + loop_count = 0; + } -private: - size_t dim{0}; - size_t loop_count{0}; - size_t exact_loop_num{1}; - const ComplexMatrix *inverse_old = NULL; - - inline void copy_matrix_U2L(ComplexMatrix &matrix); - inline void zero_matrix_L(ComplexMatrix &matrix); - -public: - size_t get_exact_loop_num() { return this->exact_loop_num; } - size_t get_loop_count() { return this->loop_count; } - Inverse_Matrix_S & set_exact_loop_num(size_t exact_loop_num_in) { this->exact_loop_num = exact_loop_num_in; return *this;} - size_t get_dim() { return this->dim; } + private: + size_t dim{0}; + size_t loop_count{0}; + size_t exact_loop_num{1}; + const ComplexMatrix* inverse_old = NULL; + + inline void copy_matrix_U2L(ComplexMatrix& matrix); + inline void zero_matrix_L(ComplexMatrix& matrix); + + public: + size_t get_exact_loop_num() + { + return this->exact_loop_num; + } + size_t get_loop_count() + { + return this->loop_count; + } + Inverse_Matrix_S& set_exact_loop_num(size_t exact_loop_num_in) + { + this->exact_loop_num = exact_loop_num_in; + return *this; + } + size_t get_dim() + { + return this->dim; + } }; #endif diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Inverse_Matrix_S_inline.h b/tools/01_NAO_generation/SIAB/src_spillage/Inverse_Matrix_S_inline.h index bc102813645..f33bf7aab1e 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Inverse_Matrix_S_inline.h +++ b/tools/01_NAO_generation/SIAB/src_spillage/Inverse_Matrix_S_inline.h @@ -1,29 +1,29 @@ #include "Inverse_Matrix_S.h" -void Inverse_Matrix_S::copy_matrix_U2L(ComplexMatrix &matrix) +void Inverse_Matrix_S::copy_matrix_U2L(ComplexMatrix& matrix) { - for( size_t iw1(0); iw1dim; ++iw1) - { - for( size_t iw2(0); iw2dim; ++iw1) + { + for (size_t iw2(0); iw2 < iw1; ++iw2) + { + matrix(iw1, iw2) = conj(matrix(iw2, iw1)); + } + } } -void Inverse_Matrix_S::zero_matrix_L(ComplexMatrix &matrix) +void Inverse_Matrix_S::zero_matrix_L(ComplexMatrix& matrix) { - for( size_t iw1(0); iw1dim; ++iw1) - { - for( size_t iw2(0); iw2dim; ++iw1) + { + for (size_t iw2(0); iw2 < iw1; ++iw2) + { + matrix(iw1, iw2) = 0; + } + } } void Inverse_Matrix_S::update() { - inverse_old = &get_inverse(); - loop_count = (loop_count+1)%exact_loop_num; + inverse_old = &get_inverse(); + loop_count = (loop_count + 1) % exact_loop_num; } \ No newline at end of file diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Metropolis.cpp b/tools/01_NAO_generation/SIAB/src_spillage/Metropolis.cpp index 466e5ecad63..b40752993a9 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Metropolis.cpp +++ b/tools/01_NAO_generation/SIAB/src_spillage/Metropolis.cpp @@ -1,159 +1,163 @@ #include "Metropolis.h" + #include "../src_parallel/parallel_common.h" -#include "Calculate_C4.h" #include "../src_tools/Simulated_Annealing.h" +#include "Calculate_C4.h" #include "Simulated_Annealing_Orbital.h" Metropolis::Metropolis() { - prob_old = 0.0; - prob_new = 0.0; - temperature = 0.0; - random_c4number = 0; - output_each_istep = 0; - kin_last = new double[1]; - kin_temp = new double[1]; + prob_old = 0.0; + prob_new = 0.0; + temperature = 0.0; + random_c4number = 0; + output_each_istep = 0; + kin_last = new double[1]; + kin_temp = new double[1]; } - Metropolis::~Metropolis() { - delete[] kin_last; - delete[] kin_temp; + delete[] kin_last; + delete[] kin_temp; } // read in Metropolis parameters. void Metropolis::init(void) { - if(test==1) TITLE("Metropolis","init"); + if (test == 1) + TITLE("Metropolis", "init"); - bool begin = false; + bool begin = false; - if(MY_RANK==0) - { - if(SCAN_BEGIN(input.ifs, "")) - { - begin = true; - } - } + if (MY_RANK == 0) + { + if (SCAN_BEGIN(input.ifs, "")) + { + begin = true; + } + } #ifdef __MPI - Parallel_Common::bcast_bool(begin); + Parallel_Common::bcast_bool(begin); #endif - ofs_running << "\n begin read in Metropolis parameter." << endl; - - if(begin) - { - if(MY_RANK==0) - { - READ_VALUE( input.ifs, this->spi_ini_temp); - READ_VALUE( input.ifs, this->spi_rate); - READ_VALUE( input.ifs, this->spi_ntemp); - READ_VALUE( input.ifs, this->spi_nsteps); - - READ_VALUE( input.ifs, this->kin_ini_temp); - READ_VALUE( input.ifs, this->kin_rate); - READ_VALUE( input.ifs, this->kin_ntemp); - READ_VALUE( input.ifs, this->kin_nsteps); - - READ_VALUE( input.ifs, this->delta_kappa); - READ_VALUE( input.ifs, this->output_each_istep); - - READ_VALUE( input.ifs, ACCEPTANCE_STEPS ); // mohan add 2009-10-31 - READ_VALUE( input.ifs, ACCEPTANCE_HIGH ); - READ_VALUE( input.ifs, ACCEPTANCE_LOW ); - - READ_VALUE( input.ifs, KINETIC_MAX ); // mohan add 2009-10-31 - READ_VALUE( input.ifs, KINETIC_DR); //mohan add 2010-04-12 - READ_VALUE( input.ifs, OPTIMIZE_METHOD); // mohan add 2010-04-14 - //SCAN_END(input.ifs, ""); - } - + ofs_running << "\n begin read in Metropolis parameter." << endl; + + if (begin) + { + if (MY_RANK == 0) + { + READ_VALUE(input.ifs, this->spi_ini_temp); + READ_VALUE(input.ifs, this->spi_rate); + READ_VALUE(input.ifs, this->spi_ntemp); + READ_VALUE(input.ifs, this->spi_nsteps); + + READ_VALUE(input.ifs, this->kin_ini_temp); + READ_VALUE(input.ifs, this->kin_rate); + READ_VALUE(input.ifs, this->kin_ntemp); + READ_VALUE(input.ifs, this->kin_nsteps); + + READ_VALUE(input.ifs, this->delta_kappa); + READ_VALUE(input.ifs, this->output_each_istep); + + READ_VALUE(input.ifs, ACCEPTANCE_STEPS); // mohan add 2009-10-31 + READ_VALUE(input.ifs, ACCEPTANCE_HIGH); + READ_VALUE(input.ifs, ACCEPTANCE_LOW); + + READ_VALUE(input.ifs, KINETIC_MAX); // mohan add 2009-10-31 + READ_VALUE(input.ifs, KINETIC_DR); // mohan add 2010-04-12 + READ_VALUE(input.ifs, OPTIMIZE_METHOD); // mohan add 2010-04-14 + // SCAN_END(input.ifs, ""); + } + #ifdef __MPI - MPI_Barrier(MPI_COMM_WORLD); + MPI_Barrier(MPI_COMM_WORLD); - Parallel_Common::bcast_double(this->spi_ini_temp); - Parallel_Common::bcast_double(this->spi_rate); - Parallel_Common::bcast_int(this->spi_ntemp); - Parallel_Common::bcast_int(this->spi_nsteps); + Parallel_Common::bcast_double(this->spi_ini_temp); + Parallel_Common::bcast_double(this->spi_rate); + Parallel_Common::bcast_int(this->spi_ntemp); + Parallel_Common::bcast_int(this->spi_nsteps); - Parallel_Common::bcast_double(this->kin_ini_temp); - Parallel_Common::bcast_double(this->kin_rate); - Parallel_Common::bcast_int(this->kin_ntemp); - Parallel_Common::bcast_int(this->kin_nsteps); + Parallel_Common::bcast_double(this->kin_ini_temp); + Parallel_Common::bcast_double(this->kin_rate); + Parallel_Common::bcast_int(this->kin_ntemp); + Parallel_Common::bcast_int(this->kin_nsteps); - Parallel_Common::bcast_double(this->delta_kappa); - Parallel_Common::bcast_int(this->output_each_istep); + Parallel_Common::bcast_double(this->delta_kappa); + Parallel_Common::bcast_int(this->output_each_istep); - Parallel_Common::bcast_int(ACCEPTANCE_STEPS ); - Parallel_Common::bcast_double(ACCEPTANCE_HIGH ); - Parallel_Common::bcast_double(ACCEPTANCE_LOW ); + Parallel_Common::bcast_int(ACCEPTANCE_STEPS); + Parallel_Common::bcast_double(ACCEPTANCE_HIGH); + Parallel_Common::bcast_double(ACCEPTANCE_LOW); - Parallel_Common::bcast_double(KINETIC_MAX ); - Parallel_Common::bcast_double(KINETIC_DR ); - Parallel_Common::bcast_int(OPTIMIZE_METHOD); + Parallel_Common::bcast_double(KINETIC_MAX); + Parallel_Common::bcast_double(KINETIC_DR); + Parallel_Common::bcast_int(OPTIMIZE_METHOD); #endif - } - this->temperature = this->spi_ini_temp; - - // check whether the parameters are in the range. - assert(spi_ini_temp >= 0.0); - assert(spi_rate >= 0.0); + this->temperature = this->spi_ini_temp; + + // check whether the parameters are in the range. + assert(spi_ini_temp >= 0.0); + assert(spi_rate >= 0.0); assert(spi_ntemp >= 0); - assert(spi_nsteps >= 0); + assert(spi_nsteps >= 0); - assert(kin_ini_temp >= 0.0); - assert(kin_rate >= 0.0); + assert(kin_ini_temp >= 0.0); + assert(kin_rate >= 0.0); assert(kin_ntemp >= 0); - assert(kin_nsteps >= 0); + assert(kin_nsteps >= 0); + + assert(delta_kappa >= 0.0); + assert(output_each_istep > 0); - assert(delta_kappa >= 0.0); - assert(output_each_istep>0); - return; } // be called in MultiZeta::init(), reset the start temperature. -void Metropolis::reset_temperature(const int &istep) +void Metropolis::reset_temperature(const int& istep) { - temperature = this->spi_ini_temp; - for(int i=0; ispi_ini_temp; + for (int i = 0; i < istep; i++) + { + temperature *= 0.8; + // I think 0.8 is more reasonable than 0.1, mohan note, 2009-08-20 + } + return; } -void Metropolis::move_various_temperature( SpillageStep &Level) +void Metropolis::move_various_temperature(SpillageStep& Level) { - timer::tick("Metropolis","move_various_t"); - - // there are two parts of this algorithm, - // first part is about spillage, - // the second part is about kinetic energy. - - Simulated_Annealing_Orbital_Spillage sap(this,Level); - sap.set_temperature_num(this->spi_ntemp) - .set_step_num(this->spi_nsteps); - try - { - sap.cal_all_temperature(); - } - catch(File_Finish_Metropolis()) { return; } - - Simulated_Annealing_Orbital_Kinetic sak(this,Level); - sak.set_temperature_num(this->kin_ntemp) - .set_step_num(this->kin_nsteps); - try - { - sak.cal_all_temperature(); - } - catch(File_Finish_Metropolis()) { return; } - - timer::tick("Metropolis","move_various_t"); - return; + timer::tick("Metropolis", "move_various_t"); + + // there are two parts of this algorithm, + // first part is about spillage, + // the second part is about kinetic energy. + + Simulated_Annealing_Orbital_Spillage sap(this, Level); + sap.set_temperature_num(this->spi_ntemp).set_step_num(this->spi_nsteps); + try + { + sap.cal_all_temperature(); + } + catch (File_Finish_Metropolis()) + { + return; + } + + Simulated_Annealing_Orbital_Kinetic sak(this, Level); + sak.set_temperature_num(this->kin_ntemp).set_step_num(this->kin_nsteps); + try + { + sak.cal_all_temperature(); + } + catch (File_Finish_Metropolis()) + { + return; + } + + timer::tick("Metropolis", "move_various_t"); + return; } \ No newline at end of file diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Metropolis.h b/tools/01_NAO_generation/SIAB/src_spillage/Metropolis.h index bfad62267e0..2aafdf1cc30 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Metropolis.h +++ b/tools/01_NAO_generation/SIAB/src_spillage/Metropolis.h @@ -1,108 +1,128 @@ #ifndef METROPOLIS_H #define METROPOLIS_H -#include "SpillageStep.h" -#include "Psi_Second.h" // use to calculaet oscillation #include "Metropolis_Info.h" +#include "Psi_Second.h" // use to calculaet oscillation +#include "SpillageStep.h" class Metropolis { - public: - Metropolis(); - ~Metropolis(); - - void init(void); - - void move_various_temperature(SpillageStep &SS); - - void move_one_step(const int &itemp, const int &istep, SpillageStep &Level, const int &nsteps_now); - - bool accept(const double &v_old, const double &v_new, const double &temperature_in); - - void reset_temperature(const int &istep); - - - // (1) start temperature - const double& get_spi_ini_temp(void) const {return spi_ini_temp;} - const double& get_kin_ini_temp(void) const {return kin_ini_temp;} - - // (2) cooling rate for each new temperature - const double& get_spi_rate(void) const { return spi_rate;} - const double& get_kin_rate(void) const { return kin_rate;} - - // (3) number of temperature - const int& get_spi_ntemp(void) const { return spi_ntemp;} - const int& get_kin_ntemp(void) const { return kin_ntemp;} - - // (4) number of steps in each temperature - const int& get_spi_nsteps(void) const { return spi_nsteps;} - const int& get_kin_nsteps(void) const { return kin_nsteps;} - - - Psi_Second Psi2; - - private: - - // (1) - double spi_ini_temp; // initial temperature for spillage. - double kin_ini_temp; // initial temperature for kinetic energy. - - // (2) - double spi_rate; // decrease rate for spillage value - double kin_rate; // decrease rate for kinetic energy - - // (3) - int spi_ntemp; // number of temperatures (for spillage) - int kin_ntemp; // number of temperatures (for kinetical energy) - - // (4) - int spi_nsteps; // number of steps per temperature for spillage - int kin_nsteps; // number of steps per temperature for kinetical energy - - - double temperature; - double kappa; // mohan add 2009-10-11 - double delta_kappa; // mohan add 2009-10-31 - int try_kappa; - - double min_spillage; // mohan add 2009-10-15 - double *kin_last; // mohan add 2009-10-15 - double *kin_temp; // mohan add 2009-10-15, temperature for each orbital. - int nchi; // mohan add 2009-10-15 - int output_each_istep; // mohan add 2009-10-16 - - double prob_new; - double prob_old; - - public: - int test; - int random_c4number; - - void preprocess_all_temperature(SpillageStep &Level); - void cout_1(size_t itemperature); - void file_finish_metropolis(size_t istep); - void update_t(void); - void change_accept_rate(size_t istep, SpillageStep &Level); - void reset_high_kinetic_energy(size_t itemperature,size_t istep); - void init_spi_states(size_t itemperature); - void cout_2(size_t itemperature, SpillageStep &Level); - void init_kin_states(); - void init_temperature_kinetic(size_t itemperature, SpillageStep &Level); - void ofs_1(size_t istep, SpillageStep &Level); - void norm_C4(size_t istep, SpillageStep &Level, size_t nsteps_now); - void small_jql(); - void trial_c4(); - double cal_spillage(SpillageStep &Level); - double cal_kinetic_energy(SpillageStep &Level); - void accept_process(SpillageStep &Level); - void reject_process(); - - Metropolis_Info info; - - private: - class Simulated_Annealing_Orbital; - class Simulated_Annealing_Orbital_Spillage; - class Simulated_Annealing_Orbital_Kinetic; + public: + Metropolis(); + ~Metropolis(); + + void init(void); + + void move_various_temperature(SpillageStep& SS); + + void move_one_step(const int& itemp, const int& istep, SpillageStep& Level, const int& nsteps_now); + + bool accept(const double& v_old, const double& v_new, const double& temperature_in); + + void reset_temperature(const int& istep); + + // (1) start temperature + const double& get_spi_ini_temp(void) const + { + return spi_ini_temp; + } + const double& get_kin_ini_temp(void) const + { + return kin_ini_temp; + } + + // (2) cooling rate for each new temperature + const double& get_spi_rate(void) const + { + return spi_rate; + } + const double& get_kin_rate(void) const + { + return kin_rate; + } + + // (3) number of temperature + const int& get_spi_ntemp(void) const + { + return spi_ntemp; + } + const int& get_kin_ntemp(void) const + { + return kin_ntemp; + } + + // (4) number of steps in each temperature + const int& get_spi_nsteps(void) const + { + return spi_nsteps; + } + const int& get_kin_nsteps(void) const + { + return kin_nsteps; + } + + Psi_Second Psi2; + + private: + // (1) + double spi_ini_temp; // initial temperature for spillage. + double kin_ini_temp; // initial temperature for kinetic energy. + + // (2) + double spi_rate; // decrease rate for spillage value + double kin_rate; // decrease rate for kinetic energy + + // (3) + int spi_ntemp; // number of temperatures (for spillage) + int kin_ntemp; // number of temperatures (for kinetical energy) + + // (4) + int spi_nsteps; // number of steps per temperature for spillage + int kin_nsteps; // number of steps per temperature for kinetical energy + + double temperature; + double kappa; // mohan add 2009-10-11 + double delta_kappa; // mohan add 2009-10-31 + int try_kappa; + + double min_spillage; // mohan add 2009-10-15 + double* kin_last; // mohan add 2009-10-15 + double* kin_temp; // mohan add 2009-10-15, temperature for each orbital. + int nchi; // mohan add 2009-10-15 + int output_each_istep; // mohan add 2009-10-16 + + double prob_new; + double prob_old; + + public: + int test; + int random_c4number; + + void preprocess_all_temperature(SpillageStep& Level); + void cout_1(size_t itemperature); + void file_finish_metropolis(size_t istep); + void update_t(void); + void change_accept_rate(size_t istep, SpillageStep& Level); + void reset_high_kinetic_energy(size_t itemperature, size_t istep); + void init_spi_states(size_t itemperature); + void cout_2(size_t itemperature, SpillageStep& Level); + void init_kin_states(); + void init_temperature_kinetic(size_t itemperature, SpillageStep& Level); + void ofs_1(size_t istep, SpillageStep& Level); + void norm_C4(size_t istep, SpillageStep& Level, size_t nsteps_now); + void small_jql(); + void trial_c4(); + double cal_spillage(SpillageStep& Level); + double cal_kinetic_energy(SpillageStep& Level); + void accept_process(SpillageStep& Level); + void reject_process(); + + Metropolis_Info info; + + private: + class Simulated_Annealing_Orbital; + class Simulated_Annealing_Orbital_Spillage; + class Simulated_Annealing_Orbital_Kinetic; }; #endif diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Metropolis_2.cpp b/tools/01_NAO_generation/SIAB/src_spillage/Metropolis_2.cpp index 05fc8a3f37c..bf3235da4c2 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Metropolis_2.cpp +++ b/tools/01_NAO_generation/SIAB/src_spillage/Metropolis_2.cpp @@ -1,502 +1,473 @@ +#include "Calculate_C4.h" #include "Metropolis.h" #include "tools.h" -#include "Calculate_C4.h" + #include #include -void Metropolis::preprocess_all_temperature(SpillageStep &Level) +void Metropolis::preprocess_all_temperature(SpillageStep& Level) { - this->nchi = Level.nchi; - - // firstly, make the initia kappa value. - // mohan add 2009-11-12 - // this 0.001 must be added to the second step - // to avoid the oscillation of first orbital. - this->kappa = 1.0 + this->delta_kappa; - - delete[] this->kin_last; - this->kin_last = new double[nchi]; - ZEROS(this->kin_last, nchi); - - delete[] this->kin_temp; - this->kin_temp = new double[nchi]; - ZEROS(this->kin_temp, nchi); + this->nchi = Level.nchi; + + // firstly, make the initia kappa value. + // mohan add 2009-11-12 + // this 0.001 must be added to the second step + // to avoid the oscillation of first orbital. + this->kappa = 1.0 + this->delta_kappa; + + delete[] this->kin_last; + this->kin_last = new double[nchi]; + ZEROS(this->kin_last, nchi); + + delete[] this->kin_temp; + this->kin_temp = new double[nchi]; + ZEROS(this->kin_temp, nchi); } void Metropolis::file_finish_metropolis(size_t istep) { - // P.S.: a method to quit this whole program. - // you only need to create a file name : FINISH - // another function: if there existep a file 'FINISH' - // in the current dir, the program quit at the first - // iteration. - // Here I choose a proper parameter: 500 steps, - // you can change on your own. - const int check_finish_file(500); - if( 0 == (istep%check_finish_file) ) - { - // mohan add 2009-08-28 - ifstream ifs("FINISH"); - if( ifs ) - { - throw File_Finish_Metropolis(); - } - } + // P.S.: a method to quit this whole program. + // you only need to create a file name : FINISH + // another function: if there existep a file 'FINISH' + // in the current dir, the program quit at the first + // iteration. + // Here I choose a proper parameter: 500 steps, + // you can change on your own. + const int check_finish_file(500); + if (0 == (istep % check_finish_file)) + { + // mohan add 2009-08-28 + ifstream ifs("FINISH"); + if (ifs) + { + throw File_Finish_Metropolis(); + } + } } void Metropolis::update_t(void) { - if(this->info.states == "Spi") - { - this->temperature *= this->spi_rate; -// cout << " New Temperature for spillage = " << temperature << endl; - } - else if(this->info.states == "Kin"||this->info.states=="Ecut") - { - // mohan add high 2009-11-05 - for(int ic=0; ic< this->nchi; ic++) - { - this->kin_temp[ic] *= this->kin_rate; - } - const double high = *max_element(this->kin_temp, this->kin_temp+this->nchi); - for(int ic=0; ic< this->nchi; ic++) - { - this->kin_temp[ic] = high; - ofs_running << setprecision(6) //<< setiosflags(ios::scientific) - << " Temperature for orbital " << ic+1 << " = " << kin_temp[ic] << endl; - } - - ofs_running << " new kappa=" << kappa << endl; - } - return; + if (this->info.states == "Spi") + { + this->temperature *= this->spi_rate; + // cout << " New Temperature for spillage = " << temperature << endl; + } + else if (this->info.states == "Kin" || this->info.states == "Ecut") + { + // mohan add high 2009-11-05 + for (int ic = 0; ic < this->nchi; ic++) + { + this->kin_temp[ic] *= this->kin_rate; + } + const double high = *max_element(this->kin_temp, this->kin_temp + this->nchi); + for (int ic = 0; ic < this->nchi; ic++) + { + this->kin_temp[ic] = high; + ofs_running << setprecision(6) //<< setiosflags(ios::scientific) + << " Temperature for orbital " << ic + 1 << " = " << kin_temp[ic] << endl; + } + + ofs_running << " new kappa=" << kappa << endl; + } + return; } -void Metropolis::change_accept_rate(size_t istep, SpillageStep &Level) +void Metropolis::change_accept_rate(size_t istep, SpillageStep& Level) { - // Here is the part about changing Metropolis algorithm's accept rate. - // we want to calculate the accept rate every 'ACCEPTANCE_STEPS' steps. - // but we don't want to calculate at istep=0; - if(istep%ACCEPTANCE_STEPS == 0) - { - if(istep>0) - { - for(int ic = 0; ic < this->nchi; ic++) - { - const int T = Level.wayc[ic].type; - const int L = Level.wayc[ic].L; - const int N = Level.wayc[ic].N; - for(int ie2=0; ie2 ACCEPTANCE_HIGH) - { - input.Coef.accept_rate(T,L,N,ie2) *= 2.0; // You can DIY 0.5 - } - else if(rate < ACCEPTANCE_LOW) - { - input.Coef.accept_rate(T,L,N,ie2) *= 0.5; // You can DIY 2.0 - } - // output information - //if(ie2%5==0) cout << "\n"; - //cout << setw(8) << rate*100 << "%"; - //cout << setw(8) << input.Coef.accept_rate(T,L,N,ie2); - } - } - } - input.Coef.accept_number.zero_out(); - } // end accept_rate -} - -void Metropolis::reset_high_kinetic_energy(size_t itemperature,size_t istep) + // Here is the part about changing Metropolis algorithm's accept rate. + // we want to calculate the accept rate every 'ACCEPTANCE_STEPS' steps. + // but we don't want to calculate at istep=0; + if (istep % ACCEPTANCE_STEPS == 0) + { + if (istep > 0) + { + for (int ic = 0; ic < this->nchi; ic++) + { + const int T = Level.wayc[ic].type; + const int L = Level.wayc[ic].L; + const int N = Level.wayc[ic].N; + for (int ie2 = 0; ie2 < Level.ne; ie2++) + { + assert(input.Coef.accept_number(T, L, N, ie2) <= ACCEPTANCE_STEPS); + + const double rate = (double)input.Coef.accept_number(T, L, N, ie2) / (double)ACCEPTANCE_STEPS; + + if (rate > ACCEPTANCE_HIGH) + { + input.Coef.accept_rate(T, L, N, ie2) *= 2.0; // You can DIY 0.5 + } + else if (rate < ACCEPTANCE_LOW) + { + input.Coef.accept_rate(T, L, N, ie2) *= 0.5; // You can DIY 2.0 + } + // output information + // if(ie2%5==0) cout << "\n"; + // cout << setw(8) << rate*100 << "%"; + // cout << setw(8) << input.Coef.accept_rate(T,L,N,ie2); + } + } + } + input.Coef.accept_number.zero_out(); + } // end accept_rate +} + +void Metropolis::reset_high_kinetic_energy(size_t itemperature, size_t istep) { - // because the initia temperature may be too high - // which caused by the high kinetic energy, so - // we need to reset the temperature. - if( 0 == itemperature && 50 == istep ) - { - for(int ic=0; icnchi; ic++) - { - this->kin_temp[ic] = this->kin_last[ic] * this->kin_ini_temp; - } - const double high = *max_element(this->kin_temp, this->kin_temp+this->nchi); - for(int ic=0; icnchi; ic++) - { - this->kin_temp[ic] = high; - ofs_running << " Notice: Change temperature of orbital " << ic+1 << " to " << kin_temp[ic] << endl; - } - } + // because the initia temperature may be too high + // which caused by the high kinetic energy, so + // we need to reset the temperature. + if (0 == itemperature && 50 == istep) + { + for (int ic = 0; ic < this->nchi; ic++) + { + this->kin_temp[ic] = this->kin_last[ic] * this->kin_ini_temp; + } + const double high = *max_element(this->kin_temp, this->kin_temp + this->nchi); + for (int ic = 0; ic < this->nchi; ic++) + { + this->kin_temp[ic] = high; + ofs_running << " Notice: Change temperature of orbital " << ic + 1 << " to " << kin_temp[ic] << endl; + } + } } void Metropolis::init_spi_states(size_t itemperature) { - ofs_running << "\n" << " ---> SpillageTemp " << setw(3) << itemperature+1 << endl; - ofs_running << " Temperature " << setw(6) << this->temperature << endl; - this->info.states = "Spi"; - ofs_running << setw(5) << "STEPS" - << setw(20) << "SPILLAGE" - << setw(10) << "UPDATES" << endl; - - if(0==itemperature) - { - cout << setw(5) << "STEP" << setw(10) << "TEMP" - << setw(20) << "SPILLAGE" << endl; - } + ofs_running << "\n" << " ---> SpillageTemp " << setw(3) << itemperature + 1 << endl; + ofs_running << " Temperature " << setw(6) << this->temperature << endl; + this->info.states = "Spi"; + ofs_running << setw(5) << "STEPS" << setw(20) << "SPILLAGE" << setw(10) << "UPDATES" << endl; + + if (0 == itemperature) + { + cout << setw(5) << "STEP" << setw(10) << "TEMP" << setw(20) << "SPILLAGE" << endl; + } } void Metropolis::init_kin_states() { - // if the kinetic energy is not small enough, we try new kappa again. - //try_new_kappa_again: - // Please DIY, mohan 2010-04-14 - if(OPTIMIZE_METHOD==1) - { - this->info.states="Kin"; - } - else if(OPTIMIZE_METHOD==2) - { - this->info.states = "Ecut"; - } - else - { - WARNING_QUIT("move_various_temperature","Check OPTIMIZE_METHOD"); - } + // if the kinetic energy is not small enough, we try new kappa again. + // try_new_kappa_again: + // Please DIY, mohan 2010-04-14 + if (OPTIMIZE_METHOD == 1) + { + this->info.states = "Kin"; + } + else if (OPTIMIZE_METHOD == 2) + { + this->info.states = "Ecut"; + } + else + { + WARNING_QUIT("move_various_temperature", "Check OPTIMIZE_METHOD"); + } } -void Metropolis::init_temperature_kinetic(size_t itemperature, SpillageStep &Level) +void Metropolis::init_temperature_kinetic(size_t itemperature, SpillageStep& Level) { - // set the initial temperature and kinetic energy value. - if(0==itemperature) - { - this->min_spillage = input.SV.cal_defined_value(0); // save the spillage value. - ofs_running << "\n\n Spillage at final temperature " << min_spillage << endl; - - for(int ic=0; icnchi; ic++) - { - const int T = Level.wayc[ic].type; - const int L = Level.wayc[ic].L; - const int N = Level.wayc[ic].N; - - double *c4_last = new double[Level.ne]; - for(int ie2=0; ie2kin_last[ic] = 0.0; - - if(this->info.states=="Kin") - { - for(int ie2=0; ie2kin_last[ic] += pow( c4_last[ie2]*(ie2+1)*PI/RCUT,2 ) - *this->Psi2.jjnorm(L,ie2); - } - } - else if(this->info.states=="Ecut") - { - this->kin_last[ic] = this->Psi2.get_ecut( c4_last, L ); - } - - // set the start temperature of kinetic part of Metropolis algorithms. - this->kin_temp[ic] = this->kin_last[ic] * this->kin_ini_temp; - - ofs_running << " Orbital " << ic+1 << " Initial E_kin " << setw(15) << this->kin_last[ic] << endl; - - delete[] c4_last; - } - } + // set the initial temperature and kinetic energy value. + if (0 == itemperature) + { + this->min_spillage = input.SV.cal_defined_value(0); // save the spillage value. + ofs_running << "\n\n Spillage at final temperature " << min_spillage << endl; + + for (int ic = 0; ic < this->nchi; ic++) + { + const int T = Level.wayc[ic].type; + const int L = Level.wayc[ic].L; + const int N = Level.wayc[ic].N; + + double* c4_last = new double[Level.ne]; + for (int ie2 = 0; ie2 < Level.ne; ie2++) + { + c4_last[ie2] = input.Coef.C4_old(T, L, N, ie2); + } + + this->kin_last[ic] = 0.0; + + if (this->info.states == "Kin") + { + for (int ie2 = 0; ie2 < Level.ne; ie2++) + { + this->kin_last[ic] += pow(c4_last[ie2] * (ie2 + 1) * PI / RCUT, 2) * this->Psi2.jjnorm(L, ie2); + } + } + else if (this->info.states == "Ecut") + { + this->kin_last[ic] = this->Psi2.get_ecut(c4_last, L); + } + + // set the start temperature of kinetic part of Metropolis algorithms. + this->kin_temp[ic] = this->kin_last[ic] * this->kin_ini_temp; + + ofs_running << " Orbital " << ic + 1 << " Initial E_kin " << setw(15) << this->kin_last[ic] << endl; + + delete[] c4_last; + } + } } -void Metropolis::norm_C4(size_t istep, SpillageStep &Level, size_t nsteps_now) +void Metropolis::norm_C4(size_t istep, SpillageStep& Level, size_t nsteps_now) { - // because we choose c4 randomly between -0.1~0.1 - // so we need to normalize the set of c4 after a few - // steps, here I choose a proper number: 50 - // you can change on you own, but please notice it can't - // be too small, because we need to update all Q and S - // matrix after this normalization, which may affect the - // speed of the program. You can't choose too large, - // because it makes the update not efficient. - if ( (istep % 50 == 0 && istep !=0 ) || istep == nsteps_now-1) // Peize Lin test -// if ( (istep % 1 == 0 && istep !=0 ) || istep == nsteps_now-1) // Peize Lin test - { - // use to describe the real space mesh. - // don't need to choose too accuray, because - // it may affect the speed of running. - // you can change on your own, too. - double dr = 0.01; // a.u. - - // norm ic for all radial wave functions. - // after this call, all c4 will change/update. - // but the spillage will not change. - // (because we just normalize the wavefunctions). - - Calculate_C4::norm_ic( - input.Coef.C4, - input.Coef.accept_rate, // Peize Lin update 2015-12-24 - input.Coef.ntype, - input.Coef.lmax, - input.Coef.nmax, - input.Coef.enumber, - input.Coef.tolerence, - input.Coef.rcut, - dr); - - // get new set of C4, so we copy the set of C4 - // to C4_old. This will also not affect the - // spillage value. - input.Coef.copy_c4_to_old(); - - // for each structure, we need to update their - // Q matrix and S matrix, which may be time - // consuming. - for (int istr(0); istr < STRNUM; ++istr) - { - Level.init_QS_matrix(istr); - } - } + // because we choose c4 randomly between -0.1~0.1 + // so we need to normalize the set of c4 after a few + // steps, here I choose a proper number: 50 + // you can change on you own, but please notice it can't + // be too small, because we need to update all Q and S + // matrix after this normalization, which may affect the + // speed of the program. You can't choose too large, + // because it makes the update not efficient. + if ((istep % 50 == 0 && istep != 0) || istep == nsteps_now - 1) // Peize Lin test + // if ( (istep % 1 == 0 && istep !=0 ) || istep == nsteps_now-1) // Peize Lin test + { + // use to describe the real space mesh. + // don't need to choose too accuray, because + // it may affect the speed of running. + // you can change on your own, too. + double dr = 0.01; // a.u. + + // norm ic for all radial wave functions. + // after this call, all c4 will change/update. + // but the spillage will not change. + // (because we just normalize the wavefunctions). + + Calculate_C4::norm_ic(input.Coef.C4, + input.Coef.accept_rate, // Peize Lin update 2015-12-24 + input.Coef.ntype, + input.Coef.lmax, + input.Coef.nmax, + input.Coef.enumber, + input.Coef.tolerence, + input.Coef.rcut, + dr); + + // get new set of C4, so we copy the set of C4 + // to C4_old. This will also not affect the + // spillage value. + input.Coef.copy_c4_to_old(); + + // for each structure, we need to update their + // Q matrix and S matrix, which may be time + // consuming. + for (int istr(0); istr < STRNUM; ++istr) + { + Level.init_QS_matrix(istr); + } + } } void Metropolis::small_jql() { - // mohan add 2009-08-26 - // this judgement means, - // If users don't want to use Jlq(ie) - // which is too oscillation(large), they - // can use only small part of Jlq. - // mohan add BLOCK_NE_MIN 2009-08-27 - if( (this->info.ieinfo.ie>=BLOCK_NE)) - throw Small_Jlq(); + // mohan add 2009-08-26 + // this judgement means, + // If users don't want to use Jlq(ie) + // which is too oscillation(large), they + // can use only small part of Jlq. + // mohan add BLOCK_NE_MIN 2009-08-27 + if ((this->info.ie < BLOCK_NE_MIN) || (this->info.ie >= BLOCK_NE)) + throw Small_Jlq(); } void Metropolis::trial_c4() { - //======================================================== - // if we want to move step, we need the set of C4 from CS - // in single zeta, the C4 can be either read in or random, - // in double zeta, the C4 should start from random, in - // most cases. - //======================================================== - const int ie = this->info.ie; - const int T = this->info.T; - const int L = this->info.L; - const int N = this->info.N; - // trial a new value of c4(T,L,N). - input.Coef.trial_c4( T, L, N, ie); + //======================================================== + // if we want to move step, we need the set of C4 from CS + // in single zeta, the C4 can be either read in or random, + // in double zeta, the C4 should start from random, in + // most cases. + //======================================================== + const int ie = this->info.ie; + const int T = this->info.T; + const int L = this->info.L; + const int N = this->info.N; + // trial a new value of c4(T,L,N). + input.Coef.trial_c4(T, L, N, ie); } -double Metropolis::cal_spillage(SpillageStep &Level) +double Metropolis::cal_spillage(SpillageStep& Level) { - const int ie = this->info.ie; - const int ic = this->info.ic; - for ( size_t istr = 0; istr < STRNUM; ++istr) - { - input.SV.value[istr] = Level.get_spillage(istr, ic, ie); - } - return input.SV.cal_defined_value(1); + const int ie = this->info.ie; + const int ic = this->info.ic; + for (size_t istr = 0; istr < STRNUM; ++istr) + { + input.SV.value[istr] = Level.get_spillage(istr, ic, ie); + } + return input.SV.cal_defined_value(1); } -double Metropolis::cal_kinetic_energy(SpillageStep &Level) +double Metropolis::cal_kinetic_energy(SpillageStep& Level) { - double *c4tmp = new double[Level.ne]; - - const int T = this->info.T; - const int L = this->info.L; - const int N = this->info.N; - for(int ie2=0; ie2info.states=="Kin") - { - // get the normalized parameters: c4 - this->Psi2.norm_c4( c4tmp, L ); - - // and then calculate the new set of kinetic energy. - for(int ie2=0; ie2Psi2.jjnorm(L,ie2);//In Rydberg Unit. - } - } - else if(this->info.states=="Ecut") - { - value2 = this->Psi2.get_ecut( c4tmp, L ); - } - delete[] c4tmp; - return value2; + double* c4tmp = new double[Level.ne]; + + const int T = this->info.T; + const int L = this->info.L; + const int N = this->info.N; + for (int ie2 = 0; ie2 < Level.ne; ie2++) + { + c4tmp[ie2] = input.Coef.C4(T, L, N, ie2); + } + + double value2 = 0.0; + if (this->info.states == "Kin") + { + // get the normalized parameters: c4 + this->Psi2.norm_c4(c4tmp, L); + + // and then calculate the new set of kinetic energy. + for (int ie2 = 0; ie2 < Level.ne; ie2++) + { + value2 += pow(c4tmp[ie2] * (ie2 + 1) * PI / RCUT, 2) * this->Psi2.jjnorm(L, ie2); // In Rydberg Unit. + } + } + else if (this->info.states == "Ecut") + { + value2 = this->Psi2.get_ecut(c4tmp, L); + } + delete[] c4tmp; + return value2; } -void Metropolis::accept_process(SpillageStep &Level) +void Metropolis::accept_process(SpillageStep& Level) { - ++this->info.update_number; - - for (int istr = 0; istr< STRNUM; ++istr) - { - Level.updateQS( istr ); // (1) update Soverlap and Qoverlap for each structure, - // S matrix and Q matrix has changed due to the change of c4 - } - - input.SV.update_value();// (2) update spillage value - - for(int istr=0; istrinfo.update_number << endl; - - // output each structure's spillage value. - // for(int is=0; is< input.str_num; is++) - // { - // cout << setw(10) << input.SV.value_old[is]*100 << "%"; - // } - } - else if(this->info.states=="Kin"||this->info.states=="Ecut") - { - //<< " ----> " << states << "Temp " << setw(15) - //<< istep+1 << "\n" ; -// << " SpiVal" -// << setprecision(6) -// << setiosflags(ios::fixed) -// << setw(15) << input.SV.cal_defined_value(0); - for(int ic=0; ic< this->nchi; ic++) - { - ofs_running << setw(4) << istep+1 << Level.wayc[ic].spd; - - ofs_running << setprecision(3) - << setiosflags(ios::scientific) - << setw(10)<< this->kin_temp[ic] - << resetiosflags(ios::scientific); - - ofs_running << setprecision(10) - << setiosflags(ios::fixed) - << setw(20)<< this->kin_last[ic] << endl - << resetiosflags(ios::fixed); - } - } - } + if (0 == (istep % this->output_each_istep)) + { + if (this->info.states == "Spi") + { + ofs_running << setw(5) << istep + 1 << setprecision(10) << setiosflags(ios::fixed) << setw(20) + << input.SV.cal_defined_value(0) << setw(10) << this->info.update_number << endl; + + // output each structure's spillage value. + // for(int is=0; is< input.str_num; is++) + // { + // cout << setw(10) << input.SV.value_old[is]*100 << "%"; + // } + } + else if (this->info.states == "Kin" || this->info.states == "Ecut") + { + //<< " ----> " << states << "Temp " << setw(15) + //<< istep+1 << "\n" ; + // << " SpiVal" + // << setprecision(6) + // << setiosflags(ios::fixed) + // << setw(15) << input.SV.cal_defined_value(0); + for (int ic = 0; ic < this->nchi; ic++) + { + ofs_running << setw(4) << istep + 1 << Level.wayc[ic].spd; + + ofs_running << setprecision(3) << setiosflags(ios::scientific) << setw(10) << this->kin_temp[ic] + << resetiosflags(ios::scientific); + + ofs_running << setprecision(10) << setiosflags(ios::fixed) << setw(20) << this->kin_last[ic] << endl + << resetiosflags(ios::fixed); + } + } + } } void Metropolis::cout_1(size_t itemperature) { - cout << setw(5) << itemperature+1 - << setiosflags(ios::showpoint) - << setiosflags(ios::scientific) - << setprecision(3) - << setw(10) - << this->temperature - << resetiosflags(ios::scientific) - << setiosflags(ios::fixed) - << setprecision(10) - << setw(19) - << input.SV.cal_defined_value(0)*100 << "%" - << resetiosflags(ios::fixed) - << resetiosflags(ios::showpoint) - << endl; + cout << setw(5) << itemperature + 1 << setiosflags(ios::showpoint) << setiosflags(ios::scientific) + << setprecision(3) << setw(10) << this->temperature << resetiosflags(ios::scientific) + << setiosflags(ios::fixed) << setprecision(10) << setw(19) << input.SV.cal_defined_value(0) * 100 << "%" + << resetiosflags(ios::fixed) << resetiosflags(ios::showpoint) << endl; } -void Metropolis::cout_2(size_t itemperature, SpillageStep &Level) +void Metropolis::cout_2(size_t itemperature, SpillageStep& Level) { - //------------------------------------------- - // it's about minimizing kinetic energy. - //------------------------------------------- - for(int ic=0; icnchi; ic++) - { - cout << setw(5) << Level.wayc[ic].spd; - cout << setiosflags(ios::showpoint) - << setiosflags(ios::scientific) - << setprecision(3) - << setw(10)<< this->kin_temp[ic] - << resetiosflags(ios::scientific); - cout << setprecision(10) - << setiosflags(ios::fixed) - << setw(20)<< this->kin_last[ic] - << resetiosflags(ios::fixed) - << resetiosflags(ios::showpoint) - << endl; - } - - for(int ic=0; icnchi; ic++) - { - if(this->kin_last[ic] > KINETIC_MAX) - { - cout << "\n kin_last[" << ic << "]=" << kin_last[ic] << " KINETIC_MAX = " << KINETIC_MAX << endl; - cout << " There maybe 2 reasons, first: The starting temperature for kinetical energy is too small." << endl; - cout << " Second: The kappa is too small." << endl; - //WARNING_QUIT("Minimizing Kinetic Energy","Kinetical Energy Too Large!"); - cout << "Kinetical Energy Too Large!"<kappa += this->delta_kappa; - //goto try_new_kappa_again; - } + //------------------------------------------- + // it's about minimizing kinetic energy. + //------------------------------------------- + for (int ic = 0; ic < this->nchi; ic++) + { + cout << setw(5) << Level.wayc[ic].spd; + cout << setiosflags(ios::showpoint) << setiosflags(ios::scientific) << setprecision(3) << setw(10) + << this->kin_temp[ic] << resetiosflags(ios::scientific); + cout << setprecision(10) << setiosflags(ios::fixed) << setw(20) << this->kin_last[ic] + << resetiosflags(ios::fixed) << resetiosflags(ios::showpoint) << endl; + } + + for (int ic = 0; ic < this->nchi; ic++) + { + if (this->kin_last[ic] > KINETIC_MAX) + { + cout << "\n kin_last[" << ic << "]=" << kin_last[ic] << " KINETIC_MAX = " << KINETIC_MAX << endl; + cout << " There maybe 2 reasons, first: The starting temperature for kinetical energy is too small." + << endl; + cout << " Second: The kappa is too small." << endl; + // WARNING_QUIT("Minimizing Kinetic Energy","Kinetical Energy Too Large!"); + cout << "Kinetical Energy Too Large!" << endl; + } + // this->kappa += this->delta_kappa; + // goto try_new_kappa_again; + } } - diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Metropolis_Info.h b/tools/01_NAO_generation/SIAB/src_spillage/Metropolis_Info.h index 0caddd3d379..5c69e8ffdbb 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Metropolis_Info.h +++ b/tools/01_NAO_generation/SIAB/src_spillage/Metropolis_Info.h @@ -3,28 +3,44 @@ struct Metropolis_Info { - int ie; - int ic; - int T; - int L; - int N; - int update_number; - string states; // mohan add 2009-10-15 + int ie; + int ic; + int T; + int L; + int N; + int update_number; + string states; // mohan add 2009-10-15 - Metropolis_Info &set_ie(int ie_in){ ie = ie_in; return *this;} - Metropolis_Info &set_ic(int ic_in){ ic = ic_in; return *this;} - Metropolis_Info &set_update_number(int update_number_in){ update_number = update_number_in; return *this;} + Metropolis_Info& set_ie(int ie_in) + { + ie = ie_in; + return *this; + } + Metropolis_Info& set_ic(int ic_in) + { + ic = ic_in; + return *this; + } + Metropolis_Info& set_update_number(int update_number_in) + { + update_number = update_number_in; + return *this; + } - void unpack_ic(SpillageStep &Level) - { - T = Level.wayc[ic].type; - L = Level.wayc[ic].L; - N = Level.wayc[ic].N; - } + void unpack_ic(SpillageStep& Level) + { + T = Level.wayc[ic].type; + L = Level.wayc[ic].L; + N = Level.wayc[ic].N; + } }; -class File_Finish_Metropolis: public exception{}; +class File_Finish_Metropolis : public exception +{ +}; -class Small_Jlq: public exception{}; +class Small_Jlq : public exception +{ +}; #endif \ No newline at end of file diff --git a/tools/01_NAO_generation/SIAB/src_spillage/MultiZeta.cpp b/tools/01_NAO_generation/SIAB/src_spillage/MultiZeta.cpp index 17f4088e35d..5e8a5b1cb05 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/MultiZeta.cpp +++ b/tools/01_NAO_generation/SIAB/src_spillage/MultiZeta.cpp @@ -1,10 +1,11 @@ #include "MultiZeta.h" -#include "SpillageStep.h" -#include "Out_Orbital.h" + #include "../src_parallel/parallel_common.h" +#include "Out_Orbital.h" +#include "SpillageStep.h" #include "tools.h" -//Peize Lin test +// Peize Lin test #include "../src_tools/complexmatrix_inline.h" #include "../src_unittest/src_tools/common_test.h" @@ -24,51 +25,50 @@ MultiZeta::~MultiZeta() // be called in main.cpp program. void MultiZeta::init(void) { - TITLE(ofs_running,"MultiZeta","init"); + TITLE(ofs_running, "MultiZeta", "init"); //============================= // nlevel : optimization levels //============================= - for (int is=0; isilevel = is+1;// 0 means level 1. + this->ilevel = is + 1; // 0 means level 1. // mohan add 2009-09-25 - //input.Coef.accumulate_num_zero(); + // input.Coef.accumulate_num_zero(); // if start temperature is different for each level // ( I don't think this is a very good idea for high level orbitals!! // mohan note:2009-8-20, it's better to make temperature from "high enough" to // "low enough" for each level) - metro.reset_temperature( is ); + metro.reset_temperature(is); this->Level[is].ilevel = is; this->Level[is].set_nchi(); this->Level[is].set_nwfc(); - this->Level[is].set_iw00();// to get Sq1q2 - this->Level[is].allocate_data( 0 ); // to get Qin - + this->Level[is].set_iw00(); // to get Sq1q2 + this->Level[is].allocate_data(0); // to get Qin // if is=0, means it's the first time, spillage0 = 1.0; - if (is==0) + if (is == 0) { - for (int istr=0; istrLevel[is].data[istr].spillage0 = 1.0; } } - else if (is>0) + else if (is > 0) { - for (int istr=0; istrLevel[is].data[istr].spillage0 = input.SV.value_old[istr]; } @@ -78,16 +78,16 @@ void MultiZeta::init(void) // (2): Orthogonal to get new Bloch space! // This subroutine is very important. //**************************************** - Orth.start( this->Level[is-1], this->Level[is] ); + Orth.start(this->Level[is - 1], this->Level[is]); } - if(USEPW) - { - PW.allocate_psi1d(is); - PW.allocate_psi3d(is); - } + if (USEPW) + { + PW.allocate_psi1d(is); + PW.allocate_psi3d(is); + } - for (int istr=0; istrLevel[is].data[istr].initQS(); this->Level[is].init_QS_matrix(istr); @@ -102,75 +102,75 @@ void MultiZeta::init(void) //******************************************************************* // main part: Use Metropolis algorithm, move different temperature, //******************************************************************* - this->metro.move_various_temperature( this->Level[is]); - + this->metro.move_various_temperature(this->Level[is]); + // save spillage value for each level input.SV.save_level(is); this->saveC(); -// Peize Lin test -for ( size_t istr = 0; istr < STRNUM; ++istr) -{ - Step_Data* Stru = &this->Level[is].data[istr]; - for (int ik = 0; ik < Stru->nks ; ik++) - { - //cout<<"ik:\t"<Strial[ik]); - //cout_matrix( Stru->inverse_S[ik].get_inverse()); - for( size_t i=0; iStrial[ik].nr; ++i) - for( size_t j=0; jStrial[ik](i,j) = conj(Stru->Strial[ik](j,i)); - ComplexMatrix inverse_tmp(Stru->inverse_S[ik].get_inverse()); - for( size_t i=0; iStrial[ik] * inverse_tmp); - } -} + // Peize Lin test + for (size_t istr = 0; istr < STRNUM; ++istr) + { + Step_Data* Stru = &this->Level[is].data[istr]; + for (int ik = 0; ik < Stru->nks; ik++) + { + // cout<<"ik:\t"<Strial[ik]); + // cout_matrix( Stru->inverse_S[ik].get_inverse()); + for (size_t i = 0; i < Stru->Strial[ik].nr; ++i) + for (size_t j = 0; j < i; ++j) + Stru->Strial[ik](i, j) = conj(Stru->Strial[ik](j, i)); + ComplexMatrix inverse_tmp(Stru->inverse_S[ik].get_inverse()); + for (size_t i = 0; i < inverse_tmp.nr; ++i) + for (size_t j = 0; j < i; ++j) + inverse_tmp(i, j) = conj(inverse_tmp(j, i)); + // cout_matrix( Stru->Strial[ik] * inverse_tmp); + } + } } - return;; + return; + ; } - - // read in paramters from BLOCK , called by MultiZeta::init(); void MultiZeta::set_levels(void) { - TITLE(ofs_running,"MultiZeta","set_levels"); + TITLE(ofs_running, "MultiZeta", "set_levels"); - bool begin = false; + bool begin = false; - if(MY_RANK==0) - { - if ( SCAN_BEGIN(input.ifs, "") ) begin = true; - } + if (MY_RANK == 0) + { + if (SCAN_BEGIN(input.ifs, "")) + begin = true; + } #ifdef __MPI - Parallel_Common::bcast_bool(begin); + Parallel_Common::bcast_bool(begin); #endif // read in parameters in BLOCK : 'OPTIMIZE' - if ( begin ) + if (begin) { // (1) read in important parameter: nlevel. // nlevel means the number of orbital shells we need. - - if(MY_RANK==0) - { - READ_VALUE(input.ifs, this->nlevel); - assert(nlevel>-1); - assert(nlevel<100); - } - -#ifdef __MPI - Parallel_Common::bcast_int(nlevel); + + if (MY_RANK == 0) + { + READ_VALUE(input.ifs, this->nlevel); + assert(nlevel > -1); + assert(nlevel < 100); + } + +#ifdef __MPI + Parallel_Common::bcast_int(nlevel); #endif - + input.nlevel = this->nlevel; // mohan add 2009-06-09 -// cout << " nlevel = " << nlevel << endl; - OUT(ofs_running,"nlevel",nlevel); + // cout << " nlevel = " << nlevel << endl; + OUT(ofs_running, "nlevel", nlevel); // (2) init 'SpillageStep' class, each level has // different parameters. @@ -179,73 +179,73 @@ void MultiZeta::set_levels(void) string useless; - if(MY_RANK==0) - { - READ_VALUE(input.ifs, useless); - } + if (MY_RANK == 0) + { + READ_VALUE(input.ifs, useless); + } #ifdef __MPI - Parallel_Common::bcast_string(useless); + Parallel_Common::bcast_string(useless); #endif // (3) read in information concerning each level. - for (int i=0; i>>>> Level=" << i+1 ; + ofs_running << "\n >>>>> Level=" << i + 1; this->Level[i].set_info(NTYPE); - for (int it=0; it> Level[i].info[it].id //id for this type - >> Level[i].info[it].na //number of atoms. - >> Level[i].info[it].state// new or skip, mohan 2009-08-27 - >> Level[i].info[it].lmax;//lmax in this level - } + + if (MY_RANK == 0) + { + input.ifs >> Level[i].info[it].id // id for this type + >> Level[i].info[it].na // number of atoms. + >> Level[i].info[it].state // new or skip, mohan 2009-08-27 + >> Level[i].info[it].lmax; // lmax in this level + } #ifdef __MPI - Parallel_Common::bcast_int(Level[i].info[it].id); - Parallel_Common::bcast_int(Level[i].info[it].na); - Parallel_Common::bcast_string(Level[i].info[it].state); - Parallel_Common::bcast_int(Level[i].info[it].lmax); + Parallel_Common::bcast_int(Level[i].info[it].id); + Parallel_Common::bcast_int(Level[i].info[it].na); + Parallel_Common::bcast_string(Level[i].info[it].state); + Parallel_Common::bcast_int(Level[i].info[it].lmax); #endif - if (!RESTART && Level[i].info[it].state=="skip") + if (!RESTART && Level[i].info[it].state == "skip") { ofs_running << "\n Level = " << i; ofs_running << "\n it = " << it; - WARNING_QUIT("MultiZeta::set_levels","skip is not available if RESTART = false."); + WARNING_QUIT("MultiZeta::set_levels", "skip is not available if RESTART = false."); } - if ( Level[i].info[it].state!="new" && Level[i].info[it].state!="skip") + if (Level[i].info[it].state != "new" && Level[i].info[it].state != "skip") { ofs_running << "\n Level = " << i; ofs_running << "\n it = " << it; ofs_running << "\n state = " << Level[i].info[it].state << endl; - WARNING_QUIT("MultiZeta::set_levels","new or skip?"); + WARNING_QUIT("MultiZeta::set_levels", "new or skip?"); } ofs_running << "\n id=" << Level[i].info[it].id; ofs_running << " na=" << Level[i].info[it].na; ofs_running << " lmax=" << Level[i].info[it].lmax; - if(Level[i].info[it].id==0) - { - ofs_running << " The id of elements should be numbers!" << endl; - ofs_running << " check the id of elements." << endl; - WARNING_QUIT("MultiZeta::set_levels","id=0"); - } + if (Level[i].info[it].id == 0) + { + ofs_running << " The id of elements should be numbers!" << endl; + ofs_running << " check the id of elements." << endl; + WARNING_QUIT("MultiZeta::set_levels", "id=0"); + } // lmax=0(s), 1(p), 2(d), 3(f), 4(g) - if (Level[i].info[it].lmax>=5 || Level[i].info[it].lmax<0) + if (Level[i].info[it].lmax >= 5 || Level[i].info[it].lmax < 0) { ofs_running << "\n lmax=" << Level[i].info[it].lmax << endl; - WARNING_QUIT("MultiZeta::set_levels"," lmax out of range."); + WARNING_QUIT("MultiZeta::set_levels", " lmax out of range."); } // 'info' stands for "Type_Information" @@ -253,7 +253,7 @@ void MultiZeta::set_levels(void) // allocate Level[i].info[it].init(); - for (int L=0; L>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> @@ -266,14 +266,14 @@ void MultiZeta::set_levels(void) input.ifs >> combine; #ifdef __MPI - Parallel_Common::bcast_string(combine); + Parallel_Common::bcast_string(combine); #endif - - //ofs_running << "\n lenght = " << combine.length() << endl; + + // ofs_running << "\n lenght = " << combine.length() << endl; // firstly, we need to find the position of the label: '/' const int index1 = combine.find_first_of("/"); - //ofs_running << "\n index1 = " << index1 << endl; + // ofs_running << "\n index1 = " << index1 << endl; // if index1==-1, means the combine didn't contain // '/' label @@ -282,14 +282,14 @@ void MultiZeta::set_levels(void) // f stands for 'full' fain = "f"; // this 'atoi' command change string to int type. - nz = std::atoi( combine.c_str() ); + nz = std::atoi(combine.c_str()); } else // mean containing "/" { - nz = std::atoi( combine.substr(0,index1+1).c_str() ); + nz = std::atoi(combine.substr(0, index1 + 1).c_str()); // substrate the front int and remain the character. // the second parameter '1' means the position is 'index+2' - fain = combine.substr(index1+1, 1); + fain = combine.substr(index1 + 1, 1); } ofs_running << "\n L=" << L; @@ -299,7 +299,7 @@ void MultiZeta::set_levels(void) // assert number of zeta < 20, you can change on your own, // typically, we only need nz=2 or 3. - assert( nz < 20 ); + assert(nz < 20); Level[i].info[it].n[L] = nz; Level[i].info[it].fa[L] = fain; @@ -310,37 +310,36 @@ void MultiZeta::set_levels(void) // (hard to understand, right? Try to figure out!) // this parameter is very useful. // nbase provide a base to count radial wave functions. - if (i==0) + if (i == 0) { Level[i].info[it].nbase[L] = 0; } - else if (i>0) + else if (i > 0) { - for (int j=0; jlmax_type = new int[NTYPE]; - for (int it=0; itlmax_type[it] = -1; } // sum information from all levels. - NMAXUSED = 0; - for (int k=0; knlevel; k++) + NMAXUSED = 0; + for (int k = 0; k < this->nlevel; k++) { - for (int it=0; it< NTYPE; it++) + for (int it = 0; it < NTYPE; it++) { - this->lmax_type[it] = std::max( this->lmax_type[it], this->Level[k].info[it].lmax ); - LMAXUSED = std::max( LMAXUSED, this->lmax_type[it] ); - NMAXUSED = std::max( NMAXUSED, this->Level[k].info[it].nmax ); + this->lmax_type[it] = std::max(this->lmax_type[it], this->Level[k].info[it].lmax); + LMAXUSED = std::max(LMAXUSED, this->lmax_type[it]); + NMAXUSED = std::max(NMAXUSED, this->Level[k].info[it].nmax); } } // for test. - if (test==2) + if (test == 2) { - for (int it=0; it *sum = new std::complex[nbands]; - for (int iw=0; iw* sum = new std::complex[nbands]; + for (int iw = 0; iw < nwfc_all; iw++) + { + for (int ie = 0; ie < ne; ie++) { - for (int ie=0; ie first_part = std::complex(0,0); - if(USEPW) - { - const int Mmu = step1.wayd[mu].m; - const int Imu = step1.wayd[mu].i; - PW.calculate_Jlq(ik,iw,ie); - first_part = PW.calculate_Jlq_Phi(ik,mu); - - // if( norm(first_part) > 0.1 ) - // { - // cout << "\n iw=" << iw << " ie=" << ie << " mu=" << mu << " first_part=" << first_part; - // BLOCK_HERE("haha"); - // } - } - else - { - for (int iemu=0; iemu first_part = std::complex(0, 0); + if (USEPW) + { + const int Mmu = step1.wayd[mu].m; + const int Imu = step1.wayd[mu].i; + PW.calculate_Jlq(ik, iw, ie); + first_part = PW.calculate_Jlq_Phi(ik, mu); + + // if( norm(first_part) > 0.1 ) + // { + // cout << "\n iw=" << iw << " ie=" << ie << " mu=" << mu << " first_part=" << + // first_part; BLOCK_HERE("haha"); + // } + } + else + { + for (int iemu = 0; iemu < ne; iemu++) { - for (int ib=0; ib sum = std::complex(0,0); - for(int mu=0; mu first_part = std::complex(0,0); - - for(int iemu=0; iemu sum = std::complex(0,0); + for(int mu=0; mu first_part = std::complex(0,0); + + for(int iemu=0; iemuversion_information( ofs ); - this->INPUTs_information( ofs ); - this->spillage_information( ofs, mz.ilevel ); - this->metropolis_information( ofs); - this->c4_information( ofs ); - } - this->mkb_information( ofs); - - if(MY_RANK==0) - { - ofs.close(); - } - - return; + TITLE(ofs_running, "Out_Orbital", "write"); + + ofstream ofs; + + if (MY_RANK == 0) + { + ofs.open("ORBITAL_RESULTS.txt"); + this->version_information(ofs); + this->INPUTs_information(ofs); + this->spillage_information(ofs, mz.ilevel); + this->metropolis_information(ofs); + this->c4_information(ofs); + } + this->mkb_information(ofs); + + if (MY_RANK == 0) + { + ofs.close(); + } + + return; } -void Out_Orbital::version_information( ofstream &ofs ) +void Out_Orbital::version_information(ofstream& ofs) { ofs << "\n "; ofs << "\n AUTHOR : Mohan Chen"; ofs << "\n StartDate : 2009-4-01"; - ofs << "\n LastModify: 2012-6-27"; + ofs << "\n LastModify: 2012-6-27"; ofs << "\n LOCATION : LQCC, Hefei, China"; ofs << "\n EMAIL : mohan@mail.ustc.edu.cn"; ofs << "\n Description : Calculate the coefficients C4 of f(r) in Spherical Bessel Basis J(qr)."; @@ -45,199 +50,193 @@ void Out_Orbital::version_information( ofstream &ofs ) ofs << "\n"; } -void Out_Orbital::INPUTs_information( ofstream &ofs ) +void Out_Orbital::INPUTs_information(ofstream& ofs) { ofs << "\n\n"; ofs << "\n" << setw(20) << ECUT << " Energy cutoff(Hartree.)."; ofs << "\n" << setw(20) << RCUT << " rcut (a.u.)"; ofs << "\n" << setw(20) << NE << " eigenvalue number( sqrt(ecut*2)*rcut/PI )."; ofs << "\n" << setw(20) << TOLERENCE << " tolerence to calculate eigenvalue."; - ofs << "\n" << setw(20) << NTYPE << " Number of atom types."; - for(int it=0; it"; } -void Out_Orbital::spillage_information( ofstream &ofs, const int &ilevel) +void Out_Orbital::spillage_information(ofstream& ofs, const int& ilevel) { - ofs << "\n\n"; - ofs << "\n" << setw(20) << STRNUM << " kinds of structures."; - - // (1) mohan add average information 2010-04-11 - ofs << "\nAverage Spillage Value"; - for(int k=0; k"; - ofs << resetiosflags(ios::scientific); - return; + ofs << "\n\n"; + ofs << "\n" << setw(20) << STRNUM << " kinds of structures."; + + // (1) mohan add average information 2010-04-11 + ofs << "\nAverage Spillage Value"; + for (int k = 0; k < ilevel; k++) + { + double average = 0.0; + for (int i = 0; i < STRNUM; i++) + { + average += input.SV.value_each_level(i, k); + } + ofs << "\n " << setw(8) << k + 1 << setiosflags(ios::scientific) << setprecision(6) << setw(20) + << average / STRNUM; + } + + // (2) SV stands for 'spillage value' + // value_old containing the newest accepted spillage value. + for (int i = 0; i < STRNUM; i++) + { + ofs << "\nStructureIndex " << i + 1; + for (int k = 0; k < ilevel; k++) + { + ofs << "\n " << setw(8) << k + 1 << setiosflags(ios::scientific) << setprecision(6) << setw(20) + << input.SV.value_each_level(i, k); + } + } + ofs << "\n"; + ofs << resetiosflags(ios::scientific); + return; } -void Out_Orbital::metropolis_information( ofstream &ofs) +void Out_Orbital::metropolis_information(ofstream& ofs) { - ofs << "\n\n"; - ofs << "\n" << setw(20) << mz.metro.get_spi_ini_temp() << " Start temperature (Kelvin) for spillage minimization."; - ofs << "\n" << setw(20) << mz.metro.get_spi_rate() << " Decreasing rate of temperature."; - ofs << "\n" << setw(20) << mz.metro.get_spi_ntemp() << " Number of different temperature (for spillage)."; - ofs << "\n" << setw(20) << mz.metro.get_spi_nsteps() << " Number of steps for each temperature (for spillage)."; - - ofs << "\n" << setw(20) << mz.metro.get_kin_ini_temp() << " Start temperature (Kelvin) for kinetical energy minimization."; - ofs << "\n" << setw(20) << mz.metro.get_kin_rate() << " Decreasing rate of temperature."; - ofs << "\n" << setw(20) << mz.metro.get_kin_ntemp() << " Number of different temperature (for kinetical)."; - ofs << "\n" << setw(20) << mz.metro.get_kin_nsteps() << " Number of steps for each temperature (for kineitcal)."; - ofs << "\n"; + ofs << "\n\n"; + ofs << "\n" << setw(20) << mz.metro.get_spi_ini_temp() << " Start temperature (Kelvin) for spillage minimization."; + ofs << "\n" << setw(20) << mz.metro.get_spi_rate() << " Decreasing rate of temperature."; + ofs << "\n" << setw(20) << mz.metro.get_spi_ntemp() << " Number of different temperature (for spillage)."; + ofs << "\n" << setw(20) << mz.metro.get_spi_nsteps() << " Number of steps for each temperature (for spillage)."; + + ofs << "\n" + << setw(20) << mz.metro.get_kin_ini_temp() << " Start temperature (Kelvin) for kinetical energy minimization."; + ofs << "\n" << setw(20) << mz.metro.get_kin_rate() << " Decreasing rate of temperature."; + ofs << "\n" << setw(20) << mz.metro.get_kin_ntemp() << " Number of different temperature (for kinetical)."; + ofs << "\n" << setw(20) << mz.metro.get_kin_nsteps() << " Number of steps for each temperature (for kineitcal)."; + ofs << "\n"; } -void Out_Orbital::c4_information( ofstream &ofs) +void Out_Orbital::c4_information(ofstream& ofs) { - ofs << "\n\n"; - ofs << "\n" << setw(20) << NCHIUSED << " Total number of radial orbitals."; - - // use multize parameters: ntype, *lmax_type, l_nchi, ne; - for(int it=0; it"; - return; + ofs << "\n\n"; + ofs << "\n" << setw(20) << NCHIUSED << " Total number of radial orbitals."; + + // use multize parameters: ntype, *lmax_type, l_nchi, ne; + for (int it = 0; it < NTYPE; it++) + { + for (int L = 0; L < mz.lmax_type[it] + 1; L++) + { + // last step, nmax is different if L is different + for (int N = 0; N < mz.l_nchi(it, L); N++) + { + ofs << "\n" << setw(20) << "Type" << setw(20) << "L" << setw(20) << "Zeta-Orbital"; + ofs << "\n" << setw(20) << it + 1 << setw(20) << L << setw(20) << N + 1; + for (int ie = 0; ie < mz.Level[0].ne; ie++) + { + if (ie % 4 == 0) + ofs << endl; + // mohan modify 2009-09-25 + // ofs << setw(20) << input.Coef.C4_accumulate(it, L, N, ie); + ofs << setw(25) << setiosflags(ios::showpoint) << setiosflags(ios::scientific) << setprecision(15) + << input.Coef.C4_old(it, L, N, ie); + } + } + } + } + ofs << "\n"; + return; } -void Out_Orbital::mkb_information( ofstream &ofs) +void Out_Orbital::mkb_information(ofstream& ofs) { - ofs_running << "\n mkb_information." << endl; - if(MY_RANK==0) - { - ofs << "\n\n"; - ofs << "\n" << setw(20) << mz.nlevel << " Total number of orbitals optimized levels." ; - } - - double sum0 = 0.0; - double *sum1 = new double[NBANDS]; - ZEROS(sum1, NBANDS); - - int bands_start; - int bands_end; - int bands_number; - - if (BANDS_CONTROL) - { - bands_start = BANDS_START; - bands_end = BANDS_END; - bands_number = bands_end - bands_start; - } - else - { - bands_start = 0; - bands_end = NBANDS; - bands_number = bands_end - bands_start; - } - - ofs << "\nBands start from " << bands_start; - ofs << "\nBands ended at " << bands_end; - ofs << "\nOptimized bands number " << bands_number; - ofs << "\nSpillage per band is " << 1.0/bands_number; - - //ofs << setioflags(ios::fixed); - for(int il=0; il"; + ofs << "\n" << setw(20) << mz.nlevel << " Total number of orbitals optimized levels."; + } + + double sum0 = 0.0; + double* sum1 = new double[NBANDS]; + ZEROS(sum1, NBANDS); + + int bands_start; + int bands_end; + int bands_number; + + if (BANDS_CONTROL) + { + bands_start = BANDS_START; + bands_end = BANDS_END; + bands_number = bands_end - bands_start; + } + else + { + bands_start = 0; + bands_end = NBANDS; + bands_number = bands_end - bands_start; + } + + ofs << "\nBands start from " << bands_start; + ofs << "\nBands ended at " << bands_end; + ofs << "\nOptimized bands number " << bands_number; + ofs << "\nSpillage per band is " << 1.0 / bands_number; + + // ofs << setioflags(ios::fixed); + for (int il = 0; il < mz.ilevel; il++) + { + if (MY_RANK == 0) + { + ofs << "\n\nFill Left Hilbert space of each band(average) by LCAO for Level " << il + 1; + ofs << "\n" + << setw(5) << "BANDS" << setw(20) << "New Fill" << setw(20) << "Total Fill" << setw(20) + << "Left Spillage"; + ofs << setprecision(10); + } + double conb = 0.0; + ofs_running << "\n NBANDS=" << NBANDS << " STRNUM=" << STRNUM << endl; + + for (int ib = bands_start; ib < bands_end; ib++) + { + double mm = 0.0; + // calculate the average fill gution from all k points. + for (int is = 0; is < STRNUM; is++) + { + + for (int ik = 0; ik < mz.Level[il].data[is].nks; ik++) + { + mm += mz.Level[il].data[is].Mkb(ik, ib) * mz.Level[il].data[is].weight[ik] / (double)bands_number + / (double)STRNUM; + } #ifdef __MPI - Parallel_Reduce::reduce_double_allpool(mm); + Parallel_Reduce::reduce_double_allpool(mm); #endif - } - sum1[ib] += mm; - - if(MY_RANK==0) - { - ofs << "\n" << setw(5) << ib+1 << setw(20) << mm << setw(20) << sum1[ib] << setw(20) << 1.0/(double)bands_number - sum1[ib]; - } - - conb += mm; - } - sum0 += conb; - ofs << "\nNew Fill Contribution = " << conb; - ofs << "\nTotal Fill Contribution = " << sum0; - - ofs << "\nLeft spillage = " << 1.0 - sum0; - } - ofs << "\n"; - delete[] sum1; - - return; + } + sum1[ib] += mm; + + if (MY_RANK == 0) + { + ofs << "\n" + << setw(5) << ib + 1 << setw(20) << mm << setw(20) << sum1[ib] << setw(20) + << 1.0 / (double)bands_number - sum1[ib]; + } + + conb += mm; + } + sum0 += conb; + ofs << "\nNew Fill Contribution = " << conb; + ofs << "\nTotal Fill Contribution = " << sum0; + + ofs << "\nLeft spillage = " << 1.0 - sum0; + } + ofs << "\n"; + delete[] sum1; + + return; } diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Out_Orbital.h b/tools/01_NAO_generation/SIAB/src_spillage/Out_Orbital.h index 228ed3857d8..b3feb7f59d4 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Out_Orbital.h +++ b/tools/01_NAO_generation/SIAB/src_spillage/Out_Orbital.h @@ -6,20 +6,19 @@ // output orbital information class Out_Orbital { - public: - Out_Orbital(); - ~Out_Orbital(); + public: + Out_Orbital(); + ~Out_Orbital(); - void write(void); + void write(void); - private: - void version_information( ofstream &ofs ); - void INPUTs_information( ofstream &ofs); - void spillage_information( ofstream &ofs, const int &ilevel); - void metropolis_information( ofstream &ofs); - void c4_information( ofstream &ofs); - void mkb_information( ofstream &ofs); - + private: + void version_information(ofstream& ofs); + void INPUTs_information(ofstream& ofs); + void spillage_information(ofstream& ofs, const int& ilevel); + void metropolis_information(ofstream& ofs); + void c4_information(ofstream& ofs); + void mkb_information(ofstream& ofs); }; #endif diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Plot_Psi.cpp b/tools/01_NAO_generation/SIAB/src_spillage/Plot_Psi.cpp index 468c14ea315..b132add98f1 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Plot_Psi.cpp +++ b/tools/01_NAO_generation/SIAB/src_spillage/Plot_Psi.cpp @@ -1,625 +1,629 @@ #include "Plot_Psi.h" -#include "tools.h" + #include "../src_tools/Gram_Schmidt_Orth.h" +#include "tools.h" Plot_Psi::Plot_Psi() { - test=0; - rcut = 0.0; - - dr = 0.01; - - meshU = 0; // uniform mesh. - meshL = 0; // Log mesh. - - rU = new double[1]; - rL = new double[1]; - - xmin = -6; - zed = 4;//charge on the atom - dx = 1.0e-2; - xmax = 10; + test = 0; + rcut = 0.0; + + dr = 0.01; + + meshU = 0; // uniform mesh. + meshL = 0; // Log mesh. + + rU = new double[1]; + rL = new double[1]; + + xmin = -6; + zed = 4; // charge on the atom + dx = 1.0e-2; + xmax = 10; } Plot_Psi::~Plot_Psi() { - delete[] rU; - delete[] rL; + delete[] rU; + delete[] rL; } void Plot_Psi::allocate(void) -{ - // (1) read in parameters. - if(SCAN_BEGIN(input.ifs, "")) +{ + // (1) read in parameters. + if (SCAN_BEGIN(input.ifs, "")) + { + READ_VALUE(input.ifs, this->dr); + assert(this->dr > 0.0); + + // mohan add 2009-11-23 + READ_VALUE(input.ifs, this->xmin); + READ_VALUE(input.ifs, this->zed); // charge? + READ_VALUE(input.ifs, this->dx); + READ_VALUE(input.ifs, this->xmax); + assert(zed > 0.0); + assert(dx > 0.0); + assert(xmax >= 0.0); + + SCAN_END(input.ifs, ""); + } + + this->ne = NE; + this->meshU = static_cast(RCUT / dr); + this->meshL = (std::log(xmax * zed) - xmin) / dx; + // mesh should be odd. + if (meshU % 2 == 0) + meshU++; + if (meshL % 2 == 0) + meshL++; + + /* + cout << "\n dr = " << dr; + cout << "\n meshU = " << meshU; + cout << "\n meshL = " << meshL; + + cout << "\n xmin = " << xmin; + cout << "\n zed = " << zed; + cout << "\n dx = " << dx; + cout << "\n xmax = " << xmax; + */ + + // other parameters. + this->rcut = RCUT; + + delete[] rU; + delete[] rL; + rU = new double[meshU]; + rL = new double[meshL]; + + for (int ir = 0; ir < meshU; ir++) + { + rU[ir] = dr * ir; + } + + for (int ir = 0; ir < meshL; ir++) { - READ_VALUE(input.ifs, this->dr); - assert( this->dr > 0.0); - - // mohan add 2009-11-23 - READ_VALUE(input.ifs, this->xmin); - READ_VALUE(input.ifs, this->zed);//charge? - READ_VALUE(input.ifs, this->dx); - READ_VALUE(input.ifs, this->xmax); - assert(zed>0.0); - assert(dx>0.0); - assert(xmax>=0.0); - - SCAN_END(input.ifs, ""); - } - - this->ne = NE; - this->meshU = static_cast(RCUT/dr); - this->meshL = ( std::log( xmax * zed ) - xmin ) / dx; - // mesh should be odd. - if(meshU%2==0) meshU++; - if(meshL%2==0) meshL++; - - /* - cout << "\n dr = " << dr; - cout << "\n meshU = " << meshU; - cout << "\n meshL = " << meshL; - - cout << "\n xmin = " << xmin; - cout << "\n zed = " << zed; - cout << "\n dx = " << dx; - cout << "\n xmax = " << xmax; - */ - - // other parameters. - this->rcut = RCUT; - - delete[] rU; - delete[] rL; - rU = new double[meshU]; - rL = new double[meshL]; - - for(int ir=0; irlmax = LMAXUSED; - this->total_nchi = NCHIUSED; - this->psiU.create( this->total_nchi, this->meshU); - this->psiL.create( this->total_nchi, this->meshL); - - return; + rL[ir] = std::exp(xmin + dx * ir) / zed; + } + + this->lmax = LMAXUSED; + this->total_nchi = NCHIUSED; + this->psiU.create(this->total_nchi, this->meshU); + this->psiL.create(this->total_nchi, this->meshL); + + return; } // be called in main.cpp void Plot_Psi::radial_wave_function(void) { - TITLE("Plot_Psi","radial_wave_function"); - timer::tick("Plot_Psi","radial_wf"); - - this->allocate(); - - // calculate all psi - int ichi = 0; - for(int it=0; it>>>>>>>>> - //The two formula are exactly same. - //<< setw(15) << abs(c4[ie])*pow(eigen1[ie],2) - //<< setw(15) << abs(c4[ie])*pow((ie+1)*PI/input.inputs.get_rcut(),2); - //<<<<<<<<<< - } - mz.metro.Psi2.ofsk << "\n Total kinetic energy = " << kin_sum; - - double *psi1 = new double[meshU]; - double *psi2 = new double[meshL]; - this->get_psi( meshU, rU, L, psi1, c4, eigen1, 0 ); - this->get_psi( meshL, rL, L, psi2, c4, eigen1, 1 ); - - for(int ir=0; irorthogonalization(meshU, rU, psiU); - this->orthogonalization(meshL, rL, psiL); - - this->print_orbitalsU_used(); - this->print_orbitalsL_used(); - - // mohan add 2010-04-13 - // use uniform grid to establish the energy cutoff of each orbital - ofstream ofs("ORBITAL_ECUT.txt"); - ofs << setprecision(12); - this->dk = 0.01; - this->meshK = 1500;// (1500*0.01)^2 * 2 = 450 Rydberg. - if(meshK%2==0)meshK++; - this->psiUK.create(total_nchi,meshK); - double *psi = new double[meshU]; - ichi = 0; - for(int it=0; itallocate(); + + // calculate all psi + int ichi = 0; + for (int it = 0; it < NTYPE; it++) + { + if (test == 2) + cout << "\n mz.lmax_type[" << it << "]=" << mz.lmax_type[it]; + for (int L = 0; L < mz.lmax_type[it] + 1; L++) + { + double* eigen1 = new double[ne]; + for (int ie = 0; ie < ne; ie++) + { + eigen1[ie] = mz.metro.Psi2.eigenvalue(L, ie); + // cout << "\n eigen1[" << ie << "]=" << eigen1[ie]; + } + + for (int N = 0; N < mz.l_nchi(it, L); N++) + { + if (test == 2) + cout << "\n T=" << it << " L=" << L << " N=" << N; + double* c4 = new double[ne]; + + for (int ie = 0; ie < ne; ie++) + { + c4[ie] = input.Coef.C4_old(it, L, N, ie); + } + + // call the function to renormalize the wave functions. + mz.metro.Psi2.norm_c4(c4, L); + + // calculate kinetical energy + mz.metro.Psi2.ofsk << "\n" << "Type=" << it + 1 << " L=" << L << " N=" << N + 1; + mz.metro.Psi2.ofsk << "\n" + << setw(5) << "nq" << setw(20) << "Coefficient" << setw(20) << "Eigen(Ry)" + << setw(20) << "Kinetic_Energy"; + double kin_sum = 0.0; + for (int ie = 0; ie < ne; ie++) + { + // double kin = abs(c4[ie])*pow((ie+1)*PI/input.inputs.get_rcut(),2);//In Rydberg Unit. + // double kin = abs(c4[ie])*pow(eigen1[ie]*PI/input.inputs.get_rcut(),2);//In Rydberg Unit. + double kin = pow(c4[ie] * eigen1[ie], 2) * mz.metro.Psi2.jjnorm(L, ie); + kin_sum += kin; + mz.metro.Psi2.ofsk << "\n" + << setw(5) << ie << setiosflags(ios::fixed) << setprecision(10) + << setiosflags(ios::showpoint) + //<< setiosflags(ios::scientific) + << setw(20) << c4[ie] << setw(20) << eigen1[ie] * eigen1[ie] << setw(20) << kin + << " (Rydberg)" << setprecision(6); + + //>>>>>>>>>> + // The two formula are exactly same. + //<< setw(15) << abs(c4[ie])*pow(eigen1[ie],2) + //<< setw(15) << abs(c4[ie])*pow((ie+1)*PI/input.inputs.get_rcut(),2); + //<<<<<<<<<< + } + mz.metro.Psi2.ofsk << "\n Total kinetic energy = " << kin_sum; + + double* psi1 = new double[meshU]; + double* psi2 = new double[meshL]; + this->get_psi(meshU, rU, L, psi1, c4, eigen1, 0); + this->get_psi(meshL, rL, L, psi2, c4, eigen1, 1); + + for (int ir = 0; ir < meshU; ir++) + { + psiU(ichi, ir) = psi1[ir]; + } + for (int ir = 0; ir < meshL; ir++) + { + psiL(ichi, ir) = psi2[ir]; + } + delete[] psi2; + delete[] psi1; + delete[] c4; + ++ichi; + } + delete[] eigen1; + } + } + + // Peize Lin add 2015-11-20 + this->orthogonalization(meshU, rU, psiU); + this->orthogonalization(meshL, rL, psiL); + + this->print_orbitalsU_used(); + this->print_orbitalsL_used(); + + // mohan add 2010-04-13 + // use uniform grid to establish the energy cutoff of each orbital + ofstream ofs("ORBITAL_ECUT.txt"); + ofs << setprecision(12); + this->dk = 0.01; + this->meshK = 1500; // (1500*0.01)^2 * 2 = 450 Rydberg. + if (meshK % 2 == 0) + meshK++; + this->psiUK.create(total_nchi, meshK); + double* psi = new double[meshU]; + ichi = 0; + for (int it = 0; it < NTYPE; it++) + { + for (int L = 0; L < mz.lmax_type[it] + 1; L++) + { + for (int N = 0; N < mz.l_nchi(it, L); N++) + { + ofs << "WaveFunctionIndex=" << ichi << endl; + ofs << "T=" << it << " L=" << L << " N=" << N << endl; + for (int ir = 0; ir < meshU; ir++) + { + psi[ir] = this->psiU(ichi, ir); + } + this->establish_ecut(ofs, ichi, psi, L); + ++ichi; + } + } + } + assert(ichi == total_nchi); + delete[] psi; + ofs.close(); + + // print the orbitals which can be viewed imediately. + this->print_orbitals_plot(); + timer::tick("Plot_Psi", "radial_wf"); + return; } -void Plot_Psi::print_orbitals_plot(void)const -{ - ofstream f1( "ORBITAL_PLOTU.dat"); - ofstream f2( "ORBITAL_PLOTL.dat"); - ofstream f3( "ORBITAL_PLOTUK.dat"); - - f1 << setprecision(12) - << setiosflags(ios::fixed); - f2 << setprecision(12) - << setiosflags(ios::fixed); - f3 << setprecision(12) - << setiosflags(ios::fixed); - for(int ir=0; ir " << setw(7) << mz.l_nchi(it,L) << endl; - } - f1 << "---------------------------------------------------------------------------"<psiU(ichi, ir); - } - ++ichi; - } - } - f1.close(); - } - return; + //==================================== + // (5) + //==================================== + // output orbital for standard use form + int ichi = 0; + // because we can't read in r = 0 term now. + // but liaochen suggest we should know psi(r) + // instead of psi(r)*r at zero point, + // I think it's correct. + int start = 0; + for (int it = 0; it < NTYPE; it++) + { + // output the wave function in uniform mesh. + stringstream ss; + ss << "ORBITAL_" << mz.Level[0].info[it].id << "U.dat"; + + ofstream f1(ss.str().c_str()); // pengfei 2014-10-13 + f1 << "---------------------------------------------------------------------------" << endl; + f1 << "Element" << setw(22) << LABEL[it] << endl; + f1 << "Energy Cutoff(Ry)" << setw(13) << ECUT << endl; + f1 << "Radius Cutoff(a.u.)" << setw(10) << RCUT << endl; + f1 << "Lmax " << setw(24) << mz.lmax_type[it] << endl; + for (int L = 0; L < mz.lmax_type[it] + 1; L++) + { + string orbital_type; + switch (L) + { + case 0: + orbital_type = "S"; + break; + case 1: + orbital_type = "P"; + break; + case 2: + orbital_type = "D"; + break; + case 3: + orbital_type = "F"; + break; + case 4: + orbital_type = "G"; + break; + } + + f1 << "Number of " << orbital_type << "orbital--> " << setw(7) << mz.l_nchi(it, L) << endl; + } + f1 << "---------------------------------------------------------------------------" << endl; + f1 << "SUMMARY" << " " << "END" << endl; + f1 << endl; + f1 << "Mesh " << setw(26) << meshU - start << endl; // mohan fix bug 2010-04-15 + f1 << "dr " << setw(29) << dr; + + for (int L = 0; L < mz.lmax_type[it] + 1; L++) + { + double* eigen1 = new double[ne]; + for (int ie = 0; ie < ne; ie++) + { + eigen1[ie] = mz.metro.Psi2.eigenvalue(L, ie); + } + for (int N = 0; N < mz.l_nchi(it, L); N++) + { + f1 << "\n" << setw(20) << "Type" << setw(20) << "L" << setw(20) << "N"; + f1 << "\n" << setw(20) << it << setw(20) << L << setw(20) << N; + f1 << setiosflags(ios::scientific) << setprecision(12); + for (int ir = start; ir < meshU; ir++) + { + if (ir % 4 == 0) + f1 << endl; + f1 << " " << this->psiU(ichi, ir); + } + ++ichi; + } + } + f1.close(); + } + return; } -void Plot_Psi::print_orbitalsL_used(void)const +void Plot_Psi::print_orbitalsL_used(void) const { - //==================================== - // (5) - //==================================== - // output orbital for standard use form - int ichi = 0; - // because we can't read in r = 0 term now. - for(int it=0; itpsiL(ichi, ir); - } - ++ichi; - } - } - f1.close(); - } - return; + //==================================== + // (5) + //==================================== + // output orbital for standard use form + int ichi = 0; + // because we can't read in r = 0 term now. + for (int it = 0; it < NTYPE; it++) + { + // output the wave function in uniform mesh. + stringstream ss; + ss << "ORBITAL_" << mz.Level[0].info[it].id << "L.dat"; + + ofstream f1(ss.str().c_str()); + + f1 << "Mesh " << meshL << endl; + f1 << "zed " << zed << endl; + f1 << "dx " << dx << endl; + f1 << "xmin " << xmin << endl; + f1 << "xmax " << xmax; + + for (int L = 0; L < mz.lmax_type[it] + 1; L++) + { + double* eigen1 = new double[ne]; + for (int ie = 0; ie < ne; ie++) + { + eigen1[ie] = mz.metro.Psi2.eigenvalue(L, ie); + } + for (int N = 0; N < mz.l_nchi(it, L); N++) + { + f1 << "\n" << setw(20) << "Type" << setw(20) << "L" << setw(20) << "N"; + f1 << "\n" << setw(20) << it << setw(20) << L << setw(20) << N << endl; + f1 << setiosflags(ios::scientific) << setprecision(12); + for (int ir = 0; ir < meshL; ir++) + { + if (ir % 4 == 0) + f1 << endl; + f1 << " " << this->psiL(ichi, ir); + } + ++ichi; + } + } + f1.close(); + } + return; } -void Plot_Psi::get_psi( - const int &mesh, - const double *r, - const int &l, - double *psi, - double *c4, - const double *eigen1, bool logmesh)const +void Plot_Psi::get_psi(const int& mesh, + const double* r, + const int& l, + double* psi, + double* c4, + const double* eigen1, + bool logmesh) const { - // (1) get psi1 - double *jle = new double[mesh]; - double *g = new double[mesh]; // smooth function - - if(SMOOTH) - { - double sigma2 = SIGMA*SIGMA; - for(int ir=0; irrcut)*(r[ir]-this->rcut)/2.0/sigma2) ); - //if(ir==mesh-1) cout << "\n g[" << ir << "]=" << g[ir] << endl; - } - } - - ZEROS(psi,mesh); - for(int ie=0; ie and - // using smooth Jlq, here we output the - // smooth functions. - if(SMOOTH) - { - for(int ir=0; irdr; - } - } - - for(int ir=0; ir0.9900)=" << pow(dk*kmesh_used,2)*2 << endl; - - while(normk*fpi< 0.999 && kmesh_used < this->meshK) - { - kmesh_used += 2; - Mathzone::Simpson_Integral(kmesh_used, functionk, rabk, normk ); - } -// ofs << "Psi(k) norm(0.999)=" << normk * fpi << endl; - ofs << "Ecut(Ry)(norm>0.9990)=" << pow(dk*kmesh_used,2)*2 << endl; - - while(normk*fpi< 0.9999 && kmesh_used < this->meshK) - { - kmesh_used += 2; - Mathzone::Simpson_Integral(kmesh_used, functionk, rabk, normk ); - } -// ofs << "Psi(k) norm(0.999)=" << normk * fpi << endl; - ofs << "Ecut(Ry)(norm>0.9999)=" << pow(dk*kmesh_used,2)*2 << endl; - - ofs << endl; - - delete[] function; - delete[] functionk; - delete[] rabk; - delete[] rab; - delete[] psik; - delete[] jj; - return; + //(1) init the parameters + double pi = 3.1415926535; + double fpi = 4 * pi; + + for (int ir = 0; ir < this->meshU; ir++) + { + psi[ir] /= sqrt(fpi); + } + + // (2) check the norm of psi + double norm = 0.0; + double* function = new double[meshU]; + double* rab = new double[meshU]; + for (int ir = 0; ir < meshU; ir++) + { + rab[ir] = dr; + function[ir] = std::pow(psi[ir] * this->rU[ir], 2); + } + Mathzone::Simpson_Integral(meshU, function, rab, norm); + ofs << "Psi(r) norm=" << norm * fpi << endl; + + // (3) get psik + double* psik = new double[meshK]; + double* jj = new double[meshU]; + for (int ik = 0; ik < meshK; ik++) + { + double q = ik * this->dk; + Mathzone::Spherical_Bessel(this->meshU, this->rU, q, L, jj); + for (int ir = 0; ir < this->meshU; ir++) + { + function[ir] = psi[ir] * this->rU[ir] * this->rU[ir] * jj[ir]; + } + Mathzone::Simpson_Integral(this->meshU, function, rab, psik[ik]); + } + double prefac_k = fpi / pow(2 * pi, 1.5); + for (int ik = 0; ik < meshK; ik++) + { + psik[ik] *= prefac_k; + this->psiUK(ichi, ik) = psik[ik]; + } + + // (4) prepare for the integral of psik. + double normk = 0.0; + double* rabk = new double[meshK]; + double* functionk = new double[meshK]; + for (int ik = 0; ik < meshK; ik++) + { + rabk[ik] = dk; + } + for (int ik = 0; ik < meshK; ik++) + { + functionk[ik] = std::pow(psik[ik] * ik * this->dk, 2); + } + Mathzone::Simpson_Integral(this->meshK, functionk, rabk, normk); + ofs << "Psi(k) norm=" << normk * fpi << endl; + + // (5) Find out the energy cut for each orbital. + normk = 0.0; + int kmesh_used = 1; + + //=================================================== + // 0.999 is the accepted value of norm, according + // to the value we establish the energy cutoff. + //=================================================== + while (normk * fpi < 0.99 && kmesh_used < this->meshK) + { + kmesh_used += 2; + Mathzone::Simpson_Integral(kmesh_used, functionk, rabk, normk); + } + // ofs << "Psi(k) norm(0.99)=" << normk * fpi << endl; + ofs << "Ecut(Ry)(norm>0.9900)=" << pow(dk * kmesh_used, 2) * 2 << endl; + + while (normk * fpi < 0.999 && kmesh_used < this->meshK) + { + kmesh_used += 2; + Mathzone::Simpson_Integral(kmesh_used, functionk, rabk, normk); + } + // ofs << "Psi(k) norm(0.999)=" << normk * fpi << endl; + ofs << "Ecut(Ry)(norm>0.9990)=" << pow(dk * kmesh_used, 2) * 2 << endl; + + while (normk * fpi < 0.9999 && kmesh_used < this->meshK) + { + kmesh_used += 2; + Mathzone::Simpson_Integral(kmesh_used, functionk, rabk, normk); + } + // ofs << "Psi(k) norm(0.999)=" << normk * fpi << endl; + ofs << "Ecut(Ry)(norm>0.9999)=" << pow(dk * kmesh_used, 2) * 2 << endl; + + ofs << endl; + + delete[] function; + delete[] functionk; + delete[] rabk; + delete[] rab; + delete[] psik; + delete[] jj; + return; } // Peize Lin add 2015-11-20 -void Plot_Psi::orthogonalization( - const int &mesh, - double *r, - const matrix &psi -) +void Plot_Psi::orthogonalization(const int& mesh, double* r, const matrix& psi) { - Gram_Schmidt_Orth gs_orth; - gs_orth.set_coordinate(1); - gs_orth.set_variable_num(mesh); - gs_orth.set_rab_element(this->dr); - gs_orth.set_r(r); - - double * tmp_index(psi.c); - - for(int it=0; it gs_orth; + gs_orth.set_coordinate(1); + gs_orth.set_variable_num(mesh); + gs_orth.set_rab_element(this->dr); + gs_orth.set_r(r); + + double* tmp_index(psi.c); + + for (int it = 0; it < NTYPE; it++) + { + for (int L = 0; L < mz.lmax_type[it] + 1; L++) + { + gs_orth.set_func_num(mz.l_nchi(it, L)); + double** psi_TL = new double*[mz.l_nchi(it, L)]; + for (int N = 0; N < mz.l_nchi(it, L); N++) + { + psi_TL[N] = tmp_index; + tmp_index += mesh; + } + gs_orth.set_func(psi_TL); + gs_orth.orth(); + delete[] psi_TL; + } + } + + return; } diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Plot_Psi.h b/tools/01_NAO_generation/SIAB/src_spillage/Plot_Psi.h index b95015b2673..95ae35f77a1 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Plot_Psi.h +++ b/tools/01_NAO_generation/SIAB/src_spillage/Plot_Psi.h @@ -5,53 +5,57 @@ class Plot_Psi { - public: - Plot_Psi(); - ~Plot_Psi(); - - void radial_wave_function(void); - - private: - - // (1) about uniform mesh. - double dr; - int meshU; // U stands for uniform. - int meshL;// L stands for log. - double* rU;// uniform mesh. - double* rL;// log mesh. - int meshK; - double dk;//classic bug, change dk to 'int' type. then see .... - - // (2) about logmesh. - double xmin; - double zed; - double dx; - double xmax; - - // (3) about the wave functions. - matrix psiU; // sum c*j(l,x) - matrix psiL; // - matrix psiUK; - - int total_nchi; - double rcut; - int lmax; - int ne; - int test; - - void allocate(void); - - // mohan add smooth, sigam and rcut 2009-08-28 - void get_psi( const int &mesh, const double *r,const int &l, - double *psi1, double *c4, const double *eigen1, bool logmesh)const; - void orthogonalization( const int &mesh, double *r, const matrix &psi ); - - void print_orbitals_plot(void)const; - void print_orbitalsU_used(void)const; - void print_orbitalsL_used(void)const; - - // mohan add 2010-04-13 - void establish_ecut( ofstream &ofs, const int &ik, double *psi, const int &L); + public: + Plot_Psi(); + ~Plot_Psi(); + + void radial_wave_function(void); + + private: + // (1) about uniform mesh. + double dr; + int meshU; // U stands for uniform. + int meshL; // L stands for log. + double* rU; // uniform mesh. + double* rL; // log mesh. + int meshK; + double dk; // classic bug, change dk to 'int' type. then see .... + + // (2) about logmesh. + double xmin; + double zed; + double dx; + double xmax; + + // (3) about the wave functions. + matrix psiU; // sum c*j(l,x) + matrix psiL; // + matrix psiUK; + + int total_nchi; + double rcut; + int lmax; + int ne; + int test; + + void allocate(void); + + // mohan add smooth, sigam and rcut 2009-08-28 + void get_psi(const int& mesh, + const double* r, + const int& l, + double* psi1, + double* c4, + const double* eigen1, + bool logmesh) const; + void orthogonalization(const int& mesh, double* r, const matrix& psi); + + void print_orbitals_plot(void) const; + void print_orbitalsU_used(void) const; + void print_orbitalsL_used(void) const; + + // mohan add 2010-04-13 + void establish_ecut(ofstream& ofs, const int& ik, double* psi, const int& L); }; #endif diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Psi_Second.cpp b/tools/01_NAO_generation/SIAB/src_spillage/Psi_Second.cpp index 2f7a5c84e77..d652870b38d 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Psi_Second.cpp +++ b/tools/01_NAO_generation/SIAB/src_spillage/Psi_Second.cpp @@ -1,332 +1,331 @@ #include "Psi_Second.h" + #include "Calculate_C4.h" #include "tools.h" Psi_Second::Psi_Second() { - test=0; - r = new double[1]; - rab = new double[1]; - psi = new double[1]; - psi_first = new double[1]; - psi_second = new double[1]; - eigen1 = new double[1]; - os_position = new double[100]; - below_max_os = new bool[1]; - count = 0; + test = 0; + r = new double[1]; + rab = new double[1]; + psi = new double[1]; + psi_first = new double[1]; + psi_second = new double[1]; + eigen1 = new double[1]; + os_position = new double[100]; + below_max_os = new bool[1]; + count = 0; } Psi_Second::~Psi_Second() { - delete[] r; - delete[] rab; - delete[] psi; - delete[] psi_first; - delete[] psi_second; - delete[] eigen1; - delete[] os_position; - delete[] below_max_os; + delete[] r; + delete[] rab; + delete[] psi; + delete[] psi_first; + delete[] psi_second; + delete[] eigen1; + delete[] os_position; + delete[] below_max_os; } -void Psi_Second::init( - ifstream &ifs, - const int &ne_in, - const int &lmax_in, - const double &rcut_in, - const int &total_nchi_in, - const double &tolerence) +void Psi_Second::init(ifstream& ifs, + const int& ne_in, + const int& lmax_in, + const double& rcut_in, + const int& total_nchi_in, + const double& tolerence) { - TITLE("Psi_Second","init"); - - this->dr = KINETIC_DR; - assert( dr > 0.0); - cout << "\n KINETIC_DR=" << dr; - -// this->start = 5; // start away from 0, which may be some strange oscillation. - this->start = 1; // mohan fix 2010-04-16 - - // number of Jlq used. - this->ne = ne_in; - // lmax needed. - this->lmax = lmax_in; - // Cutoff radius. - this->rcut = rcut_in; - // Total number of radial wave functions. - this->total_nchi = total_nchi_in; - - this->mesh = static_cast(this->rcut/this->dr); - if(this->mesh%2==0) this->mesh++; - - if(test==2) - { - cout << "\n ne=" << ne; - cout << "\n lmax=" << lmax; - cout << "\n mesh=" << mesh; - cout << "\n total_nchi=" << total_nchi; - cout << "\n dr=" << dr; - } - - // init psi matrix, generate eigenvalue matrix - // this eigenvalue will be used in many functions. - this->eigenvalue.create( this->lmax+1, this->ne); - this->jjnorm.create( this->lmax+1, this->ne); - - delete[] this->eigen1; - this->eigen1 = new double[ne]; - - delete[] this->r; - delete[] this->rab; - delete[] this->psi; - delete[] this->psi_first; - delete[] this->psi_second; - delete[] this->below_max_os; - delete[] this->os_number; - - this->r = new double[mesh]; - this->rab = new double[mesh]; - this->psi = new double[mesh]; - this->psi_first = new double[mesh]; - this->psi_second = new double[mesh]; - - this->below_max_os = new bool[total_nchi]; - this->os_number = new int[total_nchi]; - - for(int ir=0; irne, - this->eigenvalue, - this->lmax, - this->rcut, - tolerence); - - // (2) to calculate jl(qr) - this->jle = new double**[this->lmax+1]; - double* g; - if(SMOOTH) - { - g = new double[mesh]; - double sigma2 = SIGMA*SIGMA; - for(int ir=0; irrcut)*(r[ir]-this->rcut)/2.0/sigma2) ); - } - } - - for(int l=0; leigen1[ie] = this->eigenvalue(l, ie);} - - // allocate space - this->jle[l] = new double*[ne]; - - for(int ie=0; iene; ie++) - { - this->jle[l][ie] = new double[this->mesh]; - Mathzone::Spherical_Bessel( this->mesh, this->r, this->eigen1[ie], l, jle[l][ie]); - - // (3) to calculate the jjnorm - double* f = new double[mesh]; - for(int ir=0; irjjnorm(l,ie) ); - jjnorm(l,ie)*=4.0*PI; -// cout << "\n jj = " << jjnorm(l,ie); - delete[] f; - - if(SMOOTH) - for(int ir=0; irofso.open("ORBITAL_OSCILATION.txt"); - this->ofsk.open("ORBITAL_KINETIC.txt"); - - return; + TITLE("Psi_Second", "init"); + + this->dr = KINETIC_DR; + assert(dr > 0.0); + cout << "\n KINETIC_DR=" << dr; + + // this->start = 5; // start away from 0, which may be some strange oscillation. + this->start = 1; // mohan fix 2010-04-16 + + // number of Jlq used. + this->ne = ne_in; + // lmax needed. + this->lmax = lmax_in; + // Cutoff radius. + this->rcut = rcut_in; + // Total number of radial wave functions. + this->total_nchi = total_nchi_in; + + this->mesh = static_cast(this->rcut / this->dr); + if (this->mesh % 2 == 0) + this->mesh++; + + if (test == 2) + { + cout << "\n ne=" << ne; + cout << "\n lmax=" << lmax; + cout << "\n mesh=" << mesh; + cout << "\n total_nchi=" << total_nchi; + cout << "\n dr=" << dr; + } + + // init psi matrix, generate eigenvalue matrix + // this eigenvalue will be used in many functions. + this->eigenvalue.create(this->lmax + 1, this->ne); + this->jjnorm.create(this->lmax + 1, this->ne); + + delete[] this->eigen1; + this->eigen1 = new double[ne]; + + delete[] this->r; + delete[] this->rab; + delete[] this->psi; + delete[] this->psi_first; + delete[] this->psi_second; + delete[] this->below_max_os; + delete[] this->os_number; + + this->r = new double[mesh]; + this->rab = new double[mesh]; + this->psi = new double[mesh]; + this->psi_first = new double[mesh]; + this->psi_second = new double[mesh]; + + this->below_max_os = new bool[total_nchi]; + this->os_number = new int[total_nchi]; + + for (int ir = 0; ir < mesh; ir++) + { + r[ir] = dr * ir; + rab[ir] = dr; + } + + // (1) call to find the eigenvalues of Jlq. + Calculate_C4::Find_Eigenvalues(this->ne, this->eigenvalue, this->lmax, this->rcut, tolerence); + + // (2) to calculate jl(qr) + this->jle = new double**[this->lmax + 1]; + double* g; + if (SMOOTH) + { + g = new double[mesh]; + double sigma2 = SIGMA * SIGMA; + for (int ir = 0; ir < mesh; ir++) + { + g[ir] = 1.0 - std::exp(-((r[ir] - this->rcut) * (r[ir] - this->rcut) / 2.0 / sigma2)); + } + } + + for (int l = 0; l < lmax + 1; l++) + { + // get eigenvalues. + for (int ie = 0; ie < ne; ie++) + { + this->eigen1[ie] = this->eigenvalue(l, ie); + } + + // allocate space + this->jle[l] = new double*[ne]; + + for (int ie = 0; ie < this->ne; ie++) + { + this->jle[l][ie] = new double[this->mesh]; + Mathzone::Spherical_Bessel(this->mesh, this->r, this->eigen1[ie], l, jle[l][ie]); + + // (3) to calculate the jjnorm + double* f = new double[mesh]; + for (int ir = 0; ir < mesh; ir++) + { + f[ir] = std::pow(jle[l][ie][ir], 2); + } + Mathzone::Simpson_Integral(mesh, f, rab, this->jjnorm(l, ie)); + jjnorm(l, ie) *= 4.0 * PI; + // cout << "\n jj = " << jjnorm(l,ie); + delete[] f; + + if (SMOOTH) + for (int ir = 0; ir < mesh; ir++) + { + jle[l][ie][ir] *= g[ir]; + } + } + } + + if (SMOOTH) + delete[] g; + + // this->ofso.open("ORBITAL_OSCILATION.txt"); + this->ofsk.open("ORBITAL_KINETIC.txt"); + + return; } - -double Psi_Second::get_ecut( double *c4, const int &L) +double Psi_Second::get_ecut(double* c4, const int& L) { - timer::tick("Psi_Second","get_ecut"); - - static double fpi = 4*PI; - static double sfpi = sqrt(fpi); - static double prefac_k = fpi / pow(2*PI,1.5); - static double dk = 0.01; - static int meshK = 800; // 8*8*2=128Ry - static double* g; - static double* rab; - static double* rabk; - static double* functionk; - static bool init = false; - static double norm; - static double normk; - static double ecut; - - if(!init) - { - if(meshK%2==0) ++meshK; - - g = new double[mesh]; // smooth function - rab = new double[mesh]; - rabk = new double[meshK]; - functionk = new double[meshK]; - if(SMOOTH) - { - double sigma2 = SIGMA*SIGMA; - for(int ir=0; irrcut)*(r[ir]-this->rcut)/2.0/sigma2) ); - } - } - for(int ir=0; ireigen1[ie] = this->eigenvalue(L, ie);} - - ZEROS(psi,mesh); - for(int ie=0; ie0.999)=" << ecut << endl; - - delete[] jj; - delete[] psik; - delete[] function; - timer::tick("Psi_Second","get_ecut"); - return ecut; + timer::tick("Psi_Second", "get_ecut"); + + static double fpi = 4 * PI; + static double sfpi = sqrt(fpi); + static double prefac_k = fpi / pow(2 * PI, 1.5); + static double dk = 0.01; + static int meshK = 800; // 8*8*2=128Ry + static double* g; + static double* rab; + static double* rabk; + static double* functionk; + static bool init = false; + static double norm; + static double normk; + static double ecut; + + if (!init) + { + if (meshK % 2 == 0) + ++meshK; + + g = new double[mesh]; // smooth function + rab = new double[mesh]; + rabk = new double[meshK]; + functionk = new double[meshK]; + if (SMOOTH) + { + double sigma2 = SIGMA * SIGMA; + for (int ir = 0; ir < mesh; ir++) + { + g[ir] = 1.0 - std::exp(-((r[ir] - this->rcut) * (r[ir] - this->rcut) / 2.0 / sigma2)); + } + } + for (int ir = 0; ir < mesh; ir++) + { + rab[ir] = dr; + } + for (int ik = 0; ik < meshK; ik++) + { + rabk[ik] = dk; + } + init = true; + } + + // get the eigenvalues. + for (int ie = 0; ie < ne; ie++) + { + this->eigen1[ie] = this->eigenvalue(L, ie); + } + + ZEROS(psi, mesh); + for (int ie = 0; ie < ne; ie++) + { + for (int ir = 0; ir < mesh; ir++) + { + psi[ir] += c4[ie] * jle[L][ie][ir]; + } + } + + double* function = new double[mesh]; + for (int ir = 0; ir < mesh; ir++) + { + function[ir] = pow(psi[ir] * r[ir], 2); + } + + Mathzone::Simpson_Integral(mesh, function, rab, norm); + // cout << "\n norm = " << norm * fpi<< endl; + assert(norm != 0.0); + // exit(0); + + norm = sqrt(norm); + for (int ir = 0; ir < mesh; ir++) + { + psi[ir] /= norm; + } + + for (int ie = 0; ie < ne; ie++) + { + c4[ie] /= norm; + } + + // NOTICE !!!! here psi ==> psi*r*r + for (int ir = 0; ir < mesh; ir++) + { + psi[ir] = psi[ir] / sfpi * r[ir] * r[ir]; + } + + // calculate psik. + double* psik = new double[meshK]; + double* jj = new double[mesh]; + for (int ik = 0; ik < meshK; ik++) + { + double q = ik * dk; + Mathzone::Spherical_Bessel(mesh, r, q, L, jj); + for (int ir = 0; ir < mesh; ir++) + { + function[ir] = psi[ir] * jj[ir]; + } + Mathzone::Simpson_Integral(mesh, function, rab, psik[ik]); + } + + for (int ik = 0; ik < meshK; ik++) + { + psik[ik] *= prefac_k; + } + + // + for (int ik = 0; ik < meshK; ik++) + { + functionk[ik] = std::pow(psik[ik] * ik * dk, 2); + } + Mathzone::Simpson_Integral(meshK, functionk, rabk, normk); + // cout << "\n Psi(k) norm=" << normk * fpi; + + if (normk * fpi < 0.999) + { + ecut = pow(dk * meshK, 2) * 2; + } + else + { + int kmesh_used = meshK; + while (normk * fpi > 0.999 && kmesh_used >= 0) + { + kmesh_used -= 2; + Mathzone::Simpson_Integral(kmesh_used, functionk, rabk, normk); + } + ecut = pow(dk * kmesh_used, 2) * 2; + } + // cout << "Ecut(Ry)(norm>0.999)=" << ecut << endl; + + delete[] jj; + delete[] psik; + delete[] function; + timer::tick("Psi_Second", "get_ecut"); + return ecut; } - -void Psi_Second::norm_c4( double *c4, const int &L) +void Psi_Second::norm_c4(double* c4, const int& L) { - ZEROS(this->psi,mesh); - for(int ie=0; iejle, mesh); -// Mathzone::Spherical_Bessel( this->mesh, this->r, this->eigen1[ie], L, jle); - - for(int ir=0; irpsi[ir] += c4[ie] * jle[L][ie][ir];} - } - - double* function = new double[mesh]; - for(int ir=0; ir0.0); - for(int ie=0; iepsi, mesh); + for (int ie = 0; ie < ne; ie++) + { + // ZEROS(this->jle, mesh); + // Mathzone::Spherical_Bessel( this->mesh, this->r, this->eigen1[ie], L, jle); + + for (int ir = 0; ir < mesh; ir++) + { + this->psi[ir] += c4[ie] * jle[L][ie][ir]; + } + } + + double* function = new double[mesh]; + for (int ir = 0; ir < mesh; ir++) + { + function[ir] = pow(psi[ir] * r[ir], 2); + } + + double norm; + Mathzone::Simpson_Integral(mesh, function, rab, norm); + delete[] function; + + assert(norm > 0.0); + for (int ie = 0; ie < ne; ie++) + { + c4[ie] /= std::sqrt(norm); + } + + return; } - - diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Psi_Second.h b/tools/01_NAO_generation/SIAB/src_spillage/Psi_Second.h index 0772cf0b0f6..49edd212df7 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Psi_Second.h +++ b/tools/01_NAO_generation/SIAB/src_spillage/Psi_Second.h @@ -5,59 +5,58 @@ class Psi_Second { - public: - Psi_Second(); - ~Psi_Second(); - - // be called in MultiZeta.cpp. - void init(ifstream &ifs, - const int &ne_in, - const int &lmax_in, - const double &rcut_in, - const int &total_nchi_in, - const double &tolerence); - - int calculate_oscillation(const int &L, const int &N, const double *c4, const int &ic); - - int oscillation; - - ofstream ofso; - ofstream ofsk; - - bool* below_max_os; - - int *os_number; - - void norm_c4( double *c4, const int &L); - // new method, mohan 2010-04-14 - double get_ecut( double *c4, const int &L); - - matrix eigenvalue; - matrix jjnorm; //mohan add 2010-04-18 - double*** jle; - - private: - - double dr; // about mesh - int mesh; // about mesh - double* r; // about mesh - double* rab; // about mesh - double rcut; // about mesh - - double* eigen1; - - int total_nchi; - int lmax; - int ne; - int test; - int start; - int max_os; - - double* psi; // psi - double* psi_first;//first derivative - double* psi_second; //second derivative - double* os_position; - int count; + public: + Psi_Second(); + ~Psi_Second(); + + // be called in MultiZeta.cpp. + void init(ifstream& ifs, + const int& ne_in, + const int& lmax_in, + const double& rcut_in, + const int& total_nchi_in, + const double& tolerence); + + int calculate_oscillation(const int& L, const int& N, const double* c4, const int& ic); + + int oscillation; + + ofstream ofso; + ofstream ofsk; + + bool* below_max_os; + + int* os_number; + + void norm_c4(double* c4, const int& L); + // new method, mohan 2010-04-14 + double get_ecut(double* c4, const int& L); + + matrix eigenvalue; + matrix jjnorm; // mohan add 2010-04-18 + double*** jle; + + private: + double dr; // about mesh + int mesh; // about mesh + double* r; // about mesh + double* rab; // about mesh + double rcut; // about mesh + + double* eigen1; + + int total_nchi; + int lmax; + int ne; + int test; + int start; + int max_os; + + double* psi; // psi + double* psi_first; // first derivative + double* psi_second; // second derivative + double* os_position; + int count; }; #endif diff --git a/tools/01_NAO_generation/SIAB/src_spillage/ReadData.cpp b/tools/01_NAO_generation/SIAB/src_spillage/ReadData.cpp index 7fbab7cd8b4..f45ef204d87 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/ReadData.cpp +++ b/tools/01_NAO_generation/SIAB/src_spillage/ReadData.cpp @@ -1,16 +1,17 @@ #include "ReadData.h" + #include "../src_parallel/parallel_common.h" -ReadData::ReadData() -: test(0) +ReadData::ReadData() : test(0) { - weight = new double[1]; + weight = new double[1]; } ReadData::~ReadData() { - if(TEST1) cout << "\n ~ReadData()" << endl; - delete[] weight; + if (TEST1) + cout << "\n ~ReadData()" << endl; + delete[] weight; } //========================================================== @@ -20,540 +21,524 @@ ReadData::~ReadData() // The number of ie is : ne // The number of k points : nks; //========================================================== -void ReadData::OverlapQandS(const string &name) +void ReadData::OverlapQandS(const string& name) { - TITLE(ofs_running, "ReadData", "OverlapQandS"); - - ifstream ifs; - - double* weighttmp = new double[NKSTOT]; - double* carkx = new double[NKSTOT]; - double* carky = new double[NKSTOT]; - double* carkz = new double[NKSTOT]; - ZEROS(weighttmp, NKSTOT); - - if(MY_RANK==0) - { - ifs.open(name.c_str()); - - if (!ifs) - { - cout << "\n Can't find file : " << name; - WARNING_QUIT("ReadData::OverlapQandS","Can't find file."); - } - else - { - cout << " FILE : " << name << endl; - } - - double lat0; - ifs >> lat0; - assert(lat0 == LAT0); - - double e11,e12,e13,e21,e22,e23,e31,e32,e33; - - ifs >> e11 >> e12 >> e13; - ifs >> e21 >> e22 >> e23; - ifs >> e31 >> e32 >> e33; - - assert( e11 == LATVEC.e11); - assert( e12 == LATVEC.e12); - assert( e13 == LATVEC.e13); - assert( e21 == LATVEC.e21); - assert( e22 == LATVEC.e22); - assert( e23 == LATVEC.e23); - assert( e31 == LATVEC.e31); - assert( e32 == LATVEC.e32); - assert( e33 == LATVEC.e33); - - int ntype; - READ_VALUE(ifs, ntype); // 1 - assert(ntype==NTYPE); - - double x,y,z; - for(int it=0; it> x >> y >> z; - } - } - - double tmp_ecutwfc; - double tmp_ecutwfc_jlq; - double tmp_rcut; - - READ_VALUE(ifs, tmp_ecutwfc);// 4 - READ_VALUE(ifs, tmp_ecutwfc_jlq);// 5 - READ_VALUE(ifs, tmp_rcut);// 6 - - assert(tmp_ecutwfc == ECUT); - assert(tmp_ecutwfc_jlq == ECUT_JLQ); - assert(tmp_rcut == RCUT); - - // mohan add 2009-08-28 - // make sure all structures have - // same SMOOTH and SIGMA value! - bool tmp_smooth; - double tmp_sigma; - double tmp_tolerence; - - READ_VALUE(ifs, tmp_smooth); // 7 - READ_VALUE(ifs, tmp_sigma); // 8 - READ_VALUE(ifs, tmp_tolerence); // 9 - - assert(tmp_smooth == SMOOTH); - assert(tmp_sigma == SIGMA); - assert(tmp_tolerence == TOLERENCE); - - int lmaxall, nkstot, nbands, nwfcall, ne; - - READ_VALUE(ifs, lmaxall); // 10 - READ_VALUE(ifs, nkstot);// 11 - READ_VALUE(ifs, nbands);// 12 - READ_VALUE(ifs, nwfcall);// 13 - READ_VALUE(ifs, ne);// 14 - - assert(lmaxall == LMAXALL); - assert(nkstot == NKSTOT); - assert(nbands == NBANDS); - assert(nwfcall == NWFCALL); - assert(ne == NE); - - //========================================== - // if use IBZ k-points, weight is different - // for each k point - //========================================== - // 0 means don't need to search from start. - if( SCAN_BEGIN(ifs,"",0) ) - { - ofs_running << " Kx/Ky/Kz/Kweight" ; - for(int ik=0; ik> carkx[ik] - >> carky[ik] - >> carkz[ik] - >> weighttmp[ik]; - ofs_running << " " << carkx[ik] - << " " << carky[ik] - << " " << carkz[ik] - << " " << weighttmp[ik] << endl; - } - SCAN_END(ifs,""); - } - - double sum = 0; - for(int ik=0; ik> lat0; + assert(lat0 == LAT0); + + double e11, e12, e13, e21, e22, e23, e31, e32, e33; + + ifs >> e11 >> e12 >> e13; + ifs >> e21 >> e22 >> e23; + ifs >> e31 >> e32 >> e33; + + assert(e11 == LATVEC.e11); + assert(e12 == LATVEC.e12); + assert(e13 == LATVEC.e13); + assert(e21 == LATVEC.e21); + assert(e22 == LATVEC.e22); + assert(e23 == LATVEC.e23); + assert(e31 == LATVEC.e31); + assert(e32 == LATVEC.e32); + assert(e33 == LATVEC.e33); + + int ntype; + READ_VALUE(ifs, ntype); // 1 + assert(ntype == NTYPE); + + double x, y, z; + for (int it = 0; it < ntype; it++) + { + string tmp_label; + READ_VALUE(ifs, tmp_label); // 2 + assert(tmp_label == LABEL[it]); + + int na; + READ_VALUE(ifs, na); // 3 + assert(na == NA[it]); + for (int ia = 0; ia < NA[it]; ia++) + { + // mohan fix bug 2010-09-27 + // x,y,z is useless. + // And CARPOSX, CASPOSY, CASPOSZ + // only contains the first structure. + ifs >> x >> y >> z; + } + } + + double tmp_ecutwfc; + double tmp_ecutwfc_jlq; + double tmp_rcut; + + READ_VALUE(ifs, tmp_ecutwfc); // 4 + READ_VALUE(ifs, tmp_ecutwfc_jlq); // 5 + READ_VALUE(ifs, tmp_rcut); // 6 + + assert(tmp_ecutwfc == ECUT); + assert(tmp_ecutwfc_jlq == ECUT_JLQ); + assert(tmp_rcut == RCUT); + + // mohan add 2009-08-28 + // make sure all structures have + // same SMOOTH and SIGMA value! + bool tmp_smooth; + double tmp_sigma; + double tmp_tolerence; + + READ_VALUE(ifs, tmp_smooth); // 7 + READ_VALUE(ifs, tmp_sigma); // 8 + READ_VALUE(ifs, tmp_tolerence); // 9 + + assert(tmp_smooth == SMOOTH); + assert(tmp_sigma == SIGMA); + assert(tmp_tolerence == TOLERENCE); + + int lmaxall, nkstot, nbands, nwfcall, ne; + + READ_VALUE(ifs, lmaxall); // 10 + READ_VALUE(ifs, nkstot); // 11 + READ_VALUE(ifs, nbands); // 12 + READ_VALUE(ifs, nwfcall); // 13 + READ_VALUE(ifs, ne); // 14 + + assert(lmaxall == LMAXALL); + assert(nkstot == NKSTOT); + assert(nbands == NBANDS); + assert(nwfcall == NWFCALL); + assert(ne == NE); + + //========================================== + // if use IBZ k-points, weight is different + // for each k point + //========================================== + // 0 means don't need to search from start. + if (SCAN_BEGIN(ifs, "", 0)) + { + ofs_running << " Kx/Ky/Kz/Kweight"; + for (int ik = 0; ik < NKSTOT; ik++) + { + if (ik % 4 == 0) + ofs_running << endl; + ifs >> carkx[ik] >> carky[ik] >> carkz[ik] >> weighttmp[ik]; + ofs_running << " " << carkx[ik] << " " << carky[ik] << " " << carkz[ik] << " " << weighttmp[ik] << endl; + } + SCAN_END(ifs, ""); + } + + double sum = 0; + for (int ik = 0; ik < NKSTOT; ik++) + sum += weighttmp[ik]; + + // ofs_running << "\n sum kpoint weight = " << sum; + if (abs(sum - 1.0) > 1.0e-5) + { + ofs_running << "\n sum of k weight = " << sum << endl; + WARNING_QUIT("ReadData::OverlapQandS", "sum of k weight is wrong."); + } + } #ifdef __MPI - Parallel_Common::bcast_double(weighttmp, NKSTOT); - Parallel_Common::bcast_double(carkx, NKSTOT); - Parallel_Common::bcast_double(carky, NKSTOT); - Parallel_Common::bcast_double(carkz, NKSTOT); + Parallel_Common::bcast_double(weighttmp, NKSTOT); + Parallel_Common::bcast_double(carkx, NKSTOT); + Parallel_Common::bcast_double(carky, NKSTOT); + Parallel_Common::bcast_double(carkz, NKSTOT); #endif - ofs_running << " NKS=" << NKS << endl; - - - delete[] this->weight; - this->weight = new double[NKS]; - ZEROS(weight, NKS); + ofs_running << " NKS=" << NKS << endl; - CARKX = new double[NKS]; - CARKY = new double[NKS]; - CARKZ = new double[NKS]; - ZEROS(CARKX, NKS); - ZEROS(CARKY, NKS); - ZEROS(CARKZ, NKS); + delete[] this->weight; + this->weight = new double[NKS]; + ZEROS(weight, NKS); + CARKX = new double[NKS]; + CARKY = new double[NKS]; + CARKZ = new double[NKS]; + ZEROS(CARKX, NKS); + ZEROS(CARKY, NKS); + ZEROS(CARKZ, NKS); - for(int ik=0; ikweight[iknow] = weighttmp[ik]; - CARKX[iknow] = carkx[ik]; - CARKY[iknow] = carky[ik]; - CARKZ[iknow] = carkz[ik]; - //ofs_running << "\n iknow = " << iknow << " ik = " << ik << " weight = " << weight[iknow]; - } + const int iknow = ik - Pkpoints.startk_pool[MY_POOL]; + const int pool = Pkpoints.whichpool[ik]; + if (MY_POOL == pool) + { + this->weight[iknow] = weighttmp[ik]; + CARKX[iknow] = carkx[ik]; + CARKY[iknow] = carky[ik]; + CARKZ[iknow] = carkz[ik]; + // ofs_running << "\n iknow = " << iknow << " ik = " << ik << " weight = " << weight[iknow]; + } #else - weight[ik] = weighttmp[ik]; - CARKX[ik] = carkx[ik]; - CARKY[ik] = carky[ik]; - CARKZ[ik] = carkz[ik]; + weight[ik] = weighttmp[ik]; + CARKX[ik] = carkx[ik]; + CARKY[ik] = carky[ik]; + CARKZ[ik] = carkz[ik]; #endif - } - - - delete[] weighttmp; - delete[] carkx; - delete[] carky; - delete[] carkz; - - + } - this->OverlapQ(ifs); + delete[] weighttmp; + delete[] carkx; + delete[] carky; + delete[] carkz; + this->OverlapQ(ifs); + this->OverlapSq1q2(ifs); - this->OverlapSq1q2(ifs); - - - - return; + return; } -void ReadData::OverlapSq1q2(ifstream &ifs) +void ReadData::OverlapSq1q2(ifstream& ifs) { -// TITLE("ReadData", "OverlapSq1q2"); - timer::tick("ReadData", "OverlapSq1q2"); - - USEPW = true; - if(MY_RANK==0) - { - bool restart = false; - bool quit = false; - if( SCAN_BEGIN(ifs,"",restart,quit) ) - { - USEPW = false; - } - } + // TITLE("ReadData", "OverlapSq1q2"); + timer::tick("ReadData", "OverlapSq1q2"); + + USEPW = true; + if (MY_RANK == 0) + { + bool restart = false; + bool quit = false; + if (SCAN_BEGIN(ifs, "", restart, quit)) + { + USEPW = false; + } + } #ifdef __MPI - Parallel_Common::bcast_bool(USEPW); + Parallel_Common::bcast_bool(USEPW); #endif - - if(USEPW) - { - ofs_running << " USE PLANE WAVE BASIS" << endl; - } - else - { - ofs_running << " USE Q and S matrix" << endl; - } - - if(CALSPI==0) return; - - if(!USEPW) - { - this->Sq1q2 = new ComplexArray[ NKS ]; - for(int ik=0; ikSq1q2[ik].create(NWFCALL, NWFCALL, NE, NE); - } - } + if (USEPW) + { + ofs_running << " USE PLANE WAVE BASIS" << endl; + } + else + { + ofs_running << " USE Q and S matrix" << endl; + } + + if (CALSPI == 0) + return; + + if (!USEPW) + { + this->Sq1q2 = new ComplexArray[NKS]; + for (int ik = 0; ik < NKS; ik++) + { + this->Sq1q2[ik].create(NWFCALL, NWFCALL, NE, NE); + } + } #ifdef __MPI - if(!USEPW) - { - const int ndata = NWFCALL * NWFCALL * NE * NE; - double* buffera = new double[ndata]; - double* bufferb = new double[ndata]; - for(int ik=0; ik> buffera[id] >> bufferb[id]; - } - } - //ofs_running << "\n ik=" << ik; - MPI_Status ierror; - const int pool = Pkpoints.whichpool[ik]; - const int iknow = ik - Pkpoints.startk_pool[MY_POOL]; - const int startdata = 0; // this is differnt from Q. - - //ofs_running << " this k point belong to pool : " << pool; - //ofs_running << " my pool : " << MY_POOL << endl; - - if(MY_RANK==0) - { - // always ready to send data to other processors. - if(pool==0) - { - // send copies (nproc - 1). - for(int ip=1; ipSq1q2[iknow].ptr[id + startdata] = std::complex (buffera[id], bufferb[id]); - } - // ofs_running << "\n local data copy "; - } - else - { - // send copys: nproc - for(int ip=0; ipSq1q2[iknow].ptr[id + startdata] = std::complex (buffera[id], bufferb[id]); - } - // ofs_running << "\n receive data end" << endl; - } - else - { - // ofs_running << "\n do nothing"; - } - } - MPI_Barrier(MPI_COMM_WORLD); - } - delete[] buffera; - delete[] bufferb; - }// end begin + if (!USEPW) + { + const int ndata = NWFCALL * NWFCALL * NE * NE; + double* buffera = new double[ndata]; + double* bufferb = new double[ndata]; + for (int ik = 0; ik < NKSTOT; ik++) + { + ZEROS(buffera, ndata); + ZEROS(bufferb, ndata); + if (MY_RANK == 0) + { + for (int id = 0; id < ndata; id++) + { + ifs >> buffera[id] >> bufferb[id]; + } + } + // ofs_running << "\n ik=" << ik; + MPI_Status ierror; + const int pool = Pkpoints.whichpool[ik]; + const int iknow = ik - Pkpoints.startk_pool[MY_POOL]; + const int startdata = 0; // this is differnt from Q. + + // ofs_running << " this k point belong to pool : " << pool; + // ofs_running << " my pool : " << MY_POOL << endl; + + if (MY_RANK == 0) + { + // always ready to send data to other processors. + if (pool == 0) + { + // send copies (nproc - 1). + for (int ip = 1; ip < Pkpoints.nproc_pool[pool]; ip++) + { + MPI_Send(buffera, ndata, MPI_DOUBLE, ip, ik, MPI_COMM_WORLD); + MPI_Send(bufferb, ndata, MPI_DOUBLE, ip, ik, MPI_COMM_WORLD); + // ofs_running << "\n send data to processor " << ip; + } + + // local copy + for (int id = 0; id < ndata; id++) + { + this->Sq1q2[iknow].ptr[id + startdata] = std::complex(buffera[id], bufferb[id]); + } + // ofs_running << "\n local data copy "; + } + else + { + // send copys: nproc + for (int ip = 0; ip < Pkpoints.nproc_pool[pool]; ip++) + { + const int ipall = ip + Pkpoints.startpro_pool[pool]; + MPI_Send(buffera, ndata, MPI_DOUBLE, ipall, ik, MPI_COMM_WORLD); + MPI_Send(bufferb, ndata, MPI_DOUBLE, ipall, ik, MPI_COMM_WORLD); + // ofs_running << "\n send data to processor " << ipall; + } + } + } // end MY_RANK==0 + else + { + if (pool == MY_POOL) + { + // ofs_running << "\n begin receive, ndata : " << ndata << endl; + MPI_Recv(buffera, ndata, MPI_DOUBLE, 0, ik, MPI_COMM_WORLD, &ierror); + MPI_Recv(bufferb, ndata, MPI_DOUBLE, 0, ik, MPI_COMM_WORLD, &ierror); + + for (int id = 0; id < ndata; id++) + { + this->Sq1q2[iknow].ptr[id + startdata] = std::complex(buffera[id], bufferb[id]); + } + // ofs_running << "\n receive data end" << endl; + } + else + { + // ofs_running << "\n do nothing"; + } + } + MPI_Barrier(MPI_COMM_WORLD); + } + delete[] buffera; + delete[] bufferb; + } // end begin #else - assert(NKS==NKSTOT); - if(!USEPW) - { - double a,b; - for (int ik = 0; ik < NKSTOT; ik++) - { - this->Sq1q2[ik].create(NWFCALL, NWFCALL, NE, NE); - - for (int i = 0; i < Sq1q2[ik].getSize(); i++) - { - ifs >> a >> b; - Sq1q2[ik].ptr[i] = std::complex(a,b); - } - } - } + assert(NKS == NKSTOT); + if (!USEPW) + { + double a, b; + for (int ik = 0; ik < NKSTOT; ik++) + { + this->Sq1q2[ik].create(NWFCALL, NWFCALL, NE, NE); + + for (int i = 0; i < Sq1q2[ik].getSize(); i++) + { + ifs >> a >> b; + Sq1q2[ik].ptr[i] = std::complex(a, b); + } + } + } #endif -// SCAN_END(ifs,""); + // SCAN_END(ifs,""); - timer::tick("ReadData", "OverlapSq1q2"); - return; + timer::tick("ReadData", "OverlapSq1q2"); + return; } -void ReadData::OverlapQ(ifstream &ifs) +void ReadData::OverlapQ(ifstream& ifs) { - if(CALSPI==0) return; -// TITLE("ReadData", "OverlapQ"); - timer::tick("ReadData", "OverlapQ"); - - assert(NKS > 0); - assert(NBANDS > 0); - assert(NWFCALL > 0); - assert(NE > 0); - - this->Qin.create(NKS, NBANDS, NWFCALL, NE); - - ofs_running << " DIMENSION OF Q=" << endl; - OUT(ofs_running,"NKS",NKS); - OUT(ofs_running,"NBANDS",NBANDS); - OUT(ofs_running,"NWFCALL",NWFCALL); - OUT(ofs_running,"NE",NE); - - bool begin = false; - - if(MY_RANK ==0 ) - { - if( SCAN_BEGIN(ifs,"",0) ) - { - begin = true; - } - } + if (CALSPI == 0) + return; + // TITLE("ReadData", "OverlapQ"); + timer::tick("ReadData", "OverlapQ"); + + assert(NKS > 0); + assert(NBANDS > 0); + assert(NWFCALL > 0); + assert(NE > 0); + + this->Qin.create(NKS, NBANDS, NWFCALL, NE); + + ofs_running << " DIMENSION OF Q=" << endl; + OUT(ofs_running, "NKS", NKS); + OUT(ofs_running, "NBANDS", NBANDS); + OUT(ofs_running, "NWFCALL", NWFCALL); + OUT(ofs_running, "NE", NE); + + bool begin = false; + + if (MY_RANK == 0) + { + if (SCAN_BEGIN(ifs, "", 0)) + { + begin = true; + } + } #ifdef __MPI - Parallel_Common::bcast_bool(begin); + Parallel_Common::bcast_bool(begin); #endif - - - if(begin) - { - double a,b; + if (begin) + { + double a, b; #ifdef __MPI - const int ndata = NBANDS * NWFCALL * NE; - double* buffera = new double[ndata]; - double* bufferb = new double[ndata]; - for(int ik=0; ik> buffera[id] >> bufferb[id]; - } - } - - //ofs_running << "\n ik=" << ik; - - MPI_Status ierror; - const int pool = Pkpoints.whichpool[ik]; - const int iknow = ik - Pkpoints.startk_pool[MY_POOL]; - const int startdata = iknow * ndata; - - //ofs_running << "\n this k point belong to pool : " << pool; - //ofs_running << "\n my pool : " << MY_POOL; - - if(MY_RANK==0) - { - // always ready to send data to other processors. - if(pool==0) - { - // send copies (nproc - 1). - for(int ip=1; ipQin.ptr[id + startdata] = std::complex (buffera[id], bufferb[id]); - } - //ofs_running << "\n local data copy "; - } - else - { - // send copys: nproc - for(int ip=0; ipQin.ptr[id + startdata] = std::complex (buffera[id], bufferb[id]); - } - //ofs_running << "\n receive data end" << endl; - } - else - { - //ofs_running << "\n do nothing"; - } - } - - MPI_Barrier(MPI_COMM_WORLD); - }// endik - delete[] buffera; - delete[] bufferb; + const int ndata = NBANDS * NWFCALL * NE; + double* buffera = new double[ndata]; + double* bufferb = new double[ndata]; + for (int ik = 0; ik < NKSTOT; ik++) + { + ZEROS(buffera, ndata); + ZEROS(bufferb, ndata); + if (MY_RANK == 0) + { + for (int id = 0; id < ndata; id++) + { + ifs >> buffera[id] >> bufferb[id]; + } + } + + // ofs_running << "\n ik=" << ik; + + MPI_Status ierror; + const int pool = Pkpoints.whichpool[ik]; + const int iknow = ik - Pkpoints.startk_pool[MY_POOL]; + const int startdata = iknow * ndata; + + // ofs_running << "\n this k point belong to pool : " << pool; + // ofs_running << "\n my pool : " << MY_POOL; + + if (MY_RANK == 0) + { + // always ready to send data to other processors. + if (pool == 0) + { + // send copies (nproc - 1). + for (int ip = 1; ip < Pkpoints.nproc_pool[pool]; ip++) + { + MPI_Send(buffera, ndata, MPI_DOUBLE, ip, ik, MPI_COMM_WORLD); + MPI_Send(bufferb, ndata, MPI_DOUBLE, ip, ik, MPI_COMM_WORLD); + // ofs_running << "\n send data to processor " << ip; + } + + // local copy + for (int id = 0; id < ndata; id++) + { + this->Qin.ptr[id + startdata] = std::complex(buffera[id], bufferb[id]); + } + // ofs_running << "\n local data copy "; + } + else + { + // send copys: nproc + for (int ip = 0; ip < Pkpoints.nproc_pool[pool]; ip++) + { + const int ipall = ip + Pkpoints.startpro_pool[pool]; + MPI_Send(buffera, ndata, MPI_DOUBLE, ipall, ik, MPI_COMM_WORLD); + MPI_Send(bufferb, ndata, MPI_DOUBLE, ipall, ik, MPI_COMM_WORLD); + // ofs_running << "\n send data to processor " << ipall; + } + } + } // end MY_RANK==0 + else + { + if (pool == MY_POOL) + { + // ofs_running << "\n begin receive, ndata : " << ndata << endl; + MPI_Recv(buffera, ndata, MPI_DOUBLE, 0, ik, MPI_COMM_WORLD, &ierror); + MPI_Recv(bufferb, ndata, MPI_DOUBLE, 0, ik, MPI_COMM_WORLD, &ierror); + for (int id = 0; id < ndata; id++) + { + this->Qin.ptr[id + startdata] = std::complex(buffera[id], bufferb[id]); + } + // ofs_running << "\n receive data end" << endl; + } + else + { + // ofs_running << "\n do nothing"; + } + } + + MPI_Barrier(MPI_COMM_WORLD); + } // endik + delete[] buffera; + delete[] bufferb; #else - for(int i=0; i> a >> b; - this->Qin.ptr[i] = std::complex( a, b ); - } + for (int i = 0; i < Qin.getSize(); i++) + { + ifs >> a >> b; + this->Qin.ptr[i] = std::complex(a, b); + } #endif - } - + } + // SCAN_END(ifs,""); - //SCAN_END(ifs,""); - - timer::tick("ReadData", "OverlapQ"); + timer::tick("ReadData", "OverlapQ"); - return; + return; } //========================================================== // Readin inverse S matrix. // Sinv is fail to use if C4 change. //========================================================== -void ReadData::OverlapSinv(const string &name) +void ReadData::OverlapSinv(const string& name) { -// TITLE("ReadData", "OverlapSinv"); - ifstream ifs(name.c_str()); - - if (!ifs) - { - cout << "\n Can't find file : " << name; - cout << "\n But maybe we don't need it."; - return; - } - - assert(NKS > 0); - assert(NWFCALL > 0); - - this->Sinv.create(NKS, NWFCALL, NWFCALL); - - return; - - string tmp; - int row, col; - - for (int i = 0; i < NKS; i++) - { - ifs >> tmp >> row >> col; -// cout << "\n" << tmp << " row=" << row << " col=" << col << endl; - assert(row == NWFCALL); - assert(col == NWFCALL); - - for (int j = 0; j < row; j++) - { - for (int k = 0; k < col; k++) - { - double a,b; - ifs >> a >> b; - Sinv(i,j,k) = std::complex(a,b); - } - } - } - ifs.close(); - return; + // TITLE("ReadData", "OverlapSinv"); + ifstream ifs(name.c_str()); + + if (!ifs) + { + cout << "\n Can't find file : " << name; + cout << "\n But maybe we don't need it."; + return; + } + + assert(NKS > 0); + assert(NWFCALL > 0); + + this->Sinv.create(NKS, NWFCALL, NWFCALL); + + return; + + string tmp; + int row, col; + + for (int i = 0; i < NKS; i++) + { + ifs >> tmp >> row >> col; + // cout << "\n" << tmp << " row=" << row << " col=" << col << endl; + assert(row == NWFCALL); + assert(col == NWFCALL); + + for (int j = 0; j < row; j++) + { + for (int k = 0; k < col; k++) + { + double a, b; + ifs >> a >> b; + Sinv(i, j, k) = std::complex(a, b); + } + } + } + ifs.close(); + return; } diff --git a/tools/01_NAO_generation/SIAB/src_spillage/ReadData.h b/tools/01_NAO_generation/SIAB/src_spillage/ReadData.h index 4e4d356d4b0..1148c93677b 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/ReadData.h +++ b/tools/01_NAO_generation/SIAB/src_spillage/ReadData.h @@ -5,36 +5,34 @@ class ReadData { - public: - - ComplexArray Qin; - ComplexArray *Sq1q2; - ComplexArray Sinv; - - //=================================== - // read in data Q= , S= - //=================================== - void OverlapQandS( const string &name); - - //======================================================== - // inverse matrix of S directly, - // used for checking the correctness - // of spillage calculation - //======================================================== - void OverlapSinv( const string &name); - - double *weight; - - private: - void OverlapQ( ifstream &ifs); - void OverlapSq1q2( ifstream &ifs ); - - int test; - - public: - ReadData(); - ~ReadData(); - + public: + ComplexArray Qin; + ComplexArray* Sq1q2; + ComplexArray Sinv; + + //=================================== + // read in data Q= , S= + //=================================== + void OverlapQandS(const string& name); + + //======================================================== + // inverse matrix of S directly, + // used for checking the correctness + // of spillage calculation + //======================================================== + void OverlapSinv(const string& name); + + double* weight; + + private: + void OverlapQ(ifstream& ifs); + void OverlapSq1q2(ifstream& ifs); + + int test; + + public: + ReadData(); + ~ReadData(); }; #endif diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Simulated_Annealing_Orbital.h b/tools/01_NAO_generation/SIAB/src_spillage/Simulated_Annealing_Orbital.h index 6fed0be4bcb..9c8beb5c11c 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Simulated_Annealing_Orbital.h +++ b/tools/01_NAO_generation/SIAB/src_spillage/Simulated_Annealing_Orbital.h @@ -3,154 +3,167 @@ class Metropolis; -class Metropolis::Simulated_Annealing_Orbital: public Simulated_Annealing +class Metropolis::Simulated_Annealing_Orbital : public Simulated_Annealing { - public: - Simulated_Annealing_Orbital( Metropolis *metropolis_in, SpillageStep &Level_in) - :Simulated_Annealing(), - metropolis(metropolis_in), - Level(Level_in){} - protected: - virtual void preprocess_all_temperature() - { - metropolis->preprocess_all_temperature(Level); - } - virtual void preprocess_each_temperature(size_t itemperature) - { - metropolis->Psi2.ofso << "spillage_ntemp = " << itemperature << " kappa = " << metropolis->kappa << endl; - } - virtual void update_temperature(size_t itemperature) - { - metropolis->update_t(); // change temperature - } - virtual void preprocess_each_step( size_t itemperature, size_t istep) - { - metropolis->small_jql(); - } - virtual void cal_each_step( size_t itemperature, size_t istep) - { - preprocess_each_step_beyondloop(itemperature, istep); - // for each (atom_type,L,N), we need a set of C4 - // to describe the radial wave functions. - // after each radial wave function is update(ne parameters), - // we call it 'a step under a particular temperature is done' - for (int ic = 0; ic < metropolis->nchi; ic++) - { - metropolis->info.set_ic(ic).unpack_ic(Level); - const int T = Level.wayc[ic].type; // get the type index of this radial wave function. - if(Level.info[T].state=="skip") continue; // if "skip", the C4 will not changed in this 'Big step' // mohan add 2009-01-27 - for (int ie = 0; ie < Level.ne; ie++) - { - metropolis->info.set_ie(ie); - try - { - this->Simulated_Annealing::cal_each_step(itemperature, istep); - } - catch(Small_Jlq &e) { continue; } - } - } - reprocess_each_step_beyondloop(itemperature, istep); - } - virtual void preprocess_each_step_beyondloop(size_t itemperature, size_t istep){} - virtual void reprocess_each_step_beyondloop(size_t itemperature, size_t istep) - { - metropolis->ofs_1(istep, this->Level); - metropolis->norm_C4(istep, this->Level, this->step_num); - } - virtual void move_variable(size_t itemperature, size_t istep) - { - metropolis->trial_c4(); - } - virtual void accept_process( double new_function, size_t itemperature, size_t istep) - { - metropolis->accept_process(this->Level); - } - virtual void reject_process( double new_function, size_t itemperature, size_t istep) - { - metropolis->reject_process(); - } - Metropolis *metropolis; - SpillageStep &Level; + public: + Simulated_Annealing_Orbital(Metropolis* metropolis_in, SpillageStep& Level_in) + : Simulated_Annealing(), metropolis(metropolis_in), Level(Level_in) + { + } + + protected: + virtual void preprocess_all_temperature() + { + metropolis->preprocess_all_temperature(Level); + } + virtual void preprocess_each_temperature(size_t itemperature) + { + metropolis->Psi2.ofso << "spillage_ntemp = " << itemperature << " kappa = " << metropolis->kappa << endl; + } + virtual void update_temperature(size_t itemperature) + { + metropolis->update_t(); // change temperature + } + virtual void preprocess_each_step(size_t itemperature, size_t istep) + { + metropolis->small_jql(); + } + virtual void cal_each_step(size_t itemperature, size_t istep) + { + preprocess_each_step_beyondloop(itemperature, istep); + // for each (atom_type,L,N), we need a set of C4 + // to describe the radial wave functions. + // after each radial wave function is update(ne parameters), + // we call it 'a step under a particular temperature is done' + for (int ic = 0; ic < metropolis->nchi; ic++) + { + metropolis->info.set_ic(ic).unpack_ic(Level); + const int T = Level.wayc[ic].type; // get the type index of this radial wave function. + if (Level.info[T].state == "skip") + continue; // if "skip", the C4 will not changed in this 'Big step' // mohan add 2009-01-27 + for (int ie = 0; ie < Level.ne; ie++) + { + metropolis->info.set_ie(ie); + try + { + this->Simulated_Annealing::cal_each_step(itemperature, istep); + } + catch (Small_Jlq& e) + { + continue; + } + } + } + reprocess_each_step_beyondloop(itemperature, istep); + } + virtual void preprocess_each_step_beyondloop(size_t itemperature, size_t istep) + { + } + virtual void reprocess_each_step_beyondloop(size_t itemperature, size_t istep) + { + metropolis->ofs_1(istep, this->Level); + metropolis->norm_C4(istep, this->Level, this->step_num); + } + virtual void move_variable(size_t itemperature, size_t istep) + { + metropolis->trial_c4(); + } + virtual void accept_process(double new_function, size_t itemperature, size_t istep) + { + metropolis->accept_process(this->Level); + } + virtual void reject_process(double new_function, size_t itemperature, size_t istep) + { + metropolis->reject_process(); + } + Metropolis* metropolis; + SpillageStep& Level; }; -class Metropolis::Simulated_Annealing_Orbital_Spillage: public Simulated_Annealing_Orbital +class Metropolis::Simulated_Annealing_Orbital_Spillage : public Simulated_Annealing_Orbital { - public: - Simulated_Annealing_Orbital_Spillage( Metropolis *metropolis_in, SpillageStep &Level_in) - :Simulated_Annealing_Orbital(metropolis_in, Level_in){} - private: - virtual void preprocess_each_temperature(size_t itemperature) - { - metropolis->init_spi_states(itemperature); - Simulated_Annealing_Orbital::preprocess_each_temperature(itemperature); - } - virtual void reprocess_each_temperature(size_t itemperature) - { - metropolis->cout_1(itemperature); - } - virtual void preprocess_each_step_beyondloop(size_t itemperature, size_t istep) - { - metropolis->file_finish_metropolis(istep); - metropolis->change_accept_rate(istep,Level); - metropolis->info.set_update_number(0); // each step we renew every Jlq one by one. - } - virtual void move_variable(size_t itemperature, size_t istep) - { - Simulated_Annealing_Orbital::move_variable(itemperature, istep); - set_temperature( metropolis->temperature ); - } - virtual double cal_new_function(size_t itemperature, size_t istep) - { - set_old_function(input.SV.cal_defined_value(0)); - return metropolis->cal_spillage(this->Level); - } + public: + Simulated_Annealing_Orbital_Spillage(Metropolis* metropolis_in, SpillageStep& Level_in) + : Simulated_Annealing_Orbital(metropolis_in, Level_in) + { + } + + private: + virtual void preprocess_each_temperature(size_t itemperature) + { + metropolis->init_spi_states(itemperature); + Simulated_Annealing_Orbital::preprocess_each_temperature(itemperature); + } + virtual void reprocess_each_temperature(size_t itemperature) + { + metropolis->cout_1(itemperature); + } + virtual void preprocess_each_step_beyondloop(size_t itemperature, size_t istep) + { + metropolis->file_finish_metropolis(istep); + metropolis->change_accept_rate(istep, Level); + metropolis->info.set_update_number(0); // each step we renew every Jlq one by one. + } + virtual void move_variable(size_t itemperature, size_t istep) + { + Simulated_Annealing_Orbital::move_variable(itemperature, istep); + set_temperature(metropolis->temperature); + } + virtual double cal_new_function(size_t itemperature, size_t istep) + { + set_old_function(input.SV.cal_defined_value(0)); + return metropolis->cal_spillage(this->Level); + } }; -class Metropolis::Simulated_Annealing_Orbital_Kinetic: public Simulated_Annealing_Orbital +class Metropolis::Simulated_Annealing_Orbital_Kinetic : public Simulated_Annealing_Orbital { - public: - Simulated_Annealing_Orbital_Kinetic( Metropolis *metropolis_in, SpillageStep &Level_in) - :Simulated_Annealing_Orbital(metropolis_in, Level_in){} - private: - virtual void preprocess_each_temperature(size_t itemperature) - { - metropolis->init_kin_states(); - Simulated_Annealing_Orbital::preprocess_each_temperature(itemperature); - metropolis->init_temperature_kinetic(itemperature, this->Level); - } - virtual void reprocess_each_temperature(size_t itemperature) - { - metropolis->cout_1(itemperature); - metropolis->cout_2(itemperature, this->Level); - } - virtual void preprocess_each_step_beyondloop(size_t itemperature, size_t istep) - { - metropolis->file_finish_metropolis(istep); - metropolis->reset_high_kinetic_energy(itemperature, istep); - metropolis->change_accept_rate(istep, this->Level); - metropolis->info.update_number = 0; // each step we renew every Jlq one by one. - } - virtual void move_variable(size_t itemperature, size_t istep) - { - Simulated_Annealing_Orbital::move_variable(itemperature, istep); - set_temperature( metropolis->kin_temp[metropolis->info.ic] ); - } - virtual bool must_reject_before(double new_function, size_t itemperature, size_t istep) - { - // make sure the new spillage value is in the range. - // otherwise, we don't even calculate the kinetic energy. - return metropolis->cal_spillage(this->Level) >= metropolis->min_spillage * metropolis->kappa; - } - virtual double cal_new_function(size_t itemperature, size_t istep) - { - set_old_function(metropolis->kin_last[metropolis->info.ic]); - return metropolis->cal_kinetic_energy(this->Level); - } - virtual void accept_process( double new_function, size_t itemperature, size_t istep) - { - Simulated_Annealing_Orbital::accept_process( new_function, itemperature, istep); - metropolis->kin_last[metropolis->info.ic] = new_function; - } + public: + Simulated_Annealing_Orbital_Kinetic(Metropolis* metropolis_in, SpillageStep& Level_in) + : Simulated_Annealing_Orbital(metropolis_in, Level_in) + { + } + + private: + virtual void preprocess_each_temperature(size_t itemperature) + { + metropolis->init_kin_states(); + Simulated_Annealing_Orbital::preprocess_each_temperature(itemperature); + metropolis->init_temperature_kinetic(itemperature, this->Level); + } + virtual void reprocess_each_temperature(size_t itemperature) + { + metropolis->cout_1(itemperature); + metropolis->cout_2(itemperature, this->Level); + } + virtual void preprocess_each_step_beyondloop(size_t itemperature, size_t istep) + { + metropolis->file_finish_metropolis(istep); + metropolis->reset_high_kinetic_energy(itemperature, istep); + metropolis->change_accept_rate(istep, this->Level); + metropolis->info.update_number = 0; // each step we renew every Jlq one by one. + } + virtual void move_variable(size_t itemperature, size_t istep) + { + Simulated_Annealing_Orbital::move_variable(itemperature, istep); + set_temperature(metropolis->kin_temp[metropolis->info.ic]); + } + virtual bool must_reject_before(double new_function, size_t itemperature, size_t istep) + { + // make sure the new spillage value is in the range. + // otherwise, we don't even calculate the kinetic energy. + return metropolis->cal_spillage(this->Level) >= metropolis->min_spillage * metropolis->kappa; + } + virtual double cal_new_function(size_t itemperature, size_t istep) + { + set_old_function(metropolis->kin_last[metropolis->info.ic]); + return metropolis->cal_kinetic_energy(this->Level); + } + virtual void accept_process(double new_function, size_t itemperature, size_t istep) + { + Simulated_Annealing_Orbital::accept_process(new_function, itemperature, istep); + metropolis->kin_last[metropolis->info.ic] = new_function; + } }; #endif \ No newline at end of file diff --git a/tools/01_NAO_generation/SIAB/src_spillage/SpillageStep.cpp b/tools/01_NAO_generation/SIAB/src_spillage/SpillageStep.cpp index cbddc11550d..ada14cb052b 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/SpillageStep.cpp +++ b/tools/01_NAO_generation/SIAB/src_spillage/SpillageStep.cpp @@ -1,13 +1,14 @@ #include "SpillageStep.h" + #include "../src_parallel/parallel_reduce.h" -#include "tools.h" -#include "Inverse_Matrix_S.h" -#include "Inverse_Matrix_S_inline.h" #include "../src_tools/complexmatrix_inline.h" #include "../src_unittest/src_tools/common_test.h" +#include "Inverse_Matrix_S.h" +#include "Inverse_Matrix_S_inline.h" +#include "tools.h" #ifdef _OPENMP -#include +#include #endif SpillageStep::SpillageStep() @@ -26,8 +27,7 @@ SpillageStep::~SpillageStep() delete[] wayd; } - -void SpillageStep::set_info(const int &ntype_in) +void SpillageStep::set_info(const int& ntype_in) { this->ntype = ntype_in; cout << "\n Element type = " << ntype; @@ -39,14 +39,14 @@ void SpillageStep::set_info(const int &ntype_in) void SpillageStep::set_nwfc() { - TITLE("SpillageStep","set_nwfc"); - assert( ntype > 0); + TITLE("SpillageStep", "set_nwfc"); + assert(ntype > 0); //=============================== // (1) generate nwfc, nwfc2 //=============================== this->nwfc = 0; this->nwfc2 = 0; - for (int i=0; iwayd[iw2].ic = ic; - if (test==2) cout << "\n iw2=" << iw2 << " ic=" << ic; + if (test == 2) + cout << "\n iw2=" << iw2 << " ic=" << ic; ++iw2; - }// end m + } // end m } else if (info[it].fa[l] == "a") { @@ -167,29 +170,29 @@ void SpillageStep::set_nwfc() ++iw2; } ++ic; - }// end n - }// end l + } // end n + } // end l nw_this_type = iw2 - index_copy; } else { - for (int j=index_copy; jwayd[iw2].ic = this->wayd[j].ic; - if (test==2) cout << "\n iw2=" << iw2 << " ic=" << wayd[iw2].ic; + if (test == 2) + cout << "\n iw2=" << iw2 << " ic=" << wayd[iw2].ic; iw2++; } } } } - //cout << "\n iw2 = " << iw2; - //cout << "\n nwfc2 = " << nwfc2; - assert(iw2==nwfc2); + // cout << "\n iw2 = " << iw2; + // cout << "\n nwfc2 = " << nwfc2; + assert(iw2 == nwfc2); return; } - void SpillageStep::set_iw00(void) { //========================= @@ -198,30 +201,30 @@ void SpillageStep::set_iw00(void) int iw00 = 0; int iw2 = 0; - for (int it=0; it not equal."); + WARNING_QUIT("SpillageStep::set_iw00", "in not equal."); } - //cout << " type=" << it << " lmax=" << info[it].lmax << endl; - for (int i=0; i< NA[it]; i++) + // cout << " type=" << it << " lmax=" << info[it].lmax << endl; + for (int i = 0; i < NA[it]; i++) { - //for (int l=0; l<= mz.lmax_type[it]; l++) - for (int l=0; l<= LMAXALL; l++) // mohan fix bug 2010-08-06 + // for (int l=0; l<= mz.lmax_type[it]; l++) + for (int l = 0; l <= LMAXALL; l++) // mohan fix bug 2010-08-06 { - if ( l <= info[it].lmax ) + if (l <= info[it].lmax) { - for (int n=0; nwayd[iw2].iw00 = copy;// position to get overlap. + this->wayd[iw2].iw00 = copy; // position to get overlap. if (test == 2) { cout << "\n iw2=" << iw2 << " iw00=" << copy; @@ -230,42 +233,42 @@ void SpillageStep::set_iw00(void) ++iw2; } } - else if ( info[it].fa[l] == "a" ) + else if (info[it].fa[l] == "a") { this->wayd[iw2].iw00 = iw00; // just the start position to get overlap. ++iw2; } } } - iw00 += 2*l + 1; + iw00 += 2 * l + 1; } } } - if (iw00!=NWFCALL) + if (iw00 != NWFCALL) { cout << "\n iw00=" << iw00 << " NWFCALL=" << NWFCALL; - WARNING_QUIT("SpillageStep::set_iw00","iw00!=input.QS_data[0].nwfc"); + WARNING_QUIT("SpillageStep::set_iw00", "iw00!=input.QS_data[0].nwfc"); } - assert(iw2==nwfc2); + assert(iw2 == nwfc2); return; } - void SpillageStep::set_nchi() { - if (test==1)TITLE("SpillageStep","set_nchi"); - assert( ntype >0 ); + if (test == 1) + TITLE("SpillageStep", "set_nchi"); + assert(ntype > 0); //========================= // (1) generate nchi //========================= this->nchi = 0; - for (int i=0; iwayc = new Way2c4[nchi]; - int ichi=0; - for (int i=0; i 0 ); + if (test == 1) + TITLE("SpillageStep", "allocate_data"); + assert(nwfc > 0); delete[] this->data; - this->data = new Step_Data[ STRNUM ]; + this->data = new Step_Data[STRNUM]; this->ne = NE; - for ( int is=0; isnwfc, - this->nwfc2, - NWFCALL // all read in wave functions number + data[is].init(NKS, + input.QS_data[is].weight, + NBANDS, + NE, + this->nwfc, + this->nwfc2, + NWFCALL // all read in wave functions number ); if (istep == 0) { - for (int ik=0; ik 0 ) - { - for(int ik=0; ik< data[is].nks; ik++) - { - for(int ib=0; ib< data[is].nbands; ib++) - { - for(int ie=0; ie< data[is].ne; ie++) - { - data[is].Qin(ik, ib, iw_step, ie) = input.QS_data[is].Qin(ik, ib, iw, ie); - } - } - } - if(test==2)cout << "\n iw_step=" << iw_step << " iw=" << iw; - iw_step++; - } - } - iw++; - } - } - } + assert( info[it].na == input.QS_data[is].na[it] ); + for(int i=0; i< input.QS_data[is].na[it]; i++) + { + for(int l=0; l< input.QS_data[is].lmax+1; l++) + { + for(int m=0; m<2*l+1; m++) + { + // a break.. + if( l <= info[it].lmax ) + { + if( info[it].n[l] > 0 ) + { + for(int ik=0; ik< data[is].nks; ik++) + { + for(int ib=0; ib< data[is].nbands; ib++) + { + for(int ie=0; ie< data[is].ne; ie++) + { + data[is].Qin(ik, ib, iw_step, ie) = input.QS_data[is].Qin(ik, ib, iw, ie); + } + } + } + if(test==2)cout << "\n iw_step=" << iw_step << " iw=" << iw; + iw_step++; + } + } + iw++; + } + } + } } */ -// cout << "\n iw_step = " << iw_step; -// cout << "\n iw = " << iw; + // cout << "\n iw_step = " << iw_step; + // cout << "\n iw = " << iw; } return; } -void SpillageStep::init_QS_matrix( const int &istr) +void SpillageStep::init_QS_matrix(const int& istr) { - if (test==1)TITLE("SpillageStep","init_QS_matrix"); + if (test == 1) + TITLE("SpillageStep", "init_QS_matrix"); if (USEPW) { @@ -442,29 +451,26 @@ void SpillageStep::init_QS_matrix( const int &istr) const int T = this->wayd[iw].type; const int L = this->wayd[iw].L; const int N = this->wayd[iw].N; - for (int ie=0; iedata[istr].ne; ie++) + for (int ie = 0; ie < this->data[istr].ne; ie++) { //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> // if use 'full orbitals' - if ( this->wayd[iw].average == false) + if (this->wayd[iw].average == false) { - data[istr].Qoverlap[ik](iw,ib) += - input.Coef.C4( T, L, N, ie ) * - data[istr].Qin(ik, ib, jw, ie); + data[istr].Qoverlap[ik](iw, ib) += input.Coef.C4(T, L, N, ie) * data[istr].Qin(ik, ib, jw, ie); } - else// if use 'average orbitals' + else // if use 'average orbitals' { - const int nofm = 2*L+1; // number of m + const int nofm = 2 * L + 1; // number of m std::complex avalue = std::complex(0.0, 0.0); - for (int m=0; m< nofm; m++) + for (int m = 0; m < nofm; m++) { - avalue += data[istr].Qin(ik, ib, jw+m, ie); + avalue += data[istr].Qin(ik, ib, jw + m, ie); } avalue /= nofm; - data[istr].Qoverlap[ik](iw,ib) += - input.Coef.C4( T, L, N, ie ) * avalue; + data[istr].Qoverlap[ik](iw, ib) += input.Coef.C4(T, L, N, ie) * avalue; } //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> } @@ -472,7 +478,7 @@ void SpillageStep::init_QS_matrix( const int &istr) } // init S matrix - for (int iw=0; iwwayd[iw].type; const int L1 = this->wayd[iw].L; @@ -496,99 +502,95 @@ void SpillageStep::init_QS_matrix( const int &istr) } else { -// double value = 0.0; + // double value = 0.0; for (int ie = 0; ie < ne; ie++) { for (int ie2 = 0; ie2 < ne; ie2++) { //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> // case 1 - if ( !wayd[iw].average && !wayd[iw2].average) + if (!wayd[iw].average && !wayd[iw2].average) { - this->data[istr].Soverlap[ik](iw, iw2) += - input.Coef.C4( T1, L1, N1, ie ) * - input.Coef.C4( T2, L2, N2, ie2) * - input.QS_data[istr].Sq1q2[ik](iw001, iw002, ie, ie2); + this->data[istr].Soverlap[ik](iw, iw2) + += input.Coef.C4(T1, L1, N1, ie) * input.Coef.C4(T2, L2, N2, ie2) + * input.QS_data[istr].Sq1q2[ik](iw001, iw002, ie, ie2); } // case 2 else { std::complex avalue = std::complex(0.0, 0.0); // case 21 - if ( wayd[iw].average && wayd[iw2].average) + if (wayd[iw].average && wayd[iw2].average) { - for (int m1=0; m1<2*L1+1; m1++) + for (int m1 = 0; m1 < 2 * L1 + 1; m1++) { - for (int m2=0; m2<2*L2+1; m2++) + for (int m2 = 0; m2 < 2 * L2 + 1; m2++) { - avalue += input.QS_data[istr].Sq1q2[ik](iw001+m1, iw002+m2, ie, ie2); + avalue += input.QS_data[istr].Sq1q2[ik](iw001 + m1, iw002 + m2, ie, ie2); } } - avalue /= ((2*L1+1) * (2*L2+1)); + avalue /= ((2 * L1 + 1) * (2 * L2 + 1)); } // case 22 - if ( !wayd[iw].average && wayd[iw2].average) + if (!wayd[iw].average && wayd[iw2].average) { - for (int m2=0; m2<2*L2+1; m2++) + for (int m2 = 0; m2 < 2 * L2 + 1; m2++) { - avalue += input.QS_data[istr].Sq1q2[ik](iw001, iw002+m2, ie, ie2); + avalue += input.QS_data[istr].Sq1q2[ik](iw001, iw002 + m2, ie, ie2); } - avalue /= (2*L2+1); + avalue /= (2 * L2 + 1); } // case 23 - if ( wayd[iw].average && !wayd[iw2].average) + if (wayd[iw].average && !wayd[iw2].average) { - for (int m1=0; m1<2*L1+1; m1++) + for (int m1 = 0; m1 < 2 * L1 + 1; m1++) { - avalue += input.QS_data[istr].Sq1q2[ik](iw001+m1, iw002, ie, ie2); + avalue += input.QS_data[istr].Sq1q2[ik](iw001 + m1, iw002, ie, ie2); } - avalue /= (2*L1+1); + avalue /= (2 * L1 + 1); } - this->data[istr].Soverlap[ik](iw, iw2) += - input.Coef.C4( T1, L1, N1, ie ) * - input.Coef.C4( T2, L2, N2, ie2) * - avalue; + this->data[istr].Soverlap[ik](iw, iw2) + += input.Coef.C4(T1, L1, N1, ie) * input.Coef.C4(T2, L2, N2, ie2) * avalue; //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> } - }// ie2 - }// ie - }// iw - -// if(iw==3) -// { -// ofs_running << "\n iw=" << iw << " iw2=" << iw2 << " S[" << ik << "]=" << this->data[istr].Soverlap[ik](iw, iw2); -// } - }// ib + } // ie2 + } // ie + } // iw + + // if(iw==3) + // { + // ofs_running << "\n iw=" << iw << " iw2=" << iw2 << " S[" << ik << "]=" << + // this->data[istr].Soverlap[ik](iw, iw2); + // } + } // ib //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> - }//ik - + } // ik -// if(ik==0) -// { -// PRINTCM("S", data[istr].Soverlap[ik]); -// PRINTCM("Q", data[istr].Qoverlap[ik]); -// } + // if(ik==0) + // { + // PRINTCM("S", data[istr].Soverlap[ik]); + // PRINTCM("Q", data[istr].Qoverlap[ik]); + // } // because S has symmetry : S^{+} = S -// for(int iw=0; iwdata[istr].inverse_S[ik].reset_loop(); - this->data[istr].inverse_S[ik].cal_inverse(this->data[istr].Soverlap[ik]); + this->data[istr].inverse_S[ik].reset_loop(); + this->data[istr].inverse_S[ik].cal_inverse(this->data[istr].Soverlap[ik]); // notice : only half A is stored!!! this->data[istr].Sinv[ik] = this->data[istr].inverse_S[ik].get_inverse(); - //if(ik==0) - //PRINTCM("init_QS_matrix :: Sinv",data[istr].inverse.A); + // if(ik==0) + // PRINTCM("init_QS_matrix :: Sinv",data[istr].inverse.A); } return; @@ -596,10 +598,9 @@ void SpillageStep::init_QS_matrix( const int &istr) // (1) Calculate The Spillage Value // (2) For each Coefficient(type,l,n,ie) changed orbital. -double SpillageStep::get_spillage( - const int &istr, // index of structure - const int &ic, // index of orbitals: (type,l,n) - const int &ie // index of eigenvalue of Jlq +double SpillageStep::get_spillage(const int& istr, // index of structure + const int& ic, // index of orbitals: (type,l,n) + const int& ie // index of eigenvalue of Jlq ) { double spillage; // not need to initialized. @@ -629,18 +630,10 @@ double SpillageStep::get_spillage( } // (2) get new c4 coefficient from ic. - const double c4_new = input.Coef.C4( - this->wayc[ic].type, - this->wayc[ic].L, - this->wayc[ic].N, - ie); + const double c4_new = input.Coef.C4(this->wayc[ic].type, this->wayc[ic].L, this->wayc[ic].N, ie); // (3) get old c4 coefficient. - const double c4_old = input.Coef.C4_old( - this->wayc[ic].type, - this->wayc[ic].L, - this->wayc[ic].N, - ie); + const double c4_old = input.Coef.C4_old(this->wayc[ic].type, this->wayc[ic].L, this->wayc[ic].N, ie); int unstable = 0; @@ -650,22 +643,22 @@ double SpillageStep::get_spillage( PW.update_psi1d(ilevel, ic, ie, c4_new, c4_old); } - //ofs_running << "\n ic_now = " << ic << " ie_now=" << ie; - - #ifdef _OPENMP // Peize Lin add 2016-01-18 - omp_lock_t unstable_lock; - omp_init_lock(&unstable_lock); - #endif - - //#ifdef _OPENMP - //omp_set_num_threads(8); // Peize Lin test - //#endif - - #ifdef _OPENMP - #pragma omp parallel for // Peize Lin add 2016-01-18 - #endif - for (int ik = 0; ik < Stru->nks ; ik++) - { + // ofs_running << "\n ic_now = " << ic << " ie_now=" << ie; + +#ifdef _OPENMP // Peize Lin add 2016-01-18 + omp_lock_t unstable_lock; + omp_init_lock(&unstable_lock); +#endif + + // #ifdef _OPENMP + // omp_set_num_threads(8); // Peize Lin test + // #endif + +#ifdef _OPENMP +#pragma omp parallel for // Peize Lin add 2016-01-18 +#endif + for (int ik = 0; ik < Stru->nks; ik++) + { // cout << "ik=" << ik << endl; // cout << "\n ik = " << ik <<" c4_new=" << c4_new << " c4_old=" << c4_old << " nks=" << Stru->nks << endl; @@ -680,49 +673,49 @@ double SpillageStep::get_spillage( { PW.update_psi3d(ilevel, ic, ik); - timer::tick("PW_Basis","calculateS"); + timer::tick("PW_Basis", "calculateS"); // old version need to change pw_basis too - if (NPROC_IN_POOL==1) + if (NPROC_IN_POOL == 1) { - for (int mu=0; muwayd[mu ].ic; + const int ic1 = this->wayd[mu].ic; const int ic2 = this->wayd[nu].ic; - if (ic1==ic || ic2==ic) + if (ic1 == ic || ic2 == ic) { this->data[istr].Strial[ik](mu, nu) = PW.calculateS(mu, nu, ik); } } } } - else + else { - // may error when openmp? unchecked // Peize Lin note 2016-01-18 - #ifdef _OPENMP - if(0==omp_get_thread_num()) - cout<<"src_spillage/SpillageStep.cpp line1103 openmp "<wayd[mu].ic; const int ic2 = this->wayd[nu].ic; // a small trick - if (ic1==ic || ic2==ic) + if (ic1 == ic || ic2 == ic) { PW.save[count] = PW.calculateS(mu, nu, ik); PW.posmu[count] = mu; PW.posnu[count] = nu; - //this->data[istr].Strial[ik](mu, nu) = PW.calculateS(mu, nu, ik); - //ofs_running << "\n save=" << save[count] << " mu=" << mu << " nu=" << nu; + // this->data[istr].Strial[ik](mu, nu) = PW.calculateS(mu, nu, ik); + // ofs_running << "\n save=" << save[count] << " mu=" << mu << " nu=" << nu; ++count; } } @@ -730,75 +723,79 @@ double SpillageStep::get_spillage( #ifdef __MPI Parallel_Reduce::reduce_complex_double_pool(PW.save, count); #endif - for (int index=0; indexdata[istr].Strial[ik](PW.posmu[index], PW.posnu[index]) = PW.save[index]; // ofs_running << "\n mu=" << pos[index*2] << " nu=" << pos[index*2+1] << " save=" << save[index]; } } - timer::tick("PW_Basis","calculateS"); + timer::tick("PW_Basis", "calculateS"); } else { this->newS(istr, ic, ie, ik, c4_new, c4_old); } - Stru->inverse_S[ik].cal_inverse(Stru->Strial[ik]); // Peize Lin update 2015-12-05 - -// if(ik==3) -// { -// PRINTCM("S",data[istr].Strial[ik]); -// PRINTCM("Q",data[istr].Qtrial[ik]); -// PRINTCM("Sinv",Stru->inverse.A); -// } + Stru->inverse_S[ik].cal_inverse(Stru->Strial[ik]); // Peize Lin update 2015-12-05 + + // if(ik==3) + // { + // PRINTCM("S",data[istr].Strial[ik]); + // PRINTCM("Q",data[istr].Qtrial[ik]); + // PRINTCM("Sinv",Stru->inverse.A); + // } -// BLOCK_HERE("spillage"); + // BLOCK_HERE("spillage"); // assert all 'Stru' have same k points and bands. // mohan 2010-05-02 Stru->Mk[ik] = 0.00; -//>>>>> -// Keep the flexibility to calculate the spillage value -// of different number of bands for different structures - //for (int ib = 0; ib < Stru->nbands; ib++) -//<<<<< - //but now I want to calculate the spillage of given number of bands + //>>>>> + // Keep the flexibility to calculate the spillage value + // of different number of bands for different structures + // for (int ib = 0; ib < Stru->nbands; ib++) + //<<<<< + // but now I want to calculate the spillage of given number of bands //(1)from input, for all steps. (2) for each step. - for (int ib=bands_start; ibMkb_used(ik,ib) = 0.0; + Stru->Mkb_used(ik, ib) = 0.0; for (int iw = 0; iw < Stru->nwfc2; iw++) { for (int iw2 = 0; iw2 < Stru->nwfc2; iw2++) { - //Mk += (conj(Qtrial[ik](iw, ib)) * sinv(iw, iw2) * Qtrial[ik](iw2, ib)).real(); + // Mk += (conj(Qtrial[ik](iw, ib)) * sinv(iw, iw2) * Qtrial[ik](iw2, ib)).real(); if (iw <= iw2) { - Stru->Mkb_used(ik,ib) += (conj( Stru->Qtrial[ik](iw, ib)) - * Stru->inverse_S[ik].get_inverse()(iw, iw2) // Peize Lin update 2015-12-05 - * Stru->Qtrial[ik](iw2, ib)).real(); + Stru->Mkb_used(ik, ib) + += (conj(Stru->Qtrial[ik](iw, ib)) + * Stru->inverse_S[ik].get_inverse()(iw, iw2) // Peize Lin update 2015-12-05 + * Stru->Qtrial[ik](iw2, ib)) + .real(); } else { - Stru->Mkb_used(ik,ib) += (conj( Stru->Qtrial[ik](iw, ib)) - * conj( Stru->inverse_S[ik].get_inverse()(iw2, iw)) // Peize Lin update 2015-12-05 - * Stru->Qtrial[ik](iw2, ib)).real(); + Stru->Mkb_used(ik, ib) + += (conj(Stru->Qtrial[ik](iw, ib)) + * conj(Stru->inverse_S[ik].get_inverse()(iw2, iw)) // Peize Lin update 2015-12-05 + * Stru->Qtrial[ik](iw2, ib)) + .real(); } } } - Stru->Mk[ik] += Stru->Mkb_used(ik,ib); - //ofs_running << "ik=" << ik << " ib=" << ib << " Mk=" << Stru->Mk[ik] << endl; + Stru->Mk[ik] += Stru->Mkb_used(ik, ib); + // ofs_running << "ik=" << ik << " ib=" << ib << " Mk=" << Stru->Mk[ik] << endl; } Stru->Mk[ik] /= (double)bands_number; - //cout << endl; - - // in fact this is caused by numerical unstability. - #ifdef _OPENMP - omp_set_lock(&unstable_lock); // Peize Lin add 2016-01-18 - #endif - //if (Stru->Mk[ik]*Stru->weight[ik] > upbound)//mohan fix bug 2010-06-14 - if (fill > Stru->spillage0 )//mohan refix bug 2010-06-18 +// cout << endl; + +// in fact this is caused by numerical unstability. +#ifdef _OPENMP + omp_set_lock(&unstable_lock); // Peize Lin add 2016-01-18 +#endif + // if (Stru->Mk[ik]*Stru->weight[ik] > upbound)//mohan fix bug 2010-06-14 + if (fill > Stru->spillage0) // mohan refix bug 2010-06-18 { ++unstable; } @@ -810,60 +807,58 @@ double SpillageStep::get_spillage( { fill += Stru->Mk[ik] * Stru->weight[ik]; } - #ifdef _OPENMP - omp_unset_lock(&unstable_lock); // Peize Lin add 2016-01-18 - #endif - }//end ik +#ifdef _OPENMP + omp_unset_lock(&unstable_lock); // Peize Lin add 2016-01-18 +#endif + } // end ik #ifdef __MPI -// Parallel_Reduce::reduce_int_all(unstable); - if (unstable>0) + // Parallel_Reduce::reduce_int_all(unstable); + if (unstable > 0) { ofs_running << "\n unstable=" << unstable; - //return Stru->spillage0; + // return Stru->spillage0; } Parallel_Reduce::reduce_double_allpool(fill); #endif spillage = Stru->spillage0 - fill; - + if (spillage < 0.0 || spillage > 1.0 || unstable) { -// Peize Lin test -for (int ik = 0; ik < Stru->nks ; ik++) -{ - //cout<<"ik:\t"<Strial[ik]); - cout_matrix( Stru->inverse_S[ik].get_inverse()); - for( size_t i=0; iStrial[ik].nr; ++i) - for( size_t j=0; jStrial[ik](i,j) = conj(Stru->Strial[ik](j,i)); - ComplexMatrix inverse_tmp(Stru->inverse_S[ik].get_inverse()); - for( size_t i=0; iStrial[ik] * inverse_tmp); -} + // Peize Lin test + for (int ik = 0; ik < Stru->nks; ik++) + { + // cout<<"ik:\t"<Strial[ik]); + cout_matrix(Stru->inverse_S[ik].get_inverse()); + for (size_t i = 0; i < Stru->Strial[ik].nr; ++i) + for (size_t j = 0; j < i; ++j) + Stru->Strial[ik](i, j) = conj(Stru->Strial[ik](j, i)); + ComplexMatrix inverse_tmp(Stru->inverse_S[ik].get_inverse()); + for (size_t i = 0; i < inverse_tmp.nr; ++i) + for (size_t j = 0; j < i; ++j) + inverse_tmp(i, j) = conj(inverse_tmp(j, i)); + cout_matrix(Stru->Strial[ik] * inverse_tmp); + } - cout << "\n spillage0 = " << Stru->spillage0; - cout << "\n Type=" << this->wayc[ic].type - << " L=" << this->wayc[ic].L - << " N=" << this->wayc[ic].N; + cout << "\n spillage0 = " << Stru->spillage0; + cout << "\n Type=" << this->wayc[ic].type << " L=" << this->wayc[ic].L << " N=" << this->wayc[ic].N; cout << "\n bands_number = " << bands_number; - for (int ib=0; ibMkb_used(ik,ib) * Stru->weight[ik]/bands_number; + mm += Stru->Mkb_used(ik, ib) * Stru->weight[ik] / bands_number; } cout << "\n ib=" << ib << " new fill=" << mm; } double sum = 0.0; - for (int ik=0; ikMk[ik]; sum += Stru->Mk[ik] * Stru->weight[ik]; @@ -877,26 +872,28 @@ for (int ik = 0; ik < Stru->nks ; ik++) return spillage; } - -void SpillageStep::newQ( const int &istr, const int &ic, const int &ie, const int &ik, - const double &c4_now ,const double &c4_old) +void SpillageStep::newQ(const int& istr, + const int& ic, + const int& ie, + const int& ik, + const double& c4_now, + const double& c4_old) { -// TITLE("SpillageStep","newQ"); + // TITLE("SpillageStep","newQ"); const double delta = c4_now - c4_old; for (int iw = 0; iw < nwfc2; iw++) { const int jw = this->wayd[iw].iw00; - if ( this->wayd[iw].ic == ic) + if (this->wayd[iw].ic == ic) { //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> - if ( !wayd[iw].average ) + if (!wayd[iw].average) { // the change of dQ= is due to // the change of Jlq(iw)=delta(Jlq(iw)) for (int ib = 0; ib < this->data[istr].nbands; ib++) { - this->data[istr].Qtrial[ik](iw, ib) += delta - * this->data[istr].Qin(ik, ib, jw, ie); + this->data[istr].Qtrial[ik](iw, ib) += delta * this->data[istr].Qin(ik, ib, jw, ie); } } else // if use average orbital @@ -905,11 +902,11 @@ void SpillageStep::newQ( const int &istr, const int &ic, const int &ie, const in for (int ib = 0; ib < this->data[istr].nbands; ib++) { std::complex avalue = std::complex(0.0, 0.0); - for (int m=0; m<2*L+1; m++) + for (int m = 0; m < 2 * L + 1; m++) { - avalue += this->data[istr].Qin(ik, ib, jw+m, ie); + avalue += this->data[istr].Qin(ik, ib, jw + m, ie); } - avalue /= (2*L+1); + avalue /= (2 * L + 1); this->data[istr].Qtrial[ik](iw, ib) += delta * avalue; } } @@ -919,12 +916,14 @@ void SpillageStep::newQ( const int &istr, const int &ic, const int &ie, const in return; } - -void SpillageStep::newS( - const int &istr, const int &ic, const int &ie, const int &ik, - const double &c4_now ,const double &c4_old) +void SpillageStep::newS(const int& istr, + const int& ic, + const int& ie, + const int& ik, + const double& c4_now, + const double& c4_old) { -// TITLE("SpillageStep","newS"); + // TITLE("SpillageStep","newS"); // example c = c2, ie2 = ie //========================================================================================================== // iw0-s iw1p1 iw2p2 iw3p3 iw4s @@ -982,7 +981,7 @@ void SpillageStep::newS( for (int ie1 = 0; ie1 < this->data[istr].ne; ie1++) { // update for the 'ic' row - if ( ie1 != ie || ic1 !=ic ) + if (ie1 != ie || ic1 != ic) { for (int iw2 = iw; iw2 < this->data[istr].nwfc2; iw2++) { @@ -991,51 +990,49 @@ void SpillageStep::newS( { // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> // case 1 - if ( !wayd[iw].average && !wayd[iw2].average) + if (!wayd[iw].average && !wayd[iw2].average) { - this->data[istr].Strial[ik](iw, iw2) += - (c4_now * input.Coef.C4( T1, L1, N1, ie1 ) - - c4_old * input.Coef.C4_old( T1, L1, N1, ie1 ) ) - * input.QS_data[istr].Sq1q2[ik](jw, kw, ie1, ie); + this->data[istr].Strial[ik](iw, iw2) += (c4_now * input.Coef.C4(T1, L1, N1, ie1) + - c4_old * input.Coef.C4_old(T1, L1, N1, ie1)) + * input.QS_data[istr].Sq1q2[ik](jw, kw, ie1, ie); } else // case 2 { const int L2 = this->wayd[iw2].L; std::complex avalue = std::complex(0.0, 0.0); // case 2.1 - if ( wayd[iw].average && !wayd[iw2].average) + if (wayd[iw].average && !wayd[iw2].average) { - for (int m1=0; m1<2*L1+1; m1++) + for (int m1 = 0; m1 < 2 * L1 + 1; m1++) { - avalue += input.QS_data[istr].Sq1q2[ik](jw+m1, kw, ie1, ie); + avalue += input.QS_data[istr].Sq1q2[ik](jw + m1, kw, ie1, ie); } - avalue /= (2*L1+1); + avalue /= (2 * L1 + 1); } // case 2.2 - if ( !wayd[iw].average && wayd[iw2].average) + if (!wayd[iw].average && wayd[iw2].average) { - for (int m2=0; m2<2*L2+1; m2++) + for (int m2 = 0; m2 < 2 * L2 + 1; m2++) { - avalue += input.QS_data[istr].Sq1q2[ik](jw, kw+m2, ie1, ie); + avalue += input.QS_data[istr].Sq1q2[ik](jw, kw + m2, ie1, ie); } - avalue /= (2*L2+1); + avalue /= (2 * L2 + 1); } // case 2.3 - if ( wayd[iw].average && wayd[iw2].average) + if (wayd[iw].average && wayd[iw2].average) { - for (int m1=0; m1<2*L1+1; m1++) + for (int m1 = 0; m1 < 2 * L1 + 1; m1++) { - for (int m2=0; m2<2*L2+1; m2++) + for (int m2 = 0; m2 < 2 * L2 + 1; m2++) { - avalue += input.QS_data[istr].Sq1q2[ik](jw+m1, kw+m2, ie1, ie); + avalue += input.QS_data[istr].Sq1q2[ik](jw + m1, kw + m2, ie1, ie); } } - avalue /= ( (2*L1+1)*(2*L2+1) ); + avalue /= ((2 * L1 + 1) * (2 * L2 + 1)); } - this->data[istr].Strial[ik](iw, iw2) += - (c4_now * input.Coef.C4( T1, L1, N1, ie1 ) - - c4_old * input.Coef.C4_old( T1, L1, N1, ie1 ) ) - * avalue; + this->data[istr].Strial[ik](iw, iw2) += (c4_now * input.Coef.C4(T1, L1, N1, ie1) + - c4_old * input.Coef.C4_old(T1, L1, N1, ie1)) + * avalue; } // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> } @@ -1053,51 +1050,49 @@ void SpillageStep::newS( for (int ie2 = 0; ie2 < this->data[istr].ne; ie2++) { // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> - if ( !wayd[iw].average && !wayd[iw2].average) + if (!wayd[iw].average && !wayd[iw2].average) { - this->data[istr].Strial[ik](iw, iw2) += - (c4_now * input.Coef.C4( T2, L2, N2, ie2 ) - - c4_old * input.Coef.C4_old( T2, L2, N2, ie2 ) ) - * input.QS_data[istr].Sq1q2[ik](jw, kw, ie, ie2); + this->data[istr].Strial[ik](iw, iw2) += (c4_now * input.Coef.C4(T2, L2, N2, ie2) + - c4_old * input.Coef.C4_old(T2, L2, N2, ie2)) + * input.QS_data[istr].Sq1q2[ik](jw, kw, ie, ie2); } else // case 2 { const int L2 = this->wayd[iw2].L; std::complex avalue = std::complex(0.0, 0.0); // case 2.1 - if ( wayd[iw].average && !wayd[iw2].average) + if (wayd[iw].average && !wayd[iw2].average) { - for (int m1=0; m1<2*L1+1; m1++) + for (int m1 = 0; m1 < 2 * L1 + 1; m1++) { - avalue += input.QS_data[istr].Sq1q2[ik](jw+m1, kw, ie, ie2); + avalue += input.QS_data[istr].Sq1q2[ik](jw + m1, kw, ie, ie2); } - avalue /= (2*L1+1); + avalue /= (2 * L1 + 1); } // case 2.2 - if ( !wayd[iw].average && wayd[iw2].average) + if (!wayd[iw].average && wayd[iw2].average) { - for (int m2=0; m2<2*L2+1; m2++) + for (int m2 = 0; m2 < 2 * L2 + 1; m2++) { - avalue += input.QS_data[istr].Sq1q2[ik](jw, kw+m2, ie, ie2); + avalue += input.QS_data[istr].Sq1q2[ik](jw, kw + m2, ie, ie2); } - avalue /= (2*L2+1); + avalue /= (2 * L2 + 1); } // case 2.3 - if ( wayd[iw].average && wayd[iw2].average) + if (wayd[iw].average && wayd[iw2].average) { - for (int m1=0; m1<2*L1+1; m1++) + for (int m1 = 0; m1 < 2 * L1 + 1; m1++) { - for (int m2=0; m2<2*L2+1; m2++) + for (int m2 = 0; m2 < 2 * L2 + 1; m2++) { - avalue += input.QS_data[istr].Sq1q2[ik](jw+m1, kw+m2, ie, ie2); + avalue += input.QS_data[istr].Sq1q2[ik](jw + m1, kw + m2, ie, ie2); } } - avalue /= ( (2*L1+1)*(2*L2+1) ); + avalue /= ((2 * L1 + 1) * (2 * L2 + 1)); } - this->data[istr].Strial[ik](iw, iw2) += - (c4_now * input.Coef.C4( T2, L2, N2, ie2 ) - - c4_old * input.Coef.C4_old( T2, L2, N2, ie2 ) ) - * avalue; + this->data[istr].Strial[ik](iw, iw2) += (c4_now * input.Coef.C4(T2, L2, N2, ie2) + - c4_old * input.Coef.C4_old(T2, L2, N2, ie2)) + * avalue; } // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> } @@ -1110,17 +1105,14 @@ void SpillageStep::newS( return; } - - -void SpillageStep::updateQS( const int &istr ) +void SpillageStep::updateQS(const int& istr) { -// if(test==1)TITLE("SpillageStep","updateQS"); - for (int ik=0; ik< data[istr].nks; ik++) + // if(test==1)TITLE("SpillageStep","updateQS"); + for (int ik = 0; ik < data[istr].nks; ik++) { data[istr].Qoverlap[ik] = data[istr].Qtrial[ik]; data[istr].Soverlap[ik] = data[istr].Strial[ik]; - data[istr].inverse_S[ik].update(); // Peize Lin update 2015-12-05 - } + data[istr].inverse_S[ik].update(); // Peize Lin update 2015-12-05 + } return; } - diff --git a/tools/01_NAO_generation/SIAB/src_spillage/SpillageStep.h b/tools/01_NAO_generation/SIAB/src_spillage/SpillageStep.h index 6b3afcc8467..360ea0fb6d8 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/SpillageStep.h +++ b/tools/01_NAO_generation/SIAB/src_spillage/SpillageStep.h @@ -1,33 +1,33 @@ #ifndef SPILLAGE_STEP #define SPILLAGE_STEP -#include "common.h" #include "Step_Data.h" #include "Type_Information.h" +#include "common.h" // use orbital index to find the information // of each orbital. struct Way2Data // dimension : nwfc2 { - int type; - int i; // atom index - int L; - int N; - int m; // mohan add 2010-06-16 + int type; + int i; // atom index + int L; + int N; + int m; // mohan add 2010-06-16 - int iw0; // to label , which not include 'N', but not used now, - int iw00; // to - int ic; + int iw0; // to label , which not include 'N', but not used now, + int iw00; // to + int ic; - bool average; // mohan add 2009-07-12 + bool average; // mohan add 2009-07-12 }; struct Way2c4 // dimension : nchi { - int type; - int L; - int N; - char spd; //mohan add 2010-04-17, the orbital type. + int type; + int L; + int N; + char spd; // mohan add 2010-04-17, the orbital type. }; // (1) It's a memeber of MultiZeta. @@ -38,41 +38,38 @@ struct Way2c4 // dimension : nchi // (6) It also contains four calss: Type_Information; Step_Data; Way2c4; Way2Data; class SpillageStep { - public: - SpillageStep(); - ~SpillageStep(); - - Type_Information *info; // dimension : type - Step_Data *data; - Way2c4 *wayc; - Way2Data *wayd; + public: + SpillageStep(); + ~SpillageStep(); - void set_info(const int &ntype_in); - void set_nwfc(void); - void set_iw00(void); - void set_nchi(void); + Type_Information* info; // dimension : type + Step_Data* data; + Way2c4* wayc; + Way2Data* wayd; - // allocate data for current level. - void allocate_data( const int &istep); + void set_info(const int& ntype_in); + void set_nwfc(void); + void set_iw00(void); + void set_nchi(void); - void init_QS_matrix( const int &is); - double get_spillage( const int &is, const int &ic, const int &ie); - void updateQS( const int &is ); + // allocate data for current level. + void allocate_data(const int& istep); - int ntype; - int nchi; // total radial wave functions(T,L,N); - int nwfc2; - int ne; - int ilevel; + void init_QS_matrix(const int& is); + double get_spillage(const int& is, const int& ic, const int& ie); + void updateQS(const int& is); - private: + int ntype; + int nchi; // total radial wave functions(T,L,N); + int nwfc2; + int ne; + int ilevel; - void newQ( const int &is, const int &ic, const int &ie, const int &ik, - const double &c4_now, const double &c4_old); - void newS( const int &is, const int &ic, const int &ie, const int &ik, - const double &c4_now ,const double &c4_old); + private: + void newQ(const int& is, const int& ic, const int& ie, const int& ik, const double& c4_now, const double& c4_old); + void newS(const int& is, const int& ic, const int& ie, const int& ik, const double& c4_now, const double& c4_old); - int test; - int nwfc; + int test; + int nwfc; }; #endif diff --git a/tools/01_NAO_generation/SIAB/src_spillage/SpillageValue.cpp b/tools/01_NAO_generation/SIAB/src_spillage/SpillageValue.cpp index 78e7146b6ab..e5d1d705bb5 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/SpillageValue.cpp +++ b/tools/01_NAO_generation/SIAB/src_spillage/SpillageValue.cpp @@ -1,139 +1,140 @@ #include "SpillageValue.h" + #include "tools.h" SpillageValue::SpillageValue() { - dim = 0; - value = new double[1]; - value_old = new double[1]; + dim = 0; + value = new double[1]; + value_old = new double[1]; } SpillageValue::~SpillageValue() { - delete[] value; - delete[] value_old; + delete[] value; + delete[] value_old; } -void SpillageValue::allocate( const int &dim_in ) +void SpillageValue::allocate(const int& dim_in) { - this->dim = dim_in; - assert( dim < 100); - assert( dim > 0); - - delete[] value; - delete[] value_old; - value = new double[dim]; - value_old = new double[dim]; - - // assume max step < 100 - const int max_levels = 100; - value_each_level.create(dim, max_levels); - - this->reset(); - - return; + this->dim = dim_in; + assert(dim < 100); + assert(dim > 0); + + delete[] value; + delete[] value_old; + value = new double[dim]; + value_old = new double[dim]; + + // assume max step < 100 + const int max_levels = 100; + value_each_level.create(dim, max_levels); + + this->reset(); + + return; } void SpillageValue::reset(void) { - for(int i=0; ivalue; - } - else - { - v = this->value_old; - } - - double mv=0.0; - // case 1: // get the largest spillage value. - if( abs(SCHEME_VALUE)==1) - { - for(int i=0; ivalue; + } + else + { + v = this->value_old; + } + + double mv = 0.0; + // case 1: // get the largest spillage value. + if (abs(SCHEME_VALUE) == 1) + { + for (int i = 0; i < dim; i++) + { + mv = std::max(mv, v[i]); + } + return mv; + } + // case 2: // get average spillage value. + else if (abs(SCHEME_VALUE) == 2) + { + for (int i = 0; i < dim; i++) + { + mv += v[i] / dim; + } + return mv; + } + // case 3: // mixture of case 1 and 2 + else if (abs(SCHEME_VALUE) == 3) + { + double average = 0.0; + for (int i = 0; i < dim; i++) + { + average += v[i] / dim; + } + + double max_diverse = 0.0; + double avg_diverse = 0.0; + for (int i = 0; i < dim; i++) + { + // max_diverse = std::max(max_diverse, abs(average-v[i])); + avg_diverse += abs(average - v[i]) / dim; + } + + if (average < 0.1) + { + mv = average * (1.0 + avg_diverse / (average * 0.1)); + } + else + { + mv = average; + } + return mv; + } + else + { + cout << "\n Please Check SCHEME_VALUE." << endl; + WARNING_QUIT("SpillageValue", "get_value"); + } } void SpillageValue::update_value(void) { - for(int i=0; ivalue_each_level(i ,ilevel) = value_old[i]; - } - return; + // dim stands for number of structures. + for (int i = 0; i < dim; i++) + { + this->value_each_level(i, ilevel) = value_old[i]; + } + return; } diff --git a/tools/01_NAO_generation/SIAB/src_spillage/SpillageValue.h b/tools/01_NAO_generation/SIAB/src_spillage/SpillageValue.h index b0827a60f01..7cedd205400 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/SpillageValue.h +++ b/tools/01_NAO_generation/SIAB/src_spillage/SpillageValue.h @@ -6,50 +6,48 @@ class SpillageValue { - friend class Read_INPUT; - - friend class MultiZeta; - // in multizeta, 'value_old' are calculated. - - friend class Metropolis; - // in metropolis, 'value' are set. - - friend class Out_Orbital; - // in OurOrbital, 'value each level' are used. - - public: - - // 1: get new spillage value. - // 0: get old spillage value. - double cal_defined_value( bool get_new_flag); - - void update_value(); - - void out(); - - void save_level( const int &ilevel); - - private: - - // spillage value for each structure. - double *value; - - // the previous spillage value for each structure. - double *value_old; - - // spillage for all structures. - matrix value_each_level; - - // dimension is the number of structures. - // eg. 5 different dimers. - int dim; - - void allocate( const int &dim_in ); - - void reset(); - - SpillageValue(); - ~SpillageValue(); + friend class Read_INPUT; + + friend class MultiZeta; + // in multizeta, 'value_old' are calculated. + + friend class Metropolis; + // in metropolis, 'value' are set. + + friend class Out_Orbital; + // in OurOrbital, 'value each level' are used. + + public: + // 1: get new spillage value. + // 0: get old spillage value. + double cal_defined_value(bool get_new_flag); + + void update_value(); + + void out(); + + void save_level(const int& ilevel); + + private: + // spillage value for each structure. + double* value; + + // the previous spillage value for each structure. + double* value_old; + + // spillage for all structures. + matrix value_each_level; + + // dimension is the number of structures. + // eg. 5 different dimers. + int dim; + + void allocate(const int& dim_in); + + void reset(); + + SpillageValue(); + ~SpillageValue(); }; #endif diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Step_Data.cpp b/tools/01_NAO_generation/SIAB/src_spillage/Step_Data.cpp index 224f67ab796..8348e5de0ba 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Step_Data.cpp +++ b/tools/01_NAO_generation/SIAB/src_spillage/Step_Data.cpp @@ -1,84 +1,80 @@ #include "Step_Data.h" + #include "../src_tools/complexmatrix_inline.h" Step_Data::Step_Data() { -// TITLE("Step_Data","Step_Data"); - Qoverlap = new ComplexMatrix[1]; - Soverlap = new ComplexMatrix[1]; - Qtrial = new ComplexMatrix[1]; - Strial = new ComplexMatrix[1]; - weight = new double[1]; - Sinv = new ComplexMatrix[1]; - inverse_S = new Inverse_Matrix_S[1]; - Mk = new double[1]; - test = 0; + // TITLE("Step_Data","Step_Data"); + Qoverlap = new ComplexMatrix[1]; + Soverlap = new ComplexMatrix[1]; + Qtrial = new ComplexMatrix[1]; + Strial = new ComplexMatrix[1]; + weight = new double[1]; + Sinv = new ComplexMatrix[1]; + inverse_S = new Inverse_Matrix_S[1]; + Mk = new double[1]; + test = 0; } Step_Data::~Step_Data() { -// TITLE("Step_Data","~Step_Data"); - delete[] Qoverlap; - delete[] Soverlap; - delete[] Qtrial; - delete[] Strial; - delete[] weight; - delete[] Sinv; - delete[] inverse_S; - delete[] Mk; + // TITLE("Step_Data","~Step_Data"); + delete[] Qoverlap; + delete[] Soverlap; + delete[] Qtrial; + delete[] Strial; + delete[] weight; + delete[] Sinv; + delete[] inverse_S; + delete[] Mk; } -void Step_Data::init( - const int &nks_in, - const double* weight_in, - const int &nbands_in, - const int &ne_in, - const int &nwfc_in, - const int &nwfc2_in, - const int &nwfc_all_in) +void Step_Data::init(const int& nks_in, + const double* weight_in, + const int& nbands_in, + const int& ne_in, + const int& nwfc_in, + const int& nwfc2_in, + const int& nwfc_all_in) { - if(test==1)TITLE("Step_Data","init"); - this->nks = nks_in; - assert(nks > 0); - - delete[] Mk; //mohan add 2010-05-02 - this->Mk = new double[nks]; - ZEROS( Mk, nks ); - - delete[] weight; - this->weight = new double[nks]; - for(int ik=0; iknbands = nbands_in; - - assert(nbands > 0); - this->Mkb.create(nks, nbands);//mohan add 2010-05-02 - this->Mkb_used.create(nks, nbands); - this->Mkb.zero_out(); - this->Mkb_used.zero_out(); - - this->ne = ne_in; - this->nwfc = nwfc_in; - this->nwfc2 = nwfc2_in; - this->nwfc_all = nwfc_all_in; - - - if(test==1) - { - cout << "\n" << setw(10) << "nwfc" - << setw(10) << "nwfc2" - << setw(10) << "nwfc_all"; - - cout << "\n" << setw(10) << nwfc - << setw(10) << nwfc2 - << setw(10) << nwfc_all; - } - - this->allocate(); - - return; + if (test == 1) + TITLE("Step_Data", "init"); + this->nks = nks_in; + assert(nks > 0); + + delete[] Mk; // mohan add 2010-05-02 + this->Mk = new double[nks]; + ZEROS(Mk, nks); + + delete[] weight; + this->weight = new double[nks]; + for (int ik = 0; ik < nks; ik++) + { + weight[ik] = weight_in[ik]; + } + this->nbands = nbands_in; + + assert(nbands > 0); + this->Mkb.create(nks, nbands); // mohan add 2010-05-02 + this->Mkb_used.create(nks, nbands); + this->Mkb.zero_out(); + this->Mkb_used.zero_out(); + + this->ne = ne_in; + this->nwfc = nwfc_in; + this->nwfc2 = nwfc2_in; + this->nwfc_all = nwfc_all_in; + + if (test == 1) + { + cout << "\n" << setw(10) << "nwfc" << setw(10) << "nwfc2" << setw(10) << "nwfc_all"; + + cout << "\n" << setw(10) << nwfc << setw(10) << nwfc2 << setw(10) << nwfc_all; + } + + this->allocate(); + + return; } void Step_Data::allocate() @@ -90,14 +86,15 @@ void Step_Data::allocate() // be called in Multi_Zeta.cpp void Step_Data::initQS(void) { - if(test==1) TITLE("Step_Data","initQS"); - assert(nks > 0); - assert(nbands > 0); - assert(nwfc > 0); // nwfc is for , - assert(ne > 0); - - delete[] Qtrial; - delete[] Strial; + if (test == 1) + TITLE("Step_Data", "initQS"); + assert(nks > 0); + assert(nbands > 0); + assert(nwfc > 0); // nwfc is for , + assert(ne > 0); + + delete[] Qtrial; + delete[] Strial; delete[] Qoverlap; delete[] Soverlap; delete[] Sinv; @@ -110,18 +107,15 @@ void Step_Data::initQS(void) this->Sinv = new ComplexMatrix[nks]; this->inverse_S = new Inverse_Matrix_S[nks]; - for(int ik=0; ik=-1); - delete[] n; - delete[] nbase; - delete[] fa; - n = new int[lmax+1]; - nbase = new int[lmax+1]; - fa = new string[lmax+1]; - for(int i=0; i= -1); + delete[] n; + delete[] nbase; + delete[] fa; + n = new int[lmax + 1]; + nbase = new int[lmax + 1]; + fa = new string[lmax + 1]; + for (int i = 0; i < lmax + 1; i++) + { + n[i] = nbase[i] = 0; + } + return; } void Type_Information::cal_nmax() { - this->nmax = 0; - for(int l=0; lnmax = std::max( nmax, n[l]+nbase[l] ); - } -// cout << "\n nmax = " << nmax; - return; + this->nmax = 0; + for (int l = 0; l < lmax + 1; l++) + { + this->nmax = std::max(nmax, n[l] + nbase[l]); + } + // cout << "\n nmax = " << nmax; + return; } void Type_Information::cal_nw(void) { - TITLE("Type_Information","cal_nw"); - // where is this lmax from ? - // in Multi_Zeta, - // read in from file directly - // for each level and each type. + TITLE("Type_Information", "cal_nw"); + // where is this lmax from ? + // in Multi_Zeta, + // read in from file directly + // for each level and each type. - // be called in SpillageStep. + // be called in SpillageStep. - // nw2: number of total local orbitals. - // (ia, l, n, m) - // Including multi-zeta orbitals, - // which are made up from same Jl(q). + // nw2: number of total local orbitals. + // (ia, l, n, m) + // Including multi-zeta orbitals, + // which are made up from same Jl(q). this->nw2 = 0; - for(int l=0; l 1 - // f orbitals, 7 ==> 1 - } - //cout << "\n nw = " << nw; + if (fa[l] == "f") + { + nw2 += (2 * l + 1) * n[l] * na; + } + else if (fa[l] == "a") + { + nw2 += 1 * na; // 1 means average orbital used, + // for example, d orbitals 5 ==> 1 + // f orbitals, 7 ==> 1 + } + // cout << "\n nw = " << nw; } - // nw: number of overlap we must save. - // doesn't consider n[l], - // the only condition is n[l] > 0 in this atom species. - this->nw = 0; - for(int l=0; l0) - { - nw += (2*l+1)*na; - } - } + // nw: number of overlap we must save. + // doesn't consider n[l], + // the only condition is n[l] > 0 in this atom species. + this->nw = 0; + for (int l = 0; l < lmax + 1; l++) + { + // warning : in different step, n[l] value may be different. + if (n[l] > 0) + { + nw += (2 * l + 1) * na; + } + } - cout << " cal_nw::nw = " << nw; - cout << " cal_nw::nw2 = " << nw2 << endl; + cout << " cal_nw::nw = " << nw; + cout << " cal_nw::nw2 = " << nw2 << endl; return; } - - diff --git a/tools/01_NAO_generation/SIAB/src_spillage/Type_Information.h b/tools/01_NAO_generation/SIAB/src_spillage/Type_Information.h index 1a9e0d38f6e..99390a75f1d 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/Type_Information.h +++ b/tools/01_NAO_generation/SIAB/src_spillage/Type_Information.h @@ -5,37 +5,37 @@ // Information of a particular type class Type_Information { - public: + public: Type_Information(); ~Type_Information(); - int id;// index of element periodic table. + int id; // index of element periodic table. - int na;// number of atoms of this type. + int na; // number of atoms of this type. - string *fa;// mohan add 2009-07-12 - // f: full orbtial used. - // a: average orbital used. + string* fa; // mohan add 2009-07-12 + // f: full orbtial used. + // a: average orbital used. - string state; // mohan add 2009-08-27, new/skip - // state: new, start Metropolis calculation in this step. - // state: skip, skip this Metropolis calculation, - // used only in case RESTART=TRUE; - - int lmax;// lmax used. + string state; // mohan add 2009-08-27, new/skip + // state: new, start Metropolis calculation in this step. + // state: skip, skip this Metropolis calculation, + // used only in case RESTART=TRUE; - int *n;// number of radial function for each L. + int lmax; // lmax used. - int *nbase;// number of radial function of each L, sum up all the steps information before! + int* n; // number of radial function for each L. - int nmax; + int* nbase; // number of radial function of each L, sum up all the steps information before! + + int nmax; int nw; - int nw2; + int nw2; - void init();// allocate n, nbase, fa; + void init(); // allocate n, nbase, fa; void cal_nw(); - void cal_nmax(); + void cal_nmax(); }; #endif diff --git a/tools/01_NAO_generation/SIAB/src_spillage/common.cpp b/tools/01_NAO_generation/SIAB/src_spillage/common.cpp index 3ec4fd9da4d..55e83c55854 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/common.cpp +++ b/tools/01_NAO_generation/SIAB/src_spillage/common.cpp @@ -1,56 +1,59 @@ #include "common.h" -void PRINTCM(const string &s, const ComplexMatrix &m) +void PRINTCM(const string& s, const ComplexMatrix& m) { const int b1 = m.nr; const int b2 = m.nc; - cout << "\n" << s << " " << b1 << " " << b2 ; + cout << "\n" << s << " " << b1 << " " << b2; - if (b1*b2 == 0) return; + if (b1 * b2 == 0) + return; - cout << "\n norm : "; - for (int i = 0;i < b1;i++) + cout << "\n norm : "; + for (int i = 0; i < b1; i++) { - for (int j = 0;j < b2;j++) + for (int j = 0; j < b2; j++) { - if (j % 8 == 0) cout << "\n "; - double norm = ( conj( m(i, j) ) * m(i, j) ).real(); - if( abs(norm) > 1.0e-9) cout<< setw(12) << norm; - else cout< 1.0e-9) + cout << setw(12) << norm; + else + cout << setw(12) << "0"; } } - /* - cout << "\n real part : "; - for (int i = 0;i < b1;i++) + /* + cout << "\n real part : "; + for (int i = 0;i < b1;i++) { for (int j = 0;j < b2;j++) { if (j % 8 == 0) cout << "\n "; - double n = ( m(i, j) ).real(); - if( abs(n) > 1.0e-9) cout<< setw(12) << n; - else cout< 1.0e-9) cout<< setw(12) << n; + else cout< 1.0e-9) cout<< setw(12) << n; - else cout< 1.0e-9) cout<< setw(12) << n; + else cout<> SearchName; - if( SearchName == TargetName) + if (SearchName == TargetName) { find = true; -// cout << " Find the block:" << TargetName << endl; + // cout << " Find the block:" << TargetName << endl; break; } } @@ -78,34 +81,33 @@ bool SCAN_BEGIN(ifstream &ifs, const string &TargetName, const bool restart, con return find; } -void SCAN_END(ifstream &ifs, const string &TargetName) +void SCAN_END(ifstream& ifs, const string& TargetName) { string SearchName; ifs >> SearchName; - if( SearchName != TargetName) + if (SearchName != TargetName) { - cout<<"\n"<<"In SCAN_END, can't find: "<> ok; - } void QUIT(void) @@ -117,27 +119,25 @@ void QUIT(void) exit(0); } -void TITLE(const string &class_name,const string &function_name) +void TITLE(const string& class_name, const string& function_name) { -// ofs_running<<"\n\n ==> "< "< "< " << class_name << "::" << function_name << endl; + cout << " -------------------------------------" << endl; return; } -void TITLE(ofstream &ofs, const string &class_name,const string &function_name) +void TITLE(ofstream& ofs, const string& class_name, const string& function_name) { - return; - cout<<"\n -------------------------------------" << endl; - cout<<" ==> "< " << class_name << "::" << function_name << endl; + cout << " -------------------------------------" << endl; - ofs <<" -------------------------------------" << endl; - ofs <<" ==> "< " << class_name << "::" << function_name << endl; + ofs << " -------------------------------------" << endl; return; } - - diff --git a/tools/01_NAO_generation/SIAB/src_spillage/common.h b/tools/01_NAO_generation/SIAB/src_spillage/common.h index a573be54621..c793487d470 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/common.h +++ b/tools/01_NAO_generation/SIAB/src_spillage/common.h @@ -2,99 +2,94 @@ #define COMMON_H using namespace std; -#include -#include -#include +#include "../src_tools/Random.h" +#include "../src_tools/complexarray.h" +#include "../src_tools/complexmatrix.h" +#include "../src_tools/intarray.h" +#include "../src_tools/inverse_matrix.h" +#include "../src_tools/lapack_connector.h" +#include "../src_tools/mathzone.h" +#include "../src_tools/matrix.h" +#include "../src_tools/matrix3.h" +#include "../src_tools/realarray.h" +#include "../src_tools/timer.h" // added 2009-4-17 -#include +#include #include - -#include -#include -#include +#include +#include #include #include - +#include +#include +#include +#include #include -#include "../src_tools/mathzone.h" -#include "../src_tools/realarray.h" -#include "../src_tools/intarray.h" -#include "../src_tools/matrix.h" -#include "../src_tools/matrix3.h" -#include "../src_tools/complexmatrix.h" -#include "../src_tools/complexarray.h" -#include "../src_tools/lapack_connector.h" -#include "../src_tools/Random.h" -#include "../src_tools/timer.h"// added 2009-4-17 -#include "../src_tools/inverse_matrix.h" - #ifdef __MPI #include #endif template -static void READ_VALUE(ifstream &ifs, T &v) +static void READ_VALUE(ifstream& ifs, T& v) { ifs >> v; ifs.ignore(150, '\n'); return; } - -bool SCAN_BEGIN(ifstream &ifs, const string &TargetName, const bool restart=1, const bool quit=true); -void SCAN_END(ifstream &ifs, const string &TargetName); -void TITLE(const string &class_name,const string &function_name); -void TITLE(ofstream &ofs, const string &class_name,const string &function_name); -void PRINTCM(const string &s, const ComplexMatrix &m); +bool SCAN_BEGIN(ifstream& ifs, const string& TargetName, const bool restart = 1, const bool quit = true); +void SCAN_END(ifstream& ifs, const string& TargetName); +void TITLE(const string& class_name, const string& function_name); +void TITLE(ofstream& ofs, const string& class_name, const string& function_name); +void PRINTCM(const string& s, const ComplexMatrix& m); //========================================================== // GLOBAL FUNCTION : // NAME : ZEROS // set elements of u as zero which u is 1_d complex array //========================================================== -template -void ZEROS(complex *u,const int n) +template +void ZEROS(complex* u, const int n) { - assert(n!=0); - assert(u!=0); - for(int i=0;i(0.0,0.0); + u[i] = complex(0.0, 0.0); } return; } -template -void ZEROS(T *u,const int n) +template +void ZEROS(T* u, const int n) { - assert(n>=0); - for(int i=0;i= 0); + for (int i = 0; i < n; i++) { u[i] = 0; } } -void WARNING_QUIT(const string &file,const string &description); -void BLOCK_HERE( const string &description ); +void WARNING_QUIT(const string& file, const string& description); +void BLOCK_HERE(const string& description); void QUIT(); // mohan add 2010-06-12 template -void OUT(ofstream &ofs,const string &name,const T &a) +void OUT(ofstream& ofs, const string& name, const T& a) { - stringstream name2; - name2 << "[" << name << "]"; - ofs<< " " << setw(40) << name2.str() << " = " << a < -void OUT(ofstream &ofs,const string &name,const T &x, const T &y, const T &z) +void OUT(ofstream& ofs, const string& name, const T& x, const T& y, const T& z) { - stringstream name2; - name2 << "[" << name << "]"; - ofs<< " " << setw(40) <",1,0) ) - { - READ_VALUE( ifs, BANDS_CONTROL); - if (BANDS_CONTROL) - { - READ_VALUE(ifs, BANDS_START); - assert(BANDS_START > 0); - BANDS_START -= 1; - // because the band index start from 0; - READ_VALUE(ifs, BANDS_END); - // BANDS_END need not -=1; - } - - ofs_running << " ATTENTION! YOU SELECT THE BANDS: " << endl; - ofs_running << " BANDS_START=" << BANDS_START+1 << " BANDS_END=" << BANDS_END << endl; - } - - // mohan add 2009-08-26 - if ( SCAN_BEGIN(ifs,"",1,0) ) - { - READ_VALUE( ifs, BLOCK_NE); - assert(BLOCK_NE>=0); - - READ_VALUE( ifs, BLOCK_NE_MIN); - assert(BLOCK_NE >= BLOCK_NE_MIN); - // static_cast( sqrt( 2.0 * ecut )* rcut/3.1415926535897932 ); - - ofs_running << " ATTENTION! YOU SELECT BLOCK_NE=" << BLOCK_NE << endl; - ofs_running << " ALSO BLOCK_NE_MIN=" << BLOCK_NE_MIN << endl; - } - - // mohan add 2009-08-31 - if ( SCAN_BEGIN(ifs,"",1,0) ) - { - double block_ne_energy; - double rcut_in; - READ_VALUE( ifs, block_ne_energy); - READ_VALUE( ifs, rcut_in ); - BLOCK_NE = static_cast( sqrt( 2.0 * block_ne_energy )* rcut_in/3.1415926535897932 ); - assert(BLOCK_NE >= 0); - - ofs_running << "\n Attential! Now, input BLOCK_NE parameter be replaced!" << endl; - ofs_running << "\n New BLOCK_NE = " << BLOCK_NE - << " energy_cut = " << block_ne_energy - << " rcut_in = " << rcut_in << endl; - } - - // mohan add 2009-08-28 - // 1: max spillage value. - // 2: average spillage value. - // 3: add kapa value. - if ( SCAN_BEGIN(ifs,"",1,0) ) - { - READ_VALUE( ifs, SCHEME_VALUE ); - assert(abs(SCHEME_VALUE)>0); - assert(abs(SCHEME_VALUE)<4); - } - - if ( SCAN_BEGIN(ifs,"",1,0) ) - { - READ_VALUE( ifs, NKSTOT ); - assert(NKSTOT>0); + if (ifs) + { + // mohan add 2010-05-01 + if (SCAN_BEGIN(ifs, "", 1, 0)) + { + READ_VALUE(ifs, BANDS_CONTROL); + if (BANDS_CONTROL) + { + READ_VALUE(ifs, BANDS_START); + assert(BANDS_START > 0); + BANDS_START -= 1; + // because the band index start from 0; + READ_VALUE(ifs, BANDS_END); + // BANDS_END need not -=1; + } + + ofs_running << " ATTENTION! YOU SELECT THE BANDS: " << endl; + ofs_running << " BANDS_START=" << BANDS_START + 1 << " BANDS_END=" << BANDS_END << endl; + } + + // mohan add 2009-08-26 + if (SCAN_BEGIN(ifs, "", 1, 0)) + { + READ_VALUE(ifs, BLOCK_NE); + assert(BLOCK_NE >= 0); + + READ_VALUE(ifs, BLOCK_NE_MIN); + assert(BLOCK_NE >= BLOCK_NE_MIN); + // static_cast( sqrt( 2.0 * ecut )* rcut/3.1415926535897932 ); + + ofs_running << " ATTENTION! YOU SELECT BLOCK_NE=" << BLOCK_NE << endl; + ofs_running << " ALSO BLOCK_NE_MIN=" << BLOCK_NE_MIN << endl; + } + + // mohan add 2009-08-31 + if (SCAN_BEGIN(ifs, "", 1, 0)) + { + double block_ne_energy; + double rcut_in; + READ_VALUE(ifs, block_ne_energy); + READ_VALUE(ifs, rcut_in); + BLOCK_NE = static_cast(sqrt(2.0 * block_ne_energy) * rcut_in / 3.1415926535897932); + assert(BLOCK_NE >= 0); + + ofs_running << "\n Attential! Now, input BLOCK_NE parameter be replaced!" << endl; + ofs_running << "\n New BLOCK_NE = " << BLOCK_NE << " energy_cut = " << block_ne_energy + << " rcut_in = " << rcut_in << endl; + } + + // mohan add 2009-08-28 + // 1: max spillage value. + // 2: average spillage value. + // 3: add kapa value. + if (SCAN_BEGIN(ifs, "", 1, 0)) + { + READ_VALUE(ifs, SCHEME_VALUE); + assert(abs(SCHEME_VALUE) > 0); + assert(abs(SCHEME_VALUE) < 4); + } + + if (SCAN_BEGIN(ifs, "", 1, 0)) + { + READ_VALUE(ifs, NKSTOT); + assert(NKSTOT > 0); #ifdef __MPI - READ_VALUE( ifs, KPAR ); - assert(KPAR>0); -#endif - } - - }// end ifs - - else if( ifs2 ) - { - if ( SCAN_BEGIN(ifs2,"",1,0) ) - { - READ_VALUE( ifs2, WFC_FILE); - } - - if ( SCAN_BEGIN(ifs2,"",1,0) ) - { - READ_VALUE( ifs2, SMOOTH); - assert(SMOOTH == 0 || SMOOTH == 1); - READ_VALUE( ifs2, SIGMA); - assert(SIGMA >= 0.0); - READ_VALUE( ifs2, ECUT_JLQ); - assert(ECUT_JLQ > 0 && ECUT_JLQ < 1000); - // mohan set the boundaries) - READ_VALUE( ifs2, RCUT); - assert(RCUT > 0 && RCUT < 15); - READ_VALUE( ifs2, TOLERENCE); - assert(TOLERENCE> 0.0); - } - - if ( SCAN_BEGIN(ifs2,"",1,0) ) + READ_VALUE(ifs, KPAR); + assert(KPAR > 0); +#endif + } + + } // end ifs + + else if (ifs2) + { + if (SCAN_BEGIN(ifs2, "", 1, 0)) + { + READ_VALUE(ifs2, WFC_FILE); + } + + if (SCAN_BEGIN(ifs2, "", 1, 0)) + { + READ_VALUE(ifs2, SMOOTH); + assert(SMOOTH == 0 || SMOOTH == 1); + READ_VALUE(ifs2, SIGMA); + assert(SIGMA >= 0.0); + READ_VALUE(ifs2, ECUT_JLQ); + assert(ECUT_JLQ > 0 && ECUT_JLQ < 1000); + // mohan set the boundaries) + READ_VALUE(ifs2, RCUT); + assert(RCUT > 0 && RCUT < 15); + READ_VALUE(ifs2, TOLERENCE); + assert(TOLERENCE > 0.0); + } + + if (SCAN_BEGIN(ifs2, "", 1, 0)) { - READ_VALUE( ifs2, NKSTOT ); - assert(NKSTOT>0); + READ_VALUE(ifs2, NKSTOT); + assert(NKSTOT > 0); #ifdef __MPI - READ_VALUE( ifs2, KPAR ); - assert(KPAR>0); + READ_VALUE(ifs2, KPAR); + assert(KPAR > 0); #endif } - }// end ifs2 + } // end ifs2 } // end my rank0. //------------------------------------------------ - // second part + // second part //------------------------------------------------ - // 2.1 - if( ifs ) - { - if (MY_RANK==0) - { - this->read_PW_QS_1(); - } - } - else if ( ifs2 ) - { - this->read_WFC(); - } - - // 2.2 - // bcast NKSTOT. - // bcast the global information. - // need to first bcast KPAR. - this->bcast(); - - // set the parallel information - // about pool. - Pkpoints.init(); - - // init the k point information. - Pkpoints.kinfo(NKSTOT); - - // calculate the local number of k points. + // 2.1 + if (ifs) + { + if (MY_RANK == 0) + { + this->read_PW_QS_1(); + } + } + else if (ifs2) + { + this->read_WFC(); + } + + // 2.2 + // bcast NKSTOT. + // bcast the global information. + // need to first bcast KPAR. + this->bcast(); + + // set the parallel information + // about pool. + Pkpoints.init(); + + // init the k point information. + Pkpoints.kinfo(NKSTOT); + + // calculate the local number of k points. #ifdef __MPI - NKS = Pkpoints.nks_pool[MY_POOL]; + NKS = Pkpoints.nks_pool[MY_POOL]; #else - NKS = NKSTOT; + NKS = NKSTOT; #endif - OUT(ofs_running,"NKS",NKS); + OUT(ofs_running, "NKS", NKS); - // 2.3 - if( ifs ) - { - this->read_PW_QS_2(); - } - else if( ifs2 ) - { - this->read_WFC_2(); - USEPW = 1; - } + // 2.3 + if (ifs) + { + this->read_PW_QS_2(); + } + else if (ifs2) + { + this->read_WFC_2(); + USEPW = 1; + } return; } @@ -221,89 +223,89 @@ void Read_INPUT::init(void) void Read_INPUT::bcast(void) { #ifdef __MPI - Parallel_Common::bcast_bool( BANDS_CONTROL ); - Parallel_Common::bcast_int( BANDS_START ); - Parallel_Common::bcast_int( BANDS_END ); + Parallel_Common::bcast_bool(BANDS_CONTROL); + Parallel_Common::bcast_int(BANDS_START); + Parallel_Common::bcast_int(BANDS_END); - Parallel_Common::bcast_int( BLOCK_NE ); - Parallel_Common::bcast_int( BLOCK_NE_MIN ); + Parallel_Common::bcast_int(BLOCK_NE); + Parallel_Common::bcast_int(BLOCK_NE_MIN); - Parallel_Common::bcast_int( SCHEME_VALUE ); - Parallel_Common::bcast_int( KPAR ); + Parallel_Common::bcast_int(SCHEME_VALUE); + Parallel_Common::bcast_int(KPAR); - Parallel_Common::bcast_int( CALSPI ); - Parallel_Common::bcast_bool( RESTART ); - Parallel_Common::bcast_bool( this->output ); - Parallel_Common::bcast_int( STRNUM ); + Parallel_Common::bcast_int(CALSPI); + Parallel_Common::bcast_bool(RESTART); + Parallel_Common::bcast_bool(this->output); + Parallel_Common::bcast_int(STRNUM); - OUT(ofs_running,"STRNUM",STRNUM); + OUT(ofs_running, "STRNUM", STRNUM); - if (MY_RANK!=0) + if (MY_RANK != 0) { delete[] qsfile; qsfile = new string[STRNUM]; } - for (int i=0; iSV.allocate( STRNUM ); + this->SV.allocate(STRNUM); delete[] this->QS_data; - this->QS_data = new ReadData[ STRNUM ]; + this->QS_data = new ReadData[STRNUM]; - for (int i=0; iQS_data[i].OverlapQandS( qsfile[i].c_str() ); - //this->QS_data[in].OverlapSinv( sinv_file.c_str() ); // not used + this->QS_data[i].OverlapQandS(qsfile[i].c_str()); + // this->QS_data[in].OverlapSinv( sinv_file.c_str() ); // not used } return; @@ -330,105 +332,106 @@ void Read_INPUT::read_PW_QS_2(void) void Read_INPUT::read_PW_QS_1(void) { - TITLE(ofs_running, "Read_INPUT","read_PW_QS_1"); + TITLE(ofs_running, "Read_INPUT", "read_PW_QS_1"); - if ( SCAN_BEGIN(ifs,"") ) + if (SCAN_BEGIN(ifs, "")) { // <1> calculate spillage or not. - READ_VALUE( ifs, CALSPI ); + READ_VALUE(ifs, CALSPI); // <1.5> restart from file or not. // mohan add this function 2009-08-27 - READ_VALUE( ifs, RESTART); + READ_VALUE(ifs, RESTART); // <2> output file or not. - READ_VALUE( ifs, this->output ); + READ_VALUE(ifs, this->output); // <3> number of structures. - READ_VALUE( ifs, STRNUM ); - OUT(ofs_running,"STRNUM",STRNUM); - if (STRNUM<=0) + READ_VALUE(ifs, STRNUM); + OUT(ofs_running, "STRNUM", STRNUM); + if (STRNUM <= 0) { - WARNING_QUIT("Read_INPUT::read_PW_QS_1","STRNUM<=0"); + WARNING_QUIT("Read_INPUT::read_PW_QS_1", "STRNUM<=0"); } delete[] qsfile; this->qsfile = new string[STRNUM]; // <4> read in each file. - for (int in=0; in< STRNUM; in++) + for (int in = 0; in < STRNUM; in++) { READ_VALUE(ifs, qsfile[in]); ifstream check(qsfile[in].c_str()); if (!check) { ofs_running << qsfile[in] << endl; - WARNING_QUIT("Read_INPUT","the file is not exist!"); + WARNING_QUIT("Read_INPUT", "the file is not exist!"); } else { - OUT(ofs_running,"FILE",qsfile[in]); - static bool already_read = false; - if (!already_read) - { - check >> LAT0; - OUT(ofs_running,"LAT0",LAT0); - check >> LATVEC.e11 >> LATVEC.e12 >> LATVEC.e13; - check >> LATVEC.e21 >> LATVEC.e22 >> LATVEC.e23; - check >> LATVEC.e31 >> LATVEC.e32 >> LATVEC.e33; - OUT(ofs_running,"a1",LATVEC.e11,LATVEC.e12,LATVEC.e13); - OUT(ofs_running,"a2",LATVEC.e21,LATVEC.e22,LATVEC.e23); - OUT(ofs_running,"a3",LATVEC.e31,LATVEC.e32,LATVEC.e33); - READ_VALUE(check, NTYPE); - OUT(ofs_running,"NTYPE",NTYPE); - assert( NTYPE > 0 ); - NA = new int[NTYPE]; - LABEL = new string[NTYPE]; - CARPOSX = new double*[NTYPE]; - CARPOSY = new double*[NTYPE]; - CARPOSZ = new double*[NTYPE]; - for(int it=0; it> CARPOSX[it][ia] >> CARPOSY[it][ia] >> CARPOSZ[it][ia]; - OUT(ofs_running,"POS",CARPOSX[it][ia],CARPOSY[it][ia],CARPOSZ[it][ia]); - } - } - - READ_VALUE(check, ECUT); - READ_VALUE(check, ECUT_JLQ); - READ_VALUE(check, RCUT); - READ_VALUE(check, SMOOTH); - READ_VALUE(check, SIGMA); - - READ_VALUE(check, TOLERENCE); - READ_VALUE(check, LMAXALL); - READ_VALUE(check, NKSTOT); - READ_VALUE(check, NBANDS); - READ_VALUE(check, NWFCALL); - READ_VALUE(check, NE); - - OUT(ofs_running,"ECUT(Ry)",ECUT); - OUT(ofs_running,"ECUT_JLQ(Ry)",ECUT_JLQ); - OUT(ofs_running,"RCUT(Bohr)",RCUT); - OUT(ofs_running,"SMOOTH",SMOOTH); - OUT(ofs_running,"SIGMA",SIGMA); - OUT(ofs_running,"TOLERENCE",TOLERENCE); - OUT(ofs_running,"LMAXALL",LMAXALL); - OUT(ofs_running,"NKSTOT",NKSTOT); - OUT(ofs_running,"NBANDS",NBANDS); - OUT(ofs_running,"NWFCALL",NWFCALL); - OUT(ofs_running,"NE",NE); - - if(SMOOTH) assert(SIGMA!=0.0); - - already_read = true; - }//only read once. + OUT(ofs_running, "FILE", qsfile[in]); + static bool already_read = false; + if (!already_read) + { + check >> LAT0; + OUT(ofs_running, "LAT0", LAT0); + check >> LATVEC.e11 >> LATVEC.e12 >> LATVEC.e13; + check >> LATVEC.e21 >> LATVEC.e22 >> LATVEC.e23; + check >> LATVEC.e31 >> LATVEC.e32 >> LATVEC.e33; + OUT(ofs_running, "a1", LATVEC.e11, LATVEC.e12, LATVEC.e13); + OUT(ofs_running, "a2", LATVEC.e21, LATVEC.e22, LATVEC.e23); + OUT(ofs_running, "a3", LATVEC.e31, LATVEC.e32, LATVEC.e33); + READ_VALUE(check, NTYPE); + OUT(ofs_running, "NTYPE", NTYPE); + assert(NTYPE > 0); + NA = new int[NTYPE]; + LABEL = new string[NTYPE]; + CARPOSX = new double*[NTYPE]; + CARPOSY = new double*[NTYPE]; + CARPOSZ = new double*[NTYPE]; + for (int it = 0; it < NTYPE; it++) + { + READ_VALUE(check, LABEL[it]); + READ_VALUE(check, NA[it]); + OUT(ofs_running, "LABEL", LABEL[it]); + OUT(ofs_running, "NA", NA[it]); + CARPOSX[it] = new double[NA[it]]; + CARPOSY[it] = new double[NA[it]]; + CARPOSZ[it] = new double[NA[it]]; + for (int ia = 0; ia < NA[it]; ia++) + { + check >> CARPOSX[it][ia] >> CARPOSY[it][ia] >> CARPOSZ[it][ia]; + OUT(ofs_running, "POS", CARPOSX[it][ia], CARPOSY[it][ia], CARPOSZ[it][ia]); + } + } + + READ_VALUE(check, ECUT); + READ_VALUE(check, ECUT_JLQ); + READ_VALUE(check, RCUT); + READ_VALUE(check, SMOOTH); + READ_VALUE(check, SIGMA); + + READ_VALUE(check, TOLERENCE); + READ_VALUE(check, LMAXALL); + READ_VALUE(check, NKSTOT); + READ_VALUE(check, NBANDS); + READ_VALUE(check, NWFCALL); + READ_VALUE(check, NE); + + OUT(ofs_running, "ECUT(Ry)", ECUT); + OUT(ofs_running, "ECUT_JLQ(Ry)", ECUT_JLQ); + OUT(ofs_running, "RCUT(Bohr)", RCUT); + OUT(ofs_running, "SMOOTH", SMOOTH); + OUT(ofs_running, "SIGMA", SIGMA); + OUT(ofs_running, "TOLERENCE", TOLERENCE); + OUT(ofs_running, "LMAXALL", LMAXALL); + OUT(ofs_running, "NKSTOT", NKSTOT); + OUT(ofs_running, "NBANDS", NBANDS); + OUT(ofs_running, "NWFCALL", NWFCALL); + OUT(ofs_running, "NE", NE); + + if (SMOOTH) + assert(SIGMA != 0.0); + + already_read = true; + } // only read once. check.close(); } } @@ -438,77 +441,75 @@ void Read_INPUT::read_PW_QS_1(void) void Read_INPUT::read_WFC(void) { - TITLE("Read_INPUT","read_WFC"); - - ifstream ifswfc( WFC_FILE.c_str() ); - - if(!ifswfc) - { - cout << " Can't find wfc file : " << WFC_FILE << endl; - WARNING_QUIT("Read_INPUT::read_WFC","Can't find wavefunctions file."); - } - else - { - cout << " Find file: " << WFC_FILE << endl; - } - - if(MY_RANK==0) - { - - ifswfc >> LAT0; - cout << "\n LAT0 = " << LAT0; - ifswfc >> LATVEC.e11 >> LATVEC.e12 >> LATVEC.e13; - ifswfc >> LATVEC.e21 >> LATVEC.e22 >> LATVEC.e23; - ifswfc >> LATVEC.e31 >> LATVEC.e32 >> LATVEC.e33; - cout << "\n a1 = " << LATVEC.e11 << " " << LATVEC.e12 << " " << LATVEC.e13; - cout << "\n a2 = " << LATVEC.e21 << " " << LATVEC.e22 << " " << LATVEC.e23; - cout << "\n a3 = " << LATVEC.e31 << " " << LATVEC.e32 << " " << LATVEC.e33; - READ_VALUE(ifswfc, NTYPE); - cout << "\n NTYPE = " << NTYPE << endl; - assert( NTYPE > 0 ); - NA = new int[NTYPE]; - LABEL = new string[NTYPE]; - CARPOSX = new double*[NTYPE]; - CARPOSY = new double*[NTYPE]; - CARPOSZ = new double*[NTYPE]; - for(int it=0; it> LAT0; + cout << "\n LAT0 = " << LAT0; + ifswfc >> LATVEC.e11 >> LATVEC.e12 >> LATVEC.e13; + ifswfc >> LATVEC.e21 >> LATVEC.e22 >> LATVEC.e23; + ifswfc >> LATVEC.e31 >> LATVEC.e32 >> LATVEC.e33; + cout << "\n a1 = " << LATVEC.e11 << " " << LATVEC.e12 << " " << LATVEC.e13; + cout << "\n a2 = " << LATVEC.e21 << " " << LATVEC.e22 << " " << LATVEC.e23; + cout << "\n a3 = " << LATVEC.e31 << " " << LATVEC.e32 << " " << LATVEC.e33; + READ_VALUE(ifswfc, NTYPE); + cout << "\n NTYPE = " << NTYPE << endl; + assert(NTYPE > 0); + NA = new int[NTYPE]; + LABEL = new string[NTYPE]; + CARPOSX = new double*[NTYPE]; + CARPOSY = new double*[NTYPE]; + CARPOSZ = new double*[NTYPE]; + for (int it = 0; it < NTYPE; it++) + { + READ_VALUE(ifswfc, LABEL[it]); + READ_VALUE(ifswfc, NA[it]); + cout << "\n LABEL = " << LABEL[it]; + cout << "\n NA=" << NA[it]; + CARPOSX[it] = new double[NA[it]]; + CARPOSY[it] = new double[NA[it]]; + CARPOSZ[it] = new double[NA[it]]; + for (int ia = 0; ia < NA[it]; ia++) + { + ifswfc >> CARPOSX[it][ia] >> CARPOSY[it][ia] >> CARPOSZ[it][ia]; + cout << "\n " << CARPOSX[it][ia] << " " << CARPOSY[it][ia] << " " << CARPOSZ[it][ia]; + } + } + + READ_VALUE(ifswfc, ECUT); + READ_VALUE(ifswfc, NKSTOT); + READ_VALUE(ifswfc, NBANDS); + + cout << "\n ECUT = " << ECUT; + cout << "\n ECUT_JLQ = " << ECUT_JLQ; + cout << "\n RCUT = " << RCUT; + cout << "\n SMOOTH = " << SMOOTH; + cout << "\n SIGMA = " << SIGMA; + cout << "\n TOLERENCE = " << TOLERENCE; + cout << "\n NKSTOT = " << NKSTOT; + cout << "\n NBANDs = " << NBANDS; + } + + ifswfc.close(); + return; +} void Read_INPUT::read_WFC_2(void) { - TITLE("Read_INPUT","read_WFC_2"); - QUIT(); + TITLE("Read_INPUT", "read_WFC_2"); + QUIT(); } - diff --git a/tools/01_NAO_generation/SIAB/src_spillage/read_INPUT.h b/tools/01_NAO_generation/SIAB/src_spillage/read_INPUT.h index 001334396d6..e89e7053c28 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/read_INPUT.h +++ b/tools/01_NAO_generation/SIAB/src_spillage/read_INPUT.h @@ -1,48 +1,47 @@ #ifndef READ_INPUT_H #define READ_INPUT_H -#include "common.h" #include "Coefficients.h" #include "ReadData.h" #include "SpillageValue.h" +#include "common.h" // read information from file : INPUT class Read_INPUT { - public: - Read_INPUT(); - ~Read_INPUT(); - - ReadData* QS_data; - Coefficients Coef; - SpillageValue SV; + public: + Read_INPUT(); + ~Read_INPUT(); - ifstream ifs; - ifstream ifs2; - bool cal_sp; - bool output; + ReadData* QS_data; + Coefficients Coef; + SpillageValue SV; - // get this value from Multi_Zeta.cpp - int nlevel; + ifstream ifs; + ifstream ifs2; + bool cal_sp; + bool output; - void init(void); + // get this value from Multi_Zeta.cpp + int nlevel; - void read_test(void); - void read_plot_c4(void); -// void read_cal_c4(void); - void read_start_c4(const int &step); + void init(void); - private: - int test; + void read_test(void); + void read_plot_c4(void); + // void read_cal_c4(void); + void read_start_c4(const int& step); - void bcast(void); - void read_PW_QS_1(void); - void read_PW_QS_2(void); - void read_WFC(void); - void read_WFC_2(void); + private: + int test; - string *qsfile; + void bcast(void); + void read_PW_QS_1(void); + void read_PW_QS_2(void); + void read_WFC(void); + void read_WFC_2(void); + string* qsfile; }; #endif diff --git a/tools/01_NAO_generation/SIAB/src_spillage/tools.cpp b/tools/01_NAO_generation/SIAB/src_spillage/tools.cpp index 19098635486..d12ee6c425e 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/tools.cpp +++ b/tools/01_NAO_generation/SIAB/src_spillage/tools.cpp @@ -1,6 +1,6 @@ #include "tools.h" -//plane wave , mohan add 2010-06-14 +// plane wave , mohan add 2010-06-14 PW_Basis PW; MultiZeta mz; Read_INPUT input; @@ -25,11 +25,11 @@ int NE = -1; // global information int STRNUM = -1; int NTYPE = -1; -string *LABEL; -int *NA; -double **CARPOSX; -double **CARPOSY; -double **CARPOSZ; +string* LABEL; +int* NA; +double** CARPOSX; +double** CARPOSY; +double** CARPOSZ; Matrix3 LATVEC; // plane wave information @@ -37,15 +37,15 @@ bool USEPW = true; double LAT0 = -1.0; // k points -double *CARKX; -double *CARKY; -double *CARKZ; +double* CARKX; +double* CARKY; +double* CARKZ; // orbital + energy cutoff double ECUT = -1.0; double ECUT_JLQ = -1.0; double RCUT = -1.0; -bool SMOOTH = false; +bool SMOOTH = false; double SIGMA = -0.1; double TOLERENCE = 1.0e-12; int LMAXALL = -1; @@ -54,9 +54,9 @@ int NMAXUSED = -1; int NCHIUSED = -1; // band information -bool BANDS_CONTROL=false; // mohan add 2010-05-02 -int BANDS_START=0; -int BANDS_END=0; // needed to be checked carefully if not read in. +bool BANDS_CONTROL = false; // mohan add 2010-05-02 +int BANDS_START = 0; +int BANDS_END = 0; // needed to be checked carefully if not read in. // other global information int CALSPI = 0; // mohan 2010-06-17 @@ -71,16 +71,16 @@ double ACCEPTANCE_HIGH = 0.6; double ACCEPTANCE_LOW = 0.3; double KINETIC_MAX = 30.0; // (unit?????) -double KINETIC_DR = 0.01; // -int OPTIMIZE_METHOD = 1; // 1: Kin 2: Ecut +double KINETIC_DR = 0.01; // +int OPTIMIZE_METHOD = 1; // 1: Kin 2: Ecut // constant const double PI = 3.1415926535897; const double TWO_PI = 2.0 * PI; const double FOUR_PI = 4.0 * PI; -const double SQRT_INVERSE_FOUR_PI = sqrt(1.0/FOUR_PI); +const double SQRT_INVERSE_FOUR_PI = sqrt(1.0 / FOUR_PI); const std::complex IMAG_UNIT = std::complex(0, 1.0); -const std::complex NEG_IMAG_UNIT = std::complex(0,-1.0); +const std::complex NEG_IMAG_UNIT = std::complex(0, -1.0); const double PI_HALF = PI / 2.0; const double SQRT2 = 1.41421356237309504880; @@ -89,27 +89,24 @@ string WFC_FILE = "WAVEFUNC.dat"; void DESTROY(void) { - if(NTYPE>0) - { - delete[] LABEL; - delete[] NA; - for(int i=0; i0) - { - delete[] CARKX; - delete[] CARKY; - delete[] CARKZ; - } + if (NTYPE > 0) + { + delete[] LABEL; + delete[] NA; + for (int i = 0; i < NTYPE; i++) + { + delete[] CARPOSX[i]; + delete[] CARPOSY[i]; + delete[] CARPOSZ[i]; + } + delete[] CARPOSX; + delete[] CARPOSY; + delete[] CARPOSZ; + } + if (NKS > 0) + { + delete[] CARKX; + delete[] CARKY; + delete[] CARKZ; + } } - - - diff --git a/tools/01_NAO_generation/SIAB/src_spillage/tools.h b/tools/01_NAO_generation/SIAB/src_spillage/tools.h index 21486a12ae0..6f8b5111ee4 100644 --- a/tools/01_NAO_generation/SIAB/src_spillage/tools.h +++ b/tools/01_NAO_generation/SIAB/src_spillage/tools.h @@ -1,12 +1,12 @@ #ifndef TOOL_H #define TOOL_H -#include "common.h" +#include "../src_parallel/parallel_kpoints.h" // add 2010-06-12 #include "../src_pw/pw_basis.h" #include "MultiZeta.h" +#include "common.h" #include "read_INPUT.h" -#include "../src_parallel/parallel_kpoints.h" // add 2010-06-12 -//plane wave , mohan add 2010-06-14 +// plane wave , mohan add 2010-06-14 extern PW_Basis PW; extern MultiZeta mz; extern Read_INPUT input; @@ -31,11 +31,11 @@ extern int NE; // global information extern int STRNUM; extern int NTYPE; -extern string *LABEL; -extern int *NA; -extern double **CARPOSX; -extern double **CARPOSY; -extern double **CARPOSZ; +extern string* LABEL; +extern int* NA; +extern double** CARPOSX; +extern double** CARPOSY; +extern double** CARPOSZ; extern Matrix3 LATVEC; // plane wave information @@ -43,11 +43,11 @@ extern bool USEPW; extern double LAT0; // k points -extern double *CARKX; -extern double *CARKY; -extern double *CARKZ; +extern double* CARKX; +extern double* CARKY; +extern double* CARKZ; -//orbital + energy cutoff +// orbital + energy cutoff extern double ECUT; extern double ECUT_JLQ; extern double RCUT; @@ -61,12 +61,12 @@ extern int NCHIUSED; // band information extern bool BANDS_CONTROL; -extern int BANDS_START;//mohan add 2010-05-02 +extern int BANDS_START; // mohan add 2010-05-02 extern int BANDS_END; // other global information -extern int CALSPI; // mohan 2010-06-17 -extern int BLOCK_NE; // if ie < BLOCK_NE, we use this Jlq(ie), mohan add 2009-08-26 +extern int CALSPI; // mohan 2010-06-17 +extern int BLOCK_NE; // if ie < BLOCK_NE, we use this Jlq(ie), mohan add 2009-08-26 extern int BLOCK_NE_MIN; // mohan add 2009-08-27 extern bool RESTART; extern int TEST1; @@ -76,16 +76,16 @@ extern int ACCEPTANCE_STEPS; // mohan add 2009-10-31 extern double ACCEPTANCE_HIGH; extern double ACCEPTANCE_LOW; -extern double KINETIC_MAX; // mohan add 2009-10-31 -extern double KINETIC_DR; // mohan add 2010-04-12 -extern int OPTIMIZE_METHOD;// mohan add 2010-04-14 +extern double KINETIC_MAX; // mohan add 2009-10-31 +extern double KINETIC_DR; // mohan add 2010-04-12 +extern int OPTIMIZE_METHOD; // mohan add 2010-04-14 // constant -extern const double PI; // mohan add 2010-04-16 +extern const double PI; // mohan add 2010-04-16 extern const double TWO_PI; // mohan add 2010-06-14 extern const double FOUR_PI; extern const std::complex IMAG_UNIT; -extern const std::complex NEG_IMAG_UNIT;//mohan add 2010-06-14 +extern const std::complex NEG_IMAG_UNIT; // mohan add 2010-06-14 extern const double SQRT_INVERSE_FOUR_PI; extern const double PI_HALF; extern const double SQRT2;