Skip to content

Feature&Refactor: Structure symmetrization during relax #3247

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Qianruipku merged 5 commits into from
Dec 1, 2023
Merged

Feature&Refactor: Structure symmetrization during relax #3247

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
d9db0d4
put Symmetry into UnitCell
maki49 Nov 22, 2023
bc1cf7c
symmetrization for cell-relax
maki49 Nov 22, 2023
dee73a8
relax
maki49 Nov 22, 2023
59df6a1
fix a segfalt
maki49 Nov 23, 2023
1eb8181
Merge branch 'develop' into symrelax
Qianruipku Dec 1, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions source/module_cell/klist.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -621,7 +621,7 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s
ModuleBase::Vector3<double> gb01(b_optlat.e11, b_optlat.e12, b_optlat.e13);
ModuleBase::Vector3<double> gb02(b_optlat.e21, b_optlat.e22, b_optlat.e23);
ModuleBase::Vector3<double> gb03(b_optlat.e31, b_optlat.e32, b_optlat.e33);
symm.lattice_type(gb1, gb2, gb3, gb01, gb02, gb03, b_const, b0_const, bbrav, bbrav_name, ucell, false, nullptr);
symm.lattice_type(gb1, gb2, gb3, gb01, gb02, gb03, b_const, b0_const, bbrav, bbrav_name, ucell.atoms, false, nullptr);
GlobalV::ofs_running<<"(for reciprocal lattice: )"<<std::endl;
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"BRAVAIS TYPE", bbrav);
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"BRAVAIS LATTICE NAME", bbrav_name);
Expand All @@ -646,7 +646,7 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s
return;
}

symm.lattice_type(gk1, gk2, gk3, gk01, gk02, gk03, bk_const, bk0_const, bkbrav, bkbrav_name, ucell, false, nullptr);
symm.lattice_type(gk1, gk2, gk3, gk01, gk02, gk03, bk_const, bk0_const, bkbrav, bkbrav_name, ucell.atoms, false, nullptr);
GlobalV::ofs_running<<"(for k-lattice: )"<<std::endl;
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"BRAVAIS TYPE", bkbrav);
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"BRAVAIS LATTICE NAME", bkbrav_name);
Expand All @@ -656,7 +656,7 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s
ModuleBase::Matrix3 bsymop[48];
int bnop = 0;
// search again
symm.lattice_type(gb1, gb2, gb3, gb1, gb2, gb3, b_const, b0_const, bbrav, bbrav_name, ucell, false, nullptr);
symm.lattice_type(gb1, gb2, gb3, gb1, gb2, gb3, b_const, b0_const, bbrav, bbrav_name, ucell.atoms, false, nullptr);
ModuleBase::Matrix3 b_optlat_new(gb1.x, gb1.y, gb1.z, gb2.x, gb2.y, gb2.z, gb3.x, gb3.y, gb3.z);
symm.setgroup(bsymop, bnop, bbrav);
symm.gmatrix_convert(bsymop, bsymop, bnop, b_optlat_new, ucell.G);
Expand Down
2 changes: 1 addition & 1 deletion source/module_cell/klist.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
#include "module_base/global_function.h"
#include "module_base/global_variable.h"
#include "module_base/matrix3.h"
#include "module_cell/module_symmetry/symmetry.h"
#include "module_cell/unitcell.h"

class K_Vectors
{
Expand Down
2 changes: 1 addition & 1 deletion source/module_cell/module_symmetry/run_symmetry.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ void calculate()
"STRU",
ofs_running);
std::cout << "set up cell classic done." << std::endl;
symm.analy_sys(ucell, out, ofs_running);
symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running);
ofs_running.close();
// ooo.set_orb_tables();

Expand Down
140 changes: 69 additions & 71 deletions source/module_cell/module_symmetry/symmetry.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,23 +8,11 @@

namespace ModuleSymmetry
{
Symmetry::Symmetry()
{
this->epsilon = 1e-6;
this->tab = 12;
this->available = true;
}

Symmetry::~Symmetry()
{
}


int Symmetry::symm_flag = 0;
bool Symmetry::symm_autoclose = false;
bool Symmetry::pricell_loop = true;

void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
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;
Expand All @@ -47,9 +35,9 @@ void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)


// number of total atoms
this->nat = ucell.nat;
this->nat = st.nat;
// number of atom species
this->ntype = ucell.ntype;
this->ntype = st.ntype;
this->na = new int[ntype];
this->istart = new int[ntype];
this->index = new int [nat + 2];
Expand All @@ -64,9 +52,9 @@ void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
ModuleBase::GlobalFunc::ZEROS(newpos, 3*nat);
ModuleBase::GlobalFunc::ZEROS(rotpos, 3*nat);

this->a1 = ucell.a1;
this->a2 = ucell.a2;
this->a3 = ucell.a3;
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;
Expand All @@ -83,7 +71,7 @@ void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
this->itmin_start = 0;
for (int it = 0; it < ntype; ++it)
{
Atom* atom = &ucell.atoms[it];
Atom* atom = &atoms[it];
this->na[it] = atom->na;
if (it > 0) {
istart[it] = istart[it-1] + na[it-1];
Expand All @@ -104,7 +92,7 @@ void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
//a: optimized config
// find the lattice type accordiing to lattice vectors.
this->lattice_type(this->a1, this->a2, this->a3, this->s1, this->s2, this->s3,
this->cel_const, this->pre_const, this->real_brav, ilattname, ucell, true, this->newpos);
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);
Expand All @@ -118,9 +106,8 @@ void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
optlat.e21 = a2.x; optlat.e22 = a2.y; optlat.e23 = a2.z;
optlat.e31 = a3.x; optlat.e32 = a3.y; optlat.e33 = a3.z;

this->pricell(this->newpos); // pengfei Li 2018-05-14
//for( iat =0 ; iat < ucell.nat ; iat++)
// std::cout << " newpos_now = " << newpos[3*iat] << " " << newpos[3*iat+1] << " " << newpos[3*iat+2] << std::endl;
this->pricell(this->newpos, atoms); // pengfei Li 2018-05-14

test_brav = true; // output the real ibrav and point group
this->setgroup(this->symop, this->nop, this->real_brav);
this->getgroup(nrot_out, nrotk_out, ofs_running);
Expand Down Expand Up @@ -235,18 +222,20 @@ void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
}

//convert gmatrix to reciprocal space
this->gmatrix_convert_int(gmatrix, kgmatrix, nrotk, optlat, ucell.G);
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(ucell);
this->set_atom_map(atoms);

if (GlobalV::NSPIN > 1) pricell_loop = this->magmom_same_check(atoms);

if (GlobalV::NSPIN > 1) pricell_loop = this->magmom_same_check(ucell);
if (GlobalV::CALCULATION == "relax") this->all_mbl = this->is_all_movable(atoms, st);

delete[] newpos;
delete[] newpos;
delete[] na;
delete[] rotpos;
delete[] index;
Expand Down Expand Up @@ -546,9 +535,9 @@ void Symmetry::lattice_type(
double *cel_const,
double *pre_const,
int& real_brav,
std::string &bravname,
const UnitCell &ucell,
bool convert_atoms,
std::string& bravname,
const Atom* atoms,
bool convert_atoms,
double* newpos)const
{
ModuleBase::TITLE("Symmetry","lattice_type");
Expand Down Expand Up @@ -657,13 +646,13 @@ void Symmetry::lattice_type(
GlobalV::ofs_running <<" The lattice vectors have been changed (STRU_SIMPLE.cif)"<<std::endl;
GlobalV::ofs_running <<std::endl;
int at=0;
for(int it=0; it<ucell.ntype; ++it)
for (int it = 0; it < this->ntype; ++it)
{
for(int ia=0; ia<ucell.atoms[it].na; ++ia)
for (int ia = 0; ia < this->na[it]; ++ia)
{
ModuleBase::Mathzone::Cartesian_to_Direct(ucell.atoms[it].tau[ia].x,
ucell.atoms[it].tau[ia].y,
ucell.atoms[it].tau[ia].z,
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,
Expand Down Expand Up @@ -695,11 +684,11 @@ void Symmetry::lattice_type(
ofs << std::endl;
at =0;

for (int it=0; it<ucell.ntype; it++)
for (int it = 0; it < this->ntype; it++)
{
for (int ia=0; ia<ucell.atoms[it].na; ia++)
for (int ia = 0; ia < this->na[it]; ia++)
{
ofs << ucell.atoms[it].label
ofs << atoms[it].label
<< " " << newpos[3*at]
<< " " << newpos[3*at+1]
<< " " << newpos[3*at+2] << std::endl;
Expand All @@ -724,13 +713,13 @@ void Symmetry::lattice_type(
if(convert_atoms)
{
int at=0;
for(int it=0; it<ucell.ntype; ++it)
for (int it = 0; it < this->ntype; ++it)
{
for(int ia=0; ia<ucell.atoms[it].na; ++ia)
for (int ia = 0; ia < this->na[it]; ++ia)
{
ModuleBase::Mathzone::Cartesian_to_Direct(ucell.atoms[it].tau[ia].x,
ucell.atoms[it].tau[ia].y,
ucell.atoms[it].tau[ia].z,
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,
Expand Down Expand Up @@ -903,13 +892,8 @@ void Symmetry::checksym(ModuleBase::Matrix3 &s, ModuleBase::Vector3<double> &gtr
this->check_boundary(pos[j*3+1]);
this->check_boundary(pos[j*3+2]);
}
//for( int iat =0 ; iat < ucell.nat ; iat++)
//std::cout << " newpos_now1 = " << newpos[3*iat] << " " << newpos[3*iat+1] << " " << newpos[3*iat+2] << std::endl;

//order original atomic positions for current species
this->atom_ordering_new(pos + istart[it] * 3, na[it], index + istart[it]);
//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;

//Rotate atoms of current species
for (int j = istart[it]; j < istart[it] + na[it]; ++j)
Expand Down Expand Up @@ -1063,7 +1047,7 @@ void Symmetry::checksym(ModuleBase::Matrix3 &s, ModuleBase::Vector3<double> &gtr
return;
}

void Symmetry::pricell(double* pos)
void Symmetry::pricell(double* pos, const Atom* atoms)
{
bool no_diff = 0;
s_flag = 0;
Expand Down Expand Up @@ -1290,12 +1274,11 @@ void Symmetry::pricell(double* pos)
#endif

// get the optimized primitive cell
UnitCell tmp_ucell;
std::string pbravname;
ModuleBase::Vector3<double> 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,tmp_ucell, false, nullptr);
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;
Expand Down Expand Up @@ -1636,24 +1619,24 @@ for (int g_index = 0; g_index < group_index; g_index++)
ModuleBase::timer::tick("Symmetry","rhog_symmetry");
}

void Symmetry::set_atom_map(const UnitCell& ucell)
void Symmetry::set_atom_map(const Atom* atoms)
{
ModuleBase::TITLE("Symmetry", "set_atom_map");
if (this->isym_rotiat_.size() > 0) return;
if (this->isym_rotiat_.size() == this->nrotk) return;
this->isym_rotiat_.resize(this->nrotk);
for (int i = 0; i < this->nrotk; ++i)this->isym_rotiat_[i].resize(this->nat, -1);

double* pos = this->newpos;
double* rotpos = this->rotpos;
ModuleBase::GlobalFunc::ZEROS(pos, ucell.nat * 3);
ModuleBase::GlobalFunc::ZEROS(pos, this->nat * 3);
int iat = 0;
for (int it = 0; it < ucell.ntype; it++)
for (int it = 0; it < this->ntype; it++)
{
for (int ia = 0; ia < ucell.atoms[it].na; ia++)
for (int ia = 0; ia < this->na[it]; ia++)
{
pos[3 * iat] = ucell.atoms[it].taud[ia].x;
pos[3 * iat + 1] = ucell.atoms[it].taud[ia].y;
pos[3 * iat + 2] = ucell.atoms[it].taud[ia].z;
pos[3 * iat] = atoms[it].taud[ia].x;
pos[3 * iat + 1] = atoms[it].taud[ia].y;
pos[3 * iat + 2] = atoms[it].taud[ia].z;
for (int k = 0; k < 3; ++k)
{
this->check_translation(pos[iat * 3 + k], -floor(pos[iat * 3 + k]));
Expand All @@ -1662,7 +1645,7 @@ void Symmetry::set_atom_map(const UnitCell& ucell)
iat++;
}
}
for (int it = 0; it < ntype; it++)
for (int it = 0; it < this->ntype; it++)
{
for (int ia = istart[it]; ia < istart[it] + na[it]; ++ia)
{
Expand Down Expand Up @@ -1729,12 +1712,12 @@ void Symmetry::symmetrize_vec3_nat(double* v)const // pengfei 2016-12-20
return;
}

void Symmetry::symmetrize_mat3(ModuleBase::matrix& sigma, const UnitCell& ucell)const //zhengdy added 2017
void Symmetry::symmetrize_mat3(ModuleBase::matrix& sigma, const Lattice& lat)const //zhengdy added 2017
{
ModuleBase::matrix A = ucell.latvec.to_matrix();
ModuleBase::matrix AT = ucell.latvec.Transpose().to_matrix();
ModuleBase::matrix invA = ucell.GT.to_matrix();
ModuleBase::matrix invAT = ucell.G.to_matrix();
ModuleBase::matrix A = lat.latvec.to_matrix();
ModuleBase::matrix AT = lat.latvec.Transpose().to_matrix();
ModuleBase::matrix invA = lat.GT.to_matrix();
ModuleBase::matrix invAT = lat.G.to_matrix();
ModuleBase::matrix tot_sigma(3, 3, true);
sigma = A * sigma * AT;
for (int k = 0; k < nrotk; ++k)
Expand Down Expand Up @@ -2094,18 +2077,18 @@ void Symmetry::hermite_normal_form(const ModuleBase::Matrix3 &s3, ModuleBase::Ma
return;
}

bool Symmetry::magmom_same_check(const UnitCell& ucell) const
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 < ucell.atoms[it].na;++ia)
for (int ia = 1;ia < atoms[it].na;++ia)
{
if (!equal(ucell.atoms[it].m_loc_[ia].x, ucell.atoms[it].m_loc_[0].x) ||
!equal(ucell.atoms[it].m_loc_[ia].y, ucell.atoms[it].m_loc_[0].y) ||
!equal(ucell.atoms[it].m_loc_[ia].z, ucell.atoms[it].m_loc_[0].z))
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;
Expand All @@ -2115,4 +2098,19 @@ bool Symmetry::magmom_same_check(const UnitCell& ucell) const
return pricell_loop;
}

}
bool Symmetry::is_all_movable(const Atom* atoms, const Statistics& st)const
{
bool all_mbl = true;
for (int iat = 0;iat < st.nat;++iat)
{
int it = st.iat2it[iat];
int ia = st.iat2ia[iat];
if (!atoms[it].mbl[ia].x || !atoms[it].mbl[ia].y || !atoms[it].mbl[ia].z)
{
all_mbl = false;
break;
}
}
return all_mbl;
}
}
Loading