diff --git a/source/module_cell/module_symmetry/symmetry.cpp b/source/module_cell/module_symmetry/symmetry.cpp index 73009e47fd1..0f59bfe810f 100644 --- a/source/module_cell/module_symmetry/symmetry.cpp +++ b/source/module_cell/module_symmetry/symmetry.cpp @@ -25,7 +25,7 @@ int Symmetry::symm_flag=0; void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running) { - const double MAX_EPS = std::max(1e-3, epsilon); + const double MAX_EPS = std::max(1e-3, epsilon_input * 1.001); const double MULT_EPS = 2.0; if (available == false) return; ModuleBase::TITLE("Symmetry","init"); @@ -99,65 +99,114 @@ void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running) s2 = a2; s3 = a3; - //a: optimized config - // find the lattice type accordiing to lattice vectors. - this->lattice_type(this->a1, this->a2, this->a3, this->s1, this->s2, this->s3, - this->cel_const, this->pre_const, this->real_brav, ilattname, ucell, true, this->newpos); - - ofs_running<<"(for optimal symmetric configuration:)"<pricell(this->newpos); // pengfei Li 2018-05-14 - //for( iat =0 ; iat < ucell.nat ; iat++) + auto lattice_to_group = [&, this](int& nrot_out, int& nrotk_out, std::ofstream& ofs_running) -> void { + //a: optimized config + // find the lattice type accordiing to lattice vectors. + this->lattice_type(this->a1, this->a2, this->a3, this->s1, this->s2, this->s3, + this->cel_const, this->pre_const, this->real_brav, ilattname, ucell, true, this->newpos); + + ofs_running << "(for optimal symmetric configuration:)" << std::endl; + ModuleBase::GlobalFunc::OUT(ofs_running, "BRAVAIS TYPE", real_brav); + ModuleBase::GlobalFunc::OUT(ofs_running, "BRAVAIS LATTICE NAME", ilattname); + ModuleBase::GlobalFunc::OUT(ofs_running, "ibrav", real_brav); + Symm_Other::print1(real_brav, cel_const, ofs_running); + // std::cout << "a1 = " << a1.x << " " << a1.y << " " << a1.z <pricell(this->newpos); // pengfei Li 2018-05-14 + //for( iat =0 ; iat < ucell.nat ; iat++) // std::cout << " newpos_now = " << newpos[3*iat] << " " << newpos[3*iat+1] << " " << newpos[3*iat+2] << std::endl; - test_brav = true; // output the real ibrav and point group - this->setgroup(this->symop, this->nop, this->real_brav); - + test_brav = true; // output the real ibrav and point group + this->setgroup(this->symop, this->nop, this->real_brav); + this->getgroup(nrot_out, nrotk_out, ofs_running); + }; + if (GlobalV::CALCULATION == "cell-relax" && nrotk > 0) { + std::ofstream no_out; // to screen the output when trying new epsilon + if (this->nrotk > this->max_nrotk)this->max_nrotk = this->nrotk; int tmp_nrot, tmp_nrotk; - this->getgroup(tmp_nrot, tmp_nrotk, ofs_running); + lattice_to_group(tmp_nrot, tmp_nrotk, ofs_running); //some different method to enlarge symmetry_prec - bool eps_changed = false; + bool eps_enlarged = false; auto eps_mult = [this](double mult) {epsilon *= mult;}; auto eps_to = [this](double new_eps) {epsilon = new_eps;}; + + // store the symmetry_prec and nrotk for each try + std::vector precs_try; + std::vector nrotks_try; + // store the initial result + precs_try.push_back(epsilon); + nrotks_try.push_back(tmp_nrotk); //enlarge epsilon and regenerate pointgroup - while (tmp_nrotk < this->nrotk && epsilon < MAX_EPS) + while (tmp_nrotk < this->max_nrotk && epsilon < MAX_EPS) { eps_mult(MULT_EPS); - eps_changed = true; - this->getgroup(tmp_nrot, tmp_nrotk, ofs_running); + eps_enlarged = true; + // lattice_to_group(tmp_nrot, tmp_nrotk, no_out); + lattice_to_group(tmp_nrot, tmp_nrotk, no_out); + precs_try.push_back(epsilon); + nrotks_try.push_back(tmp_nrotk); } - if (epsilon > MAX_EPS) + if (tmp_nrotk > this->nrotk) { - ofs_running << "ERROR: Symmetry cannot be kept due to the lost of accuracy with atom position during cell-relax." << std::endl; - ofs_running << "Please set `symmetry` to 0 or -1 in INPUT file. " << std::endl; - ModuleBase::QUIT(); + this->nrotk = tmp_nrotk; + ofs_running << "Find new symmtry operations during cell-relax." << std::endl; + if (this->nrotk > this->max_nrotk)this->max_nrotk = this->nrotk; } - - if (eps_changed) + if (eps_enlarged) { - ofs_running << "WARNING: current `symmetry_prec` is too small to give the right number of symmtry operations." << std::endl; - ofs_running << " Changed `symmetry_prec` to " << epsilon <<"." << std::endl; + if (epsilon > MAX_EPS) + { + ofs_running << "WARNING: Symmetry cannot be kept due to the lost of accuracy with atom position during cell-relax." << std::endl; + ofs_running << "Continue cell-relax with a lower symmetry. " << std::endl; + // find the smallest epsilon that gives the current number of symmetry operations + int valid_index = nrotks_try.size() - 1; + while (valid_index > 0 && tmp_nrotk <= nrotks_try[valid_index - 1])--valid_index; + eps_to(precs_try[valid_index]); + if (valid_index > 0) ofs_running << " Enlarging `symmetry_prec` to " << epsilon << " ..." << std::endl; + else eps_enlarged = false; + // regenerate pointgroup after change epsilon (may not give the same result) + lattice_to_group(tmp_nrot, tmp_nrotk, ofs_running); + this->nrotk = tmp_nrotk; + } + else ofs_running << " Enlarging `symmetry_prec` to " << epsilon << " ..." << std::endl; } - if (tmp_nrotk > this->nrotk) - { + if (!eps_enlarged && epsilon > epsilon_input * 1.001) // not "else" here. "eps_enlarged" can be set to false in the above "if" + { // try a smaller symmetry_prec until the number of symmetry operations decreases + precs_try.erase(precs_try.begin() + 1, precs_try.end()); + nrotks_try.erase(nrotks_try.begin() + 1, nrotks_try.end()); + double eps_current = epsilon; // record the current symmetry_prec + do { + eps_mult(1 / MULT_EPS); + lattice_to_group(tmp_nrot, tmp_nrotk, no_out); + precs_try.push_back(epsilon); + nrotks_try.push_back(tmp_nrotk); + } while (tmp_nrotk >= nrotks_try[0] && epsilon > epsilon_input * 1.001 && precs_try.size() < 5); + int valid_index = (tmp_nrotk < nrotks_try[0]) ? nrotks_try.size() - 2 : nrotks_try.size() - 1; +#ifdef __DEBUG + assert(valid_index >= 0); + assert(nrotks_try[valid_index] >= nrotks_try[0]); +#endif + epsilon = precs_try[valid_index]; + // regenerate pointgroup after change epsilon + lattice_to_group(tmp_nrot, tmp_nrotk, ofs_running); this->nrotk = tmp_nrotk; - ofs_running << "Find new symmtry operations during cell-relax." << std::endl; + if (valid_index > 0)//epsilon is set smaller + ofs_running << " Narrowing `symmetry_prec` from " << eps_current << " to " << epsilon << " ..." << std::endl; } - assert(tmp_nrotk == this->nrotk); + // final number of symmetry operations +#ifdef __DEBUG + ofs_running << "symmetry_prec(epsilon) in current ion step: " << this->epsilon << std::endl; + ofs_running << "number of symmetry operations in current ion step: " << this->nrotk << std::endl; +#endif } else - this->getgroup(this->nrot, this->nrotk, ofs_running); + lattice_to_group(this->nrot, this->nrotk, ofs_running); this->pointgroup(this->nrot, this->pgnumber, this->pgname, this->gmatrix, ofs_running); ModuleBase::GlobalFunc::OUT(ofs_running,"POINT GROUP", this->pgname); diff --git a/source/module_cell/module_symmetry/symmetry.h b/source/module_cell/module_symmetry/symmetry.h index 5959c34f8de..5881f0f8188 100644 --- a/source/module_cell/module_symmetry/symmetry.h +++ b/source/module_cell/module_symmetry/symmetry.h @@ -62,8 +62,9 @@ class Symmetry : public Symmetry_Basic int nop; //the number of point group operations of the pure bravais lattice without basis int s_flag; //whether the current matrix is one of all space group operations int nrot; //the number of pure point group rotations - int nrotk=-1; //the number of all space group operations - int pgnumber; //the serial number of point group + int nrotk = -1; //the number of all space group operations + int max_nrotk = -1; ///< record the maximum number of symmetry operations during cell-relax + int pgnumber; //the serial number of point group int spgnumber; //the serial number of point group in space group std::string pgname; //the Schoenflies name of the point group R in {R|0} std::string spgname; //the Schoenflies name of the point group R in the space group {R|t} diff --git a/source/module_cell/module_symmetry/symmetry_basic.h b/source/module_cell/module_symmetry/symmetry_basic.h index fbc94cb3a4f..1bfa581ba4c 100644 --- a/source/module_cell/module_symmetry/symmetry_basic.h +++ b/source/module_cell/module_symmetry/symmetry_basic.h @@ -17,8 +17,9 @@ class Symmetry_Basic Symmetry_Basic(); ~Symmetry_Basic(); - double epsilon; - + double epsilon; + double epsilon_input; ///< the input value of symmetry_prec, should not be changed + // control accuray bool equal(const double &m, const double &n)const; void check_boundary(double &x)const; diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index 0b89c89cda9..5ec2964bfea 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -15,8 +15,8 @@ namespace ModuleESolver pw_big->setbxyz(INPUT.bx, INPUT.by, INPUT.bz); sf.set(INPUT.nbspline); - this->symm.epsilon = INPUT.symmetry_prec; - } + this->symm.epsilon = this->symm.epsilon_input = INPUT.symmetry_prec; +} ESolver_FP::~ESolver_FP() { if (this->psi != nullptr)