Skip to content

Commit

Permalink
update atom mapping in symmetry
Browse files Browse the repository at this point in the history
  • Loading branch information
@maki49
maki49 committed Aug 12, 2023
1 parent 8f3391b commit 4a89e0f
Show file tree
Hide file tree
Showing 6 changed files with 79 additions and 25 deletions.
34 changes: 25 additions & 9 deletions source/module_cell/module_symmetry/symmetry.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,10 @@ void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
this->index = new int [nat + 2];
ModuleBase::GlobalFunc::ZEROS(na, ntype);
ModuleBase::GlobalFunc::ZEROS(istart, ntype);
ModuleBase::GlobalFunc::ZEROS(index, nat+2);
ModuleBase::GlobalFunc::ZEROS(index, nat + 2);
#ifdef __EXX
this->invf1.resize(this->nat, 0);
#endif

// atom positions
// used in checksym.
Expand Down Expand Up @@ -852,6 +855,16 @@ void Symmetry::getgroup(int &nrot, int &nrotk, std::ofstream &ofs_running)
gtrans[i].z = 0;
}
}
#ifdef __EXX
// output the inverse map of each symmetry operation
ofs_running << "ISYM: new atom index after symmetry operation" << std::endl;
for (int isym = 0;isym < nrotk;++isym)
{
ofs_running << isym << ": ";
for (int iat = 0;iat < this->nat;++iat) ofs_running << isym_rotiat_iat[isym][iat] << " ";
ofs_running << std::endl;
}
#endif
return;
}

Expand All @@ -865,9 +878,6 @@ void Symmetry::checksym(ModuleBase::Matrix3 &s, ModuleBase::Vector3<double> &gtr
bool no_diff = 0;
ModuleBase::Vector3<double> trans(2.0, 2.0, 2.0);
s_flag = 0;
#ifdef __EXX
std::vector<int> f1(this->nat, 0); //ordering-map before rotatation
#endif
for (int it = 0; it < ntype; it++)
{
//------------------------------------
Expand All @@ -885,9 +895,6 @@ void Symmetry::checksym(ModuleBase::Matrix3 &s, ModuleBase::Vector3<double> &gtr

//order original atomic positions for current species
this->atom_ordering_new(pos + istart[it] * 3, na[it], index + istart[it]);
#ifdef __EXX
if (this->firstsort) for (int ia = istart[it]; ia < istart[it] + na[it]; ++ia) f1[ia] = index[ia]; // save the ordering-map before rotation
#endif
//for( int iat =0 ; iat < ucell.nat ; iat++)
//std::cout << " newpos_now2 = " << newpos[3*iat] << " " << newpos[3*iat+1] << " " << newpos[3*iat+2] << std::endl;

Expand Down Expand Up @@ -1044,8 +1051,12 @@ void Symmetry::checksym(ModuleBase::Matrix3 &s, ModuleBase::Vector3<double> &gtr
#ifdef __EXX
// now this->index stores the ordering-map after symmetry operation: f2
// we can calculate the atom-index-map of group operation: g=f2^{-1}f1
std::vector<int> invf2 = Symmetry::invmap(this->index, this->nat);
this->isym_iat_rotiat.push_back(Symmetry::mapmul(f1.data(), invf2.data(), this->nat));
// but here we calculate g^{-1}=f1^{-1}f2
for (int it = 0; it < ntype; it++)
for (int ia = istart[it]; ia < na[it] + istart[it]; ia++)
this->index[ia] += istart[it];
std::vector<int> f2 = Symmetry::invmap(this->index, this->nat);
this->isym_rotiat_iat.push_back(Symmetry::mapmul(f2.data(), invf1.data(), this->nat));
#endif
}
return;
Expand All @@ -1072,6 +1083,11 @@ void Symmetry::pricell(double* pos)

//order original atomic positions for current species
this->atom_ordering_new(pos + istart[it] * 3, na[it], index + istart[it]);
#ifdef __EXX
// save the ordering-map before rotation
// index is the structure after sort, which is the inverse map of sort
for (int ia = istart[it]; ia < istart[it] + na[it]; ++ia) invf1[ia] = istart[it] + index[ia];
#endif
//copy pos to rotpos
for (int j = istart[it]; j < istart[it] + na[it]; ++j)
{
Expand Down
7 changes: 4 additions & 3 deletions source/module_cell/module_symmetry/symmetry.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,9 +119,10 @@ class Symmetry : public Symmetry_Basic
void hermite_normal_form(const ModuleBase::Matrix3& s, ModuleBase::Matrix3& H, ModuleBase::Matrix3& b) const;

#ifdef __EXX
bool firstsort = true;
/// @brief size: atom index map corresponding to each symmetry operation. size: [nrotk][nat]
std::vector<std::vector<int>> isym_iat_rotiat;
/// inverse ordering-map before rotatation
std::vector<int> invf1;
/// size: atom index map corresponding to each symmetry operation. size: [nrotk][nat]
std::vector<std::vector<int>> isym_rotiat_iat;
static std::vector<int> invmap(const int* map, const int& size)
{
std::vector<int> invf(size);
Expand Down
16 changes: 11 additions & 5 deletions source/module_cell/module_symmetry/symmetry_basic.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1059,15 +1059,21 @@ void Symmetry_Basic::atom_ordering_new(double *posi, const int natom, int *subin
while(ix_right<natom && equal(tmpx[ix_right],tmpx[i])) ++ix_right;
int nxequal=ix_right-i;
if(nxequal>1) //need a new sort
{
subindex[0] = 0;
{
std::vector<int>tmpindex(nxequal, 0);
for(int j=0; j<nxequal; ++j)
{
weighted_func[j]=1/epsilon*tmpy[i+j]+tmpz[i+j];
}
ModuleBase::heapsort(nxequal, weighted_func, subindex);
this->order_atoms(&posi[i*3], nxequal, subindex);
}
ModuleBase::heapsort(nxequal, weighted_func, tmpindex.data());
this->order_atoms(&posi[i * 3], nxequal, tmpindex.data());
//rearange subindex using tmpindex
for (int j = 0; j < nxequal; ++j)
{
tmpindex[j] = subindex[i + tmpindex[j]];
subindex[i + j] = tmpindex[j];
}
}
i=ix_right;
}

Expand Down
27 changes: 25 additions & 2 deletions source/module_ri/exx_symmetry.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ namespace ExxSym
const std::vector<std::complex<double>> sloc_ikibz,
const int nbasis,
const Parallel_2D& p2d,
const std::vector<std::vector<int>>& isym_iat_rotiat,
const std::vector<std::vector<int>>& isym_rotiat_iat,
const std::map<int, ModuleBase::Vector3<double>>& kstar_ibz,
const UnitCell& ucell,
const bool col_inside)
Expand All @@ -32,7 +32,7 @@ namespace ExxSym
{
int isym = isym_kvecd.first;
//get invmap : iat1 to iat0
std::vector<int> rotiat_iat = ModuleSymmetry::Symmetry::invmap(isym_iat_rotiat[isym].data(), ucell.nat);
std::vector<int> rotiat_iat = ModuleSymmetry::Symmetry::invmap(isym_rotiat_iat[isym].data(), ucell.nat);
// 1 symmetry operation - may more than one kvec_d ?? check it !!!
// ModuleBase::Vector3<double> kvds = isym_kvecd.second;
// set rearanged sloc_ik
Expand Down Expand Up @@ -107,6 +107,29 @@ namespace ExxSym
return psi_full;
}

psi::Psi<std::complex<double>, psi::DEVICE_CPU> restore_psik(
const int& nkstot_full,
const psi::Psi<std::complex<double>, psi::DEVICE_CPU>& psi_ibz,
const std::vector<std::vector<std::vector<std::complex<double>>>>& invSkrot_Sgk,
const int& nbasis,
const int& nbands,
const Parallel_Orbitals& pv)
{
ModuleBase::TITLE("ExxSym", "restore_psik");
psi::Psi<std::complex<double>, psi::DEVICE_CPU> psi_full(nkstot_full, pv.ncol_bands, pv.get_row_size());
int ikfull_start = 0;
for (int ik_ibz = 0;ik_ibz < psi_ibz.get_nk();++ik_ibz)
{
#ifdef __MPI
restore_psik_scalapack(ik_ibz, ikfull_start, psi_ibz, invSkrot_Sgk[ik_ibz], nbasis, nbands, pv, &psi_full);
#else
restore_psik_lapack(ik_ibz, ikfull_start, psi_ibz, invSkrot_Sgk[ik_ibz], nbasis, nbands, &psi_full);
#endif
ikfull_start += invSkrot_Sgk[ik_ibz].size();
}
return psi_full;
}

void restore_psik_lapack(
const int& ikibz,
const int& ikfull_start,
Expand Down
12 changes: 10 additions & 2 deletions source/module_ri/exx_symmetry.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ namespace ExxSym
/// @param sloc_ikibz [in] local overlap matrices of current ibz-kpoint: S(gk)
/// @param nbasis [in] global number of basis
/// @param p2d [in]2d parallel info
/// @param isym_iat_rotiat [in] atom index map corresponding to each symmetry operation
/// @param isym_rotiat_iat [in] inversion of atom index map corresponding to each symmetry operation
/// @param kstar_ibz [in] symmetry-equal k points to current ibz-kpont: [isym][kvec_d]
/// @param ucell [in] unitcell
/// @param col_inside [in] whether the matrix is column-major (major means memory continuity)
Expand All @@ -20,7 +20,7 @@ namespace ExxSym
const std::vector<std::complex<double>> sloc_ikibz,
const int nbasis,
const Parallel_2D& p2d,
const std::vector<std::vector<int>>& isym_iat_rotiat,
const std::vector<std::vector<int>>& isym_rotiat_iat,
const std::map<int, ModuleBase::Vector3<double>>& kstar_ibz,
const UnitCell& ucell,
const bool col_inside);
Expand All @@ -42,6 +42,14 @@ namespace ExxSym
const int& nbasis,
const int& nbands,
const Parallel_Orbitals& pv);
///@param invSkrot_Sgk [in] $\tilde{S}^{-1}(k)S(gk)$ for all the ibz-kpoints gk
psi::Psi<std::complex<double>, psi::DEVICE_CPU> restore_psik(
const int& nkstot_full,
const psi::Psi<std::complex<double>, psi::DEVICE_CPU>& psi_ibz,
const std::vector<std::vector<std::vector<std::complex<double>>>>& invSkrot_Sgk,
const int& nbasis,
const int& nbands,
const Parallel_Orbitals& pv);

#ifdef __MPI
/// @brief calculate $\tilde{S}^{-1}(k)S(gk)$ for one ibz-kpoint
Expand Down
8 changes: 4 additions & 4 deletions source/module_ri/test/exx_symmetry_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -171,7 +171,7 @@ TEST_F(SymExxTest, cal_Sk_rot)
{
// case
std::vector<std::pair<int, int>> atoms = { {2, 1}, {3, 3} };
std::vector<std::vector<int>> isym_iat_rotiat = { {0, 1, 2, 3, 4}, {0,1,3,4,2}, {0,1,4,2,3}, {1,0,2,3,4}, {1,0,3,4,2}, {1,0,4,2,3} };
std::vector<std::vector<int>> isym_rotiat_iat = { {0, 1, 2, 3, 4}, {0,1,3,4,2}, {0,1,4,2,3}, {1,0,2,3,4}, {1,0,3,4,2}, {1,0,4,2,3} };
ModuleBase::Vector3<double> kvd(0, 0, 0); //only one ibzkpt and its coordinate is important in this function
std::map<int, ModuleBase::Vector3<double>> kstar_ibz = { {0, kvd}, {1, kvd}, {2, kvd}, {3, kvd}, {4, kvd}, {5, kvd} };

Expand All @@ -190,17 +190,17 @@ TEST_F(SymExxTest, cal_Sk_rot)
this->copy_from_global(sfull_gk.data(), sloc_gk.data(), nbasis, nbasis, pv.get_row_size(), pv.get_col_size(), false, false);

// run (row-major)
std::vector<std::vector<std::complex<double>>> sloc_ks = ExxSym::cal_Sk_rot(sloc_gk, nbasis, pv, isym_iat_rotiat, kstar_ibz, ucell, false);
std::vector<std::vector<std::complex<double>>> sloc_ks = ExxSym::cal_Sk_rot(sloc_gk, nbasis, pv, isym_rotiat_iat, kstar_ibz, ucell, false);
// check
for (int isym = 0;isym < isym_iat_rotiat.size();++isym)
for (int isym = 0;isym < isym_rotiat_iat.size();++isym)
{
std::vector<std::complex<double>> sfull_ik = ExxSym::get_full_smat(sloc_ks[isym], nbasis, pv, false);
for (int i = 0;i < pv.get_row_size();i++)
{
int iwt0 = pv.local2global_row(i);
int iat0 = ucell.iwt2iat[iwt0];
int iw = ucell.iwt2iw[iwt0];
int iat1 = isym_iat_rotiat[isym][iat0];
int iat1 = isym_rotiat_iat[isym][iat0];
int iwt1 = ucell.iat2iwt[iat1] + iw;
EXPECT_EQ(sfull_gk[iwt0 * nbasis + iwt0], sfull_ik[iwt1 * nbasis + iwt0]);
}
Expand Down

0 comments on commit 4a89e0f

Please sign in to comment.