From 2fad65295495d6946560ff8142d0a3e79297f5cf Mon Sep 17 00:00:00 2001 From: mohanchen Date: Tue, 25 Mar 2025 14:18:04 +0800 Subject: [PATCH 01/12] update symmetry output format --- .../module_cell/module_symmetry/symmetry.cpp | 120 +++++---- .../module_symmetry/symmetry_basic.cpp | 242 ++---------------- tests/integrate/801_PW_LT_sc/INPUT | 6 +- 3 files changed, 91 insertions(+), 277 deletions(-) diff --git a/source/module_cell/module_symmetry/symmetry.cpp b/source/module_cell/module_symmetry/symmetry.cpp index e9b9624b224..e422de6e655 100644 --- a/source/module_cell/module_symmetry/symmetry.cpp +++ b/source/module_cell/module_symmetry/symmetry.cpp @@ -20,18 +20,18 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, const double MAX_EPS = std::max(1e-3, epsilon_input * 1.001); const double MULT_EPS = 2.0; - ModuleBase::TITLE("Symmetry","init"); + ModuleBase::TITLE("Symmetry","analy_sys"); ModuleBase::timer::tick("Symmetry","analy_sys"); ofs_running << "\n\n\n\n"; ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; ofs_running << " | |" << std::endl; - ofs_running << " | Doing symmetry analysis: |" << std::endl; + ofs_running << " | Performing symmetry analysis: |" << std::endl; ofs_running << " | We calculate the norm of 3 vectors and the angles between them, |" << std::endl; ofs_running << " | the type of Bravais lattice is given. We can judge if the unticell |" << std::endl; ofs_running << " | is a primitive cell. Finally we give the point group operation for |" << std::endl; - ofs_running << " | this unitcell. We use the point group operations to do symmetry |" << std::endl; - ofs_running << " | analysis on given k-point mesh and the charge density. |" << std::endl; + ofs_running << " | this unitcell. We use the point group operations to perform |" << std::endl; + ofs_running << " | symmetry analysis on given k-point mesh and the charge density. |" << std::endl; ofs_running << " | |" << std::endl; ofs_running << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; ofs_running << "\n\n\n\n"; @@ -43,9 +43,13 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, this->nat = st.nat; // number of atom species this->ntype = st.ntype; + + assert(ntype>0); + this->na = new int[ntype]; this->istart = new int[ntype]; // start number of atom. this->index = new int [nat + 2]; // index of atoms + ModuleBase::GlobalFunc::ZEROS(na, ntype); ModuleBase::GlobalFunc::ZEROS(istart, ntype); ModuleBase::GlobalFunc::ZEROS(index, nat+2); @@ -91,52 +95,55 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, s3 = a3; - auto lattice_to_group = [&, this](int& nrot_out, int& nrotk_out, std::ofstream& ofs_running) -> void { - // a: the optimized lattice vectors, output - // s: the input lattice vectors, input - // find the real_brav 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, atoms, 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); - - optlat.e11 = a1.x; optlat.e12 = a1.y; optlat.e13 = a1.z; - optlat.e21 = a2.x; optlat.e22 = a2.y; optlat.e23 = a2.z; - optlat.e31 = a3.x; optlat.e32 = a3.y; optlat.e33 = a3.z; - - // count the number of primitive cells in the supercell - this->pricell(this->newpos, atoms); - - test_brav = true; // output the real ibrav and point group - - // list all possible point group operations - this->setgroup(this->symop, this->nop, this->real_brav); - - // special case for AFM analysis - // which should be loop over all atoms, f.e only loop over spin-up atoms - // -------------------------------- - // AFM analysis Start - if (PARAM.inp.nspin > 1) { - pricell_loop = this->magmom_same_check(atoms); - } + auto lattice_to_group = [&, this](int& nrot_out, int& nrotk_out, std::ofstream& ofs_running) -> void + { + // a: the optimized lattice vectors, output + // s: the input lattice vectors, input + // find the real_brav 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, atoms, 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); + + optlat.e11 = a1.x; optlat.e12 = a1.y; optlat.e13 = a1.z; + optlat.e21 = a2.x; optlat.e22 = a2.y; optlat.e23 = a2.z; + optlat.e31 = a3.x; optlat.e32 = a3.y; optlat.e33 = a3.z; + + // count the number of primitive cells in the supercell + this->pricell(this->newpos, atoms); + + test_brav = true; // output the real ibrav and point group + + // list all possible point group operations + this->setgroup(this->symop, this->nop, this->real_brav); + + // special case for AFM analysis + // which should be loop over all atoms, f.e only loop over spin-up atoms + // -------------------------------- + // AFM analysis Start + if (PARAM.inp.nspin > 1) + { + pricell_loop = this->magmom_same_check(atoms); + } - if (!pricell_loop && PARAM.inp.nspin == 2) - { - this->analyze_magnetic_group(atoms, st, nrot_out, nrotk_out); - } - else - { - // get the real symmetry operations according to the input structure - // nrot_out: the number of pure point group rotations - // nrotk_out: the number of all space group operations - this->getgroup(nrot_out, nrotk_out, ofs_running, this->nop, this->symop, this->gmatrix, this->gtrans, - this->newpos, this->rotpos, this->index, this->ntype, this->itmin_type, this->itmin_start, this->istart, this->na); - } - }; + if (!pricell_loop && PARAM.inp.nspin == 2) + { + this->analyze_magnetic_group(atoms, st, nrot_out, nrotk_out); + } + else + { + // get the real symmetry operations according to the input structure + // nrot_out: the number of pure point group rotations + // nrotk_out: the number of all space group operations + this->getgroup(nrot_out, nrotk_out, ofs_running, this->nop, this->symop, + this->gmatrix, this->gtrans, this->newpos, this->rotpos, this->index, + this->ntype, this->itmin_type, this->itmin_start, this->istart, this->na); + } + }; // -------------------------------- // 2. analyze the symmetry @@ -181,17 +188,18 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, if (tmp_nrotk > this->nrotk) { 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; - } + ofs_running << " Find new symmtry operations during cell-relax." << std::endl; + if (this->nrotk > this->max_nrotk) + { + this->max_nrotk = this->nrotk; + } } if (eps_enlarged) { 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; + 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 @@ -1220,7 +1228,7 @@ void Symmetry::pricell(double* pos, const Atom* atoms) int ntrans=ptrans.size(); if (ntrans <= 1) { - GlobalV::ofs_running<<"Original cell was already a primitive cell."<p1=this->a1; this->p2=this->a2; this->p3=this->a3; diff --git a/source/module_cell/module_symmetry/symmetry_basic.cpp b/source/module_cell/module_symmetry/symmetry_basic.cpp index f0c7cd93144..dffa59f69d7 100644 --- a/source/module_cell/module_symmetry/symmetry_basic.cpp +++ b/source/module_cell/module_symmetry/symmetry_basic.cpp @@ -1,10 +1,7 @@ -//========================================================== -// AUTHOR : Zhengpan , mohan , spshu -// DATE : 2007-9 -//========================================================== #include "symmetry.h" #include "module_base/mymath.h" #include "module_parameter/parameter.h" + bool ModuleSymmetry::test_brav = 0; namespace ModuleSymmetry @@ -111,162 +108,7 @@ void Symmetry_Basic::order_atoms(double* pos, const int& nat, const int* index) return; } -/* -// atom ordering for each atom type. -void Symmetry_Basic::atom_ordering(double *posi, const int natom, int *subindex) -{ - //order the atomic positions inside a supercell by a unique ordering scheme - subindex[0] = 0; - - if(natom == 1) - { - //if there is only one atom, it is not necessary to order - return; - } - - //order the x-dimension of the atomic position, - //only get a permutation table, data in momery are not changed - //this->heapsort_pos(natom, posi, subindex); - - double *tmppos = new double[natom]; - for(int i=0; iorder_y(&posi[oldpos*3], oldpos, newpos); - oldv = vnew; - oldpos = ia; - } - else if(ia==natom-1) - { - // in this case, the new one is the last one, - // also the same as the old one, - // so we need to add 1 - this->order_y(&posi[oldpos*3], oldpos, newpos+1); - } - } - delete[] tmppos; - return; -} - -void Symmetry_Basic::order_y(double *pos1, const int &oldpos1, const int &newpos1) -{ -// ModuleBase::TITLE("Symmetry_Basic","order_y"); - // how many atoms need to be reordered according to same x value. - const int nat1 = newpos1 - oldpos1; -// std::cout << "\n nat1=" << nat1 << std::endl; - if(nat1 == 1) return; - - double* tmp1 = new double[nat1]; - int* index1 = new int[nat1]; - ModuleBase::GlobalFunc::ZEROS(index1, nat1); - - for(int ia=0; iaorder_atoms(pos2,nat2,index2); - -// for(int ia=0; ia 48) { @@ -447,8 +261,6 @@ void Symmetry_Basic::matrigen(ModuleBase::Matrix3 *symgen, const int ngen, Modul return; } symop[now - 1] = temp2; - // std::cout<<"\n symop[now]: "< 24) @@ -732,7 +532,10 @@ int Symmetry_Basic::subgroup(const int& nrot, const int& ninv, const int& nc2, c } return 1;//C_1 } -bool Symmetry_Basic::pointgroup(const int& nrot, int& pgnumber, std::string& pgname, const ModuleBase::Matrix3* gmatrix, std::ofstream& ofs_running)const + + +bool Symmetry_Basic::pointgroup(const int& nrot, int& pgnumber, + std::string& pgname, const ModuleBase::Matrix3* gmatrix, std::ofstream& ofs_running)const { //------------------------------------------------------------------------- //return the name of the point group @@ -751,8 +554,10 @@ bool Symmetry_Basic::pointgroup(const int& nrot, int& pgnumber, std::string& pgn //because the number of their elements are exclusive if (PARAM.inp.test_symmetry) ModuleBase::TITLE("Symmetry_Basic", "pointgroup"); - std::vector pgdict = { "none", "C_1", "S_2", "C_2", "C_1h", "C_2h", "D_2", "C_2v", "D_2h", "C_3", "S_6", "D_3", "C_3v", "D_3d", "C_4", "S_4", "C_4h", "D_4", "C_4v", "D_2d", "D_4h", - "C_6", "C_3h", "C_6h", "D_6", "C_6v", "D_3h", "D_6h", "T", "T_h", "O", "T_d", "O_h" }; + std::vector pgdict = { "none", "C_1", "S_2", "C_2", "C_1h", "C_2h", + "D_2", "C_2v", "D_2h", "C_3", "S_6", "D_3", "C_3v", "D_3d", "C_4", "S_4", "C_4h", + "D_4", "C_4v", "D_2d", "D_4h", "C_6", "C_3h", "C_6h", "D_6", "C_6v", "D_3h", "D_6h", + "T", "T_h", "O", "T_d", "O_h" }; if(nrot == 1) { @@ -805,7 +610,6 @@ bool Symmetry_Basic::pointgroup(const int& nrot, int& pgnumber, std::string& pgn int ns4 = 0; int ns6 = 0; //mohan add 2012-01-15 -// GlobalV::ofs_running << " " << std::setw(5) << "NROT" << std::setw(15) << "TRACE" << std::setw(15) << "DET" << std::endl; for(int i = 0; i < nrot; ++i) { //calculate the trace of a matrix @@ -813,8 +617,6 @@ bool Symmetry_Basic::pointgroup(const int& nrot, int& pgnumber, std::string& pgn //calculate the determinant of a matrix det = int(gmatrix[i].Det()); -// GlobalV::ofs_running << " " << std::setw(5) << i+1 << std::setw(15) << trace << std::setw(15) << det << std::endl; - if(trace == 3) { continue; //found unity operator (trivial) @@ -836,7 +638,7 @@ bool Symmetry_Basic::pointgroup(const int& nrot, int& pgnumber, std::string& pgn else if(trace == -2 && det == -1) ++ns3; //mohan add 2012-01-15 } - if(test_brav) + if(test_brav) { ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"C2",nc2); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"C3",nc3); @@ -1109,7 +911,7 @@ void Symmetry_Basic::atom_ordering_new(double *posi, const int natom, int *subin double z_max = *max_element(tmpz.begin(),tmpz.end()); double z_min = *min_element(tmpz.begin(),tmpz.end()); - double* weighted_func = new double[natom]; + double* weighted_func = new double[natom]; //the first time: f(x, y, z) for(int i=0; i1) //need a new sort { @@ -1153,4 +959,4 @@ void Symmetry_Basic::test_atom_ordering(double *posi, const int natom, int *subi //an interface to test a protected function this->atom_ordering_new(posi, natom, subindex); } -} \ No newline at end of file +} diff --git a/tests/integrate/801_PW_LT_sc/INPUT b/tests/integrate/801_PW_LT_sc/INPUT index 92b92f06145..f24f5a38334 100644 --- a/tests/integrate/801_PW_LT_sc/INPUT +++ b/tests/integrate/801_PW_LT_sc/INPUT @@ -1,12 +1,12 @@ INPUT_PARAMETERS #Parameters (System) -suffix autotest +suffix autotest nbands 4 stru_file STRU kpoint_file KPT -pseudo_dir ../../PP_ORB -calculation scf +pseudo_dir ../../PP_ORB +calculation scf latname sc #Parameters (PW) ecutwfc 25.0 # Rydberg From 6d8f7f50b1ac972d64e101e7af9a04d2d38b6292 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Tue, 25 Mar 2025 15:08:56 +0800 Subject: [PATCH 02/12] update symmetry --- .../module_symmetry/symm_check.cpp | 159 ++ .../module_symmetry/symm_getgroup.cpp | 105 + .../module_symmetry/symm_lattice.cpp | 451 +++++ .../module_symmetry/symm_magnetic.cpp | 94 + .../module_cell/module_symmetry/symm_other.h | 2 + .../module_symmetry/symm_pricell.cpp | 319 +++ .../module_cell/module_symmetry/symm_rho.cpp | 299 +++ .../module_cell/module_symmetry/symmetry.cpp | 1734 ++--------------- source/module_cell/module_symmetry/symmetry.h | 70 +- 9 files changed, 1618 insertions(+), 1615 deletions(-) create mode 100644 source/module_cell/module_symmetry/symm_check.cpp create mode 100644 source/module_cell/module_symmetry/symm_getgroup.cpp create mode 100644 source/module_cell/module_symmetry/symm_lattice.cpp create mode 100644 source/module_cell/module_symmetry/symm_magnetic.cpp create mode 100644 source/module_cell/module_symmetry/symm_pricell.cpp create mode 100644 source/module_cell/module_symmetry/symm_rho.cpp diff --git a/source/module_cell/module_symmetry/symm_check.cpp b/source/module_cell/module_symmetry/symm_check.cpp new file mode 100644 index 00000000000..5cb72e3d4cc --- /dev/null +++ b/source/module_cell/module_symmetry/symm_check.cpp @@ -0,0 +1,159 @@ + +bool Symmetry::checksym(const ModuleBase::Matrix3 &s, + ModuleBase::Vector3& gtrans, + double* pos, double* rotpos, int* index, + const int ntype, const int itmin_type, const int itmin_start, + int* istart, int* na)const +{ + //---------------------------------------------- + // checks whether a point group symmetry element + // is a valid symmetry operation on a supercell + //---------------------------------------------- + // the start atom index. + bool no_diff = false; + ModuleBase::Vector3 trans(2.0, 2.0, 2.0); + bool s_flag = false; + + for (int it = 0; it < ntype; it++) + { + //------------------------------------ + // impose periodic boundary condition + // 0.5 -> -0.5 + //------------------------------------ + for (int j = istart[it]; j < istart[it] + na[it]; ++j) + { + this->check_boundary(pos[j*3+0]); + this->check_boundary(pos[j*3+1]); + this->check_boundary(pos[j*3+2]); + } + //order original atomic positions for current species + this->atom_ordering_new(pos + istart[it] * 3, na[it], index + istart[it]); + + //Rotate atoms of current species + for (int j = istart[it]; j < istart[it] + na[it]; ++j) + { + const int xx=j*3; + const int yy=j*3+1; + const int zz=j*3+2; + + + rotpos[xx] = pos[xx] * s.e11 + + pos[yy] * s.e21 + + pos[zz] * s.e31; + + rotpos[yy] = pos[xx] * s.e12 + + pos[yy] * s.e22 + + pos[zz] * s.e32; + + rotpos[zz] = pos[xx] * s.e13 + + pos[yy] * s.e23 + + pos[zz] * s.e33; + + rotpos[xx] = fmod(rotpos[xx] + 100.5,1) - 0.5; + rotpos[yy] = fmod(rotpos[yy] + 100.5,1) - 0.5; + rotpos[zz] = fmod(rotpos[zz] + 100.5,1) - 0.5; + this->check_boundary(rotpos[xx]); + this->check_boundary(rotpos[yy]); + this->check_boundary(rotpos[zz]); + } + //order rotated atomic positions for current species + this->atom_ordering_new(rotpos + istart[it] * 3, na[it], index + istart[it]); + } + + ModuleBase::Vector3 diff; + + //--------------------------------------------------------- + // itmin_start = the start atom positions of species itmin + //--------------------------------------------------------- + // (s)tart (p)osition of atom (t)ype which has (min)inal number. + ModuleBase::Vector3 sptmin(rotpos[itmin_start * 3], rotpos[itmin_start * 3 + 1], rotpos[itmin_start * 3 + 2]); + + for (int i = itmin_start; i < itmin_start + na[itmin_type]; ++i) + { + //set up the current test std::vector "gtrans" + //and "gtrans" could possibly contain trivial translations: + gtrans.x = this->get_translation_vector( sptmin.x, pos[i*3+0]); + gtrans.y = this->get_translation_vector( sptmin.y, pos[i*3+1]); + gtrans.z = this->get_translation_vector( sptmin.z, pos[i*3+2]); + + //If we had already detected some translation, + //we must only look at the vectors with coordinates smaller than those + //of the previously detected std::vector (find the smallest) + if (gtrans.x > trans.x + epsilon || + gtrans.y > trans.y + epsilon || + gtrans.z > trans.z + epsilon + ) + { + continue; + } + + //translate all the atomic coordinates BACK by "gtrans" + for (int it = 0; it < ntype; it++) + { + for (int ia = istart[it]; ia < na[it] + istart[it]; ia++) + { + this->check_translation( rotpos[ia*3+0], gtrans.x ); + this->check_translation( rotpos[ia*3+1], gtrans.y ); + this->check_translation( rotpos[ia*3+2], gtrans.z ); + + this->check_boundary( rotpos[ia*3+0] ); + this->check_boundary( rotpos[ia*3+1] ); + this->check_boundary( rotpos[ia*3+2] ); + } + //order translated atomic positions for current species + this->atom_ordering_new(rotpos + istart[it] * 3, na[it], index + istart[it]); + } + + no_diff = true; + //compare the two lattices 'one-by-one' whether they are identical + for (int it = 0; it < ntype; it++) + { + for (int ia = istart[it]; ia < na[it] + istart[it]; ia++) + { + //take the difference of the rotated and the original coordinates + diff.x = this->check_diff( pos[ia*3+0], rotpos[ia*3+0]); + diff.y = this->check_diff( pos[ia*3+1], rotpos[ia*3+1]); + diff.z = this->check_diff( pos[ia*3+2], rotpos[ia*3+2]); + //only if all "diff" are zero vectors, flag will remain "1" + if ( no_diff == false|| + !equal(diff.x,0.0)|| + !equal(diff.y,0.0)|| + !equal(diff.z,0.0) + ) + { + no_diff = false; + } + } + } + + //the current test is successful + if (no_diff == true) + { + s_flag = true; + //save the detected translation std::vector temporarily + trans.x = gtrans.x; + trans.y = gtrans.y; + trans.z = gtrans.z; + } + + //restore the original rotated coordinates by subtracting "gtrans" + for (int it = 0; it < ntype; it++) + { + for (int ia = istart[it]; ia < na[it] + istart[it]; ia++) + { + rotpos[ia*3+0] -= gtrans.x; + rotpos[ia*3+1] -= gtrans.y; + rotpos[ia*3+2] -= gtrans.z; + } + } + } + + if (s_flag == 1) + { + gtrans.x = trans.x; + gtrans.y = trans.y; + gtrans.z = trans.z; + } + return s_flag; +} + diff --git a/source/module_cell/module_symmetry/symm_getgroup.cpp b/source/module_cell/module_symmetry/symm_getgroup.cpp new file mode 100644 index 00000000000..977566af9d7 --- /dev/null +++ b/source/module_cell/module_symmetry/symm_getgroup.cpp @@ -0,0 +1,105 @@ +void Symmetry::getgroup(int& nrot, int& nrotk, std::ofstream& ofs_running, + const int& nop, const ModuleBase::Matrix3* symop, ModuleBase::Matrix3* gmatrix, + ModuleBase::Vector3* gtrans, double* pos, double* rotpos, + int* index, const int ntype, const int itmin_type, + const int itmin_start, int* istart, int* na)const +{ + ModuleBase::TITLE("Symmetry", "getgroup"); + + //-------------------------------------------------------------------------------- + //return all possible space group operators that reproduce a lattice with basis + //out of a (maximum) pool of point group operations that is compatible with + //the symmetry of the pure translation lattice without any basic. + //-------------------------------------------------------------------------------- + + ModuleBase::Matrix3 zero(0,0,0,0,0,0,0,0,0); + ModuleBase::Matrix3 help[48]; + ModuleBase::Vector3 temp[48]; + + nrot = 0; + nrotk = 0; + + //------------------------------------------------------------------------- + //pass through the pool of (possibly allowed) symmetry operations and + //check each operation whether it can reproduce the lattice with basis + //------------------------------------------------------------------------- + //std::cout << "nop = " <checksym(symop[i], gtrans[i], pos, rotpos, index, ntype, itmin_type, itmin_start, istart, na); + if (s_flag == 1) + { + //------------------------------ + // this is a symmetry operation + // with no translation vectors + // so ,this is pure point group + // operations + //------------------------------ + if ( equal(gtrans[i].x,0.0) && + equal(gtrans[i].y,0.0) && + equal(gtrans[i].z,0.0)) + { + ++nrot; + gmatrix[nrot - 1] = symop[i]; + gtrans[nrot - 1].x = 0; + gtrans[nrot - 1].y = 0; + gtrans[nrot - 1].z = 0; + } + //------------------------------ + // this is a symmetry operation + // with translation vectors + // so ,this is space group + // operations + //------------------------------ + else + { + ++nrotk; + help[nrotk - 1] = symop[i]; + temp[nrotk - 1].x = gtrans[i].x; + temp[nrotk - 1].y = gtrans[i].y; + temp[nrotk - 1].z = gtrans[i].z; + } + } + } + + //----------------------------------------------------- + //If there are operations with nontrivial translations + //then store them together in the momory + //----------------------------------------------------- + if (nrotk > 0) + { + for (int i = 0; i < nrotk; ++i) + { + gmatrix[nrot + i] = help[i]; + gtrans[nrot + i].x = temp[i].x; + gtrans[nrot + i].y = temp[i].y; + gtrans[nrot + i].z = temp[i].z; + } + } + + //----------------------------------------------------- + //total number of space group operations + //----------------------------------------------------- + nrotk += nrot; + + if(test_brav) + { + ModuleBase::GlobalFunc::OUT(ofs_running,"PURE POINT GROUP OPERATIONS",nrot); + ModuleBase::GlobalFunc::OUT(ofs_running,"SPACE GROUP OPERATIONS",nrotk); + } + + //----------------------------------------------------- + //fill the rest of matrices and vectors with zeros + //----------------------------------------------------- + if (nrotk < 48) + { + for (int i = nrotk; i < 48; ++i) + { + gmatrix[i] = zero; + gtrans[i].x = 0; + gtrans[i].y = 0; + gtrans[i].z = 0; + } + } + return; +} diff --git a/source/module_cell/module_symmetry/symm_lattice.cpp b/source/module_cell/module_symmetry/symm_lattice.cpp new file mode 100644 index 00000000000..bef6387b49d --- /dev/null +++ b/source/module_cell/module_symmetry/symm_lattice.cpp @@ -0,0 +1,451 @@ +//--------------------------------------------------- +// The lattice will be transformed to a 'standard +// cystallographic setting', the relation between +// 'origin' and 'transformed' lattice vectors will +// be givin in matrix form +//--------------------------------------------------- +int Symmetry::standard_lat( + ModuleBase::Vector3 &a, + ModuleBase::Vector3 &b, + ModuleBase::Vector3 &c, + double *cel_const) const +{ + static bool first = true; + // there are only 14 types of Bravais lattice. + int type = 15; + //---------------------------------------------------- + // used to calculte the volume to judge whether + // the lattice vectors corrispond the right-hand-sense + //---------------------------------------------------- + double volume = 0; + //the lattice vectors have not been changed + + const double aa = a * a; + const double bb = b * b; + const double cc = c * c; + const double ab = a * b; //std::vector: a * b * cos(alpha) + const double bc = b * c; //std::vector: b * c * cos(beta) + const double ca = c * a; //std::vector: c * a * cos(gamma) + double norm_a = a.norm(); + double norm_b = b.norm(); + double norm_c = c.norm(); + double gamma = ab /( norm_a * norm_b ); // cos(gamma) + double alpha = bc /( norm_b * norm_c ); // cos(alpha) + double beta = ca /( norm_a * norm_c ); // cos(beta) + double amb = sqrt( aa + bb - 2 * ab ); //amb = |a - b| + double bmc = sqrt( bb + cc - 2 * bc ); + double cma = sqrt( cc + aa - 2 * ca ); + double apb = sqrt( aa + bb + 2 * ab ); //amb = |a + b| + double bpc = sqrt( bb + cc + 2 * bc ); + double cpa = sqrt( cc + aa + 2 * ca ); + double apbmc = sqrt( aa + bb + cc + 2 * ab - 2 * bc - 2 * ca ); //apbmc = |a + b - c| + double bpcma = sqrt( bb + cc + aa + 2 * bc - 2 * ca - 2 * ab ); + double cpamb = sqrt( cc + aa + bb + 2 * ca - 2 * ab - 2 * bc ); + double abc = ab + bc + ca; + + if (first) + { + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"NORM_A",norm_a); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"NORM_B",norm_b); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"NORM_C",norm_c); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"ALPHA (DEGREE)", acos(alpha)/ModuleBase::PI*180.0 ); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"BETA (DEGREE)" ,acos(beta)/ModuleBase::PI*180.0 ); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"GAMMA (DEGREE)" ,acos(gamma)/ModuleBase::PI*180.0 ); + first = false; + } + + Symm_Other::right_hand_sense(a, b, c); + ModuleBase::GlobalFunc::ZEROS(cel_const, 6); + const double small = PARAM.inp.symmetry_prec; + + //--------------------------- + // 1. alpha == beta == gamma + //--------------------------- + if( equal(alpha, gamma) && equal(alpha, beta) ) + { + //-------------- + // a == b == c + //-------------- + if( equal(norm_a, norm_b) && equal(norm_b, norm_c)) + { + //--------------------------------------- + // alpha == beta == gamma == 90 degree + //--------------------------------------- + if ( equal(alpha,0.0) ) + { + type=1; + cel_const[0]=norm_a; + } + //---------------------------------------- + // cos(alpha) = -1.0/3.0 + //---------------------------------------- + else if( equal(alpha, -1.0/3.0) ) + { + type=2; + cel_const[0]=norm_a*2.0/sqrt(3.0); + } + //---------------------------------------- + // cos(alpha) = 0.5 + //---------------------------------------- + else if( equal(alpha, 0.5) ) + { + type=3; + cel_const[0]=norm_a*sqrt(2.0); + } + //---------------------------------------- + // cos(alpha) = all the others + //---------------------------------------- + else + { + type=7; + cel_const[0]=norm_a; + cel_const[3]=alpha; + } + } + // Crystal classes with inequal length of lattice vectors but also with + // A1*A2=A1*A3=A2*A3: + // Orthogonal axes: + else if(equal(gamma,0.0)) + { + // Two axes with equal lengths means simple tetragonal: (IBRAV=5) + // Adjustment: 'c-axis' shall be the special axis. + if (equal(norm_a, norm_b)) + { + type=5; + cel_const[0]=norm_a; + cel_const[2]=norm_c/norm_a; + // No axes with equal lengths means simple orthorhombic (IBRAV=8): + // Adjustment: Sort the axis by increasing lengths: + } + else if(((norm_c-norm_b)>small) && ((norm_b-norm_a)>small) ) + { + type=8; + cel_const[0]=norm_a; + cel_const[1]=norm_b/norm_a; + cel_const[2]=norm_c/norm_a; + } + // Crystal classes with A1*A3=A2*A3=/A1*A2: + } + }//end alpha=beta=gamma + //----------------------- + // TWO EQUAL ANGLES + // alpha == beta != gamma (gamma is special) + //------------------------ + else if (equal(alpha-beta, 0)) + { + //--------------------------------------------------------- + // alpha = beta = 90 degree + // One axis orthogonal with respect to the other two axes: + //--------------------------------------------------------- + if (equal(alpha, 0.0)) + { + //----------------------------------------------- + // a == b + // Equal length of the two nonorthogonal axes: + //----------------------------------------------- + if (equal(norm_a, norm_b)) + { + // Cosine(alpha) equal to -1/2 means hexagonal: (IBRAV=4) + // Adjustment: 'c-axis' shall be the special axis. + if ( equal(gamma, -0.5)) //gamma = 120 degree + { + type=4; + cel_const[0]=norm_a; + cel_const[2]=norm_c/norm_a; + // Other angles mean base-centered orthorhombic: (IBRAV=11) + // Adjustment: Cosine between A1 and A2 shall be lower than zero, the + // 'c-axis' shall be the special axis. + } + else if(gamma<(-1.0*small)) //gamma > 90 degree + { + type=11; + cel_const[0]=apb; + cel_const[1]=amb/apb; + cel_const[2]=norm_c/apb; + cel_const[5]=gamma; + } + // Different length of the two axes means simple monoclinic (IBRAV=12): + // Adjustment: Cosine(gamma) should be lower than zero, special axis + // shall be the 'b-axis'(!!!) and |A1|<|A3|: + } + //---------- + // a!=b!=c + //---------- + else if( gamma<(-1.0*small) && (norm_a-norm_b)>small) + { + type=12; + cel_const[0]=norm_b; + cel_const[1]=norm_c/norm_b; + cel_const[2]=norm_a/norm_b; + cel_const[4]=gamma; + //adjust: a->c, b->a, c->b + ModuleBase::Vector3 tmp=c; + c=a; + a=b; + b=tmp; + } + }//end gamma) + cel_const[4]=(a+b)*c/apb/norm_c; + } + } + } //end alpha==beta + //------------------------------- + // three angles are not equal + //------------------------------- + else + { + // Crystal classes with A1*A2=/A1*A3=/A2*A3 + // |A1|=|A2|=|A3| means body-centered orthorhombic (IBRAV=9): + // Further additional criterions are: (A1+A2), (A1+A3) and (A2+A3) are + // orthogonal to one another and (adjustment//): |A1+A2|>|A1+A3|>|A2+A3| + if (equal(norm_a, norm_b) && + equal(norm_b, norm_c) && + ((cpa-bpc)>small) && + ((apb-cpa)>small) && + equal(norm_c*norm_c+abc, 0)) + { + type=9; + cel_const[0]=bpc; + cel_const[1]=cpa/bpc; + cel_const[2]=apb/bpc; + } + // |A1|=|A2-A3| and |A2|=|A1-A3| and |A3|=|A1-A2| means face-centered + // orthorhombic (IBRAV=10): + // Adjustment: |A1+A2-A3|>|A1+A3-A2|>|A2+A3-A1| + else if(equal(amb, norm_c) && + equal(cma, norm_b) && + equal(bmc, norm_a) && + ((apbmc-cpamb)>small) && + ((cpamb-bpcma)>small)) + { + type=10; + cel_const[0]=bpcma; + cel_const[1]=cpamb/bpcma; + cel_const[2]=apbmc/bpcma; + } + // Now there exists only one further possibility - triclinic (IBRAV=14): + // Adjustment: All three cosines shall be greater than zero and ordered: + else if((gamma>beta) && (beta>alpha) && (alpha>small)) + { + type=14; + cel_const[0]=norm_a; + cel_const[1]=norm_b/norm_a; + cel_const[2]=norm_c/norm_a; + cel_const[3]=alpha; + cel_const[4]=beta; + cel_const[5]=gamma; + } + } + + return type; +} + +//--------------------------------------------------- +// The lattice will be transformed to a 'standard +// cystallographic setting', the relation between +// 'origin' and 'transformed' lattice vectors will +// be givin in matrix form +// must be called before symmetry analysis +// only need to called once for each ion step +//--------------------------------------------------- +void Symmetry::lattice_type( + ModuleBase::Vector3 &v1, + ModuleBase::Vector3 &v2, + ModuleBase::Vector3 &v3, + ModuleBase::Vector3 &v01, + ModuleBase::Vector3 &v02, + ModuleBase::Vector3 &v03, + double *cel_const, + double *pre_const, + int& real_brav, + std::string& bravname, + const Atom* atoms, + bool convert_atoms, + double* newpos)const +{ + ModuleBase::TITLE("Symmetry","lattice_type"); + + //---------------------------------------------- + // (1) adjustement of the basis to right hand + // sense by inversion of all three lattice + // vectors if necessary + //---------------------------------------------- + const bool right = Symm_Other::right_hand_sense(v1, v2, v3); + + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"right hand lattice",right); + + //------------------------------------------------- + // (2) save and copy the original lattice vectors. + //------------------------------------------------- + v01 = v1; + v02 = v2; + v03 = v3; + + //-------------------------------------------- + // (3) calculate the 'pre_const' + //-------------------------------------------- + ModuleBase::GlobalFunc::ZEROS(pre_const, 6); + + int pre_brav = standard_lat(v1, v2, v3, cel_const); + + for ( int i = 0; i < 6; ++i) + { + pre_const[i] = cel_const[i]; + } + + // find the shortest basis vectors of the lattice + this->get_shortest_latvec(v1, v2, v3); + + Symm_Other::right_hand_sense(v1, v2, v3); + + real_brav = 15; + double temp_const[6]; + + //then we should find the best lattice vectors to make much easier the determination of the lattice symmetry + //the method is to contrast the combination of the shortest vectors and determine their symmmetry + + ModuleBase::Vector3 w1, w2, w3; + ModuleBase::Vector3 q1, q2, q3; + this->get_optlat(v1, v2, v3, w1, w2, w3, real_brav, cel_const, temp_const); + + //now, the highest symmetry of the combination of the shortest vectors has been found + //then we compare it with the original symmetry + + bool change_flag=false; + for (int i = 0; i < 6; ++i) { + if(!equal(cel_const[i], pre_const[i])) + {change_flag=true; break; + } + } + + if ( real_brav < pre_brav || change_flag ) + { + //if the symmetry of the new vectors is higher, store the new ones + for (int i = 0; i < 6; ++i) + { + cel_const[i] = temp_const[i]; + } + q1 = w1; + q2 = w2; + q3 = w3; + if(convert_atoms) + { + GlobalV::ofs_running <ntype; ++it) + { + for (int ia = 0; ia < this->na[it]; ++ia) + { + ModuleBase::Mathzone::Cartesian_to_Direct(atoms[it].tau[ia].x, + atoms[it].tau[ia].y, + atoms[it].tau[ia].z, + q1.x, q1.y, q1.z, + q2.x, q2.y, q2.z, + q3.x, q3.y, q3.z, + newpos[3*at],newpos[3*at+1],newpos[3*at+2]); + + for(int k=0; k<3; ++k) + { + this->check_translation( newpos[at*3+k], -floor(newpos[at*3+k])); + this->check_boundary( newpos[at*3+k] ); + } + ++at; + } + } + } + // return the optimized lattice in v1, v2, v3 + v1=q1; + v2=q2; + v3=q3; + } + else + { + //else, store the original ones + for (int i = 0; i < 6; ++i) + { + cel_const[i] = pre_const[i]; + } + //newpos also need to be set + if(convert_atoms) + { + int at=0; + for (int it = 0; it < this->ntype; ++it) + { + for (int ia = 0; ia < this->na[it]; ++ia) + { + ModuleBase::Mathzone::Cartesian_to_Direct(atoms[it].tau[ia].x, + atoms[it].tau[ia].y, + atoms[it].tau[ia].z, + v1.x, v1.y, v1.z, + v2.x, v2.y, v2.z, + v3.x, v3.y, v3.z, + newpos[3*at],newpos[3*at+1],newpos[3*at+2]); + for(int k=0; k<3; ++k) + { + this->check_translation( newpos[at*3+k], -floor(newpos[at*3+k])); + this->check_boundary( newpos[at*3+k] ); + } + ++at; + } + } + } + } + + /* + bool flag3; + if (pre_brav == temp_brav) + { + flag3 = 0; + if (!equal(temp_const[0], pre_const[0]) || + !equal(temp_const[1], pre_const[1]) || + !equal(temp_const[2], pre_const[2]) || + !equal(temp_const[3], pre_const[3]) || + !equal(temp_const[4], pre_const[4]) || + !equal(temp_const[5], pre_const[5]) + ) + { + flag3 = 1; + } + if (flag3==0) { + //particularly, if the symmetry of origin and new are exactly the same, we choose the original ones + //Hey! the original vectors have been changed!!! + v1 = s1; + v2 = s2; + v3 = s3; + change=0; + GlobalV::ofs_running<<" The lattice vectors have been set back!"<> mag_type_atoms; + for (int it = 0;it < ntype;++it) + { + for (int ia = 0; ia < atoms[it].na; ++ia) + { + bool find = false; + for (auto& mt : mag_type_atoms) + { + const int mag_iat = *mt.begin(); + const int mag_it = st.iat2it[mag_iat]; + const int mag_ia = st.iat2ia[mag_iat]; + if (it == mag_it && this->equal(atoms[it].mag[ia], atoms[mag_it].mag[mag_ia])) + { + mt.insert(st.itia2iat(it, ia)); + find = true; + break; + } + } + if (!find) + { + mag_type_atoms.push_back(std::set({ st.itia2iat(it,ia) })); + } + } + } + + // 2. get the start index, number of atoms and positions for each mag_type + std::vector mag_istart(mag_type_atoms.size()); + std::vector mag_na(mag_type_atoms.size()); + std::vector mag_pos; + int mag_itmin_type = 0; + int mag_itmin_start = 0; + for (int mag_it = 0;mag_it < mag_type_atoms.size(); ++mag_it) + { + mag_na[mag_it] = mag_type_atoms.at(mag_it).size(); + if (mag_it > 0) + { + mag_istart[mag_it] = mag_istart[mag_it - 1] + mag_na[mag_it - 1]; + } + if (mag_na[mag_it] < mag_na[itmin_type]) + { + mag_itmin_type = mag_it; + mag_itmin_start = mag_istart[mag_it]; + } + for (auto& mag_iat : mag_type_atoms.at(mag_it)) + { + // this->newpos have been ordered by original structure(ntype, na), it cannot be directly used here. + // we need to reset the calculate again the coordinate of the new structure. + const ModuleBase::Vector3 direct_tmp = atoms[st.iat2it[mag_iat]].tau[st.iat2ia[mag_iat]] * this->optlat.Inverse(); + std::array direct = { direct_tmp.x, direct_tmp.y, direct_tmp.z }; + for (int i = 0; i < 3; ++i) + { + this->check_translation(direct[i], -floor(direct[i])); + this->check_boundary(direct[i]); + mag_pos.push_back(direct[i]); + } + } + } + + // 3. analyze the effective structure + this->getgroup(nrot_out, nrotk_out, GlobalV::ofs_running, + this->nop, this->symop, this->gmatrix, + this->gtrans, mag_pos.data(), this->rotpos, + this->index, mag_type_atoms.size(), mag_itmin_type, + mag_itmin_start, mag_istart.data(), mag_na.data()); + +} + +bool Symmetry::magmom_same_check(const Atom* atoms)const +{ + ModuleBase::TITLE("Symmetry", "magmom_same_check"); + bool pricell_loop = true; + for (int it = 0;it < ntype;++it) + { + if (pricell_loop) { + for (int ia = 1;ia < atoms[it].na;++ia) + { + if (!equal(atoms[it].m_loc_[ia].x, atoms[it].m_loc_[0].x) || + !equal(atoms[it].m_loc_[ia].y, atoms[it].m_loc_[0].y) || + !equal(atoms[it].m_loc_[ia].z, atoms[it].m_loc_[0].z)) + { + pricell_loop = false; + break; + } + } + } + } + return pricell_loop; +} + diff --git a/source/module_cell/module_symmetry/symm_other.h b/source/module_cell/module_symmetry/symm_other.h index 5e571a29a75..8107d7f6320 100644 --- a/source/module_cell/module_symmetry/symm_other.h +++ b/source/module_cell/module_symmetry/symm_other.h @@ -8,7 +8,9 @@ namespace ModuleSymmetry namespace Symm_Other { void print1(const int &ibrav, const double *cel_const, std::ofstream &ofs_running); + bool right_hand_sense(ModuleBase::Vector3 &v1,ModuleBase::Vector3 &v2,ModuleBase::Vector3 &v3); + double celvol(const ModuleBase::Vector3 &a, const ModuleBase::Vector3 &b, const ModuleBase::Vector3 &c); diff --git a/source/module_cell/module_symmetry/symm_pricell.cpp b/source/module_cell/module_symmetry/symm_pricell.cpp new file mode 100644 index 00000000000..14f8f8d56cd --- /dev/null +++ b/source/module_cell/module_symmetry/symm_pricell.cpp @@ -0,0 +1,319 @@ +void Symmetry::pricell(double* pos, const Atom* atoms) +{ + bool no_diff = false; + ptrans.clear(); + + for (int it = 0; it < ntype; it++) + { + //------------------------------------ + // impose periodic boundary condition + // 0.5 -> -0.5 + //------------------------------------ + for (int j = istart[it]; j < istart[it] + na[it]; ++j) + { + this->check_boundary(pos[j*3+0]); + this->check_boundary(pos[j*3+1]); + this->check_boundary(pos[j*3+2]); + } + + //order original atomic positions for current species + this->atom_ordering_new(pos + istart[it] * 3, na[it], index + istart[it]); + //copy pos to rotpos + for (int j = istart[it]; j < istart[it] + na[it]; ++j) + { + const int xx=j*3; + const int yy=j*3+1; + const int zz=j*3+2; + rotpos[xx] = pos[xx]; + rotpos[yy] = pos[yy]; + rotpos[zz] = pos[zz]; + } + } + + ModuleBase::Vector3 diff; + double tmp_ptrans[3]; + + //--------------------------------------------------------- + // itmin_start = the start atom positions of species itmin + //--------------------------------------------------------- + // (s)tart (p)osition of atom (t)ype which has (min)inal number. + ModuleBase::Vector3 sptmin(pos[itmin_start * 3], pos[itmin_start * 3 + 1], pos[itmin_start * 3 + 2]); + + for (int i = itmin_start; i < itmin_start + na[itmin_type]; ++i) + { + //set up the current test std::vector "gtrans" + //and "gtrans" could possibly contain trivial translations: + tmp_ptrans[0] = this->get_translation_vector( pos[i*3+0], sptmin.x); + tmp_ptrans[1] = this->get_translation_vector( pos[i*3+1], sptmin.y); + tmp_ptrans[2] = this->get_translation_vector( pos[i*3+2], sptmin.z); + //translate all the atomic coordinates by "gtrans" + for (int it = 0; it < ntype; it++) + { + for (int ia = istart[it]; ia < na[it] + istart[it]; ia++) + { + this->check_translation( rotpos[ia*3+0], tmp_ptrans[0] ); + this->check_translation( rotpos[ia*3+1], tmp_ptrans[1] ); + this->check_translation( rotpos[ia*3+2], tmp_ptrans[2] ); + + this->check_boundary( rotpos[ia*3+0] ); + this->check_boundary( rotpos[ia*3+1] ); + this->check_boundary( rotpos[ia*3+2] ); + } + //order translated atomic positions for current species + this->atom_ordering_new(rotpos + istart[it] * 3, na[it], index + istart[it]); + } + + no_diff = true; + //compare the two lattices 'one-by-one' whether they are identical + for (int it = 0; it < ntype; it++) + { + for (int ia = istart[it]; ia < na[it] + istart[it]; ia++) + { + //take the difference of the rotated and the original coordinates + diff.x = this->check_diff( pos[ia*3+0], rotpos[ia*3+0]); + diff.y = this->check_diff( pos[ia*3+1], rotpos[ia*3+1]); + diff.z = this->check_diff( pos[ia*3+2], rotpos[ia*3+2]); + //only if all "diff" are zero vectors, flag will remain "1" + if (!equal(diff.x,0.0)|| + !equal(diff.y,0.0)|| + !equal(diff.z,0.0)) + { + no_diff = false; + break; + } + } + if (!no_diff) { + break; + } + } + + //the current test is successful + if (no_diff) { + ptrans.push_back(ModuleBase::Vector3(tmp_ptrans[0], + tmp_ptrans[1], + tmp_ptrans[2])); + } + //restore the original rotated coordinates by subtracting "ptrans" + for (int it = 0; it < ntype; it++) + { + for (int ia = istart[it]; ia < na[it] + istart[it]; ia++) + { + rotpos[ia*3+0] -= tmp_ptrans[0]; + rotpos[ia*3+1] -= tmp_ptrans[1]; + rotpos[ia*3+2] -= tmp_ptrans[2]; + } + } + } + int ntrans=ptrans.size(); + if (ntrans <= 1) + { + GlobalV::ofs_running<<" Original cell was already a primitive cell."<p1=this->a1; + this->p2=this->a2; + this->p3=this->a3; + this->pbrav=this->real_brav; + this->ncell=1; + for (int i = 0; i < 6; ++i) { + this->pcel_const[i] = this->cel_const[i]; + } + return; + } + + //sort ptrans: + double* ptrans_array = new double[ntrans*3]; + for(int i=0;iatom_ordering_new(ptrans_array, ntrans, index); + // std::cout<<"final ptrans:"< b1, b2, b3; + int iplane=0, jplane=0, kplane=0; + //1. kplane for b3 + while (kplane < ntrans + && std::abs(ptrans[kplane].z - ptrans[0].z) < this->epsilon) { + ++kplane; + } + if (kplane == ntrans) { + kplane = 0; // a3-direction have no smaller pricell + } + b3=kplane>0 ? + ModuleBase::Vector3(ptrans[kplane].x, ptrans[kplane].y, ptrans[kplane].z) : + ModuleBase::Vector3(0, 0, 1); + //2. jplane for b2 (not collinear with b3) + jplane=kplane+1; + while (jplane < ntrans + && (std::abs(ptrans[jplane].y - ptrans[0].y) < this->epsilon + || equal((ptrans[jplane] ^ b3).norm(), 0))) { + ++jplane; + } + if (jplane == ntrans) { + jplane = kplane; // a2-direction have no smaller pricell + } + b2=jplane>kplane ? + ModuleBase::Vector3(ptrans[jplane].x, ptrans[jplane].y, ptrans[jplane].z) : + ModuleBase::Vector3(0, 1, 0); + //3. iplane for b1 (not coplane with ) + iplane=jplane+1; + while (iplane < ntrans + && (std::abs(ptrans[iplane].x - ptrans[0].x) < this->epsilon + || equal(ptrans[iplane] * (b2 ^ b3), 0))) { + ++iplane; + } + b1=(iplane>jplane && iplane(ptrans[iplane].x, ptrans[iplane].y, ptrans[iplane].z) : + ModuleBase::Vector3(1, 0, 0); //a1-direction have no smaller pricell + + + ModuleBase::Matrix3 coeff(b1.x, b1.y, b1.z, b2.x, b2.y, b2.z, b3.x, b3.y, b3.z); + this->plat=coeff*this->optlat; + + //deal with collineation caused by default b1, b2, b3 + if(equal(plat.Det(), 0)) + { + if(kplane==0) //try a new b3 + { + std::cout<<"try a new b3"<kplane) // use default b2 + { + coeff.e31=0; + coeff.e32=1; + coeff.e33=0; + } + else //use default b1 + { + coeff.e31=1; + coeff.e32=0; + coeff.e33=0; + } + } + else if(jplane<=kplane) + { + coeff.e21=0; + coeff.e22=0; + coeff.e23=1; + } + else + { + coeff.e11=0; + coeff.e12=0; + coeff.e13=1; + } + this->plat=coeff*this->optlat; + assert(!equal(plat.Det(), 0)); + } + + this->p1.x=plat.e11; + this->p1.y=plat.e12; + this->p1.z=plat.e13; + this->p2.x=plat.e21; + this->p2.y=plat.e22; + this->p2.z=plat.e23; + this->p3.x=plat.e31; + this->p3.y=plat.e32; + this->p3.z=plat.e33; + +#ifdef __DEBUG + GlobalV::ofs_running<<"lattice vectors of primitive cell (initial):"< p01=p1, p02=p2, p03=p3; + double pcel_pre_const[6]; + for (int i = 0; i < 6; ++i) { + pcel_pre_const[i] = pcel_const[i]; + } + this->lattice_type(p1, p2, p3, p01, p02, p03, pcel_const, pcel_pre_const, pbrav, pbravname, atoms, false, nullptr); + + this->plat.e11=p1.x; + this->plat.e12=p1.y; + this->plat.e13=p1.z; + this->plat.e21=p2.x; + this->plat.e22=p2.y; + this->plat.e23=p2.z; + this->plat.e31=p3.x; + this->plat.e32=p3.y; + this->plat.e33=p3.z; + +#ifdef __DEBUG + GlobalV::ofs_running<<"lattice vectors of primitive cell (optimized):"<pbrav, this->pcel_const, GlobalV::ofs_running); + + //count the number of pricells + GlobalV::ofs_running<<"optimized lattice volume: "<optlat.Det()<plat.Det()<optlat.Det()/this->plat.Det()); + this->ncell=floor(ncell_double+0.5); + + auto reset_pcell = [this]() -> void { + std::cout << " Now regard the structure as a primitive cell." << std::endl; + this->ncell = 1; + this->ptrans = std::vector >(1, ModuleBase::Vector3(0, 0, 0)); + GlobalV::ofs_running << "WARNING: Original cell may have more than one primitive cells, \ + but we have to treat it as a primitive cell. Use a larger `symmetry_prec`to avoid this warning." << std::endl; + }; + if (this->ncell != ntrans) + { + std::cout << " WARNING: PRICELL: NCELL != NTRANS !" << std::endl; + std::cout << " NCELL=" << ncell << ", NTRANS=" << ntrans << std::endl; + std::cout << " Suggest solution: Use a larger `symmetry_prec`. " << std::endl; + reset_pcell(); + return; + } + if(std::abs(ncell_double-double(this->ncell)) > this->epsilon*100) + { + std::cout << " WARNING: THE NUMBER OF PRIMITIVE CELL IS NOT AN INTEGER !" << std::endl; + std::cout << " NCELL(double)=" << ncell_double << ", NTRANS=" << ncell << std::endl; + std::cout << " Suggest solution: Use a larger `symmetry_prec`. " << std::endl; + reset_pcell(); + return; + } + GlobalV::ofs_running<<"Original cell was built up by "<ncell<<" primitive cells."<gtrans_convert(ptrans.data(), ptrans.data(), ntrans, this->optlat, inputlat ); + + //how many pcell in supercell + int n1=0; + int n2=0; + int n3=0; + + ModuleBase::Matrix3 nummat0=this->optlat*this->plat.Inverse(); + ModuleBase::Matrix3 nummat, transmat; + hermite_normal_form(nummat0, nummat, transmat); + n1=floor (nummat.e11 + epsilon); + n2=floor (nummat.e22 + epsilon); + n3=floor (nummat.e33 + epsilon); + if(n1*n2*n3 != this->ncell) + { + std::cout << " WARNING: Number of cells and number of vectors did not agree."; + std::cout<<"Try to change symmetry_prec in INPUT." << std::endl; + reset_pcell(); + } + return; +} + diff --git a/source/module_cell/module_symmetry/symm_rho.cpp b/source/module_cell/module_symmetry/symm_rho.cpp new file mode 100644 index 00000000000..fb7ea296e35 --- /dev/null +++ b/source/module_cell/module_symmetry/symm_rho.cpp @@ -0,0 +1,299 @@ +void Symmetry::rho_symmetry( double *rho, + const int &nr1, const int &nr2, const int &nr3) +{ + ModuleBase::timer::tick("Symmetry","rho_symmetry"); + + // allocate flag for each FFT grid. + bool* symflag = new bool[nr1 * nr2 * nr3]; + for (int i=0; i0 ); + assert(nrotk <=48 ); + int *ri = new int[nrotk]; + int *rj = new int[nrotk]; + int *rk = new int[nrotk]; + + int ci = 0; + for (int i = 0; i< nr1; ++i) + { + for (int j = 0; j< nr2; ++j) + { + for (int k = 0; k< nr3; ++k) + { + if (!symflag[i * nr2 * nr3 + j * nr3 + k]) + { + double sum = 0; + + for (int isym = 0; isym < nrotk; ++isym) + { + this->rotate(gmatrix[isym], gtrans[isym], i, j, k, nr1, nr2, nr3, ri[isym], rj[isym], rk[isym]); + const int index = ri[isym] * nr2 * nr3 + rj[isym] * nr3 + rk[isym]; + sum += rho[ index ]; + } + sum /= nrotk; + + for (int isym = 0; isym < nrotk; ++isym) + { + const int index = ri[isym] * nr2 * nr3 + rj[isym] * nr3 + rk[isym]; + rho[index] = sum; + symflag[index] = true; + } + } + } + } + } + + delete[] symflag; + delete[] ri; + delete[] rj; + delete[] rk; + ModuleBase::timer::tick("Symmetry","rho_symmetry"); +} + +void Symmetry::rhog_symmetry(std::complex *rhogtot, + int* ixyz2ipw, const int &nx, const int &ny, const int &nz, + const int &fftnx, const int &fftny, const int &fftnz) +{ + ModuleBase::timer::tick("Symmetry","rhog_symmetry"); + // ---------------------------------------------------------------------- + // the current way is to cluster the FFT grid points into groups in advance. + // and use OpenMP to realize parallel calculation, one thread works in one group. + // ---------------------------------------------------------------------- + + const int nxyz = fftnx*fftny*fftnz; + assert(nxyz>0); + + // allocate flag for each FFT grid. + // which group the grid belongs to + int* symflag = new int[nxyz]; + + // which rotration operation the grid corresponds to + int(*isymflag)[48] = new int[nxyz][48]; + + // group information + int(*table_xyz)[48] = new int[nxyz][48]; + + // how many symmetry operations have been covered + int* count_xyz = new int[nxyz]; + + for (int i = 0; i < nxyz; i++) + { + symflag[i] = -1; + } + int group_index = 0; + + assert(nrotk >0 ); + assert(nrotk <=48 ); + + //map the gmatrix to inv + std::vectorinvmap(this->nrotk, -1); + this->gmatrix_invmap(kgmatrix, nrotk, invmap.data()); + + // ------------------------------------------------------------------------ + // This code defines a lambda function called "rotate_recip" that takes + // a 3x3 matrix and a 3D vector as input. It performs a rotation operation + // on the vector using the matrix and returns the rotated vector. + // Specifically, it calculates the new coordinates of the vector after + // the rotation and applies periodic boundary conditions to ensure that + // the coordinates are within the FFT-grid dimensions. + // The rotated vector is returned by modifying the input vector. + // ------------------------------------------------------------------------ + //rotate function (different from real space, without scaling gmatrix) + auto rotate_recip = [&] (ModuleBase::Matrix3& g, ModuleBase::Vector3& g0, int& ii, int& jj, int& kk) + { + ii = int(g.e11 * g0.x + g.e21 * g0.y + g.e31 * g0.z) ; + if (ii < 0) + { + ii += 10 * nx; + } + ii = ii%nx; + jj = int(g.e12 * g0.x + g.e22 * g0.y + g.e32 * g0.z) ; + if (jj < 0) + { + jj += 10 * ny; + } + jj = jj%ny; + kk = int(g.e13 * g0.x + g.e23 * g0.y + g.e33 * g0.z); + if (kk < 0) + { + kk += 10 * nz; + } + kk = kk%nz; + return; + }; + + // ------------------------------------------------------------------------ + // Trying to group fft grids first. + // It iterates over each FFT-grid point and checks if it is within the + // PW-sphere. If it is, put all the FFT-grid points connected by the + // rotation operation into one group( the index is stored in int(*table_xyz)). + // The code marks the point as processed to avoid redundant calculations + // by using int* symflag. + // ------------------------------------------------------------------------ + + ModuleBase::timer::tick("Symmetry","group_fft_grids"); + for (int i = 0; i< fftnx; ++i) + { + //tmp variable + ModuleBase::Vector3 tmp_gdirect0(0, 0, 0); + tmp_gdirect0.x=(i>int(nx/2)+1)?(i-nx):i; + for (int j = 0; j< fftny; ++j) + { + tmp_gdirect0.y=(j>int(ny/2)+1)?(j-ny):j; + for (int k = 0; k< fftnz; ++k) + { + int ixyz0=(i*fftny+j)*fftnz+k; + if (symflag[ixyz0] == -1) + { + int ipw0=ixyz2ipw[ixyz0]; + //if a fft-grid is not in pw-sphere, just do not consider it. + if (ipw0 == -1) { + continue; + } + tmp_gdirect0.z=(k>int(nz/2)+1)?(k-nz):k; + int rot_count=0; + for (int isym = 0; isym < nrotk; ++isym) + { + if (invmap[isym] < 0 || invmap[isym] > nrotk) { continue; } + //tmp variables + int ii, jj, kk=0; + rotate_recip(kgmatrix[invmap[isym]], tmp_gdirect0, ii, jj, kk); + if(ii>=fftnx || jj>=fftny || kk>= fftnz) + { + if(!PARAM.globalv.gamma_only_pw) + { + std::cout << " ROTATE OUT OF FFT-GRID IN RHOG_SYMMETRY !" << std::endl; + ModuleBase::QUIT(); + } + // for gamma_only_pw, just do not consider this rotation. + continue; + } + int ixyz=(ii*fftny+jj)*fftnz+kk; + //fft-grid index to (ip, ig) + int ipw=ixyz2ipw[ixyz]; + if(ipw==-1) //not in pw-sphere + { + continue; //else, just skip it + } + symflag[ixyz] = group_index; + isymflag[group_index][rot_count] = invmap[isym]; + table_xyz[group_index][rot_count] = ixyz; + ++rot_count; + assert(rot_count <= nrotk); + count_xyz[group_index] = rot_count; + } + group_index++; + } + } + } + } + ModuleBase::timer::tick("Symmetry","group_fft_grids"); + + // ------------------------------------------------------------------- + // This code performs symmetry operations on the reciprocal space + // charge density using FFT-grids. It iterates over each FFT-grid + // point in a particular group, applies a phase factor and sums the + // charge density over the symmetry operations, and then divides by + // the number of symmetry operations. Finally, it updates the charge + // density for each FFT-grid point using the calculated sum. + // ------------------------------------------------------------------- + +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (int g_index = 0; g_index < group_index; g_index++) + { + // record the index and gphase but not the final gdirect for each symm-opt + int *ipw_record = new int[nrotk]; + int *ixyz_record = new int[nrotk]; + std::complex* gphase_record = new std::complex [nrotk]; + std::complex sum(0, 0); + int rot_count=0; + + for (int c_index = 0; c_index < count_xyz[g_index]; ++c_index) + { + int ixyz0 = table_xyz[g_index][c_index]; + int ipw0 = ixyz2ipw[ixyz0]; + + if (symflag[ixyz0] == g_index) + { + // note : do not use PBC after rotation. + // we need a real gdirect to get the correspoding rhogtot. + int k = ixyz0%fftnz; + int j = ((ixyz0-k)/fftnz)%fftny; + int i = ((ixyz0-k)/fftnz-j)/fftny; + + //fft-grid index to gdirect + ModuleBase::Vector3 tmp_gdirect_double(0.0, 0.0, 0.0); + tmp_gdirect_double.x=static_cast((i>int(nx/2)+1)?(i-nx):i); + tmp_gdirect_double.y=static_cast((j>int(ny/2)+1)?(j-ny):j); + tmp_gdirect_double.z=static_cast((k>int(nz/2)+1)?(k-nz):k); + + //calculate phase factor + tmp_gdirect_double = tmp_gdirect_double * ModuleBase::TWO_PI; + + double cos_arg = 0.0, sin_arg = 0.0; + double arg_gtrans = tmp_gdirect_double * gtrans[isymflag[g_index][c_index]]; + + std::complex phase_gtrans (ModuleBase::libm::cos(arg_gtrans), + ModuleBase::libm::sin(arg_gtrans)); + + // for each pricell in supercell: + for (int ipt = 0;ipt < ((ModuleSymmetry::Symmetry::pricell_loop) ? this->ncell : 1);++ipt) + { + double arg = tmp_gdirect_double * ptrans[ipt]; + double tmp_cos = 0.0, tmp_sin = 0.0; + ModuleBase::libm::sincos(arg, &tmp_sin, &tmp_cos); + cos_arg += tmp_cos; + sin_arg += tmp_sin; + } + + // add nothing to sum, so don't consider this isym into rot_count + cos_arg/=static_cast(ncell); + sin_arg/=static_cast(ncell); + + //deal with double-zero + if (equal(cos_arg, 0.0) && equal(sin_arg, 0.0)) + { + continue; + } + + std::complex gphase(cos_arg, sin_arg); + gphase = phase_gtrans * gphase; + + //deal with small difference from 1 + if (equal(gphase.real(), 1.0) && equal(gphase.imag(), 0)) + { + gphase = std::complex(1.0, 0.0); + } + + gphase_record[rot_count]=gphase; + sum += rhogtot[ipw0]*gphase; + //record + ipw_record[rot_count]=ipw0; + ixyz_record[rot_count]=ixyz0; + ++rot_count; + //assert(rot_count<=nrotk); + }//end if section + }//end c_index loop + sum /= rot_count; + for (int isym = 0; isym < rot_count; ++isym) + { + rhogtot[ipw_record[isym]] = sum/gphase_record[isym]; + } + + //Clean the records variables for each fft grid point + delete[] ipw_record; + delete[] ixyz_record; + delete[] gphase_record; + }//end g_index loop + + delete[] symflag; + delete[] isymflag; + delete[] table_xyz; + delete[] count_xyz; + ModuleBase::timer::tick("Symmetry","rhog_symmetry"); +} diff --git a/source/module_cell/module_symmetry/symmetry.cpp b/source/module_cell/module_symmetry/symmetry.cpp index e422de6e655..0f4bfb41dc7 100644 --- a/source/module_cell/module_symmetry/symmetry.cpp +++ b/source/module_cell/module_symmetry/symmetry.cpp @@ -321,1407 +321,6 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, return; } -//--------------------------------------------------- -// The lattice will be transformed to a 'standard -// cystallographic setting', the relation between -// 'origin' and 'transformed' lattice vectors will -// be givin in matrix form -//--------------------------------------------------- -int Symmetry::standard_lat( - ModuleBase::Vector3 &a, - ModuleBase::Vector3 &b, - ModuleBase::Vector3 &c, - double *cel_const) const -{ - static bool first = true; - // there are only 14 types of Bravais lattice. - int type = 15; - //---------------------------------------------------- - // used to calculte the volume to judge whether - // the lattice vectors corrispond the right-hand-sense - //---------------------------------------------------- - double volume = 0; - //the lattice vectors have not been changed - - const double aa = a * a; - const double bb = b * b; - const double cc = c * c; - const double ab = a * b; //std::vector: a * b * cos(alpha) - const double bc = b * c; //std::vector: b * c * cos(beta) - const double ca = c * a; //std::vector: c * a * cos(gamma) - double norm_a = a.norm(); - double norm_b = b.norm(); - double norm_c = c.norm(); - double gamma = ab /( norm_a * norm_b ); // cos(gamma) - double alpha = bc /( norm_b * norm_c ); // cos(alpha) - double beta = ca /( norm_a * norm_c ); // cos(beta) - double amb = sqrt( aa + bb - 2 * ab ); //amb = |a - b| - double bmc = sqrt( bb + cc - 2 * bc ); - double cma = sqrt( cc + aa - 2 * ca ); - double apb = sqrt( aa + bb + 2 * ab ); //amb = |a + b| - double bpc = sqrt( bb + cc + 2 * bc ); - double cpa = sqrt( cc + aa + 2 * ca ); - double apbmc = sqrt( aa + bb + cc + 2 * ab - 2 * bc - 2 * ca ); //apbmc = |a + b - c| - double bpcma = sqrt( bb + cc + aa + 2 * bc - 2 * ca - 2 * ab ); - double cpamb = sqrt( cc + aa + bb + 2 * ca - 2 * ab - 2 * bc ); - double abc = ab + bc + ca; - - if (first) - { - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"NORM_A",norm_a); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"NORM_B",norm_b); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"NORM_C",norm_c); -// OUT(GlobalV::ofs_running,"alpha = ", alpha ); -// OUT(GlobalV::ofs_running,"beta = " ,beta ); -// OUT(GlobalV::ofs_running,"gamma = " ,gamma ); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"ALPHA (DEGREE)", acos(alpha)/ModuleBase::PI*180.0 ); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"BETA (DEGREE)" ,acos(beta)/ModuleBase::PI*180.0 ); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"GAMMA (DEGREE)" ,acos(gamma)/ModuleBase::PI*180.0 ); - first = false; - } - - -// std::cout << " a=" << norm_a << std::endl; -// std::cout << " b=" << norm_b << std::endl; -// std::cout << " c=" << norm_c << std::endl; -// std::cout << " alpha=" << alpha << std::endl; -// std::cout << " beta=" << beta << std::endl; -// std::cout << " gamma=" << gamma << std::endl; - - - Symm_Other::right_hand_sense(a, b, c); - ModuleBase::GlobalFunc::ZEROS(cel_const, 6); - const double small = PARAM.inp.symmetry_prec; - - //--------------------------- - // 1. alpha == beta == gamma - //--------------------------- - if( equal(alpha, gamma) && equal(alpha, beta) ) - { - //-------------- - // a == b == c - //-------------- - if( equal(norm_a, norm_b) && equal(norm_b, norm_c)) - { - //--------------------------------------- - // alpha == beta == gamma == 90 degree - //--------------------------------------- - if ( equal(alpha,0.0) ) - { - type=1; - cel_const[0]=norm_a; - } - //---------------------------------------- - // cos(alpha) = -1.0/3.0 - //---------------------------------------- - else if( equal(alpha, -1.0/3.0) ) - { - type=2; - cel_const[0]=norm_a*2.0/sqrt(3.0); - } - //---------------------------------------- - // cos(alpha) = 0.5 - //---------------------------------------- - else if( equal(alpha, 0.5) ) - { - type=3; - cel_const[0]=norm_a*sqrt(2.0); - } - //---------------------------------------- - // cos(alpha) = all the others - //---------------------------------------- - else - { - type=7; - cel_const[0]=norm_a; - cel_const[3]=alpha; - } - } - // Crystal classes with inequal length of lattice vectors but also with - // A1*A2=A1*A3=A2*A3: - // Orthogonal axes: - else if(equal(gamma,0.0)) - { - // Two axes with equal lengths means simple tetragonal: (IBRAV=5) - // Adjustment: 'c-axis' shall be the special axis. - if (equal(norm_a, norm_b)) - { - type=5; - cel_const[0]=norm_a; - cel_const[2]=norm_c/norm_a; - // No axes with equal lengths means simple orthorhombic (IBRAV=8): - // Adjustment: Sort the axis by increasing lengths: - } - else if(((norm_c-norm_b)>small) && ((norm_b-norm_a)>small) ) - { - type=8; - cel_const[0]=norm_a; - cel_const[1]=norm_b/norm_a; - cel_const[2]=norm_c/norm_a; - } - // Crystal classes with A1*A3=A2*A3=/A1*A2: - } - }//end alpha=beta=gamma - //----------------------- - // TWO EQUAL ANGLES - // alpha == beta != gamma (gamma is special) - //------------------------ - else if (equal(alpha-beta, 0)) - { - //--------------------------------------------------------- - // alpha = beta = 90 degree - // One axis orthogonal with respect to the other two axes: - //--------------------------------------------------------- - if (equal(alpha, 0.0)) - { - //----------------------------------------------- - // a == b - // Equal length of the two nonorthogonal axes: - //----------------------------------------------- - if (equal(norm_a, norm_b)) - { - // Cosine(alpha) equal to -1/2 means hexagonal: (IBRAV=4) - // Adjustment: 'c-axis' shall be the special axis. - if ( equal(gamma, -0.5)) //gamma = 120 degree - { - type=4; - cel_const[0]=norm_a; - cel_const[2]=norm_c/norm_a; - // Other angles mean base-centered orthorhombic: (IBRAV=11) - // Adjustment: Cosine between A1 and A2 shall be lower than zero, the - // 'c-axis' shall be the special axis. - } - else if(gamma<(-1.0*small)) //gamma > 90 degree - { - type=11; - cel_const[0]=apb; - cel_const[1]=amb/apb; - cel_const[2]=norm_c/apb; - cel_const[5]=gamma; - } - // Different length of the two axes means simple monoclinic (IBRAV=12): - // Adjustment: Cosine(gamma) should be lower than zero, special axis - // shall be the 'b-axis'(!!!) and |A1|<|A3|: - } - //---------- - // a!=b!=c - //---------- - else if( gamma<(-1.0*small) && (norm_a-norm_b)>small) - { - type=12; - cel_const[0]=norm_b; - cel_const[1]=norm_c/norm_b; - cel_const[2]=norm_a/norm_b; - cel_const[4]=gamma; - //adjust: a->c, b->a, c->b - ModuleBase::Vector3 tmp=c; - c=a; - a=b; - b=tmp; - } - }//end gamma) - cel_const[4]=(a+b)*c/apb/norm_c; - } - } - } //end alpha==beta - //------------------------------- - // three angles are not equal - //------------------------------- - else - { - // Crystal classes with A1*A2=/A1*A3=/A2*A3 - // |A1|=|A2|=|A3| means body-centered orthorhombic (IBRAV=9): - // Further additional criterions are: (A1+A2), (A1+A3) and (A2+A3) are - // orthogonal to one another and (adjustment//): |A1+A2|>|A1+A3|>|A2+A3| - if (equal(norm_a, norm_b) && - equal(norm_b, norm_c) && - ((cpa-bpc)>small) && - ((apb-cpa)>small) && - equal(norm_c*norm_c+abc, 0)) - { - type=9; - cel_const[0]=bpc; - cel_const[1]=cpa/bpc; - cel_const[2]=apb/bpc; - } - // |A1|=|A2-A3| and |A2|=|A1-A3| and |A3|=|A1-A2| means face-centered - // orthorhombic (IBRAV=10): - // Adjustment: |A1+A2-A3|>|A1+A3-A2|>|A2+A3-A1| - else if(equal(amb, norm_c) && - equal(cma, norm_b) && - equal(bmc, norm_a) && - ((apbmc-cpamb)>small) && - ((cpamb-bpcma)>small)) - { - type=10; - cel_const[0]=bpcma; - cel_const[1]=cpamb/bpcma; - cel_const[2]=apbmc/bpcma; - } - // Now there exists only one further possibility - triclinic (IBRAV=14): - // Adjustment: All three cosines shall be greater than zero and ordered: - else if((gamma>beta) && (beta>alpha) && (alpha>small)) - { - type=14; - cel_const[0]=norm_a; - cel_const[1]=norm_b/norm_a; - cel_const[2]=norm_c/norm_a; - cel_const[3]=alpha; - cel_const[4]=beta; - cel_const[5]=gamma; - } - } - - return type; -} - -//--------------------------------------------------- -// The lattice will be transformed to a 'standard -// cystallographic setting', the relation between -// 'origin' and 'transformed' lattice vectors will -// be givin in matrix form -// must be called before symmetry analysis -// only need to called once for each ion step -//--------------------------------------------------- -void Symmetry::lattice_type( - ModuleBase::Vector3 &v1, - ModuleBase::Vector3 &v2, - ModuleBase::Vector3 &v3, - ModuleBase::Vector3 &v01, - ModuleBase::Vector3 &v02, - ModuleBase::Vector3 &v03, - double *cel_const, - double *pre_const, - int& real_brav, - std::string& bravname, - const Atom* atoms, - bool convert_atoms, - double* newpos)const -{ - ModuleBase::TITLE("Symmetry","lattice_type"); -// std::cout << "v1 = " << v1.x << " " << v1.y << " " << v1.z <get_shortest_latvec(v1, v2, v3); -// std::cout << "a1 = " << v1.x << " " << v1.y << " " << v1.z < w1, w2, w3; - ModuleBase::Vector3 q1, q2, q3; - this->get_optlat(v1, v2, v3, w1, w2, w3, real_brav, cel_const, temp_const); -// std::cout << "a1 = " << v1.x << " " << v1.y << " " << v1.z <ntype; ++it) - { - for (int ia = 0; ia < this->na[it]; ++ia) - { - ModuleBase::Mathzone::Cartesian_to_Direct(atoms[it].tau[ia].x, - atoms[it].tau[ia].y, - atoms[it].tau[ia].z, - q1.x, q1.y, q1.z, - q2.x, q2.y, q2.z, - q3.x, q3.y, q3.z, - newpos[3*at],newpos[3*at+1],newpos[3*at+2]); - - // std::cout << " newpos_before = " << newpos[3*iat] << " " << newpos[3*iat+1] << " " << newpos[3*iat+2] << std::endl; - // GlobalV::ofs_running << " newpos_before = " << newpos[3*iat] << " " << newpos[3*iat+1] << " " << newpos[3*iat+2] << std::endl; - for(int k=0; k<3; ++k) - { - this->check_translation( newpos[at*3+k], -floor(newpos[at*3+k])); - this->check_boundary( newpos[at*3+k] ); - } - // std::cout << " newpos = " << newpos[3*iat] << " " << newpos[3*iat+1] << " " << newpos[3*iat+2] << std::endl; - - // GlobalV::ofs_running << " newpos = " << newpos[3*iat] << " " << newpos[3*iat+1] << " " << newpos[3*iat+2] << std::endl; - ++at; - } - } - /* - std::stringstream ss; - ss << PARAM.globalv.global_out_dir << "STRU_SIMPLE.cif"; - - std::ofstream ofs( ss.str().c_str() ); - ofs << "Lattice vector : " << std::endl; - ofs << q1.x <<" "<ntype; it++) - { - for (int ia = 0; ia < this->na[it]; ia++) - { - ofs << atoms[it].label - << " " << newpos[3*at] - << " " << newpos[3*at+1] - << " " << newpos[3*at+2] << std::endl; - at++; - } - } - ofs.close(); - */ - } - // return the optimized lattice in v1, v2, v3 - v1=q1; - v2=q2; - v3=q3; - } - else - { - //else, store the original ones - for (int i = 0; i < 6; ++i) - { - cel_const[i] = pre_const[i]; - } - //newpos also need to be set - if(convert_atoms) - { - int at=0; - for (int it = 0; it < this->ntype; ++it) - { - for (int ia = 0; ia < this->na[it]; ++ia) - { - ModuleBase::Mathzone::Cartesian_to_Direct(atoms[it].tau[ia].x, - atoms[it].tau[ia].y, - atoms[it].tau[ia].z, - v1.x, v1.y, v1.z, - v2.x, v2.y, v2.z, - v3.x, v3.y, v3.z, - newpos[3*at],newpos[3*at+1],newpos[3*at+2]); - for(int k=0; k<3; ++k) - { - this->check_translation( newpos[at*3+k], -floor(newpos[at*3+k])); - this->check_boundary( newpos[at*3+k] ); - } - ++at; - } - } - } - } - - /* - bool flag3; - if (pre_brav == temp_brav) - { - flag3 = 0; - if (!equal(temp_const[0], pre_const[0]) || - !equal(temp_const[1], pre_const[1]) || - !equal(temp_const[2], pre_const[2]) || - !equal(temp_const[3], pre_const[3]) || - !equal(temp_const[4], pre_const[4]) || - !equal(temp_const[5], pre_const[5]) - ) - { - flag3 = 1; - } - if (flag3==0) { - //particularly, if the symmetry of origin and new are exactly the same, we choose the original ones - //Hey! the original vectors have been changed!!! - v1 = s1; - v2 = s2; - v3 = s3; - change=0; - GlobalV::ofs_running<<" The lattice vectors have been set back!"<* gtrans, - double* pos, double* rotpos, int* index, const int ntype, const int itmin_type, const int itmin_start, int* istart, int* na)const -{ - ModuleBase::TITLE("Symmetry", "getgroup"); - - //-------------------------------------------------------------------------------- - //return all possible space group operators that reproduce a lattice with basis - //out of a (maximum) pool of point group operations that is compatible with - //the symmetry of the pure translation lattice without any basic. - //-------------------------------------------------------------------------------- - - ModuleBase::Matrix3 zero(0,0,0,0,0,0,0,0,0); - ModuleBase::Matrix3 help[48]; - ModuleBase::Vector3 temp[48]; - - nrot = 0; - nrotk = 0; - - //------------------------------------------------------------------------- - //pass through the pool of (possibly allowed) symmetry operations and - //check each operation whether it can reproduce the lattice with basis - //------------------------------------------------------------------------- - //std::cout << "nop = " <checksym(symop[i], gtrans[i], pos, rotpos, index, ntype, itmin_type, itmin_start, istart, na); - if (s_flag == 1) - { - //------------------------------ - // this is a symmetry operation - // with no translation vectors - // so ,this is pure point group - // operations - //------------------------------ - if ( equal(gtrans[i].x,0.0) && - equal(gtrans[i].y,0.0) && - equal(gtrans[i].z,0.0)) - { - ++nrot; - gmatrix[nrot - 1] = symop[i]; - gtrans[nrot - 1].x = 0; - gtrans[nrot - 1].y = 0; - gtrans[nrot - 1].z = 0; - } - //------------------------------ - // this is a symmetry operation - // with translation vectors - // so ,this is space group - // operations - //------------------------------ - else - { - ++nrotk; - help[nrotk - 1] = symop[i]; - temp[nrotk - 1].x = gtrans[i].x; - temp[nrotk - 1].y = gtrans[i].y; - temp[nrotk - 1].z = gtrans[i].z; - } - } - } - - //----------------------------------------------------- - //If there are operations with nontrivial translations - //then store them together in the momory - //----------------------------------------------------- - if (nrotk > 0) - { - for (int i = 0; i < nrotk; ++i) - { - gmatrix[nrot + i] = help[i]; - gtrans[nrot + i].x = temp[i].x; - gtrans[nrot + i].y = temp[i].y; - gtrans[nrot + i].z = temp[i].z; - } - } - - //----------------------------------------------------- - //total number of space group operations - //----------------------------------------------------- - nrotk += nrot; - - if(test_brav) - { - ModuleBase::GlobalFunc::OUT(ofs_running,"PURE POINT GROUP OPERATIONS",nrot); - ModuleBase::GlobalFunc::OUT(ofs_running,"SPACE GROUP OPERATIONS",nrotk); - } - - //----------------------------------------------------- - //fill the rest of matrices and vectors with zeros - //----------------------------------------------------- - if (nrotk < 48) - { - for (int i = nrotk; i < 48; ++i) - { - gmatrix[i] = zero; - gtrans[i].x = 0; - gtrans[i].y = 0; - gtrans[i].z = 0; - } - } - return; -} - -bool Symmetry::checksym(const ModuleBase::Matrix3& s, ModuleBase::Vector3& gtrans, - double* pos, double* rotpos, int* index, const int ntype, const int itmin_type, const int itmin_start, int* istart, int* na)const -{ - //---------------------------------------------- - // checks whether a point group symmetry element - // is a valid symmetry operation on a supercell - //---------------------------------------------- - // the start atom index. - bool no_diff = false; - ModuleBase::Vector3 trans(2.0, 2.0, 2.0); - bool s_flag = false; - - for (int it = 0; it < ntype; it++) - { - //------------------------------------ - // impose periodic boundary condition - // 0.5 -> -0.5 - //------------------------------------ - for (int j = istart[it]; j < istart[it] + na[it]; ++j) - { - this->check_boundary(pos[j*3+0]); - this->check_boundary(pos[j*3+1]); - this->check_boundary(pos[j*3+2]); - } - //order original atomic positions for current species - this->atom_ordering_new(pos + istart[it] * 3, na[it], index + istart[it]); - - //Rotate atoms of current species - for (int j = istart[it]; j < istart[it] + na[it]; ++j) - { - const int xx=j*3; - const int yy=j*3+1; - const int zz=j*3+2; - - - rotpos[xx] = pos[xx] * s.e11 - + pos[yy] * s.e21 - + pos[zz] * s.e31; - - rotpos[yy] = pos[xx] * s.e12 - + pos[yy] * s.e22 - + pos[zz] * s.e32; - - rotpos[zz] = pos[xx] * s.e13 - + pos[yy] * s.e23 - + pos[zz] * s.e33; - - // std::cout << "pos = " << pos[xx] <<" "<check_boundary(rotpos[xx]); - this->check_boundary(rotpos[yy]); - this->check_boundary(rotpos[zz]); - } - //order rotated atomic positions for current species - this->atom_ordering_new(rotpos + istart[it] * 3, na[it], index + istart[it]); - } - - ModuleBase::Vector3 diff; - - //--------------------------------------------------------- - // itmin_start = the start atom positions of species itmin - //--------------------------------------------------------- - // (s)tart (p)osition of atom (t)ype which has (min)inal number. - ModuleBase::Vector3 sptmin(rotpos[itmin_start * 3], rotpos[itmin_start * 3 + 1], rotpos[itmin_start * 3 + 2]); - - for (int i = itmin_start; i < itmin_start + na[itmin_type]; ++i) - { - //set up the current test std::vector "gtrans" - //and "gtrans" could possibly contain trivial translations: - gtrans.x = this->get_translation_vector( sptmin.x, pos[i*3+0]); - gtrans.y = this->get_translation_vector( sptmin.y, pos[i*3+1]); - gtrans.z = this->get_translation_vector( sptmin.z, pos[i*3+2]); - - //If we had already detected some translation, - //we must only look at the vectors with coordinates smaller than those - //of the previously detected std::vector (find the smallest) - if (gtrans.x > trans.x + epsilon || - gtrans.y > trans.y + epsilon || - gtrans.z > trans.z + epsilon - ) - { - continue; - } - - //translate all the atomic coordinates BACK by "gtrans" - for (int it = 0; it < ntype; it++) - { - for (int ia = istart[it]; ia < na[it] + istart[it]; ia++) - { - this->check_translation( rotpos[ia*3+0], gtrans.x ); - this->check_translation( rotpos[ia*3+1], gtrans.y ); - this->check_translation( rotpos[ia*3+2], gtrans.z ); - - this->check_boundary( rotpos[ia*3+0] ); - this->check_boundary( rotpos[ia*3+1] ); - this->check_boundary( rotpos[ia*3+2] ); - } - //order translated atomic positions for current species - this->atom_ordering_new(rotpos + istart[it] * 3, na[it], index + istart[it]); - } - - no_diff = true; - //compare the two lattices 'one-by-one' whether they are identical - for (int it = 0; it < ntype; it++) - { - for (int ia = istart[it]; ia < na[it] + istart[it]; ia++) - { - //take the difference of the rotated and the original coordinates - diff.x = this->check_diff( pos[ia*3+0], rotpos[ia*3+0]); - diff.y = this->check_diff( pos[ia*3+1], rotpos[ia*3+1]); - diff.z = this->check_diff( pos[ia*3+2], rotpos[ia*3+2]); - //only if all "diff" are zero vectors, flag will remain "1" - if ( no_diff == false|| - !equal(diff.x,0.0)|| - !equal(diff.y,0.0)|| - !equal(diff.z,0.0) - ) - { - no_diff = false; - } - } - } - - - //BLOCK_HERE("check symm"); - - //the current test is successful - if (no_diff == true) - { - s_flag = true; - //save the detected translation std::vector temporarily - trans.x = gtrans.x; - trans.y = gtrans.y; - trans.z = gtrans.z; - } - - //restore the original rotated coordinates by subtracting "gtrans" - for (int it = 0; it < ntype; it++) - { - for (int ia = istart[it]; ia < na[it] + istart[it]; ia++) - { - rotpos[ia*3+0] -= gtrans.x; - rotpos[ia*3+1] -= gtrans.y; - rotpos[ia*3+2] -= gtrans.z; - } - } - } - - if (s_flag == 1) - { - gtrans.x = trans.x; - gtrans.y = trans.y; - gtrans.z = trans.z; - } - return s_flag; -} - -void Symmetry::pricell(double* pos, const Atom* atoms) -{ - bool no_diff = false; - ptrans.clear(); - - for (int it = 0; it < ntype; it++) - { - //------------------------------------ - // impose periodic boundary condition - // 0.5 -> -0.5 - //------------------------------------ - for (int j = istart[it]; j < istart[it] + na[it]; ++j) - { - this->check_boundary(pos[j*3+0]); - this->check_boundary(pos[j*3+1]); - this->check_boundary(pos[j*3+2]); - } - - //order original atomic positions for current species - this->atom_ordering_new(pos + istart[it] * 3, na[it], index + istart[it]); - //copy pos to rotpos - for (int j = istart[it]; j < istart[it] + na[it]; ++j) - { - const int xx=j*3; - const int yy=j*3+1; - const int zz=j*3+2; - rotpos[xx] = pos[xx]; - rotpos[yy] = pos[yy]; - rotpos[zz] = pos[zz]; - } - } - - ModuleBase::Vector3 diff; - double tmp_ptrans[3]; - - //--------------------------------------------------------- - // itmin_start = the start atom positions of species itmin - //--------------------------------------------------------- - // (s)tart (p)osition of atom (t)ype which has (min)inal number. - ModuleBase::Vector3 sptmin(pos[itmin_start * 3], pos[itmin_start * 3 + 1], pos[itmin_start * 3 + 2]); - - for (int i = itmin_start; i < itmin_start + na[itmin_type]; ++i) - { - //set up the current test std::vector "gtrans" - //and "gtrans" could possibly contain trivial translations: - tmp_ptrans[0] = this->get_translation_vector( pos[i*3+0], sptmin.x); - tmp_ptrans[1] = this->get_translation_vector( pos[i*3+1], sptmin.y); - tmp_ptrans[2] = this->get_translation_vector( pos[i*3+2], sptmin.z); - //translate all the atomic coordinates by "gtrans" - for (int it = 0; it < ntype; it++) - { - for (int ia = istart[it]; ia < na[it] + istart[it]; ia++) - { - this->check_translation( rotpos[ia*3+0], tmp_ptrans[0] ); - this->check_translation( rotpos[ia*3+1], tmp_ptrans[1] ); - this->check_translation( rotpos[ia*3+2], tmp_ptrans[2] ); - - this->check_boundary( rotpos[ia*3+0] ); - this->check_boundary( rotpos[ia*3+1] ); - this->check_boundary( rotpos[ia*3+2] ); - } - //order translated atomic positions for current species - this->atom_ordering_new(rotpos + istart[it] * 3, na[it], index + istart[it]); - } - - no_diff = true; - //compare the two lattices 'one-by-one' whether they are identical - for (int it = 0; it < ntype; it++) - { - for (int ia = istart[it]; ia < na[it] + istart[it]; ia++) - { - //take the difference of the rotated and the original coordinates - diff.x = this->check_diff( pos[ia*3+0], rotpos[ia*3+0]); - diff.y = this->check_diff( pos[ia*3+1], rotpos[ia*3+1]); - diff.z = this->check_diff( pos[ia*3+2], rotpos[ia*3+2]); - //only if all "diff" are zero vectors, flag will remain "1" - if (!equal(diff.x,0.0)|| - !equal(diff.y,0.0)|| - !equal(diff.z,0.0)) - { - no_diff = false; - break; - } - } - if (!no_diff) { - break; - } - } - - //the current test is successful - if (no_diff) { - ptrans.push_back(ModuleBase::Vector3(tmp_ptrans[0], - tmp_ptrans[1], - tmp_ptrans[2])); - } - //restore the original rotated coordinates by subtracting "ptrans" - for (int it = 0; it < ntype; it++) - { - for (int ia = istart[it]; ia < na[it] + istart[it]; ia++) - { - rotpos[ia*3+0] -= tmp_ptrans[0]; - rotpos[ia*3+1] -= tmp_ptrans[1]; - rotpos[ia*3+2] -= tmp_ptrans[2]; - } - } - } - int ntrans=ptrans.size(); - if (ntrans <= 1) - { - GlobalV::ofs_running<<" Original cell was already a primitive cell."<p1=this->a1; - this->p2=this->a2; - this->p3=this->a3; - this->pbrav=this->real_brav; - this->ncell=1; - for (int i = 0; i < 6; ++i) { - this->pcel_const[i] = this->cel_const[i]; - } - return; - } - - //sort ptrans: - double* ptrans_array = new double[ntrans*3]; - for(int i=0;iatom_ordering_new(ptrans_array, ntrans, index); - // std::cout<<"final ptrans:"< b1, b2, b3; - int iplane=0, jplane=0, kplane=0; - //1. kplane for b3 - while (kplane < ntrans - && std::abs(ptrans[kplane].z - ptrans[0].z) < this->epsilon) { - ++kplane; - } - if (kplane == ntrans) { - kplane = 0; // a3-direction have no smaller pricell - } - b3=kplane>0 ? - ModuleBase::Vector3(ptrans[kplane].x, ptrans[kplane].y, ptrans[kplane].z) : - ModuleBase::Vector3(0, 0, 1); - //2. jplane for b2 (not collinear with b3) - jplane=kplane+1; - while (jplane < ntrans - && (std::abs(ptrans[jplane].y - ptrans[0].y) < this->epsilon - || equal((ptrans[jplane] ^ b3).norm(), 0))) { - ++jplane; - } - if (jplane == ntrans) { - jplane = kplane; // a2-direction have no smaller pricell - } - b2=jplane>kplane ? - ModuleBase::Vector3(ptrans[jplane].x, ptrans[jplane].y, ptrans[jplane].z) : - ModuleBase::Vector3(0, 1, 0); - //3. iplane for b1 (not coplane with ) - iplane=jplane+1; - while (iplane < ntrans - && (std::abs(ptrans[iplane].x - ptrans[0].x) < this->epsilon - || equal(ptrans[iplane] * (b2 ^ b3), 0))) { - ++iplane; - } - b1=(iplane>jplane && iplane(ptrans[iplane].x, ptrans[iplane].y, ptrans[iplane].z) : - ModuleBase::Vector3(1, 0, 0); //a1-direction have no smaller pricell - - // std::cout<<"iplane="<plat=coeff*this->optlat; - - //deal with collineation caused by default b1, b2, b3 - if(equal(plat.Det(), 0)) - { - if(kplane==0) //try a new b3 - { - std::cout<<"try a new b3"<kplane) // use default b2 - { - coeff.e31=0; - coeff.e32=1; - coeff.e33=0; - } - else //use default b1 - { - coeff.e31=1; - coeff.e32=0; - coeff.e33=0; - } - } - else if(jplane<=kplane) - { - //std::cout<<"try a new b2"<plat=coeff*this->optlat; - assert(!equal(plat.Det(), 0)); - } - - this->p1.x=plat.e11; - this->p1.y=plat.e12; - this->p1.z=plat.e13; - this->p2.x=plat.e21; - this->p2.y=plat.e22; - this->p2.z=plat.e23; - this->p3.x=plat.e31; - this->p3.y=plat.e32; - this->p3.z=plat.e33; - -#ifdef __DEBUG - GlobalV::ofs_running<<"lattice vectors of primitive cell (initial):"< p01=p1, p02=p2, p03=p3; - double pcel_pre_const[6]; - for (int i = 0; i < 6; ++i) { - pcel_pre_const[i] = pcel_const[i]; - } - this->lattice_type(p1, p2, p3, p01, p02, p03, pcel_const, pcel_pre_const, pbrav, pbravname, atoms, false, nullptr); - - this->plat.e11=p1.x; - this->plat.e12=p1.y; - this->plat.e13=p1.z; - this->plat.e21=p2.x; - this->plat.e22=p2.y; - this->plat.e23=p2.z; - this->plat.e31=p3.x; - this->plat.e32=p3.y; - this->plat.e33=p3.z; - -#ifdef __DEBUG - GlobalV::ofs_running<<"lattice vectors of primitive cell (optimized):"<pbrav, this->pcel_const, GlobalV::ofs_running); - - //count the number of pricells - GlobalV::ofs_running<<"optimized lattice volume: "<optlat.Det()<plat.Det()<optlat.Det()/this->plat.Det()); - this->ncell=floor(ncell_double+0.5); - - auto reset_pcell = [this]() -> void { - std::cout << " Now regard the structure as a primitive cell." << std::endl; - this->ncell = 1; - this->ptrans = std::vector >(1, ModuleBase::Vector3(0, 0, 0)); - GlobalV::ofs_running << "WARNING: Original cell may have more than one primitive cells, \ - but we have to treat it as a primitive cell. Use a larger `symmetry_prec`to avoid this warning." << std::endl; - }; - if (this->ncell != ntrans) - { - std::cout << " WARNING: PRICELL: NCELL != NTRANS !" << std::endl; - std::cout << " NCELL=" << ncell << ", NTRANS=" << ntrans << std::endl; - std::cout << " Suggest solution: Use a larger `symmetry_prec`. " << std::endl; - reset_pcell(); - return; - } - if(std::abs(ncell_double-double(this->ncell)) > this->epsilon*100) - { - std::cout << " WARNING: THE NUMBER OF PRIMITIVE CELL IS NOT AN INTEGER !" << std::endl; - std::cout << " NCELL(double)=" << ncell_double << ", NTRANS=" << ncell << std::endl; - std::cout << " Suggest solution: Use a larger `symmetry_prec`. " << std::endl; - reset_pcell(); - return; - } - GlobalV::ofs_running<<"Original cell was built up by "<ncell<<" primitive cells."<gtrans_convert(ptrans.data(), ptrans.data(), ntrans, this->optlat, inputlat ); - - //how many pcell in supercell - int n1, n2, n3; - ModuleBase::Matrix3 nummat0=this->optlat*this->plat.Inverse(); - ModuleBase::Matrix3 nummat, transmat; - hermite_normal_form(nummat0, nummat, transmat); - n1=floor (nummat.e11 + epsilon); - n2=floor (nummat.e22 + epsilon); - n3=floor (nummat.e33 + epsilon); - if(n1*n2*n3 != this->ncell) - { - std::cout << " WARNING: Number of cells and number of vectors did not agree."; - std::cout<<"Try to change symmetry_prec in INPUT." << std::endl; - reset_pcell(); - } - return; -} - - -//modified by shu on 2010.01.20 -void Symmetry::rho_symmetry( double *rho, - const int &nr1, const int &nr2, const int &nr3) -{ - ModuleBase::timer::tick("Symmetry","rho_symmetry"); - - // allocate flag for each FFT grid. - bool* symflag = new bool[nr1 * nr2 * nr3]; - for (int i=0; i0 ); - assert(nrotk <=48 ); - int *ri = new int[nrotk]; - int *rj = new int[nrotk]; - int *rk = new int[nrotk]; - - int ci = 0; - for (int i = 0; i< nr1; ++i) - { - for (int j = 0; j< nr2; ++j) - { - for (int k = 0; k< nr3; ++k) - { - if (!symflag[i * nr2 * nr3 + j * nr3 + k]) - { - double sum = 0; - - for (int isym = 0; isym < nrotk; ++isym) - { - this->rotate(gmatrix[isym], gtrans[isym], i, j, k, nr1, nr2, nr3, ri[isym], rj[isym], rk[isym]); - const int index = ri[isym] * nr2 * nr3 + rj[isym] * nr3 + rk[isym]; - sum += rho[ index ]; - } - sum /= nrotk; - - for (int isym = 0; isym < nrotk; ++isym) - { - const int index = ri[isym] * nr2 * nr3 + rj[isym] * nr3 + rk[isym]; - rho[index] = sum; - symflag[index] = true; - } - } - } - } - } - - delete[] symflag; - delete[] ri; - delete[] rj; - delete[] rk; - ModuleBase::timer::tick("Symmetry","rho_symmetry"); -} - - -void Symmetry::rhog_symmetry(std::complex *rhogtot, - int* ixyz2ipw, const int &nx, const int &ny, const int &nz, - const int &fftnx, const int &fftny, const int &fftnz) -{ - ModuleBase::timer::tick("Symmetry","rhog_symmetry"); -// ---------------------------------------------------------------------- -// the current way is to cluster the FFT grid points into groups in advance. -// and use OpenMP to realize parallel calculation, one thread works in one group. -// ---------------------------------------------------------------------- - - // allocate flag for each FFT grid. - int* symflag = new int[fftnx*fftny*fftnz];// which group the grid belongs to - int(*isymflag)[48] = new int[fftnx*fftny*fftnz][48];//which rotration operation the grid corresponds to - int(*table_xyz)[48] = new int[fftnx * fftny * fftnz][48];// group information - int* count_xyz = new int[fftnx * fftny * fftnz];// how many symmetry operations has been covered - for (int i = 0; i < fftnx * fftny * fftnz; i++) - { - symflag[i] = -1; - } - int group_index = 0; - - assert(nrotk >0 ); - assert(nrotk <=48 ); - - //map the gmatrix to inv - std::vectorinvmap(this->nrotk, -1); - this->gmatrix_invmap(kgmatrix, nrotk, invmap.data()); - -// --------------------------------------------------- -/* This code defines a lambda function called "rotate_recip" that takes - a 3x3 matrix and a 3D vector as input. It performs a rotation operation - on the vector using the matrix and returns the rotated vector. - Specifically, it calculates the new coordinates of the vector after - the rotation and applies periodic boundary conditions to ensure that - the coordinates are within the FFT-grid dimensions. - The rotated vector is returned by modifying the input vector. -*/ -// --------------------------------------------------- - //rotate function (different from real space, without scaling gmatrix) - auto rotate_recip = [&] (ModuleBase::Matrix3& g, ModuleBase::Vector3& g0, int& ii, int& jj, int& kk) - { - ii = int(g.e11 * g0.x + g.e21 * g0.y + g.e31 * g0.z) ; - if (ii < 0) - { - ii += 10 * nx; - } - ii = ii%nx; - jj = int(g.e12 * g0.x + g.e22 * g0.y + g.e32 * g0.z) ; - if (jj < 0) - { - jj += 10 * ny; - } - jj = jj%ny; - kk = int(g.e13 * g0.x + g.e23 * g0.y + g.e33 * g0.z); - if (kk < 0) - { - kk += 10 * nz; - } - kk = kk%nz; - return; - }; - -// ------------------------------------------- -/* - Trying to group fft grids first. - It iterates over each FFT-grid point and checks if it is within the PW-sphere. If it is, put all the FFT-grid points connected by the - rotation operation into one group( the index is stored in int(*table_xyz)). - The code marks the point as processed to avoid redundant calculations - by using int* symflag. -*/ -// ------------------------------------------- - ModuleBase::timer::tick("Symmetry","group_fft_grids"); - for (int i = 0; i< fftnx; ++i) - { - //tmp variable - ModuleBase::Vector3 tmp_gdirect0(0, 0, 0); - tmp_gdirect0.x=(i>int(nx/2)+1)?(i-nx):i; - for (int j = 0; j< fftny; ++j) - { - tmp_gdirect0.y=(j>int(ny/2)+1)?(j-ny):j; - for (int k = 0; k< fftnz; ++k) - { - int ixyz0=(i*fftny+j)*fftnz+k; - if (symflag[ixyz0] == -1) - { - int ipw0=ixyz2ipw[ixyz0]; - //if a fft-grid is not in pw-sphere, just do not consider it. - if (ipw0 == -1) { - continue; - } - tmp_gdirect0.z=(k>int(nz/2)+1)?(k-nz):k; - int rot_count=0; - for (int isym = 0; isym < nrotk; ++isym) - { - if (invmap[isym] < 0 || invmap[isym] > nrotk) { continue; } - //tmp variables - int ii, jj, kk=0; - rotate_recip(kgmatrix[invmap[isym]], tmp_gdirect0, ii, jj, kk); - if(ii>=fftnx || jj>=fftny || kk>= fftnz) - { - if(!PARAM.globalv.gamma_only_pw) - { - std::cout << " ROTATE OUT OF FFT-GRID IN RHOG_SYMMETRY !" << std::endl; - ModuleBase::QUIT(); - } - // for gamma_only_pw, just do not consider this rotation. - continue; - } - int ixyz=(ii*fftny+jj)*fftnz+kk; - //fft-grid index to (ip, ig) - int ipw=ixyz2ipw[ixyz]; - if(ipw==-1) //not in pw-sphere - { - continue; //else, just skip it - } - symflag[ixyz] = group_index; - isymflag[group_index][rot_count] = invmap[isym]; - table_xyz[group_index][rot_count] = ixyz; - ++rot_count; - assert(rot_count <= nrotk); - count_xyz[group_index] = rot_count; - } - group_index++; - } - } - } - } - ModuleBase::timer::tick("Symmetry","group_fft_grids"); - -// ------------------------------------------- -/* This code performs symmetry operations on the reciprocal space - charge density using FFT-grids. It iterates over each FFT-grid - point in a particular group, applies a phase factor and sums the - charge density over the symmetry operations, and then divides by - the number of symmetry operations. Finally, it updates the charge - density for each FFT-grid point using the calculated sum. -*/ -// ------------------------------------------- -#ifdef _OPENMP -#pragma omp parallel for schedule(static) -#endif -for (int g_index = 0; g_index < group_index; g_index++) -{ - // record the index and gphase but not the final gdirect for each symm-opt - int *ipw_record = new int[nrotk]; - int *ixyz_record = new int[nrotk]; - std::complex* gphase_record = new std::complex [nrotk]; - std::complex sum(0, 0); - int rot_count=0; - - for (int c_index = 0; c_index < count_xyz[g_index]; ++c_index) - { - int ixyz0 = table_xyz[g_index][c_index]; - int ipw0 = ixyz2ipw[ixyz0]; - if (symflag[ixyz0] == g_index) - { - // note : do not use PBC after rotation. - // we need a real gdirect to get the correspoding rhogtot. - int k = ixyz0%fftnz; - int j = ((ixyz0-k)/fftnz)%fftny; - int i = ((ixyz0-k)/fftnz-j)/fftny; - //fft-grid index to gdirect - ModuleBase::Vector3 tmp_gdirect_double(0.0, 0.0, 0.0); - tmp_gdirect_double.x=static_cast((i>int(nx/2)+1)?(i-nx):i); - tmp_gdirect_double.y=static_cast((j>int(ny/2)+1)?(j-ny):j); - tmp_gdirect_double.z=static_cast((k>int(nz/2)+1)?(k-nz):k); - //calculate phase factor - tmp_gdirect_double = tmp_gdirect_double * ModuleBase::TWO_PI; - double cos_arg = 0.0, sin_arg = 0.0; - double arg_gtrans = tmp_gdirect_double * gtrans[isymflag[g_index][c_index]]; - std::complex phase_gtrans (ModuleBase::libm::cos(arg_gtrans), ModuleBase::libm::sin(arg_gtrans)); - // for each pricell in supercell: - for (int ipt = 0;ipt < ((ModuleSymmetry::Symmetry::pricell_loop) ? this->ncell : 1);++ipt) - { - double arg = tmp_gdirect_double * ptrans[ipt]; - double tmp_cos = 0.0, tmp_sin = 0.0; - ModuleBase::libm::sincos(arg, &tmp_sin, &tmp_cos); - cos_arg += tmp_cos; - sin_arg += tmp_sin; - } - // add nothing to sum, so don't consider this isym into rot_count - cos_arg/=static_cast(ncell); - sin_arg/=static_cast(ncell); - //deal with double-zero - if (equal(cos_arg, 0.0) && equal(sin_arg, 0.0)) { - continue; - } - std::complex gphase(cos_arg, sin_arg); - gphase = phase_gtrans * gphase; - //deal with small difference from 1 - if (equal(gphase.real(), 1.0) && equal(gphase.imag(), 0)) { - gphase = std::complex(1.0, 0.0); - } - gphase_record[rot_count]=gphase; - sum += rhogtot[ipw0]*gphase; - //record - ipw_record[rot_count]=ipw0; - ixyz_record[rot_count]=ixyz0; - ++rot_count; - //assert(rot_count<=nrotk); - }//end if section - }//end c_index loop - sum /= rot_count; - for (int isym = 0; isym < rot_count; ++isym) - { - rhogtot[ipw_record[isym]] = sum/gphase_record[isym]; - } - //g_index++; - //Clean the records variables for each fft grid point - delete[] ipw_record; - delete[] ixyz_record; - delete[] gphase_record; -}//end g_index loop - delete[] symflag; - delete[] isymflag; - delete[] table_xyz; - delete[] count_xyz; - ModuleBase::timer::tick("Symmetry","rhog_symmetry"); -} void Symmetry::set_atom_map(const Atom* atoms) { @@ -1760,9 +359,15 @@ void Symmetry::set_atom_map(const Atom* atoms) const int xx = ia * 3; const int yy = ia * 3 + 1; const int zz = ia * 3 + 2; for (int k = 0;k < this->nrotk;++k) { - rotpos[xx] = pos[xx] * gmatrix[k].e11 + pos[yy] * gmatrix[k].e21 + pos[zz] * gmatrix[k].e31 + gtrans[k].x; - rotpos[yy] = pos[xx] * gmatrix[k].e12 + pos[yy] * gmatrix[k].e22 + pos[zz] * gmatrix[k].e32 + gtrans[k].y; - rotpos[zz] = pos[xx] * gmatrix[k].e13 + pos[yy] * gmatrix[k].e23 + pos[zz] * gmatrix[k].e33 + gtrans[k].z; + rotpos[xx] = pos[xx] * gmatrix[k].e11 + + pos[yy] * gmatrix[k].e21 + + pos[zz] * gmatrix[k].e31 + gtrans[k].x; + rotpos[yy] = pos[xx] * gmatrix[k].e12 + + pos[yy] * gmatrix[k].e22 + + pos[zz] * gmatrix[k].e32 + gtrans[k].y; + rotpos[zz] = pos[xx] * gmatrix[k].e13 + + pos[yy] * gmatrix[k].e23 + + pos[zz] * gmatrix[k].e33 + gtrans[k].z; check_translation(rotpos[xx], -floor(rotpos[xx])); check_boundary(rotpos[xx]); @@ -1805,9 +410,9 @@ void Symmetry::symmetrize_vec3_nat(double* v)const // pengfei 2016-12-20 if (l < 0) { continue; } - vtot[l * 3] = vtot[l * 3] + v[jx] * gmatrix[k].e11 + v[jy] * gmatrix[k].e21 + v[jz] * gmatrix[k].e31; - vtot[l * 3 + 1] = vtot[l * 3 + 1] + v[jx] * gmatrix[k].e12 + v[jy] * gmatrix[k].e22 + v[jz] * gmatrix[k].e32; - vtot[l * 3 + 2] = vtot[l * 3 + 2] + v[jx] * gmatrix[k].e13 + v[jy] * gmatrix[k].e23 + v[jz] * gmatrix[k].e33; + vtot[l*3] = vtot[l*3] + v[jx] * gmatrix[k].e11 + v[jy] * gmatrix[k].e21 + v[jz] * gmatrix[k].e31; + vtot[l*3+1] = vtot[l*3+1] + v[jx] * gmatrix[k].e12 + v[jy] * gmatrix[k].e22 + v[jz] * gmatrix[k].e32; + vtot[l*3+2] = vtot[l*3+2] + v[jx] * gmatrix[k].e13 + v[jy] * gmatrix[k].e23 + v[jz] * gmatrix[k].e33; n[l]++; } } @@ -1859,6 +464,7 @@ void Symmetry::gmatrix_convert_int(const ModuleBase::Matrix3* sa, ModuleBase::Ma sb[i].e33=round(sb[i].e33); } } + void Symmetry::gmatrix_convert(const ModuleBase::Matrix3* sa, ModuleBase::Matrix3* sb, const int n, const ModuleBase::Matrix3 &a, const ModuleBase::Matrix3 &b)const { @@ -1869,6 +475,7 @@ void Symmetry::gmatrix_convert(const ModuleBase::Matrix3* sa, ModuleBase::Matrix sb[i]=b*ai*sa[i]*a*bi; } } + void Symmetry::gtrans_convert(const ModuleBase::Vector3* va, ModuleBase::Vector3* vb, const int n, const ModuleBase::Matrix3 &a, const ModuleBase::Matrix3 &b)const { @@ -1878,6 +485,7 @@ void Symmetry::gtrans_convert(const ModuleBase::Vector3* va, ModuleBase: vb[i]=va[i]*a*bi; } } + void Symmetry::gmatrix_invmap(const ModuleBase::Matrix3* s, const int n, int* invmap) const { ModuleBase::Matrix3 eig(1, 0, 0, 0, 1, 0, 0, 0, 1); @@ -1992,38 +600,8 @@ void Symmetry::get_optlat(ModuleBase::Vector3 &v1, ModuleBase::Vector3 &v1, ModuleBase::Vector3 h(1, 1) - epsilon) { - for(int r=0;r<3;++r) {h(r, 0) -= h(r, 1); b(r, 0) -= b(r, 1); - } - } - while (h(1, 0) < -epsilon) { - for(int r=0;r<3;++r) {h(r, 0) += h(r, 1); b(r, 0) += b(r, 1); - } - } + while (h(1, 0) > h(1, 1) - epsilon) + { + for(int r=0;r<3;++r) + { + h(r, 0) -= h(r, 1); + b(r, 0) -= b(r, 1); + } + } + while (h(1, 0) < -epsilon) + { + for(int r=0;r<3;++r) + { + h(r, 0) += h(r, 1); + b(r, 0) += b(r, 1); + } + } for(int j=0;j<2;++j) { - while (h(2, j) > h(2, 2) - epsilon) { - for(int r=0;r<3;++r) {h(r, j) -= h(r, 2); b(r, j) -= b(r, 2); - } + while (h(2, j) > h(2, 2) - epsilon) + { + for(int r=0;r<3;++r) + { + h(r, j) -= h(r, 2); + b(r, j) -= b(r, 2); + } } - while (h(2, j) < -epsilon) { - for(int r=0;r<3;++r) {h(r, j) += h(r, 2); b(r, j) += b(r, 2); - } + while (h(2, j) < -epsilon) + { + for(int r=0;r<3;++r) + { + h(r, j) += h(r, 2); + b(r, j) += b(r, 2); + } } } @@ -2221,35 +842,17 @@ void Symmetry::hermite_normal_form(const ModuleBase::Matrix3 &s3, ModuleBase::Ma //check s*b=h ModuleBase::matrix check_zeros = s3.to_matrix() * b - h; #ifdef __DEBUG - for (int i = 0;i < 3;++i) - for(int j=0;j<3;++j) - assert(near_equal(check_zeros(i, j), 0)); + for (int i = 0;i < 3;++i) + { + for(int j=0;j<3;++j) + { + assert(near_equal(check_zeros(i, j), 0)); + } + } #endif return; } -bool Symmetry::magmom_same_check(const Atom* atoms)const -{ - ModuleBase::TITLE("Symmetry", "magmom_same_check"); - bool pricell_loop = true; - for (int it = 0;it < ntype;++it) - { - if (pricell_loop) { - for (int ia = 1;ia < atoms[it].na;++ia) - { - if (!equal(atoms[it].m_loc_[ia].x, atoms[it].m_loc_[0].x) || - !equal(atoms[it].m_loc_[ia].y, atoms[it].m_loc_[0].y) || - !equal(atoms[it].m_loc_[ia].z, atoms[it].m_loc_[0].z)) - { - pricell_loop = false; - break; - } - } - } - } - return pricell_loop; -} - bool Symmetry::is_all_movable(const Atom* atoms, const Statistics& st)const { bool all_mbl = true; @@ -2265,70 +868,3 @@ bool Symmetry::is_all_movable(const Atom* atoms, const Statistics& st)const } return all_mbl; } - -void Symmetry::analyze_magnetic_group(const Atom* atoms, const Statistics& st, int& nrot_out, int& nrotk_out) -{ - // 1. classify atoms with different magmom - // (use symmetry_prec to judge if two magmoms are the same) - std::vector> mag_type_atoms; - for (int it = 0;it < ntype;++it) - { - for (int ia = 0; ia < atoms[it].na; ++ia) - { - bool find = false; - for (auto& mt : mag_type_atoms) - { - const int mag_iat = *mt.begin(); - const int mag_it = st.iat2it[mag_iat]; - const int mag_ia = st.iat2ia[mag_iat]; - if (it == mag_it && this->equal(atoms[it].mag[ia], atoms[mag_it].mag[mag_ia])) - { - mt.insert(st.itia2iat(it, ia)); - find = true; - break; - } - } - if (!find) - { - mag_type_atoms.push_back(std::set({ st.itia2iat(it,ia) })); - } - } - } - - // 2. get the start index, number of atoms and positions for each mag_type - std::vector mag_istart(mag_type_atoms.size()); - std::vector mag_na(mag_type_atoms.size()); - std::vector mag_pos; - int mag_itmin_type = 0; - int mag_itmin_start = 0; - for (int mag_it = 0;mag_it < mag_type_atoms.size(); ++mag_it) - { - mag_na[mag_it] = mag_type_atoms.at(mag_it).size(); - if (mag_it > 0) - { - mag_istart[mag_it] = mag_istart[mag_it - 1] + mag_na[mag_it - 1]; - } - if (mag_na[mag_it] < mag_na[itmin_type]) - { - mag_itmin_type = mag_it; - mag_itmin_start = mag_istart[mag_it]; - } - for (auto& mag_iat : mag_type_atoms.at(mag_it)) - { - // this->newpos have been ordered by original structure(ntype, na), it cannot be directly used here. - // we need to reset the calculate again the coordinate of the new structure. - const ModuleBase::Vector3 direct_tmp = atoms[st.iat2it[mag_iat]].tau[st.iat2ia[mag_iat]] * this->optlat.Inverse(); - std::array direct = { direct_tmp.x, direct_tmp.y, direct_tmp.z }; - for (int i = 0; i < 3; ++i) - { - this->check_translation(direct[i], -floor(direct[i])); - this->check_boundary(direct[i]); - mag_pos.push_back(direct[i]); - } - } - } - // 3. analyze the effective structure - this->getgroup(nrot_out, nrotk_out, GlobalV::ofs_running, this->nop, this->symop, this->gmatrix, this->gtrans, - mag_pos.data(), this->rotpos, this->index, mag_type_atoms.size(), mag_itmin_type, mag_itmin_start, mag_istart.data(), mag_na.data()); -} -} diff --git a/source/module_cell/module_symmetry/symmetry.h b/source/module_cell/module_symmetry/symmetry.h index c3ef81607d4..0d9133edbec 100644 --- a/source/module_cell/module_symmetry/symmetry.h +++ b/source/module_cell/module_symmetry/symmetry.h @@ -7,12 +7,16 @@ namespace ModuleSymmetry { + class Symmetry : public Symmetry_Basic { + public: - Symmetry() { - this->epsilon = 1e-6; - }; + + Symmetry() + { + this->epsilon = 1e-6; + }; ~Symmetry() {}; //symmetry flag for levels @@ -82,19 +86,42 @@ class Symmetry : public Symmetry_Basic bool all_mbl = true; ///< whether all the atoms are movable in all the directions - int standard_lat(ModuleBase::Vector3& a, ModuleBase::Vector3& b, ModuleBase::Vector3& c, double* celconst)const; - - void lattice_type(ModuleBase::Vector3 &v1,ModuleBase::Vector3 &v2,ModuleBase::Vector3 &v3, - ModuleBase::Vector3& v01, ModuleBase::Vector3& v02, ModuleBase::Vector3& v03, - double* cel_const, double* pre_const, int& real_brav, std::string& bravname, const Atom* atoms, - bool convert_atoms, double* newpos = nullptr)const; + int standard_lat(ModuleBase::Vector3& a, + ModuleBase::Vector3& b, + ModuleBase::Vector3& c, + double* celconst)const; + + void lattice_type(ModuleBase::Vector3 &v1, + ModuleBase::Vector3 &v2, + ModuleBase::Vector3 &v3, + ModuleBase::Vector3 &v01, + ModuleBase::Vector3 &v02, + ModuleBase::Vector3 &v03, + double* cel_const, + double* pre_const, + int& real_brav, + std::string& bravname, + const Atom* atoms, + bool convert_atoms, + double* newpos = nullptr)const; + + void getgroup(int& nrot, + int& nrotk, + std::ofstream& ofs_running, + const int& nop, + const ModuleBase::Matrix3* symop, + ModuleBase::Matrix3* gmatrix, + ModuleBase::Vector3* gtrans, + double* pos, double* rotpos, int* index, + const int ntype, const int itmin_type, const int itmin_start, + int* istart, int* na)const; + + bool checksym(const ModuleBase::Matrix3 &s, + ModuleBase::Vector3& gtrans, + double* pos, double* rotpos, int* index, + const int itmin_type, const int ntype, const int itmin_start, + int* istart, int* na)const; - void getgroup(int& nrot, int& nrotk, std::ofstream& ofs_running, const int& nop, - const ModuleBase::Matrix3* symop, ModuleBase::Matrix3* gmatrix, ModuleBase::Vector3* gtrans, - double* pos, double* rotpos, int* index, const int ntype, const int itmin_type, const int itmin_start, int* istart, int* na)const; - bool checksym(const ModuleBase::Matrix3& s, ModuleBase::Vector3& gtrans, - double* pos, double* rotpos, int* index, const int itmin_type, - const int ntype, const int itmin_start, int* istart, int* na)const; /// @brief primitive cell analysis void pricell(double* pos, const Atom* atoms); @@ -102,28 +129,37 @@ class Symmetry : public Symmetry_Basic /// Symmetrize the charge density, the forces, and the stress /// ----------------------- void rho_symmetry(double *rho, const int &nr1, const int &nr2, const int &nr3); + void rhog_symmetry(std::complex *rhogtot, int* ixyz2ipw, const int &nx, const int &ny, const int &nz, const int & fftnx, const int &fftny, const int &fftnz); + /// symmetrize a vector3 with nat elements, which can be forces or variation of atom positions in relax void symmetrize_vec3_nat(double* v)const; // force + /// symmetrize a 3*3 tensor, which can be stress or variation of unitcell in cell-relax void symmetrize_mat3(ModuleBase::matrix& sigma, const Lattice& lat)const; // stress //convert n rotation-matrices from sa on basis {a1, a2, a3} to sb on basis {b1, b2, b3} void gmatrix_convert(const ModuleBase::Matrix3* sa, ModuleBase::Matrix3* sb, const int n, const ModuleBase::Matrix3 &a, const ModuleBase::Matrix3 &b)const; + void gmatrix_convert_int(const ModuleBase::Matrix3* sa, ModuleBase::Matrix3* sb, const int n, const ModuleBase::Matrix3 &a, const ModuleBase::Matrix3 &b)const; + //convert n translation-vectors from va on basis {a1, a2, a3} to vb on basis {b1, b2, b3} void gtrans_convert(const ModuleBase::Vector3* va, ModuleBase::Vector3* vb, const int n, const ModuleBase::Matrix3 &a, const ModuleBase::Matrix3 &b)const; + void gmatrix_invmap(const ModuleBase::Matrix3* s, const int n, int* invmap) const; + void hermite_normal_form(const ModuleBase::Matrix3 &s, ModuleBase::Matrix3 &H, ModuleBase::Matrix3 &b) const; + int get_rotated_atom(int isym, int iat)const { if (!this->isym_rotiat_.empty()) { return this->isym_rotiat_[isym][iat]; } else { return -1; } } + private: /// atom-map for each symmetry operation: isym_rotiat[isym][iat]=rotiat @@ -138,12 +174,14 @@ class Symmetry : public Symmetry_Basic // to be called in lattice_type void get_shortest_latvec(ModuleBase::Vector3 &a1, ModuleBase::Vector3 &a2, ModuleBase::Vector3 &a3)const; + void get_optlat(ModuleBase::Vector3 &v1, ModuleBase::Vector3 &v2, ModuleBase::Vector3 &v3, ModuleBase::Vector3 &w1, ModuleBase::Vector3 &w2, ModuleBase::Vector3 &w3, int& real_brav, double* cel_const, double* tmp_const)const; - /// Loop the magmom of each atoms in its type when NSPIN>1. If not all the same, primitive cells should not be looped in rhog_symmetry. + /// Loop the magmom of each atoms in its type when NSPIN>1. + /// If not all the same, primitive cells should not be looped in rhog_symmetry. bool magmom_same_check(const Atom* atoms)const; /// Analyze magnetic group without time-reversal symmetry From e3f687badbaf3a6e89d06becaae20d62ffbae8e2 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Tue, 25 Mar 2025 15:16:37 +0800 Subject: [PATCH 03/12] update symmetry --- .../module_symmetry/symm_hermite.cpp | 219 ++++++++++++++++++ .../module_cell/module_symmetry/symmetry.cpp | 219 +----------------- 2 files changed, 220 insertions(+), 218 deletions(-) create mode 100644 source/module_cell/module_symmetry/symm_hermite.cpp diff --git a/source/module_cell/module_symmetry/symm_hermite.cpp b/source/module_cell/module_symmetry/symm_hermite.cpp new file mode 100644 index 00000000000..f74192a70a1 --- /dev/null +++ b/source/module_cell/module_symmetry/symm_hermite.cpp @@ -0,0 +1,219 @@ + +void Symmetry::hermite_normal_form(const ModuleBase::Matrix3 &s3, + ModuleBase::Matrix3 &h3, + ModuleBase::Matrix3 &b3) const +{ + ModuleBase::TITLE("Symmetry","hermite_normal_form"); + // check the non-singularity and integer elements of s +#ifdef __DEBUG + assert(!equal(s3.Det(), 0.0)); +#endif + + auto near_equal = [this](double x, double y) + { + return fabs(x - y) < 10 * epsilon; + }; + + ModuleBase::matrix s = s3.to_matrix(); + + for (int i = 0; i < 3; ++i) + { + for (int j = 0;j < 3;++j) + { + double sij_round = std::round(s(i, j)); +#ifdef __DEBUG + assert(near_equal(s(i, j), sij_round)); +#endif + s(i, j) = sij_round; + } + } + + // convert Matrix3 to matrix + ModuleBase::matrix h=s, b(3, 3, true); + b(0, 0)=1; + b(1, 1)=1; + b(2, 2)=1; + + // transform H into lower triangular form + auto max_min_index = [&h, this](int row, int &i1_to_max, int &i2_to_min) + { + if(fabs(h(row, i1_to_max)) < fabs(h(row, i2_to_min)) - epsilon) + { + int tmp = i2_to_min; + i2_to_min = i1_to_max; + i1_to_max = tmp; + } + return; + }; + + auto max_min_index_row1 = [&max_min_index, &h, this](int &imax, int &imin) + { + int imid=1; + imax=0; imin=2; + max_min_index(0, imid, imin); + max_min_index(0, imax, imid); + max_min_index(0, imid, imin); + if (equal(h(0, imin), 0)) + { + imin = imid; + } + else if (equal(h(0, imax), 0)) + { + imax = imid; + } + return; + }; + + auto swap_col = [&h, &b](int c1, int c2) + { + double tmp=0.0; + for(int r=0;r<3;++r) + { + tmp = h(r, c2); + h(r, c2)=h(r, c1); + h(r, c1)=tmp; + tmp = b(r, c2); + b(r, c2)=b(r, c1); + b(r, c1)=tmp; + } + return; + }; + + // row 1 + int imax=0; + int imin=0; + + while(int(equal(h(0, 0), 0)) + int(equal(h(0, 1), 0)) + int(equal(h(0, 2), 0)) < 2) + { + max_min_index_row1(imax, imin); + double f = floor((fabs(h(0, imax) )+ epsilon)/fabs(h(0, imin))); + if (h(0, imax) * h(0, imin) < -epsilon) + { + f *= -1; + } + for(int r=0;r<3;++r) + { + h(r, imax) -= f*h(r, imin); + b(r, imax) -= f*b(r, imin); + } + } + + if (equal(h(0, 0), 0)) + { + equal(h(0, 1), 0) ? swap_col(0, 2) : swap_col(0, 1); + } + + if (h(0, 0) < -epsilon) + { + for (int r = 0; r < 3; ++r) + { + h(r, 0) *= -1; + b(r, 0) *= -1; + } + } + + //row 2 + if (equal(h(1, 1), 0)) + { + swap_col(1, 2); + } + + while(!equal(h(1, 2), 0)) + { + imax=1, imin=2; + max_min_index(1, imax, imin); + double f = floor((fabs(h(1, imax) )+ epsilon)/fabs(h(1, imin))); + + if (h(1, imax) * h(1, imin) < -epsilon) + { + f *= -1; + } + + for(int r=0;r<3;++r) + { + h(r, imax) -= f*h(r, imin); + b(r, imax) -= f*b(r, imin); + } + + if (equal(h(1, 1), 0)) + { + swap_col(1, 2); + } + } + if (h(1, 1) < -epsilon) + { + for (int r = 0; r < 3; ++r) + { + h(r, 1) *= -1; + b(r, 1) *= -1; + } + } + + //row3 + if (h(2, 2) < -epsilon) + { + for (int r = 0; r < 3; ++r) + { + h(r, 2) *= -1; + b(r, 2) *= -1; + } + } + + // deal with off-diagonal elements + while (h(1, 0) > h(1, 1) - epsilon) + { + for(int r=0;r<3;++r) + { + h(r, 0) -= h(r, 1); + b(r, 0) -= b(r, 1); + } + } + while (h(1, 0) < -epsilon) + { + for(int r=0;r<3;++r) + { + h(r, 0) += h(r, 1); + b(r, 0) += b(r, 1); + } + } + for(int j=0;j<2;++j) + { + while (h(2, j) > h(2, 2) - epsilon) + { + for(int r=0;r<3;++r) + { + h(r, j) -= h(r, 2); + b(r, j) -= b(r, 2); + } + } + while (h(2, j) < -epsilon) + { + for(int r=0;r<3;++r) + { + h(r, j) += h(r, 2); + b(r, j) += b(r, 2); + } + } + } + + //convert matrix to Matrix3 + h3.e11=h(0, 0); h3.e12=h(0, 1); h3.e13=h(0, 2); + h3.e21=h(1, 0); h3.e22=h(1, 1); h3.e23=h(1, 2); + h3.e31=h(2, 0); h3.e32=h(2, 1); h3.e33=h(2, 2); + b3.e11=b(0, 0); b3.e12=b(0, 1); b3.e13=b(0, 2); + b3.e21=b(1, 0); b3.e22=b(1, 1); b3.e23=b(1, 2); + b3.e31=b(2, 0); b3.e32=b(2, 1); b3.e33=b(2, 2); + + //check s*b=h + ModuleBase::matrix check_zeros = s3.to_matrix() * b - h; +#ifdef __DEBUG + for (int i = 0;i < 3;++i) + { + for(int j=0;j<3;++j) + { + assert(near_equal(check_zeros(i, j), 0)); + } + } +#endif + return; +} diff --git a/source/module_cell/module_symmetry/symmetry.cpp b/source/module_cell/module_symmetry/symmetry.cpp index 0f4bfb41dc7..916faaa158e 100644 --- a/source/module_cell/module_symmetry/symmetry.cpp +++ b/source/module_cell/module_symmetry/symmetry.cpp @@ -566,7 +566,7 @@ void Symmetry::get_optlat(ModuleBase::Vector3 &v1, ModuleBase::Vector3 &v1, ModuleBase::Vector3 h(1, 1) - epsilon) - { - for(int r=0;r<3;++r) - { - h(r, 0) -= h(r, 1); - b(r, 0) -= b(r, 1); - } - } - while (h(1, 0) < -epsilon) - { - for(int r=0;r<3;++r) - { - h(r, 0) += h(r, 1); - b(r, 0) += b(r, 1); - } - } - for(int j=0;j<2;++j) - { - while (h(2, j) > h(2, 2) - epsilon) - { - for(int r=0;r<3;++r) - { - h(r, j) -= h(r, 2); - b(r, j) -= b(r, 2); - } - } - while (h(2, j) < -epsilon) - { - for(int r=0;r<3;++r) - { - h(r, j) += h(r, 2); - b(r, j) += b(r, 2); - } - } - } - - //convert matrix to Matrix3 - h3.e11=h(0, 0); h3.e12=h(0, 1); h3.e13=h(0, 2); - h3.e21=h(1, 0); h3.e22=h(1, 1); h3.e23=h(1, 2); - h3.e31=h(2, 0); h3.e32=h(2, 1); h3.e33=h(2, 2); - b3.e11=b(0, 0); b3.e12=b(0, 1); b3.e13=b(0, 2); - b3.e21=b(1, 0); b3.e22=b(1, 1); b3.e23=b(1, 2); - b3.e31=b(2, 0); b3.e32=b(2, 1); b3.e33=b(2, 2); - - //check s*b=h - ModuleBase::matrix check_zeros = s3.to_matrix() * b - h; -#ifdef __DEBUG - for (int i = 0;i < 3;++i) - { - for(int j=0;j<3;++j) - { - assert(near_equal(check_zeros(i, j), 0)); - } - } -#endif - return; -} - bool Symmetry::is_all_movable(const Atom* atoms, const Statistics& st)const { bool all_mbl = true; From 089bdb41778b4925ff7a00f866f8370319fcd6ea Mon Sep 17 00:00:00 2001 From: mohanchen Date: Tue, 25 Mar 2025 15:28:35 +0800 Subject: [PATCH 04/12] refactor symmetry.cpp --- .../module_symmetry/symm_analysis.cpp | 307 +++++++++++++++++ .../module_cell/module_symmetry/symmetry.cpp | 314 +----------------- 2 files changed, 312 insertions(+), 309 deletions(-) create mode 100644 source/module_cell/module_symmetry/symm_analysis.cpp diff --git a/source/module_cell/module_symmetry/symm_analysis.cpp b/source/module_cell/module_symmetry/symm_analysis.cpp new file mode 100644 index 00000000000..3c336d91b34 --- /dev/null +++ b/source/module_cell/module_symmetry/symm_analysis.cpp @@ -0,0 +1,307 @@ + +void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, std::ofstream& ofs_running) +{ + const double MAX_EPS = std::max(1e-3, epsilon_input * 1.001); + const double MULT_EPS = 2.0; + + ModuleBase::TITLE("Symmetry","analy_sys"); + ModuleBase::timer::tick("Symmetry","analy_sys"); + + ofs_running << "\n\n\n\n"; + ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; + ofs_running << " | |" << std::endl; + ofs_running << " | Performing symmetry analysis: |" << std::endl; + ofs_running << " | We calculate the norm of 3 vectors and the angles between them, |" << std::endl; + ofs_running << " | the type of Bravais lattice is given. We can judge if the unticell |" << std::endl; + ofs_running << " | is a primitive cell. Finally we give the point group operation for |" << std::endl; + ofs_running << " | this unitcell. We use the point group operations to perform |" << std::endl; + ofs_running << " | symmetry analysis on given k-point mesh and the charge density. |" << std::endl; + ofs_running << " | |" << std::endl; + ofs_running << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; + ofs_running << "\n\n\n\n"; + + // -------------------------------- + // 1. copy data and allocate memory + // -------------------------------- + // number of total atoms + this->nat = st.nat; + // number of atom species + this->ntype = st.ntype; + + assert(ntype>0); + + this->na = new int[ntype]; + this->istart = new int[ntype]; // start number of atom. + this->index = new int [nat + 2]; // index of atoms + + ModuleBase::GlobalFunc::ZEROS(na, ntype); + ModuleBase::GlobalFunc::ZEROS(istart, ntype); + ModuleBase::GlobalFunc::ZEROS(index, nat+2); + + // atom positions + // used in checksym. + newpos = new double[3*nat]; // positions of atoms before rotation + rotpos = new double[3*nat]; // positions of atoms after rotation + ModuleBase::GlobalFunc::ZEROS(newpos, 3*nat); + ModuleBase::GlobalFunc::ZEROS(rotpos, 3*nat); + + this->a1 = lat.a1; + this->a2 = lat.a2; + this->a3 = lat.a3; + + ModuleBase::Matrix3 latvec1; + latvec1.e11 = a1.x; latvec1.e12 = a1.y; latvec1.e13 = a1.z; + latvec1.e21 = a2.x; latvec1.e22 = a2.y; latvec1.e23 = a2.z; + latvec1.e31 = a3.x; latvec1.e32 = a3.y; latvec1.e33 = a3.z; + + output::printM3(ofs_running,"LATTICE VECTORS: (CARTESIAN COORDINATE: IN UNIT OF A0)",latvec1); + + istart[0] = 0; + this->itmin_type = 0; + this->itmin_start = 0; + for (int it = 0; it < ntype; ++it) + { + Atom* atom = &atoms[it]; + this->na[it] = atom->na; + if (it > 0) { + istart[it] = istart[it-1] + na[it-1]; + } + //std::cout << "\n istart = " << istart[it]; + if (na[it] < na[itmin_type]) + { + this->itmin_type = it; + this->itmin_start = istart[it]; + } + } + //s: input config + s1 = a1; + s2 = a2; + s3 = a3; + + + auto lattice_to_group = [&, this](int& nrot_out, int& nrotk_out, std::ofstream& ofs_running) -> void + { + // a: the optimized lattice vectors, output + // s: the input lattice vectors, input + // find the real_brav 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, atoms, 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); + + optlat.e11 = a1.x; optlat.e12 = a1.y; optlat.e13 = a1.z; + optlat.e21 = a2.x; optlat.e22 = a2.y; optlat.e23 = a2.z; + optlat.e31 = a3.x; optlat.e32 = a3.y; optlat.e33 = a3.z; + + // count the number of primitive cells in the supercell + this->pricell(this->newpos, atoms); + + test_brav = true; // output the real ibrav and point group + + // list all possible point group operations + this->setgroup(this->symop, this->nop, this->real_brav); + + // special case for AFM analysis + // which should be loop over all atoms, f.e only loop over spin-up atoms + // -------------------------------- + // AFM analysis Start + if (PARAM.inp.nspin > 1) + { + pricell_loop = this->magmom_same_check(atoms); + } + + if (!pricell_loop && PARAM.inp.nspin == 2) + { + this->analyze_magnetic_group(atoms, st, nrot_out, nrotk_out); + } + else + { + // get the real symmetry operations according to the input structure + // nrot_out: the number of pure point group rotations + // nrotk_out: the number of all space group operations + this->getgroup(nrot_out, nrotk_out, ofs_running, this->nop, this->symop, + this->gmatrix, this->gtrans, this->newpos, this->rotpos, this->index, + this->ntype, this->itmin_type, this->itmin_start, this->istart, this->na); + } + }; + + // -------------------------------- + // 2. analyze the symmetry + // -------------------------------- + // 2.1 skip the symmetry analysis if the symmetry has been analyzed + if (PARAM.inp.calculation == "cell-relax" && nrotk > 0) + { + std::ofstream no_out; // to screen the output when trying new epsilon + + // For the cases where cell-relax cause the number of symmetry operations to increase + if (this->nrotk > this->max_nrotk) { + this->max_nrotk = this->nrotk; + } + + int tmp_nrot, tmp_nrotk; + lattice_to_group(tmp_nrot, tmp_nrotk, ofs_running); // get the real symmetry operations + + // Actually, the analysis of symmetry has been done now + // Following implementation is find the best epsilon to keep the symmetry + // some different method to enlarge symmetry_prec + 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 + // Try to find the symmetry operations by increasing epsilon + while (tmp_nrotk < this->max_nrotk && epsilon < MAX_EPS) + { + eps_mult(MULT_EPS); + 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 (tmp_nrotk > this->nrotk) + { + 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_enlarged) + { + 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 (!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; + if (valid_index > 0) { // epsilon is set smaller + ofs_running << " Narrowing `symmetry_prec` from " << eps_current + << " to " << epsilon << " ..." << std::endl; + } + } + } else { + lattice_to_group(this->nrot, this->nrotk, ofs_running); + } + // Symmetry analysis End! + //------------------------------------------- + + // 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 + //---------------------------------- + // 3. output to running.log + //---------------------------------- + // output the point group + bool valid_group = this->pointgroup(this->nrot, this->pgnumber, this->pgname, this->gmatrix, ofs_running); + ModuleBase::GlobalFunc::OUT(ofs_running,"POINT GROUP", this->pgname); + // output the space group + valid_group = this->pointgroup(this->nrotk, this->spgnumber, this->spgname, this->gmatrix, ofs_running); + ModuleBase::GlobalFunc::OUT(ofs_running, "POINT GROUP IN SPACE GROUP", this->spgname); + + //----------------------------- + // 4. For the case where point group is not complete due to symmetry_prec + //----------------------------- + if (!valid_group) + { // select the operations that have the inverse + std::vectorinvmap(this->nrotk, -1); + this->gmatrix_invmap(this->gmatrix, this->nrotk, invmap.data()); + int nrotk_new = 0; + for (int isym = 0;isym < this->nrotk;++isym) + { + if (invmap[isym] != -1) + { + if(nrotk_new < isym) + { + this->gmatrix[nrotk_new] = this->gmatrix[isym]; + this->gtrans[nrotk_new] = this->gtrans[isym]; + } + ++nrotk_new; + } + } + this->nrotk = nrotk_new; + } + + // convert gmatrix to reciprocal space + this->gmatrix_convert_int(gmatrix, kgmatrix, nrotk, optlat, lat.G); + + // convert the symmetry operations from the basis of optimal symmetric configuration + // to the basis of input configuration + this->gmatrix_convert_int(gmatrix, gmatrix, nrotk, optlat, latvec1); + this->gtrans_convert(gtrans, gtrans, nrotk, optlat, latvec1); + + this->set_atom_map(atoms); // find the atom mapping according to the symmetry operations + + // Do this here for debug + if (PARAM.inp.calculation == "relax") + { + this->all_mbl = this->is_all_movable(atoms, st); + if (!this->all_mbl) + { + std::cout << "WARNING: Symmetry cannot be kept when not all atoms are movable.\n "; + std::cout << "Continue with symmetry=0 ... \n"; + ModuleSymmetry::Symmetry::symm_flag = 0; + } + } + + delete[] newpos; + delete[] na; + delete[] rotpos; + delete[] index; + delete[] istart; + ModuleBase::timer::tick("Symmetry","analy_sys"); + return; +} + diff --git a/source/module_cell/module_symmetry/symmetry.cpp b/source/module_cell/module_symmetry/symmetry.cpp index 916faaa158e..445ed526137 100644 --- a/source/module_cell/module_symmetry/symmetry.cpp +++ b/source/module_cell/module_symmetry/symmetry.cpp @@ -15,313 +15,6 @@ int Symmetry::symm_flag = 0; bool Symmetry::symm_autoclose = false; bool Symmetry::pricell_loop = true; -void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, std::ofstream& ofs_running) -{ - const double MAX_EPS = std::max(1e-3, epsilon_input * 1.001); - const double MULT_EPS = 2.0; - - ModuleBase::TITLE("Symmetry","analy_sys"); - ModuleBase::timer::tick("Symmetry","analy_sys"); - - ofs_running << "\n\n\n\n"; - ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; - ofs_running << " | |" << std::endl; - ofs_running << " | Performing symmetry analysis: |" << std::endl; - ofs_running << " | We calculate the norm of 3 vectors and the angles between them, |" << std::endl; - ofs_running << " | the type of Bravais lattice is given. We can judge if the unticell |" << std::endl; - ofs_running << " | is a primitive cell. Finally we give the point group operation for |" << std::endl; - ofs_running << " | this unitcell. We use the point group operations to perform |" << std::endl; - ofs_running << " | symmetry analysis on given k-point mesh and the charge density. |" << std::endl; - ofs_running << " | |" << std::endl; - ofs_running << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; - ofs_running << "\n\n\n\n"; - - // -------------------------------- - // 1. copy data and allocate memory - // -------------------------------- - // number of total atoms - this->nat = st.nat; - // number of atom species - this->ntype = st.ntype; - - assert(ntype>0); - - this->na = new int[ntype]; - this->istart = new int[ntype]; // start number of atom. - this->index = new int [nat + 2]; // index of atoms - - ModuleBase::GlobalFunc::ZEROS(na, ntype); - ModuleBase::GlobalFunc::ZEROS(istart, ntype); - ModuleBase::GlobalFunc::ZEROS(index, nat+2); - - // atom positions - // used in checksym. - newpos = new double[3*nat]; // positions of atoms before rotation - rotpos = new double[3*nat]; // positions of atoms after rotation - ModuleBase::GlobalFunc::ZEROS(newpos, 3*nat); - ModuleBase::GlobalFunc::ZEROS(rotpos, 3*nat); - - this->a1 = lat.a1; - this->a2 = lat.a2; - this->a3 = lat.a3; - - ModuleBase::Matrix3 latvec1; - latvec1.e11 = a1.x; latvec1.e12 = a1.y; latvec1.e13 = a1.z; - latvec1.e21 = a2.x; latvec1.e22 = a2.y; latvec1.e23 = a2.z; - latvec1.e31 = a3.x; latvec1.e32 = a3.y; latvec1.e33 = a3.z; - - output::printM3(ofs_running,"LATTICE VECTORS: (CARTESIAN COORDINATE: IN UNIT OF A0)",latvec1); - - istart[0] = 0; - this->itmin_type = 0; - this->itmin_start = 0; - for (int it = 0; it < ntype; ++it) - { - Atom* atom = &atoms[it]; - this->na[it] = atom->na; - if (it > 0) { - istart[it] = istart[it-1] + na[it-1]; - } - //std::cout << "\n istart = " << istart[it]; - if (na[it] < na[itmin_type]) - { - this->itmin_type = it; - this->itmin_start = istart[it]; - } - } - //s: input config - s1 = a1; - s2 = a2; - s3 = a3; - - - auto lattice_to_group = [&, this](int& nrot_out, int& nrotk_out, std::ofstream& ofs_running) -> void - { - // a: the optimized lattice vectors, output - // s: the input lattice vectors, input - // find the real_brav 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, atoms, 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); - - optlat.e11 = a1.x; optlat.e12 = a1.y; optlat.e13 = a1.z; - optlat.e21 = a2.x; optlat.e22 = a2.y; optlat.e23 = a2.z; - optlat.e31 = a3.x; optlat.e32 = a3.y; optlat.e33 = a3.z; - - // count the number of primitive cells in the supercell - this->pricell(this->newpos, atoms); - - test_brav = true; // output the real ibrav and point group - - // list all possible point group operations - this->setgroup(this->symop, this->nop, this->real_brav); - - // special case for AFM analysis - // which should be loop over all atoms, f.e only loop over spin-up atoms - // -------------------------------- - // AFM analysis Start - if (PARAM.inp.nspin > 1) - { - pricell_loop = this->magmom_same_check(atoms); - } - - if (!pricell_loop && PARAM.inp.nspin == 2) - { - this->analyze_magnetic_group(atoms, st, nrot_out, nrotk_out); - } - else - { - // get the real symmetry operations according to the input structure - // nrot_out: the number of pure point group rotations - // nrotk_out: the number of all space group operations - this->getgroup(nrot_out, nrotk_out, ofs_running, this->nop, this->symop, - this->gmatrix, this->gtrans, this->newpos, this->rotpos, this->index, - this->ntype, this->itmin_type, this->itmin_start, this->istart, this->na); - } - }; - - // -------------------------------- - // 2. analyze the symmetry - // -------------------------------- - // 2.1 skip the symmetry analysis if the symmetry has been analyzed - if (PARAM.inp.calculation == "cell-relax" && nrotk > 0) - { - std::ofstream no_out; // to screen the output when trying new epsilon - - // For the cases where cell-relax cause the number of symmetry operations to increase - if (this->nrotk > this->max_nrotk) { - this->max_nrotk = this->nrotk; - } - - int tmp_nrot, tmp_nrotk; - lattice_to_group(tmp_nrot, tmp_nrotk, ofs_running); // get the real symmetry operations - - // Actually, the analysis of symmetry has been done now - // Following implementation is find the best epsilon to keep the symmetry - // some different method to enlarge symmetry_prec - 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 - // Try to find the symmetry operations by increasing epsilon - while (tmp_nrotk < this->max_nrotk && epsilon < MAX_EPS) - { - eps_mult(MULT_EPS); - 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 (tmp_nrotk > this->nrotk) - { - 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_enlarged) - { - 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 (!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; - if (valid_index > 0) { // epsilon is set smaller - ofs_running << " Narrowing `symmetry_prec` from " << eps_current - << " to " << epsilon << " ..." << std::endl; - } - } - } else { - lattice_to_group(this->nrot, this->nrotk, ofs_running); - } - // Symmetry analysis End! - //------------------------------------------- - - // 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 - //---------------------------------- - // 3. output to running.log - //---------------------------------- - // output the point group - bool valid_group = this->pointgroup(this->nrot, this->pgnumber, this->pgname, this->gmatrix, ofs_running); - ModuleBase::GlobalFunc::OUT(ofs_running,"POINT GROUP", this->pgname); - // output the space group - valid_group = this->pointgroup(this->nrotk, this->spgnumber, this->spgname, this->gmatrix, ofs_running); - ModuleBase::GlobalFunc::OUT(ofs_running, "POINT GROUP IN SPACE GROUP", this->spgname); - - //----------------------------- - // 4. For the case where point group is not complete due to symmetry_prec - //----------------------------- - if (!valid_group) - { // select the operations that have the inverse - std::vectorinvmap(this->nrotk, -1); - this->gmatrix_invmap(this->gmatrix, this->nrotk, invmap.data()); - int nrotk_new = 0; - for (int isym = 0;isym < this->nrotk;++isym) - { - if (invmap[isym] != -1) - { - if(nrotk_new < isym) - { - this->gmatrix[nrotk_new] = this->gmatrix[isym]; - this->gtrans[nrotk_new] = this->gtrans[isym]; - } - ++nrotk_new; - } - } - this->nrotk = nrotk_new; - } - - // convert gmatrix to reciprocal space - this->gmatrix_convert_int(gmatrix, kgmatrix, nrotk, optlat, lat.G); - - // convert the symmetry operations from the basis of optimal symmetric configuration - // to the basis of input configuration - this->gmatrix_convert_int(gmatrix, gmatrix, nrotk, optlat, latvec1); - this->gtrans_convert(gtrans, gtrans, nrotk, optlat, latvec1); - - this->set_atom_map(atoms); // find the atom mapping according to the symmetry operations - - // Do this here for debug - if (PARAM.inp.calculation == "relax") - { - this->all_mbl = this->is_all_movable(atoms, st); - if (!this->all_mbl) - { - std::cout << "WARNING: Symmetry cannot be kept when not all atoms are movable.\n "; - std::cout << "Continue with symmetry=0 ... \n"; - ModuleSymmetry::Symmetry::symm_flag = 0; - } - } - - delete[] newpos; - delete[] na; - delete[] rotpos; - delete[] index; - delete[] istart; - ModuleBase::timer::tick("Symmetry","analy_sys"); - return; -} - - void Symmetry::set_atom_map(const Atom* atoms) { ModuleBase::TITLE("Symmetry", "set_atom_map"); @@ -356,8 +49,11 @@ void Symmetry::set_atom_map(const Atom* atoms) { for (int ia = istart[it]; ia < istart[it] + na[it]; ++ia) { - const int xx = ia * 3; const int yy = ia * 3 + 1; const int zz = ia * 3 + 2; - for (int k = 0;k < this->nrotk;++k) + const int xx = ia * 3; + const int yy = ia * 3 + 1; + const int zz = ia * 3 + 2; + + for (int k = 0;k < this->nrotk;++k) { rotpos[xx] = pos[xx] * gmatrix[k].e11 + pos[yy] * gmatrix[k].e21 From 0f09fb46f7d243142337d87747f3a6ff0fe5cbd0 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Tue, 25 Mar 2025 15:34:34 +0800 Subject: [PATCH 05/12] update Makefile.Objects --- source/Makefile.Objects | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 1f450083a79..19948c6a222 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -461,6 +461,14 @@ OBJS_SURCHEM=surchem.o\ OBJS_SYMMETRY=symm_other.o\ symmetry_basic.o\ + symm_analysis.o\ + symm_check.o\ + symm_getgroup.o\ + symm_hermite.o\ + symm_lattice.o\ + symm_magnetic.o\ + symm_pricell.o\ + symm_rho.o\ symmetry.o\ OBJS_XC=xc_functional.o\ From d6f04520b103ed67a9fd1eba26741cfb32d86dd6 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Tue, 25 Mar 2025 15:36:00 +0800 Subject: [PATCH 06/12] update cmakelists --- source/module_cell/module_symmetry/CMakeLists.txt | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/source/module_cell/module_symmetry/CMakeLists.txt b/source/module_cell/module_symmetry/CMakeLists.txt index 33b990a0648..8e1c81fb32e 100644 --- a/source/module_cell/module_symmetry/CMakeLists.txt +++ b/source/module_cell/module_symmetry/CMakeLists.txt @@ -3,6 +3,14 @@ add_library( OBJECT symm_other.cpp symmetry_basic.cpp + symm_analysis.cpp + symm_check.cpp + symm_getgroup.cpp + symm_hermite.cpp + symm_lattice.cpp + symm_magnetic.cpp + symm_pricell.cpp + symm_rho.cpp symmetry.cpp ) From 37e79f2b85c3daa38f32b6ae51b93024365f22d5 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Tue, 25 Mar 2025 15:50:44 +0800 Subject: [PATCH 07/12] update symmetry files --- source/module_cell/module_symmetry/symm_analysis.cpp | 3 +++ source/module_cell/module_symmetry/symm_check.cpp | 2 ++ source/module_cell/module_symmetry/symm_getgroup.cpp | 3 +++ source/module_cell/module_symmetry/symm_hermite.cpp | 2 ++ source/module_cell/module_symmetry/symm_lattice.cpp | 3 +++ source/module_cell/module_symmetry/symm_magnetic.cpp | 4 ++++ source/module_cell/module_symmetry/symm_pricell.cpp | 3 +++ source/module_cell/module_symmetry/symm_rho.cpp | 5 +++++ source/module_cell/module_symmetry/symmetry.cpp | 11 +++-------- source/module_cell/module_symmetry/symmetry.h | 6 ++++++ 10 files changed, 34 insertions(+), 8 deletions(-) diff --git a/source/module_cell/module_symmetry/symm_analysis.cpp b/source/module_cell/module_symmetry/symm_analysis.cpp index 3c336d91b34..be1af5e8230 100644 --- a/source/module_cell/module_symmetry/symm_analysis.cpp +++ b/source/module_cell/module_symmetry/symm_analysis.cpp @@ -1,3 +1,6 @@ +#include "symmetry.h" + +using namespace ModuleSymmetry; void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, std::ofstream& ofs_running) { diff --git a/source/module_cell/module_symmetry/symm_check.cpp b/source/module_cell/module_symmetry/symm_check.cpp index 5cb72e3d4cc..bd51cf9c546 100644 --- a/source/module_cell/module_symmetry/symm_check.cpp +++ b/source/module_cell/module_symmetry/symm_check.cpp @@ -1,3 +1,5 @@ +#include "symmetry.h" +using namespace ModuleSymmetry; bool Symmetry::checksym(const ModuleBase::Matrix3 &s, ModuleBase::Vector3& gtrans, diff --git a/source/module_cell/module_symmetry/symm_getgroup.cpp b/source/module_cell/module_symmetry/symm_getgroup.cpp index 977566af9d7..7c7085e6ed3 100644 --- a/source/module_cell/module_symmetry/symm_getgroup.cpp +++ b/source/module_cell/module_symmetry/symm_getgroup.cpp @@ -1,3 +1,6 @@ +#include "symmetry.h" +using namespace ModuleSymmetry; + void Symmetry::getgroup(int& nrot, int& nrotk, std::ofstream& ofs_running, const int& nop, const ModuleBase::Matrix3* symop, ModuleBase::Matrix3* gmatrix, ModuleBase::Vector3* gtrans, double* pos, double* rotpos, diff --git a/source/module_cell/module_symmetry/symm_hermite.cpp b/source/module_cell/module_symmetry/symm_hermite.cpp index f74192a70a1..9eef65405e0 100644 --- a/source/module_cell/module_symmetry/symm_hermite.cpp +++ b/source/module_cell/module_symmetry/symm_hermite.cpp @@ -1,3 +1,5 @@ +#include "symmetry.h" +using namespace ModuleSymmetry; void Symmetry::hermite_normal_form(const ModuleBase::Matrix3 &s3, ModuleBase::Matrix3 &h3, diff --git a/source/module_cell/module_symmetry/symm_lattice.cpp b/source/module_cell/module_symmetry/symm_lattice.cpp index bef6387b49d..2945ed0d94f 100644 --- a/source/module_cell/module_symmetry/symm_lattice.cpp +++ b/source/module_cell/module_symmetry/symm_lattice.cpp @@ -1,3 +1,6 @@ +#include "symmetry.h" +using namespace ModuleSymmetry; + //--------------------------------------------------- // The lattice will be transformed to a 'standard // cystallographic setting', the relation between diff --git a/source/module_cell/module_symmetry/symm_magnetic.cpp b/source/module_cell/module_symmetry/symm_magnetic.cpp index ca107c2bcf4..f28c1d963f3 100644 --- a/source/module_cell/module_symmetry/symm_magnetic.cpp +++ b/source/module_cell/module_symmetry/symm_magnetic.cpp @@ -1,3 +1,7 @@ +#include "symmetry.h" +using namespace ModuleSymmetry; + +#include void Symmetry::analyze_magnetic_group(const Atom* atoms, const Statistics& st, int& nrot_out, int& nrotk_out) { diff --git a/source/module_cell/module_symmetry/symm_pricell.cpp b/source/module_cell/module_symmetry/symm_pricell.cpp index 14f8f8d56cd..f70325dcab1 100644 --- a/source/module_cell/module_symmetry/symm_pricell.cpp +++ b/source/module_cell/module_symmetry/symm_pricell.cpp @@ -1,3 +1,6 @@ +#include "symmetry.h" +using namespace ModuleSymmetry; + void Symmetry::pricell(double* pos, const Atom* atoms) { bool no_diff = false; diff --git a/source/module_cell/module_symmetry/symm_rho.cpp b/source/module_cell/module_symmetry/symm_rho.cpp index fb7ea296e35..2fbc1de4987 100644 --- a/source/module_cell/module_symmetry/symm_rho.cpp +++ b/source/module_cell/module_symmetry/symm_rho.cpp @@ -1,3 +1,8 @@ +#include "symmetry.h" +using namespace ModuleSymmetry; + +#include "module_base/libm/libm.h" + void Symmetry::rho_symmetry( double *rho, const int &nr1, const int &nr2, const int &nr3) { diff --git a/source/module_cell/module_symmetry/symmetry.cpp b/source/module_cell/module_symmetry/symmetry.cpp index 445ed526137..d65bb4ffe4d 100644 --- a/source/module_cell/module_symmetry/symmetry.cpp +++ b/source/module_cell/module_symmetry/symmetry.cpp @@ -1,14 +1,7 @@ #include #include -#include -#include "module_parameter/parameter.h" #include "symmetry.h" -#include "module_parameter/parameter.h" -#include "module_base/libm/libm.h" -#include "module_base/mathzone.h" -#include "module_base/constants.h" -#include "module_base/timer.h" -#include "module_io/output.h" + namespace ModuleSymmetry { int Symmetry::symm_flag = 0; @@ -347,3 +340,5 @@ bool Symmetry::is_all_movable(const Atom* atoms, const Statistics& st)const } return all_mbl; } + +} diff --git a/source/module_cell/module_symmetry/symmetry.h b/source/module_cell/module_symmetry/symmetry.h index 0d9133edbec..bd83b017846 100644 --- a/source/module_cell/module_symmetry/symmetry.h +++ b/source/module_cell/module_symmetry/symmetry.h @@ -1,8 +1,14 @@ #ifndef SYMMETRY_H #define SYMMETRY_H + #include "module_cell/unitcell_data.h" #include "module_cell/atom_spec.h" +#include "module_base/timer.h" +#include "module_parameter/parameter.h" +#include "module_base/mathzone.h" +#include "module_base/constants.h" +#include "module_io/output.h" #include "symmetry_basic.h" namespace ModuleSymmetry From cea665bf6989574748305107d6d49d4a023a2c74 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Tue, 25 Mar 2025 15:55:54 +0800 Subject: [PATCH 08/12] update outputs --- source/module_cell/klist.cpp | 8 ++++---- source/module_cell/module_symmetry/symm_lattice.cpp | 2 +- source/module_cell/module_symmetry/symm_pricell.cpp | 2 +- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/source/module_cell/klist.cpp b/source/module_cell/klist.cpp index 7648ee3f5f4..52859f694dc 100644 --- a/source/module_cell/klist.cpp +++ b/source/module_cell/klist.cpp @@ -635,9 +635,9 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry& symm, recip_brav_name, ucell.atoms, false, - nullptr); - GlobalV::ofs_running << "(for reciprocal lattice: )" << std::endl; - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "BRAVAIS TYPE", recip_brav_type); + nullptr); + GlobalV::ofs_running << "\n For reciprocal-space lattice:" << std::endl; + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "BRAVAIS TYPE", recip_brav_type); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "BRAVAIS LATTICE NAME", recip_brav_name); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "ibrav", recip_brav_type); @@ -680,7 +680,7 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry& symm, ucell.atoms, false, nullptr); - GlobalV::ofs_running << "(for k-lattice: )" << std::endl; + GlobalV::ofs_running << "\n For k vectors:" << std::endl; ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "BRAVAIS TYPE", k_brav_type); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "BRAVAIS LATTICE NAME", k_brav_name); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "ibrav", k_brav_type); diff --git a/source/module_cell/module_symmetry/symm_lattice.cpp b/source/module_cell/module_symmetry/symm_lattice.cpp index 2945ed0d94f..680ed3eb764 100644 --- a/source/module_cell/module_symmetry/symm_lattice.cpp +++ b/source/module_cell/module_symmetry/symm_lattice.cpp @@ -303,7 +303,7 @@ void Symmetry::lattice_type( //---------------------------------------------- const bool right = Symm_Other::right_hand_sense(v1, v2, v3); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"right hand lattice",right); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"RIGHT HAND LATTICE",right); //------------------------------------------------- // (2) save and copy the original lattice vectors. diff --git a/source/module_cell/module_symmetry/symm_pricell.cpp b/source/module_cell/module_symmetry/symm_pricell.cpp index f70325dcab1..c4dfad8f06e 100644 --- a/source/module_cell/module_symmetry/symm_pricell.cpp +++ b/source/module_cell/module_symmetry/symm_pricell.cpp @@ -110,7 +110,7 @@ void Symmetry::pricell(double* pos, const Atom* atoms) int ntrans=ptrans.size(); if (ntrans <= 1) { - GlobalV::ofs_running<<" Original cell was already a primitive cell."<p1=this->a1; this->p2=this->a2; this->p3=this->a3; From 2a92999b0c2ce1bf3708ee9c8438f5df7ea1de5e Mon Sep 17 00:00:00 2001 From: mohanchen Date: Tue, 25 Mar 2025 16:02:01 +0800 Subject: [PATCH 09/12] update k point --- source/module_cell/klist.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/source/module_cell/klist.cpp b/source/module_cell/klist.cpp index 52859f694dc..0b36eceea64 100644 --- a/source/module_cell/klist.cpp +++ b/source/module_cell/klist.cpp @@ -989,7 +989,7 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry& symm, ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nkstot_ibz", nkstot_ibz); table.clear(); - table += "K-POINTS REDUCTION ACCORDING TO SYMMETRY\n"; + table += "\n K-POINTS REDUCTION ACCORDING TO SYMMETRY:\n"; table += FmtCore::format("%8s%12s%12s%12s%8s%8s\n", "IBZ", "DIRECT_X", "DIRECT_Y", "DIRECT_Z", "WEIGHT", "ibz2bz"); for (int ik = 0; ik < nkstot_ibz; ik++) { @@ -1094,7 +1094,7 @@ void K_Vectors::set_both_kvec(const ModuleBase::Matrix3& G, const ModuleBase::Ma kd_done = true; } std::string table; - table += "K-POINTS DIRECT COORDINATES\n"; + table += " K-POINTS DIRECT COORDINATES\n"; table += FmtCore::format("%8s%12s%12s%12s%8s\n", "KPOINTS", "DIRECT_X", "DIRECT_Y", "DIRECT_Z", "WEIGHT"); for (int i = 0; i < nkstot; i++) { @@ -1331,7 +1331,7 @@ void K_Vectors::print_klists(std::ofstream& ofs) ModuleBase::WARNING_QUIT("print_klists", "nkstot < nks"); } std::string table; - table += "K-POINTS CARTESIAN COORDINATES\n"; + table += " K-POINTS CARTESIAN COORDINATES\n"; table += FmtCore::format("%8s%12s%12s%12s%8s\n", "KPOINTS", "CARTESIAN_X", "CARTESIAN_Y", "CARTESIAN_Z", "WEIGHT"); for (int i = 0; i < nks; i++) { @@ -1345,7 +1345,7 @@ void K_Vectors::print_klists(std::ofstream& ofs) GlobalV::ofs_running << "\n" << table << std::endl; table.clear(); - table += "K-POINTS DIRECT COORDINATES\n"; + table += " K-POINTS DIRECT COORDINATES\n"; table += FmtCore::format("%8s%12s%12s%12s%8s\n", "KPOINTS", "DIRECT_X", "DIRECT_Y", "DIRECT_Z", "WEIGHT"); for (int i = 0; i < nks; i++) { From 19c9f72d34e516b1f56fc51287d5df501f797025 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Tue, 25 Mar 2025 16:13:11 +0800 Subject: [PATCH 10/12] update output of surfchem --- .../module_surchem/H_correction_pw.cpp | 60 +++++++++++-------- .../module_surchem/cal_vcav.cpp | 30 +++++----- .../module_surchem/minimize_cg.cpp | 11 +--- 3 files changed, 51 insertions(+), 50 deletions(-) diff --git a/source/module_hamilt_general/module_surchem/H_correction_pw.cpp b/source/module_hamilt_general/module_surchem/H_correction_pw.cpp index 3d97bc99952..b9f6a39548e 100644 --- a/source/module_hamilt_general/module_surchem/H_correction_pw.cpp +++ b/source/module_hamilt_general/module_surchem/H_correction_pw.cpp @@ -14,41 +14,49 @@ ModuleBase::matrix surchem::v_correction(const UnitCell& cell, const double* vlocal, Structure_Factor* sf) { - ModuleBase::TITLE("surchem", "v_correction"); - ModuleBase::timer::tick("surchem", "v_correction"); - - double* Porter = new double[rho_basis->nrxx]; - for (int i = 0; i < rho_basis->nrxx; i++) - Porter[i] = 0.0; + ModuleBase::TITLE("surchem", "v_cor"); + ModuleBase::timer::tick("surchem", "v_cor"); + + assert(rho_basis->nrxx>0); + + double* porter = new double[rho_basis->nrxx]; + for (int i = 0; i < rho_basis->nrxx; i++) + { + porter[i] = 0.0; + } const int nspin0 = (nspin == 2) ? 2 : 1; - for (int is = 0; is < nspin0; is++) - for (int ir = 0; ir < rho_basis->nrxx; ir++) - Porter[ir] += rho[is][ir]; + for (int is = 0; is < nspin0; is++) + { + for (int ir = 0; ir < rho_basis->nrxx; ir++) + { + porter[ir] += rho[is][ir]; + } + } - complex* Porter_g = new complex[rho_basis->npw]; - ModuleBase::GlobalFunc::ZEROS(Porter_g, rho_basis->npw); + complex* porter_g = new complex[rho_basis->npw]; + ModuleBase::GlobalFunc::ZEROS(porter_g, rho_basis->npw); - rho_basis->real2recip(Porter, Porter_g); + rho_basis->real2recip(porter, porter_g); - complex* N = new complex[rho_basis->npw]; - complex* TOTN = new complex[rho_basis->npw]; - complex* PS_TOTN = new complex[rho_basis->npw]; + complex* n = new complex[rho_basis->npw]; + complex* total_n = new complex[rho_basis->npw]; + complex* ps_totn = new complex[rho_basis->npw]; - cal_totn(cell, rho_basis, Porter_g, N, TOTN, vlocal); + cal_totn(cell, rho_basis, porter_g, n, total_n, vlocal); - cal_pseudo(cell, pgrid, rho_basis, Porter_g, PS_TOTN, sf); + cal_pseudo(cell, pgrid, rho_basis, porter_g, ps_totn, sf); ModuleBase::matrix v(nspin, rho_basis->nrxx); - v += cal_vel(cell, rho_basis, TOTN, PS_TOTN, nspin); - v += cal_vcav(cell, rho_basis, PS_TOTN, nspin); + v += cal_vel(cell, rho_basis, total_n, ps_totn, nspin); + v += cal_vcav(cell, rho_basis, ps_totn, nspin); - delete[] Porter; - delete[] Porter_g; - delete[] N; - delete[] PS_TOTN; - delete[] TOTN; + delete[] porter; + delete[] porter_g; + delete[] n; + delete[] ps_totn; + delete[] total_n; - ModuleBase::timer::tick("surchem", "v_correction"); + ModuleBase::timer::tick("surchem", "v_cor"); return v; -} \ No newline at end of file +} diff --git a/source/module_hamilt_general/module_surchem/cal_vcav.cpp b/source/module_hamilt_general/module_surchem/cal_vcav.cpp index 3f7e768f4b6..86d6923524e 100644 --- a/source/module_hamilt_general/module_surchem/cal_vcav.cpp +++ b/source/module_hamilt_general/module_surchem/cal_vcav.cpp @@ -39,12 +39,12 @@ void lapl_rho(const double& tpiba2, // calculates first derivative of the shape function in realspace // exp(-(log(n/n_c))^2 /(2 sigma^2)) /(sigma * sqrt(2*pi) )/n -void shape_gradn(const complex* PS_TOTN, const ModulePW::PW_Basis* rho_basis, double* eprime) +void shape_gradn(const complex* ps_totn, const ModulePW::PW_Basis* rho_basis, double* eprime) { - double *PS_TOTN_real = new double[rho_basis->nrxx]; - ModuleBase::GlobalFunc::ZEROS(PS_TOTN_real, rho_basis->nrxx); - rho_basis->recip2real(PS_TOTN, PS_TOTN_real); + double *ps_totn_real = new double[rho_basis->nrxx]; + ModuleBase::GlobalFunc::ZEROS(ps_totn_real, rho_basis->nrxx); + rho_basis->recip2real(ps_totn, ps_totn_real); double epr_c = 1.0 / sqrt(ModuleBase::TWO_PI) / PARAM.inp.sigma_k; double epr_z = 0; @@ -52,16 +52,16 @@ void shape_gradn(const complex* PS_TOTN, const ModulePW::PW_Basis* rho_b for (int ir = 0; ir < rho_basis->nrxx; ir++) { - epr_z = log(std::max(PS_TOTN_real[ir], min) / PARAM.inp.nc_k) / sqrt(2) / PARAM.inp.sigma_k; - eprime[ir] = epr_c * exp(-pow(epr_z, 2)) / std::max(PS_TOTN_real[ir], min); + epr_z = log(std::max(ps_totn_real[ir], min) / PARAM.inp.nc_k) / sqrt(2) / PARAM.inp.sigma_k; + eprime[ir] = epr_c * exp(-pow(epr_z, 2)) / std::max(ps_totn_real[ir], min); } - delete[] PS_TOTN_real; + delete[] ps_totn_real; } void surchem::createcavity(const UnitCell& ucell, const ModulePW::PW_Basis* rho_basis, - const complex* PS_TOTN, + const complex* ps_totn, double* vwork) { ModuleBase::Vector3 *nablan = new ModuleBase::Vector3[rho_basis->nrxx]; @@ -75,7 +75,7 @@ void surchem::createcavity(const UnitCell& ucell, ModuleBase::GlobalFunc::ZEROS(lapn, rho_basis->nrxx); // nabla n - XC_Functional::grad_rho(PS_TOTN, nablan, rho_basis, ucell.tpiba); + XC_Functional::grad_rho(ps_totn, nablan, rho_basis, ucell.tpiba); // |\nabla n |^2 = nablan_2 for (int ir = 0; ir < rho_basis->nrxx; ir++) @@ -84,7 +84,7 @@ void surchem::createcavity(const UnitCell& ucell, } // Laplacian of n - lapl_rho(ucell.tpiba2,PS_TOTN, lapn, rho_basis); + lapl_rho(ucell.tpiba2,ps_totn, lapn, rho_basis); //------------------------------------------------------------- // add -Lap(n)/|\nabla n| to vwork and copy \sqrt(|\nabla n|^2) @@ -105,7 +105,7 @@ void surchem::createcavity(const UnitCell& ucell, // gamma * A = exp(-(log(n/n_c))^2 /(2 sigma^2)) /(sigma * sqrt(2*pi) ) //------------------------------------------------------------- double *term1 = new double[rho_basis->nrxx]; - shape_gradn(PS_TOTN, rho_basis, term1); + shape_gradn(ps_totn, rho_basis, term1); //------------------------------------------------------------- // quantum surface area, integral of (gamma*A / n) * |\nabla n| @@ -164,7 +164,7 @@ void surchem::createcavity(const UnitCell& ucell, ModuleBase::matrix surchem::cal_vcav(const UnitCell& ucell, const ModulePW::PW_Basis* rho_basis, - complex* PS_TOTN, + complex* ps_totn, int nspin) { ModuleBase::TITLE("surchem", "cal_vcav"); @@ -173,9 +173,9 @@ ModuleBase::matrix surchem::cal_vcav(const UnitCell& ucell, double *tmp_Vcav = new double[rho_basis->nrxx]; ModuleBase::GlobalFunc::ZEROS(tmp_Vcav, rho_basis->nrxx); - createcavity(ucell, rho_basis, PS_TOTN, tmp_Vcav); + createcavity(ucell, rho_basis, ps_totn, tmp_Vcav); - ModuleBase::GlobalFunc::ZEROS(Vcav.c, nspin * rho_basis->nrxx); + ModuleBase::GlobalFunc::ZEROS(Vcav.c, nspin * rho_basis->nrxx); if (nspin == 4) { for (int ir = 0; ir < rho_basis->nrxx; ir++) @@ -197,4 +197,4 @@ ModuleBase::matrix surchem::cal_vcav(const UnitCell& ucell, delete[] tmp_Vcav; ModuleBase::timer::tick("surchem", "cal_vcav"); return Vcav; -} \ No newline at end of file +} diff --git a/source/module_hamilt_general/module_surchem/minimize_cg.cpp b/source/module_hamilt_general/module_surchem/minimize_cg.cpp index 8f88cff09eb..35a7cc6a322 100644 --- a/source/module_hamilt_general/module_surchem/minimize_cg.cpp +++ b/source/module_hamilt_general/module_surchem/minimize_cg.cpp @@ -171,13 +171,10 @@ void surchem::Leps2(const UnitCell& ucell, complex* phi_work, complex* lp) { - // cout<<"leps2!"< *grad_phi = new ModuleBase::Vector3[rho_basis->nrxx]; XC_Functional::grad_rho(phi, grad_phi, rho_basis, ucell.tpiba); - // for (int i = 0; i < 10; i++) { - // grad_phi[i].print(); - // } + for (int ir = 0; ir < rho_basis->nrxx; ir++) { grad_phi[ir].x *= epsilon[ir]; @@ -235,17 +232,13 @@ void surchem::Leps2(const UnitCell& ucell, lp_real[ir] += tmp_vector3[ir].z; } - // for (int i = 0; i < 10; i++) { - // cout << lp_real << [i] << endl; - // } rho_basis->real2recip(lp_real.data(), lp); - // cout<<"lp: "<().swap(lp_real); std::vector().swap(grad_grad_phi); + delete[] grad_grad_phi_G; delete[] tmp_vector3; } From 687661831cc3b91d2bc9070889cfefa246626a0d Mon Sep 17 00:00:00 2001 From: mohanchen Date: Tue, 25 Mar 2025 16:43:12 +0800 Subject: [PATCH 11/12] update documents about KPT.info --- docs/quick_start/output.md | 4 +-- source/module_cell/klist.cpp | 2 +- .../module_charge/charge_init.cpp | 28 ++++++------------- tests/integrate/107_PW_CHG_mismatch/INPUT | 15 +++++----- 4 files changed, 19 insertions(+), 30 deletions(-) diff --git a/docs/quick_start/output.md b/docs/quick_start/output.md index 6a6e150d8e3..d00a6647d6f 100644 --- a/docs/quick_start/output.md +++ b/docs/quick_start/output.md @@ -14,7 +14,7 @@ For a complete list of input parameters, please consult this [instruction](../ad `running_scf.log` contains information on nearly all function calls made during the execution of ABACUS. -## *KPT* +## *KPT.info* This file contains the information of all generated k-points, as well as the list of k-points actually used for calculations after considering symmetry. @@ -39,4 +39,4 @@ ABACUS generates a `.cif` format structure file based on the input file `STRU`, ## *warning.log* -The file contains all the warning messages generated during the ABACUS run. \ No newline at end of file +The file contains all the warning messages generated during the ABACUS run. diff --git a/source/module_cell/klist.cpp b/source/module_cell/klist.cpp index 0b36eceea64..3c771d3dc6d 100644 --- a/source/module_cell/klist.cpp +++ b/source/module_cell/klist.cpp @@ -134,7 +134,7 @@ void K_Vectors::set(const UnitCell& ucell, { // output kpoints file std::stringstream skpt; - skpt << PARAM.globalv.global_readin_dir << "kpoints"; + skpt << PARAM.globalv.global_out_dir << "KPT.info"; //mohan modified 20250325 std::ofstream ofkpt(skpt.str().c_str()); // clear kpoints ofkpt << skpt2 << skpt1; ofkpt.close(); diff --git a/source/module_elecstate/module_charge/charge_init.cpp b/source/module_elecstate/module_charge/charge_init.cpp index 930a0c2b7aa..6e671f887aa 100644 --- a/source/module_elecstate/module_charge/charge_init.cpp +++ b/source/module_elecstate/module_charge/charge_init.cpp @@ -40,7 +40,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout, bool read_kin_error = false; if (PARAM.inp.init_chg == "file" || PARAM.inp.init_chg == "auto") { - GlobalV::ofs_running << " try to read charge from file" << std::endl; + GlobalV::ofs_running << " Read electron density from file" << std::endl; // try to read charge from binary file first, which is the same as QE // liuyu 2023-12-05 @@ -48,7 +48,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout, binary << PARAM.globalv.global_readin_dir << PARAM.inp.suffix + "-CHARGE-DENSITY.restart"; if (ModuleIO::read_rhog(binary.str(), rhopw, rhog)) { - GlobalV::ofs_running << " Read in the charge density: " << binary.str() << std::endl; + GlobalV::ofs_running << " Read electron density from file: " << binary.str() << std::endl; for (int is = 0; is < PARAM.inp.nspin; ++is) { rhopw->recip2real(rhog[is], rho[is]); @@ -67,24 +67,24 @@ void Charge::init_rho(elecstate::efermi& eferm_iout, this->rho[is], ucell.nat)) { - GlobalV::ofs_running << " Read in the charge density: " << ssc.str() << std::endl; + GlobalV::ofs_running << " Read electron density from file: " << ssc.str() << std::endl; } else if (is > 0) // nspin=2 or 4 { if (is == 1) // failed at the second spin { - std::cout << "Incomplete charge density file!" << std::endl; + std::cout << "Incomplete electron density file!" << std::endl; read_error = true; break; } else if (is == 2) // read 2 files when nspin=4 { - GlobalV::ofs_running << " Didn't read in the charge density but would rearrange it later. " + GlobalV::ofs_running << " Didn't read in the electron density but would rearrange it later. " << std::endl; } else if (is == 3) // read 2 files when nspin=4 { - GlobalV::ofs_running << " rearrange charge density " << std::endl; + GlobalV::ofs_running << " rearrange electron density " << std::endl; for (int ir = 0; ir < this->rhopw->nrxx; ir++) { this->rho[3][ir] = this->rho[0][ir] - this->rho[1][ir]; @@ -229,7 +229,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout, this->rho[is], ucell.nat)) { - GlobalV::ofs_running << " Read in the charge density: " << ssc.str() << std::endl; + GlobalV::ofs_running << " Read in the electron density: " << ssc.str() << std::endl; } } } @@ -262,17 +262,6 @@ void Charge::set_rho_core(const UnitCell& ucell, ModuleBase::TITLE("Charge","set_rho_core"); ModuleBase::timer::tick("Charge","set_rho_core"); - // double eps = 1.e-10; - //---------------------------------------------------------- - // LOCAL VARIABLES : - // counter on mesh points - // counter on atomic types - // counter on g vectors - //---------------------------------------------------------- - // int ir = 0; - // int it = 0; - // int ig = 0; - bool bl = false; for (int it = 0; it Date: Tue, 25 Mar 2025 17:11:55 +0800 Subject: [PATCH 12/12] update symmetry files --- .../module_cell/module_symmetry/symm_analysis.cpp | 1 + .../module_cell/module_symmetry/symm_lattice.cpp | 2 ++ source/module_cell/module_symmetry/symm_rho.cpp | 1 + source/module_cell/module_symmetry/symmetry.h | 2 -- tests/integrate/107_PW_OBOD_MemSaver/INPUT | 15 +++++++-------- tests/integrate/107_PW_OB_outputbands/INPUT | 1 - 6 files changed, 11 insertions(+), 11 deletions(-) diff --git a/source/module_cell/module_symmetry/symm_analysis.cpp b/source/module_cell/module_symmetry/symm_analysis.cpp index be1af5e8230..43d200a4df3 100644 --- a/source/module_cell/module_symmetry/symm_analysis.cpp +++ b/source/module_cell/module_symmetry/symm_analysis.cpp @@ -1,4 +1,5 @@ #include "symmetry.h" +#include "module_parameter/parameter.h" using namespace ModuleSymmetry; diff --git a/source/module_cell/module_symmetry/symm_lattice.cpp b/source/module_cell/module_symmetry/symm_lattice.cpp index 680ed3eb764..7ab5388b5ed 100644 --- a/source/module_cell/module_symmetry/symm_lattice.cpp +++ b/source/module_cell/module_symmetry/symm_lattice.cpp @@ -1,6 +1,8 @@ #include "symmetry.h" using namespace ModuleSymmetry; +#include "module_parameter/parameter.h" + //--------------------------------------------------- // The lattice will be transformed to a 'standard // cystallographic setting', the relation between diff --git a/source/module_cell/module_symmetry/symm_rho.cpp b/source/module_cell/module_symmetry/symm_rho.cpp index 2fbc1de4987..858054df233 100644 --- a/source/module_cell/module_symmetry/symm_rho.cpp +++ b/source/module_cell/module_symmetry/symm_rho.cpp @@ -2,6 +2,7 @@ using namespace ModuleSymmetry; #include "module_base/libm/libm.h" +#include "module_parameter/parameter.h" void Symmetry::rho_symmetry( double *rho, const int &nr1, const int &nr2, const int &nr3) diff --git a/source/module_cell/module_symmetry/symmetry.h b/source/module_cell/module_symmetry/symmetry.h index bd83b017846..748ca49d58a 100644 --- a/source/module_cell/module_symmetry/symmetry.h +++ b/source/module_cell/module_symmetry/symmetry.h @@ -1,11 +1,9 @@ #ifndef SYMMETRY_H #define SYMMETRY_H - #include "module_cell/unitcell_data.h" #include "module_cell/atom_spec.h" #include "module_base/timer.h" -#include "module_parameter/parameter.h" #include "module_base/mathzone.h" #include "module_base/constants.h" #include "module_io/output.h" diff --git a/tests/integrate/107_PW_OBOD_MemSaver/INPUT b/tests/integrate/107_PW_OBOD_MemSaver/INPUT index 7779daee4c3..6588402957d 100644 --- a/tests/integrate/107_PW_OBOD_MemSaver/INPUT +++ b/tests/integrate/107_PW_OBOD_MemSaver/INPUT @@ -5,25 +5,24 @@ calculation nscf nbands 6 symmetry 0 -pseudo_dir ../../PP_ORB +pseudo_dir ../../PP_ORB read_file_dir ./ #Parameters (2.Iteration) ecutwfc 20 -scf_thr 1e-8 -scf_nmax 100 +scf_thr 1e-8 +scf_nmax 100 -#out_chg 1 init_chg file -out_band 1 -out_dos 1 +out_band 1 +out_dos 1 #Parameters (3.Basis) basis_type pw #Parameters (4.Smearing) -smearing_method gauss -smearing_sigma 0.002 +smearing_method gauss +smearing_sigma 0.002 #Parameters (5.Mixing) mixing_type plain diff --git a/tests/integrate/107_PW_OB_outputbands/INPUT b/tests/integrate/107_PW_OB_outputbands/INPUT index ded78ffae6e..d42e2e194c1 100644 --- a/tests/integrate/107_PW_OB_outputbands/INPUT +++ b/tests/integrate/107_PW_OB_outputbands/INPUT @@ -14,7 +14,6 @@ ecutwfc 20 scf_thr 1e-8 scf_nmax 100 -#out_chg 1 init_chg file init_wfc random out_band 1