diff --git a/AGENTS.md b/AGENTS.md index 8f157fecfdb..ea0fcd28fc7 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -115,11 +115,6 @@ rules. Read the complete governance document before making or reviewing changes: - Do not relax existing tests or references merely to make a failure pass. Update references only when the intended behavior changed and the PR explains why. -- When mocking `UnitCell` in a test fixture, do not `delete[] iat2it` or - `iat2ia` in `TearDown`: they are owned by `UnitCell`'s internal `Statistics` - member, whose destructor releases them. Deleting them again causes a double - free. Mirror the ownership pattern of existing fixtures such as - `source/source_lcao/module_dftu/test/dftu_lcao_test.cpp`. ## Review And Exception Flow diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 60c74cad612..9238934d249 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -632,7 +632,7 @@ OBJS_SYMMETRY=symm_other.o\ symm_rho.o\ little_group.o\ symmetry.o\ - irreducible_sector.o\ + irred_sec.o\ irred_sec_bvk.o\ symm_rotation_k.o\ diff --git a/source/source_cell/cal_wfc.cpp b/source/source_cell/cal_wfc.cpp index b6307e5fc9e..5044ea33a47 100644 --- a/source/source_cell/cal_wfc.cpp +++ b/source/source_cell/cal_wfc.cpp @@ -53,10 +53,8 @@ namespace unitcell // This assertion ensures consistency between the two calculation paths. assert(nlocal_tmp > 0); assert(nlocal_tmp == nlocal); - delete[] ucell.iwt2iat; - delete[] ucell.iwt2iw; - ucell.iwt2iat = new int[nlocal_tmp]; - ucell.iwt2iw = new int[nlocal_tmp]; + ucell.iwt2iat.resize(nlocal_tmp); + ucell.iwt2iw.resize(nlocal_tmp); ucell.itia2iat.create(ntype, ucell.namax); ucell.set_iat2iwt(npol); diff --git a/source/source_cell/module_symmetry/CMakeLists.txt b/source/source_cell/module_symmetry/CMakeLists.txt index 5b2a03cd967..cd6eec7056c 100644 --- a/source/source_cell/module_symmetry/CMakeLists.txt +++ b/source/source_cell/module_symmetry/CMakeLists.txt @@ -17,7 +17,7 @@ add_library( symmetry.cpp symm_rot_spin.cpp little_group.cpp - irreducible_sector.cpp + irred_sec.cpp irred_sec_bvk.cpp symm_rotation_k.cpp ) diff --git a/source/source_cell/module_symmetry/irreducible_sector.cpp b/source/source_cell/module_symmetry/irred_sec.cpp similarity index 85% rename from source/source_cell/module_symmetry/irreducible_sector.cpp rename to source/source_cell/module_symmetry/irred_sec.cpp index 2b9ca66d7bf..b719a53c33a 100644 --- a/source/source_cell/module_symmetry/irreducible_sector.cpp +++ b/source/source_cell/module_symmetry/irred_sec.cpp @@ -1,4 +1,11 @@ -#include "source_cell/module_symmetry/irreducible_sector.h" +#include "source_cell/module_symmetry/irred_sec.h" +#include "source_base/global_variable.h" + +#include +#include +#include +#include + namespace ModuleSymmetry { // Raw-index dispatch shared by the real-space sector helpers, matching the convention used @@ -16,7 +23,7 @@ namespace ModuleSymmetry : symm.get_rotated_atom_anti(isym - symm.nrotk, iat); } - TC Irreducible_Sector::rotate_R(const Symmetry& symm, + TC Irred_Sector::rotate_R(const Symmetry& symm, const int isym, const int iat1, const int iat2, const TC& R, const char gauge) const { auto round2int = [symm](const double x) -> int { return x > 0 ? static_cast(x + symm.epsilon) : static_cast(x - symm.epsilon); }; @@ -28,7 +35,7 @@ namespace ModuleSymmetry : R_double * gmat + this->return_lattice_[iat2][isym] - this->return_lattice_[iat1][isym]; return { round2int(Rrot_double.x), round2int(Rrot_double.y), round2int(Rrot_double.z) }; } - TapR Irreducible_Sector::rotate_apR_by_formula(const Symmetry& symm, + TapR Irred_Sector::rotate_apR_by_formula(const Symmetry& symm, const int isym, const TapR& apR, const char gauge) const { const Tap& aprot = { sector_rotated_atom(symm, isym, apR.first.first), @@ -36,19 +43,6 @@ namespace ModuleSymmetry return { aprot, this->rotate_R(symm, isym, apR.first.first, apR.first.second, apR.second, gauge) }; } - TCdouble Irreducible_Sector::get_aRb_direct(const Atom* atoms, const Statistics& st, - const int iat1, const int iat2, const TCdouble& R, const char gauge)const - { - return TCdouble(atoms[st.iat2it[iat1]].taud[st.iat2ia[iat1]] - atoms[st.iat2it[iat2]].taud[st.iat2ia[iat2]]) + (gauge == 'L' ? R : TCdouble(-R)); - } - - TCdouble Irreducible_Sector::get_aRb_direct(const Atom* atoms, const Statistics& st, - const int iat1, const int iat2, const TC& R, const char gauge) const - { - const TCdouble R_double(static_cast(R[0]), static_cast(R[1]), static_cast(R[2])); - return get_aRb_direct(atoms, st, iat1, iat2, R_double); - } - inline void output_return_lattice(const std::vector>& return_lattice) { std::cout << "return lattice:" << std::endl; @@ -64,7 +58,7 @@ namespace ModuleSymmetry // Perfoming {R|t} to atom position r in the R=0 lattice, we get Rr+t, which may get out of R=0 lattice, // whose image in R=0 lattice is r'=Rr+t-O. This function is to get O for each atom and each symmetry operation. // the range of direct position is [-0.5, 0.5). - TCdouble Irreducible_Sector::get_return_lattice(const Symmetry& symm, + TCdouble Irred_Sector::get_return_lattice(const Symmetry& symm, const ModuleBase::Matrix3& gmatd, const TCdouble gtransd, const TCdouble& posd_a1, const TCdouble& posd_a2)const { @@ -84,27 +78,28 @@ namespace ModuleSymmetry return TCdouble(std::round(return_lattice_double.x), std::round(return_lattice_double.y), std::round(return_lattice_double.z)); } - void Irreducible_Sector::cal_return_lattice_all(const Symmetry& symm, const Atom* atoms, const Statistics& st) + void Irred_Sector::cal_return_lattice_all(const Symmetry& symm, const Atom* atoms, + const int nat, const std::vector& iat2it, const std::vector& iat2ia) { - ModuleBase::TITLE("Symmetry_rotation", "cal_return_lattice_all"); + ModuleBase::TITLE("Irred_Sector", "cal_return_lattice_all"); // Columns [0, nrotk) are the unitary operations; columns [nrotk, nrotk+nrotk_anti) are the // spatial parts of the antiunitary elements Theta*g of the Shubnikov group (nspin=4 magnetic), // so that Symmetry_rotation can address both with one raw index. - this->return_lattice_.resize(st.nat, std::vector(symm.nrotk + symm.nrotk_anti)); - for (int iat1 = 0;iat1 < st.nat;++iat1) + this->return_lattice_.resize(nat, std::vector(symm.nrotk + symm.nrotk_anti)); + for (int iat1 = 0;iat1 < nat;++iat1) { - int it = st.iat2it[iat1]; - int ia1 = st.iat2ia[iat1]; + int it = iat2it[iat1]; + int ia1 = iat2ia[iat1]; for (int isym = 0;isym < symm.nrotk;++isym) { int iat2 = symm.get_rotated_atom(isym, iat1); - int ia2 = st.iat2ia[iat2]; + int ia2 = iat2ia[iat2]; this->return_lattice_[iat1][isym] = get_return_lattice(symm, symm.gmatrix[isym], symm.gtrans[isym], atoms[it].taud[ia1], atoms[it].taud[ia2]); } for (int j = 0;j < symm.nrotk_anti;++j) { int iat2 = symm.get_rotated_atom_anti(j, iat1); - int ia2 = st.iat2ia[iat2]; + int ia2 = iat2ia[iat2]; this->return_lattice_[iat1][symm.nrotk + j] = get_return_lattice(symm, symm.gmatrix_anti[j], symm.gtrans_anti[j], atoms[it].taud[ia1], atoms[it].taud[ia2]); } } @@ -112,7 +107,7 @@ namespace ModuleSymmetry // output_return_lattice(this->return_lattice_); } - void Irreducible_Sector::output_full_map_to_irreducible_sector(const int nat) + void Irred_Sector::output_full_map_to_irreducible_sector(const int nat) { std::cout << "Final map to irreducible sector: " << std::endl; for (auto& apR_isym_irapR : this->full_map_to_irreducible_sector_) @@ -127,7 +122,7 @@ namespace ModuleSymmetry } } - void Irreducible_Sector::output_sector_star() + void Irred_Sector::output_sector_star() { std::cout << "Found " << this->sector_stars_.size() << " irreducible sector stars:" << std::endl; // for (auto& irs_star : this->sector_stars_) @@ -150,7 +145,7 @@ namespace ModuleSymmetry std::cout << std::endl; } } - void Irreducible_Sector::write_irreducible_sector(const std::string& output_dir) + void Irred_Sector::write_irreducible_sector(const std::string& output_dir) { if(GlobalV::MY_RANK == 0 && !output_dir.empty()) { @@ -164,20 +159,23 @@ namespace ModuleSymmetry } } - void Irreducible_Sector::find_irreducible_sector(const Symmetry& symm, const Atom* atoms, const Statistics& st, const std::vector& Rs, const TC& period, const Lattice& lat, const std::string& output_dir) + void Irred_Sector::find_irred_sector(const Symmetry& symm, const Atom* atoms, + const int nat, const int ntype, + const std::vector& iat2it, const std::vector& iat2ia, + const std::vector& Rs, const TC& period, const Lattice& lat, const std::string& output_dir) { this->full_map_to_irreducible_sector_.clear(); this->irreducible_sector_.clear(); this->sector_stars_.clear(); - if (this->return_lattice_.empty()) this->cal_return_lattice_all(symm, atoms, st); + if (this->return_lattice_.empty()) this->cal_return_lattice_all(symm, atoms, nat, iat2it, iat2ia); // if (this->atompair_stars_.empty()) this->find_irreducible_atom_pairs(symm); // contruct {atom pair, R} set // constider different number of Rs for different atom pairs later. std::map> apR_all; - for (int iat1 = 0;iat1 < st.nat; iat1++) - for (int iat2 = 0; iat2 < st.nat; iat2++) + for (int iat1 = 0;iat1 < nat; iat1++) + for (int iat2 = 0; iat2 < nat; iat2++) for (auto& R : Rs) apR_all[{iat1, iat2}].insert(R); @@ -200,7 +198,7 @@ namespace ModuleSymmetry // get symmetry of BvK supercell if (this->isymbvk_to_isym_.empty()) - this->gen_symmetry_BvK(symm, atoms, lat, st, period); + this->gen_symm_bvk(symm, atoms, lat, nat, ntype, period); assert(!this->isymbvk_to_isym_.empty()); // std::vector in_2d_plain; // const bool judge_2d = (symm.real_brav == 4); @@ -255,7 +253,7 @@ namespace ModuleSymmetry for (auto& sector : this->sector_stars_) total_apR_in_star += sector.second.size(); assert(total_apR_in_star == this->full_map_to_irreducible_sector_.size()); - // this->output_full_map_to_irreducible_sector(st.nat); + // this->output_full_map_to_irreducible_sector(nat); // this->output_sector_star(); this->write_irreducible_sector(output_dir); } diff --git a/source/source_cell/module_symmetry/irreducible_sector.h b/source/source_cell/module_symmetry/irred_sec.h similarity index 87% rename from source/source_cell/module_symmetry/irreducible_sector.h rename to source/source_cell/module_symmetry/irred_sec.h index 73a9248130e..0f68e88103d 100644 --- a/source/source_cell/module_symmetry/irreducible_sector.h +++ b/source/source_cell/module_symmetry/irred_sec.h @@ -1,5 +1,5 @@ -#ifndef ABACUS_SOURCE_CELL_MODULE_SYMMETRY_IRREDUCIBLE_SECTOR_H -#define ABACUS_SOURCE_CELL_MODULE_SYMMETRY_IRREDUCIBLE_SECTOR_H +#ifndef ABACUS_SOURCE_CELL_MODULE_SYMMETRY_IRRED_SECTOR_H +#define ABACUS_SOURCE_CELL_MODULE_SYMMETRY_IRRED_SECTOR_H #include #include @@ -16,7 +16,7 @@ namespace ModuleSymmetry using TapR = std::pair; using TCdouble = ModuleBase::Vector3; - class Irreducible_Sector + class Irred_Sector { public: struct ap_less_func @@ -61,9 +61,11 @@ namespace ModuleSymmetry /// output directory setting); kept as an explicit argument rather than /// reading the global config directly so this LibRI-free class has no /// module_parameter link dependency. - void find_irreducible_sector(const Symmetry& symm, const Atom* atoms, const Statistics& st, + void find_irred_sector(const Symmetry& symm, const Atom* atoms, + const int nat, const int ntype, + const std::vector& iat2it, const std::vector& iat2ia, const std::vector& Rs, const TC& period, const Lattice& lat, const std::string& output_dir = ""); - const std::map>& get_irreducible_sector()const { return this->irreducible_sector_; } + const std::map>& get_irred_sector()const { return this->irreducible_sector_; } // const std::map>> convirt_irreducible_sector() {}; //-------------------------------------------------------------------------------- @@ -86,7 +88,8 @@ namespace ModuleSymmetry return this->return_lattice_[iat][isym]; } protected: - void cal_return_lattice_all(const Symmetry& symm, const Atom* atoms, const Statistics& st); + void cal_return_lattice_all(const Symmetry& symm, const Atom* atoms, + const int nat, const std::vector& iat2it, const std::vector& iat2ia); //-------------------------------------------------------------------------------- /// The sub functions to find irreducible sector: {abR} @@ -95,9 +98,6 @@ namespace ModuleSymmetry /// gauge='L': R'=R+O_1-O_2; gauge='R': R'=R+O_2-O_1 TC rotate_R(const Symmetry& symm, const int isym, const int iat1, const int iat2, const TC& R, const char gauge = 'R')const; TapR rotate_apR_by_formula(const Symmetry& symm, const int isym, const TapR& iapR, const char gauge = 'R')const; - /// gauge='L': tau_a + R - tau_b; gauge='R': tau_a - tau_b - R (direct) - TCdouble get_aRb_direct(const Atom* atoms, const Statistics& st, const int iat1, const int iat2, const TC& R, const char gauge = 'R')const; - TCdouble get_aRb_direct(const Atom* atoms, const Statistics& st, const int iat1, const int iat2, const TCdouble& R, const char gauge = 'R')const; ModuleBase::Matrix3 direct_to_cartesian(const ModuleBase::Matrix3& d, const ModuleBase::Matrix3& latvec)const; @@ -115,7 +115,8 @@ namespace ModuleSymmetry //-------------------------------------------------------------------------------- /// The sub functions judge special symmetry - void gen_symmetry_BvK(const Symmetry& symm, const Atom* atoms, const Lattice& lat, const Statistics& st, const TC bvk_period); + void gen_symm_bvk(const Symmetry& symm, const Atom* atoms, const Lattice& lat, + const int nat, const int ntype, const TC bvk_period); /// whether in 2D plain or not for each symmetry operation // std::vector in_plain(const ModuleSymmetry::Symmetry& symm, const ModuleBase::Matrix3& latvec)const; //-------------------------------------------------------------------------------- diff --git a/source/source_cell/module_symmetry/irred_sec_bvk.cpp b/source/source_cell/module_symmetry/irred_sec_bvk.cpp index 5d506cab434..714e9cf5063 100644 --- a/source/source_cell/module_symmetry/irred_sec_bvk.cpp +++ b/source/source_cell/module_symmetry/irred_sec_bvk.cpp @@ -1,7 +1,8 @@ -#include "irreducible_sector.h" +#include "irred_sec.h" +#include "source_base/global_variable.h" namespace ModuleSymmetry { - ModuleBase::Matrix3 Irreducible_Sector::direct_to_cartesian(const ModuleBase::Matrix3& d, const ModuleBase::Matrix3& latvec)const + ModuleBase::Matrix3 Irred_Sector::direct_to_cartesian(const ModuleBase::Matrix3& d, const ModuleBase::Matrix3& latvec)const { return latvec.Inverse() * d * latvec; } @@ -43,7 +44,7 @@ namespace ModuleSymmetry // Unmatched stays -1. That is legitimate for nspin=4 magnetic: the unit-cell set is the // Shubnikov group H (union) A, generally a PROPER subset of the crystallographic group // (operations that merely tilt the moment belong to neither), so a BvK operation may - // have no counterpart. The consumer in find_irreducible_sector skips negative entries; + // have no counterpart. The consumer in find_irred_sector skips negative entries; // it must never use one as an index. } return isymbvk2isym; @@ -55,9 +56,10 @@ namespace ModuleSymmetry int c = a % b; return (c == 0) ? b : gcd(b, c); } - void Irreducible_Sector::gen_symmetry_BvK(const ModuleSymmetry::Symmetry& symm, const Atom* atoms, const Lattice& lat, const Statistics& st, const TC bvk_period) + void Irred_Sector::gen_symm_bvk(const ModuleSymmetry::Symmetry& symm, const Atom* atoms, const Lattice& lat, + const int nat, const int ntype, const TC bvk_period) { - ModuleBase::TITLE("Irreducible_Sector", "gen_symmetry_BvK"); + ModuleBase::TITLE("Irred_Sector", "gen_symm_bvk"); auto set_matrix3 = [](const ModuleBase::Vector3& a1, const ModuleBase::Vector3& a2, const ModuleBase::Vector3& a3) -> ModuleBase::Matrix3 {return ModuleBase::Matrix3(a1.x, a1.y, a1.z, a2.x, a2.y, a2.z, a3.x, a3.y, a3.z);}; auto set_bvk_same_as_ucell = [&symm, this]()->void @@ -78,11 +80,11 @@ namespace ModuleSymmetry // extern lattice to minimal BvK lattice, and set direct coordinates in min BvK lattice int bvk_gcd = gcd(bvk_period[0], gcd(bvk_period[1], bvk_period[2])); const TC bvk_min_period = TC({ bvk_period[0] / bvk_gcd, bvk_period[1] / bvk_gcd, bvk_period[2] / bvk_gcd }); - const int bvk_nat = st.nat * bvk_min_period[0] * bvk_min_period[1] * bvk_min_period[2]; - std::vector bvk_na(st.ntype); - std::vector bvk_istart(st.ntype, 0); + const int bvk_nat = nat * bvk_min_period[0] * bvk_min_period[1] * bvk_min_period[2]; + std::vector bvk_na(ntype); + std::vector bvk_istart(ntype, 0); int bvk_itmin_start = 0, bvk_itmin_type = 0; - for (int it = 0;it < st.ntype;++it) + for (int it = 0;it < ntype;++it) { bvk_na[it] = atoms[it].na * bvk_min_period[0] * bvk_min_period[1] * bvk_min_period[2]; if (it > 0) { bvk_istart[it] = bvk_istart[it - 1] + bvk_na[it - 1]; @@ -103,7 +105,7 @@ namespace ModuleSymmetry s3 = a3 = lat.a3 * static_cast(bvk_min_period[2]); ModuleBase::Matrix3 bvk_min_lat = set_matrix3(s1, s2, s3); int at = 0; - for (int it = 0; it < st.ntype; ++it) { + for (int it = 0; it < ntype; ++it) { for (int c1 = 0;c1 < bvk_min_period[0];++c1) { for (int c2 = 0;c2 < bvk_min_period[1];++c2) { for (int c3 = 0;c3 < bvk_min_period[2];++c3) { @@ -184,7 +186,7 @@ namespace ModuleSymmetry return; } - // std::vector Irreducible_Sector::in_plain(const ModuleSymmetry::Symmetry& symm, const ModuleBase::Matrix3& latvec)const + // std::vector Irred_Sector::in_plain(const ModuleSymmetry::Symmetry& symm, const ModuleBase::Matrix3& latvec)const // { // // get euler angel of the cartesian gmatrix in optimal lattice // std::vector gmatc(symm.nrotk); diff --git a/source/source_cell/module_symmetry/symm_analysis.cpp b/source/source_cell/module_symmetry/symm_analysis.cpp index 9f4361a06d3..330aba413af 100644 --- a/source/source_cell/module_symmetry/symm_analysis.cpp +++ b/source/source_cell/module_symmetry/symm_analysis.cpp @@ -3,7 +3,12 @@ using namespace ModuleSymmetry; -void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, std::ofstream& ofs_running, +void Symmetry::analy_sys(const Lattice& lat, Atom* atoms, + const int nat, const int ntype, + const std::vector& iat2it, + const std::vector& iat2ia, + const ModuleBase::IntArray& itia2iat, + std::ofstream& ofs_running, const double symmetry_prec, const int nspin, const std::string& calculation, const int* cal_symm_repr) { @@ -30,9 +35,9 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, // 1. copy data and allocate memory // -------------------------------- // number of total atoms - this->nat = st.nat; + this->nat = nat; // number of atom species - this->ntype = st.ntype; + this->ntype = ntype; assert(ntype>0); @@ -116,7 +121,7 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, if (!pricell_loop && nspin == 2) { - this->analyze_magnetic_group(atoms, st, nrot_out, nrotk_out); + this->analyze_magnetic_group(atoms, iat2it, iat2ia, itia2iat, nrot_out, nrotk_out); } else { @@ -288,13 +293,13 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, // the magnetization (pseudovector), so they are not applied in k-reduction / density symmetrization. if (nspin == 4) { - this->analyze_magnetic_group_nspin4(atoms, st, latvec1); + this->analyze_magnetic_group_nspin4(atoms, iat2it, iat2ia, latvec1); } // Do this here for debug if (calculation == "relax") { - this->all_mbl = this->is_all_movable(atoms, st); + this->all_mbl = this->is_all_movable(atoms, iat2it, iat2ia); if (!this->all_mbl) { std::cout << "WARNING: Symmetry cannot be kept when not all atoms are movable.\n "; diff --git a/source/source_cell/module_symmetry/symm_magnetic.cpp b/source/source_cell/module_symmetry/symm_magnetic.cpp index c236a0fdad9..c0f74ea8d0c 100644 --- a/source/source_cell/module_symmetry/symm_magnetic.cpp +++ b/source/source_cell/module_symmetry/symm_magnetic.cpp @@ -7,7 +7,11 @@ using namespace ModuleSymmetry; #include #include -void Symmetry::analyze_magnetic_group(const Atom* atoms, const Statistics& st, int& nrot_out, int& nrotk_out) +void Symmetry::analyze_magnetic_group(const Atom* atoms, + const std::vector& iat2it, + const std::vector& iat2ia, + const ModuleBase::IntArray& itia2iat, + int& nrot_out, int& nrotk_out) { // 1. classify atoms with different magmom // (use symmetry_prec to judge if two magmoms are the same) @@ -20,18 +24,18 @@ void Symmetry::analyze_magnetic_group(const Atom* atoms, const Statistics& st, i 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]; + const int mag_it = iat2it[mag_iat]; + const int mag_ia = 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)); + mt.insert(itia2iat(it, ia)); find = true; break; } } if (!find) { - mag_type_atoms.push_back(std::set({ st.itia2iat(it,ia) })); + mag_type_atoms.push_back(std::set({ itia2iat(it, ia) })); } } } @@ -58,7 +62,7 @@ void Symmetry::analyze_magnetic_group(const Atom* atoms, const Statistics& st, i { // 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(); + const ModuleBase::Vector3 direct_tmp = atoms[iat2it[mag_iat]].tau[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) { @@ -78,7 +82,10 @@ void Symmetry::analyze_magnetic_group(const Atom* atoms, const Statistics& st, i } -void Symmetry::analyze_magnetic_group_nspin4(const Atom* atoms, const Statistics& st, const ModuleBase::Matrix3& latvec) +void Symmetry::analyze_magnetic_group_nspin4(const Atom* atoms, + const std::vector& iat2it, + const std::vector& iat2ia, + const ModuleBase::Matrix3& latvec) { // Restrict the space group to the unitary magnetic subgroup (nspin=4 / SOC): // operation g survives if it preserves the magnetic configuration as a pseudovector, @@ -97,7 +104,7 @@ void Symmetry::analyze_magnetic_group_nspin4(const Atom* atoms, const Statistics bool has_moment = false; for (int iat = 0; iat < this->nat && !has_moment; ++iat) { - const ModuleBase::Vector3& m = atoms[st.iat2it[iat]].m_loc_[st.iat2ia[iat]]; + const ModuleBase::Vector3& m = atoms[iat2it[iat]].m_loc_[iat2ia[iat]]; if (!this->equal(m.x, 0.0) || !this->equal(m.y, 0.0) || !this->equal(m.z, 0.0)) { has_moment = true; } } std::vector anti; // operations that REVERSE the moment: Theta*g is a symmetry @@ -109,11 +116,11 @@ void Symmetry::analyze_magnetic_group_nspin4(const Atom* atoms, const Statistics bool ok = true; for (int iat = 0; iat < this->nat && ok; ++iat) { - const ModuleBase::Vector3& m = atoms[st.iat2it[iat]].m_loc_[st.iat2ia[iat]]; + const ModuleBase::Vector3& m = atoms[iat2it[iat]].m_loc_[iat2ia[iat]]; // pseudovector-rotated moment W*m (column-vector convention: m'^i = W_ij m^j) const ModuleBase::Vector3& mrot = W * m; const int jat = this->get_rotated_atom(isym, iat); - const ModuleBase::Vector3& mj = atoms[st.iat2it[jat]].m_loc_[st.iat2ia[jat]]; + const ModuleBase::Vector3& mj = atoms[iat2it[jat]].m_loc_[iat2ia[jat]]; if (!this->equal(mrot.x, mj.x) || !this->equal(mrot.y, mj.y) || !this->equal(mrot.z, mj.z)) { ok = false; } } if (ok) @@ -129,10 +136,10 @@ void Symmetry::analyze_magnetic_group_nspin4(const Atom* atoms, const Statistics bool anti_ok = true; for (int iat = 0; iat < this->nat && anti_ok; ++iat) { - const ModuleBase::Vector3& m = atoms[st.iat2it[iat]].m_loc_[st.iat2ia[iat]]; + const ModuleBase::Vector3& m = atoms[iat2it[iat]].m_loc_[iat2ia[iat]]; const ModuleBase::Vector3& mrot = W * m; const int jat = this->get_rotated_atom(isym, iat); - const ModuleBase::Vector3& mj = atoms[st.iat2it[jat]].m_loc_[st.iat2ia[jat]]; + const ModuleBase::Vector3& mj = atoms[iat2it[jat]].m_loc_[iat2ia[jat]]; if (!this->equal(mrot.x, -mj.x) || !this->equal(mrot.y, -mj.y) || !this->equal(mrot.z, -mj.z)) { anti_ok = false; } } if (anti_ok) { anti.push_back(isym); } diff --git a/source/source_cell/module_symmetry/symm_rotation_k.cpp b/source/source_cell/module_symmetry/symm_rotation_k.cpp index 71b145a2147..ac145a67e16 100644 --- a/source/source_cell/module_symmetry/symm_rotation_k.cpp +++ b/source/source_cell/module_symmetry/symm_rotation_k.cpp @@ -12,7 +12,7 @@ namespace ModuleSymmetry { void Symmetry_rotation_k::reset_symmetry() { - this->irs_ = Irreducible_Sector(); + this->irs_ = Irred_Sector(); this->Ms_.clear(); this->little_groups_.clear(); this->spin_U_.clear(); @@ -33,29 +33,31 @@ namespace ModuleSymmetry return cells; } - void Symmetry_rotation_k::cal_Ms(const K_Vectors& kv, - const UnitCell& ucell, const Parallel_2D& pv, const int nspin) + void Symmetry_rotation_k::cal_Ms(const K_Vectors& kv, const Symmetry& symm, const Atom* atoms, + const ModuleBase::Matrix3& latvec, const int lmax, + const int nat, const std::vector& iat2it, const std::vector& iat2ia, + const Parallel_2D& pv, const int nspin) { ModuleBase::TITLE("Symmetry_rotation_k", "cal_Ms"); ModuleBase::timer::start("Symmetry_rotation_k", "cal_Ms"); this->nspin_ = nspin; - this->nsym_ = ucell.symm.nrotk; - this->nanti_ = ucell.symm.nrotk_anti; - this->magnetic_nspin4_ = ucell.symm.magnetic_nspin4; - this->eps_ = ucell.symm.epsilon; + this->nsym_ = symm.nrotk; + this->nanti_ = symm.nrotk_anti; + this->magnetic_nspin4_ = symm.magnetic_nspin4; + this->eps_ = symm.epsilon; if (this->irs_.invmap_.empty()) { - this->irs_.invmap_.resize(ucell.symm.nrotk); - ucell.symm.gmatrix_invmap(ucell.symm.gmatrix, ucell.symm.nrotk, this->irs_.invmap_.data()); + this->irs_.invmap_.resize(symm.nrotk); + symm.gmatrix_invmap(symm.gmatrix, symm.nrotk, this->irs_.invmap_.data()); } // 1. calculate the rotation matrix in real spherical harmonics representation for each symmetry operation: [T_l (isym)]_mm' const int nop_tot = this->nsym_ + this->nanti_; std::vector gmatc(nop_tot); - for (int i = 0;i < nsym_;++i) { gmatc[i] = this->irs_.direct_to_cartesian(ucell.symm.gmatrix[i], ucell.latvec); } + for (int i = 0;i < nsym_;++i) { gmatc[i] = this->irs_.direct_to_cartesian(symm.gmatrix[i], latvec); } for (int j = 0;j < this->nanti_;++j) - { gmatc[nsym_ + j] = this->irs_.direct_to_cartesian(ucell.symm.gmatrix_anti[j], ucell.latvec); } - this->cal_rotmat_Slm(gmatc.data(), std::max(this->abfs_Lmax_, ucell.lmax), nop_tot); + { gmatc[nsym_ + j] = this->irs_.direct_to_cartesian(symm.gmatrix_anti[j], latvec); } + this->cal_rotmat_Slm(gmatc.data(), std::max(this->abfs_Lmax_, lmax), nop_tot); // 1.5 (nspin=4) the SU(2) spin-1/2 rotation U(isym) for each symmetry operation. The AO // rotation matrix M becomes the spinor operator T(isym) (x) U(isym) so that the same @@ -111,7 +113,7 @@ namespace ModuleSymmetry } for (int op = 0; op < nsym_; ++op) { - const ModuleBase::Vector3 delta = kvec_d_ibz_global[ik_ibz] * ucell.symm.kgmatrix[op] - kvec_d_ibz_global[ik_ibz]; + const ModuleBase::Vector3 delta = kvec_d_ibz_global[ik_ibz] * symm.kgmatrix[op] - kvec_d_ibz_global[ik_ibz]; if (std::abs(delta.x - std::round(delta.x)) < this->eps_ && std::abs(delta.y - std::round(delta.y)) < this->eps_ && std::abs(delta.z - std::round(delta.z)) < this->eps_) @@ -123,7 +125,7 @@ namespace ModuleSymmetry for (const int op : needed) { this->Ms_[ik_ibz][op] = this->contruct_2d_rot_mat_ao( - ucell.symm, ucell.atoms, ucell.st, kvec_d_ibz_global[ik_ibz], op, pv, spin_U[op]); + symm, atoms, nat, iat2it, iat2ia, kvec_d_ibz_global[ik_ibz], op, pv, spin_U[op]); } } @@ -399,7 +401,8 @@ namespace ModuleSymmetry // 2d-block parallized rotation matrix in AO-representation, denoted as M. // finally we will use D(k)=M(R, k)^\dagger*D(Rk)*M(R, k) to D(k) from D(Rk) in cal_Ms. - std::vector> Symmetry_rotation_k::contruct_2d_rot_mat_ao(const Symmetry& symm, const Atom* atoms, const Statistics& cell_st, + std::vector> Symmetry_rotation_k::contruct_2d_rot_mat_ao(const Symmetry& symm, const Atom* atoms, + const int nat, const std::vector& iat2it, const std::vector& iat2ia, const TCdouble& kvec_d_ibz, int isym, const Parallel_2D& pv, const SpinRotation::Su2& spin_U) const { const bool soc = (this->nspin_ == 4); @@ -413,12 +416,12 @@ namespace ModuleSymmetry return (is < nrotk_u) ? symm.get_rotated_atom(is, iat) : symm.get_rotated_atom_anti(is - nrotk_u, iat); }; - for (int iat1 = 0;iat1 < cell_st.nat;++iat1) + for (int iat1 = 0;iat1 < nat;++iat1) { - int it = cell_st.iat2it[iat1]; // it1=it2 - int ia1 = cell_st.iat2ia[iat1]; + int it = iat2it[iat1]; // it1=it2 + int ia1 = iat2ia[iat1]; int iat2 = rotated_atom(isym, iat1); //iat2=rot(iat1) - int ia2 = cell_st.iat2ia[iat2]; + int ia2 = iat2ia[iat2]; // cal phase factor from return lattice: exp(-ik_ibz*O) double arg = -2 * ModuleBase::PI * kvec_d_ibz * this->irs_.return_lattice_[iat1][isym]; std::complexphase_factor = std::complex(std::cos(arg), std::sin(arg)); diff --git a/source/source_cell/module_symmetry/symm_rotation_k.h b/source/source_cell/module_symmetry/symm_rotation_k.h index 40efdc18ba9..ee45b40b34b 100644 --- a/source/source_cell/module_symmetry/symm_rotation_k.h +++ b/source/source_cell/module_symmetry/symm_rotation_k.h @@ -1,6 +1,6 @@ #ifndef SYMM_ROTATION_K_H #define SYMM_ROTATION_K_H -#include "irreducible_sector.h" +#include "irred_sec.h" #include "source_base/parallel_2d.h" #include "source_base/complexmatrix.h" #include "source_cell/module_symmetry/symm_rot_spin.h" @@ -29,7 +29,7 @@ namespace ModuleSymmetry //-------------------------------------------------------------------------------- // getters - const std::map>& get_irreducible_sector()const { return this->irs_.get_irreducible_sector(); } + const std::map>& get_irred_sector()const { return this->irs_.get_irred_sector(); } TCdouble get_return_lattice(const Symmetry& symm, const ModuleBase::Matrix3& gmatd, const TCdouble gtransd, const TCdouble& posd_a1, const TCdouble& posd_a2)const @@ -45,10 +45,13 @@ namespace ModuleSymmetry const int& abfs_Lmax = this->abfs_Lmax_; //-------------------------------------------------------------------------------- // setters - void find_irreducible_sector(const Symmetry& symm, const Atom* atoms, const Statistics& st, + void find_irred_sector(const Symmetry& symm, const Atom* atoms, + const int nat, const int ntype, + const std::vector& iat2it, const std::vector& iat2ia, const std::vector& Rs, const TC& period, const Lattice& lat, const std::string& output_dir = "") { - this->irs_.find_irreducible_sector(symm, atoms, st, Rs, period, lat, output_dir); + this->irs_.find_irred_sector(symm, atoms, nat, ntype, iat2it, iat2ia, + Rs, period, lat, output_dir); } void set_abfs_Lmax(const int l) { this->abfs_Lmax_ = l; } //-------------------------------------------------------------------------------- @@ -57,12 +60,21 @@ namespace ModuleSymmetry /// The top-level calculation interface of this class. calculate the rotation matrix in AO representation: M /// only need once call in each ion step (decided by the configuration) /// @param kstars equal k points to each ibz-kpont, corresponding to a certain symmetry operations. + /// @param symm analyzed symmetry of the cell (gmatrix, gmatrix_anti, kgmatrix, ...) + /// @param atoms all atoms + /// @param latvec Cartesian lattice vectors, as rows of a 3x3 matrix + /// @param lmax maximum angular momentum of the atomic orbitals + /// @param nat total number of atoms + /// @param iat2it map from global atom index iat to its species index it + /// @param iat2ia map from global atom index iat to its index ia within the species /// @param nspin stored as a member so restore_dm()/contruct_2d_rot_mat_ao() do not each /// need to read the global nspin config setting (keeps this LibRI-free class /// free of a module_parameter link dependency; every existing caller already /// has nspin in scope). - void cal_Ms(const K_Vectors& kv, - const UnitCell& ucell, const Parallel_2D& pv, const int nspin); + void cal_Ms(const K_Vectors& kv, const Symmetry& symm, const Atom* atoms, + const ModuleBase::Matrix3& latvec, const int lmax, + const int nat, const std::vector& iat2it, const std::vector& iat2ia, + const Parallel_2D& pv, const int nspin); /// Use calculated M matrix to recover D(k) from D(k_ibz): D(k) = M(R, k)^\dagger D(k_ibz) M(R, k) /// the link "ik_ibz-isym-ik" can be found in kstars: k_bz = gmat[isym](k) @@ -108,7 +120,8 @@ namespace ModuleSymmetry /// 2d-block parallized rotation matrix in AO-representation, denoted as M. /// finally we will use D(k)=M(R, k)^\dagger*D(Rk)*M(R, k) to recover D(k) from D(Rk). - std::vector> contruct_2d_rot_mat_ao(const Symmetry& symm, const Atom* atoms, const Statistics& cell_st, + std::vector> contruct_2d_rot_mat_ao(const Symmetry& symm, const Atom* atoms, + const int nat, const std::vector& iat2it, const std::vector& iat2ia, const TCdouble& kvec_d_ibz, int isym, const Parallel_2D& pv, const SpinRotation::Su2& spin_U /*= SpinRotation::Su2{ 1.0, 0.0, 0.0, 1.0 }*/) const; @@ -174,7 +187,7 @@ namespace ModuleSymmetry std::vector spin_U_; /// irreducible sector - Irreducible_Sector irs_; + Irred_Sector irs_; }; } #endif // SYMM_ROTATION_K_H diff --git a/source/source_cell/module_symmetry/symmetry.cpp b/source/source_cell/module_symmetry/symmetry.cpp index a291e21a405..c0f0bf3189e 100644 --- a/source/source_cell/module_symmetry/symmetry.cpp +++ b/source/source_cell/module_symmetry/symmetry.cpp @@ -364,13 +364,14 @@ void Symmetry::get_optlat(ModuleBase::Vector3 &v1, ModuleBase::Vector3& iat2it, + const std::vector& iat2ia) const { bool all_mbl = true; - for (int iat = 0;iat < st.nat;++iat) + for (int iat = 0;iat < this->nat;++iat) { - int it = st.iat2it[iat]; - int ia = st.iat2ia[iat]; + int it = iat2it[iat]; + int ia = iat2ia[iat]; if (!atoms[it].mbl[ia].x || !atoms[it].mbl[ia].y || !atoms[it].mbl[ia].z) { all_mbl = false; diff --git a/source/source_cell/module_symmetry/symmetry.h b/source/source_cell/module_symmetry/symmetry.h index be36a233213..7d9714ef36e 100644 --- a/source/source_cell/module_symmetry/symmetry.h +++ b/source/source_cell/module_symmetry/symmetry.h @@ -41,15 +41,24 @@ class Symmetry : public Symmetry_Basic /// @brief analyze the symmetry of the system /// @param lat structure of lattice - /// @param st /// @param atoms all atoms - /// @param ofs_running + /// @param nat total number of atoms + /// @param ntype number of atom species + /// @param iat2it map from global atom index iat to its species index it + /// @param iat2ia map from global atom index iat to its index ia within the species + /// @param itia2iat map from (it, ia) to the global atom index iat + /// @param ofs_running /// @param symmetry_prec precision for symmetry analysis /// @param nspin number of spin components /// @param calculation calculation type (scf, relax, cell-relax, etc.) /// @param cal_symm_repr control for symmetry representation output [0]=flag, [1]=precision /// get the symmetry information of the system, gmatries (rotation 3*3 matrixs), gtrans (transfer a collections vector3), etc. - void analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, std::ofstream& ofs_running, + void analy_sys(const Lattice& lat, Atom* atoms, + const int nat, const int ntype, + const std::vector& iat2it, + const std::vector& iat2ia, + const ModuleBase::IntArray& itia2iat, + std::ofstream& ofs_running, const double symmetry_prec, const int nspin, const std::string& calculation, const int* cal_symm_repr); @@ -351,7 +360,10 @@ class Symmetry : public Symmetry_Basic void set_atom_map(const Atom* atoms); /// @brief check if all the atoms are movable /// delta_pos symmetrization in relax is only meaningful when all the atoms are movable in all the directions. - bool is_all_movable(const Atom* atoms, const Statistics& st)const; + /// @param iat2it map from global atom index iat to its species index it + /// @param iat2ia map from global atom index iat to its index ia within the species + bool is_all_movable(const Atom* atoms, const std::vector& iat2it, + const std::vector& iat2ia) const; // to be called in lattice_type void get_shortest_latvec(ModuleBase::Vector3 &a1, @@ -366,10 +378,14 @@ class Symmetry : public Symmetry_Basic /// 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 + /// Analyze magnetic group without time-reversal symmetry /// (because currently the charge density symmetrization does not support it) /// Method: treat atoms with different magmom as atoms of different type - void analyze_magnetic_group(const Atom* atoms, const Statistics& st, int& nrot_out, int& nrotk_out); + void analyze_magnetic_group(const Atom* atoms, + const std::vector& iat2it, + const std::vector& iat2ia, + const ModuleBase::IntArray& itia2iat, + int& nrot_out, int& nrotk_out); /// (nspin=4 / SOC) Restrict the already-built space group to the unitary magnetic /// subgroup: keep operation g only if it preserves the magnetization as a pseudovector, @@ -377,7 +393,10 @@ class Symmetry : public Symmetry_Basic /// that reverse the moment (which are only symmetries when combined with time reversal) /// from being applied in k-reduction and density symmetrization. /// Non-magnetic (m_i=0) keeps all operations. - void analyze_magnetic_group_nspin4(const Atom* atoms, const Statistics& st, const ModuleBase::Matrix3& latvec); + void analyze_magnetic_group_nspin4(const Atom* atoms, + const std::vector& iat2it, + const std::vector& iat2ia, + const ModuleBase::Matrix3& latvec); }; /** diff --git a/source/source_cell/module_symmetry/test/symm_test_analysis.cpp b/source/source_cell/module_symmetry/test/symm_test_analysis.cpp index 24dcbbd26d8..389342700bd 100644 --- a/source/source_cell/module_symmetry/test/symm_test_analysis.cpp +++ b/source/source_cell/module_symmetry/test/symm_test_analysis.cpp @@ -55,8 +55,9 @@ TEST_F(SymmetryTest, AtomMapRefreshWithUnchangedOperationCount) ModuleSymmetry::Symmetry symm; const int representation[2] = {0, 0}; const std::string calculation = "cell-relax"; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, - 1e-6, 1, calculation, representation); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, calculation, representation); const int original_count = symm.nrotk; EXPECT_EQ(original_count, 48); @@ -70,8 +71,9 @@ TEST_F(SymmetryTest, AtomMapRefreshWithUnchangedOperationCount) { ucell.atoms[0].tau[atom] = ucell.atoms[0].taud[atom] * ucell.latvec; } - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, - 1e-6, 1, calculation, representation); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, calculation, representation); EXPECT_EQ(symm.nrotk, original_count); for (int operation = 0; operation < symm.nrotk; ++operation) { @@ -104,7 +106,9 @@ TEST_F(SymmetryTest, AnalySys) ModuleSymmetry::Symmetry symm; construct_ucell(stru_lib[stru]); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); //1. ibrav std::string ref_point_group = stru_lib[stru].point_group; @@ -308,7 +312,9 @@ TEST_F(SymmetryTest, SG_Pricell) symm.epsilon = 1e-5; construct_ucell(supercell_lib[stru]); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); std::string ref_point_group = supercell_lib[stru].point_group; std::string cal_point_group = symm.pgname; diff --git a/source/source_cell/module_symmetry/test/symm_test_symtrz.cpp b/source/source_cell/module_symmetry/test/symm_test_symtrz.cpp index 79f32a9a9cd..845a446a880 100644 --- a/source/source_cell/module_symmetry/test/symm_test_symtrz.cpp +++ b/source/source_cell/module_symmetry/test/symm_test_symtrz.cpp @@ -71,7 +71,9 @@ TEST_F(SymmetryTest, ForceSymmetry) ModuleSymmetry::Symmetry symm; construct_ucell(supercell_lib[stru]); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); ModuleBase::matrix force(ucell.nat, 3, true); //generate random number for force and restrict to [-100,100) @@ -103,7 +105,9 @@ TEST_F(SymmetryTest, StressSymmetry) ModuleSymmetry::Symmetry symm; construct_ucell(supercell_lib[stru]); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); ModuleBase::matrix stress(3, 3, true); //generate random number for stress and restrict to [-1e5,1e5) diff --git a/source/source_cell/test/klist_test.cpp b/source/source_cell/test/klist_test.cpp index 62aaef39837..bdc98c45600 100644 --- a/source/source_cell/test/klist_test.cpp +++ b/source/source_cell/test/klist_test.cpp @@ -844,7 +844,9 @@ TEST_F(KlistTest, IbzKpoint) construct_ucell(stru_lib[0]); ofs_running.open("tmp_klist_3"); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); std::string k_file = "./support/KPT1"; kv->set_spin_mult(1); kv->read_kpoints_for_testing(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); @@ -871,7 +873,9 @@ TEST_F(KlistTest, IbzKpointIsMP) construct_ucell(stru_lib[0]); ofs_running.open("tmp_klist_4"); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); std::string k_file = "./support/KPT1"; kv->set_spin_mult(1); kv->read_kpoints_for_testing(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); @@ -899,7 +903,9 @@ TEST_F(KlistTest, IbzKpointCustomWeights) construct_ucell(stru_lib[0]); ofs_running.open("tmp_klist_custom_weights"); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); // Test 1: Non-MP k-points with uniform weights (KPT4) { diff --git a/source/source_cell/test/klist_test_para.cpp b/source/source_cell/test/klist_test_para.cpp index 7f8f76d524c..1e20b685970 100644 --- a/source/source_cell/test/klist_test_para.cpp +++ b/source/source_cell/test/klist_test_para.cpp @@ -178,7 +178,9 @@ TEST_F(KlistParaTest, Set) GlobalV::ofs_running.open("tmp_klist_5"); } const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + GlobalV::ofs_running, 1e-6, 1, "scf", cal_symm_repr); // read KPT std::string k_file = "./support/KPT1"; // note: do NOT pre-set kv->spin_mult here; set() takes the physical @@ -302,7 +304,9 @@ TEST_F(KlistParaTest, SetAfterVC) GlobalV::ofs_running.open("tmp_klist_6"); } const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + GlobalV::ofs_running, 1e-6, 1, "scf", cal_symm_repr); // read KPT std::string k_file = "./support/KPT1"; // note: do NOT pre-set kv->spin_mult here; set() takes the physical diff --git a/source/source_cell/test/little_group_test.cpp b/source/source_cell/test/little_group_test.cpp index 153dc9f501d..1f0c6c6b2b3 100644 --- a/source/source_cell/test/little_group_test.cpp +++ b/source/source_cell/test/little_group_test.cpp @@ -120,7 +120,9 @@ class LittleGroupTest : public testing::Test } std::ofstream ofs_running("tmp_little_group"); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); } void TearDown() override diff --git a/source/source_cell/test/qlist_test.cpp b/source/source_cell/test/qlist_test.cpp index 637898edf44..1344a692ebc 100644 --- a/source/source_cell/test/qlist_test.cpp +++ b/source/source_cell/test/qlist_test.cpp @@ -173,7 +173,9 @@ TEST_F(QListTest, GenerateMeshFullSymmetry) ofs_running.open("tmp_qlist_1"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {8, 8, 8}, true); @@ -211,7 +213,9 @@ TEST_F(QListTest, GenerateMeshSmallGrid) ofs_running.open("tmp_qlist_2"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {2, 2, 2}, true); @@ -235,7 +239,9 @@ TEST_F(QListTest, GammaOnlyGrid) ofs_running.open("tmp_qlist_3"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {1, 1, 1}, true); @@ -255,7 +261,9 @@ TEST_F(QListTest, IrrepPlaceholder) ofs_running.open("tmp_qlist_4"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {2, 2, 2}, true); @@ -282,7 +290,9 @@ TEST_F(QListTest, CartesianCoordinatesComputed) ofs_running.open("tmp_qlist_cart"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {2, 2, 2}, true); @@ -312,7 +322,9 @@ TEST_F(QListTest, UseIrrepsSwitch) ofs_running.open("tmp_qlist_irreps"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); // use_irreps = false: the q mesh is still reduced, but no irrep data qlist.generate_mesh(ucell, symm, {2, 2, 2}, false); @@ -331,7 +343,9 @@ TEST_F(QListTest, PrintQlists) ofs_running.open("tmp_qlist_print"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {1, 1, 1}, false); diff --git a/source/source_cell/test/test_unitcell_index.cpp b/source/source_cell/test/test_unitcell_index.cpp index 478fc66b04e..792096a86c1 100644 --- a/source/source_cell/test/test_unitcell_index.cpp +++ b/source/source_cell/test/test_unitcell_index.cpp @@ -22,8 +22,7 @@ class TestUnitCellIndex : public ::testing::Test void TearDown() override { // Do NOT delete[] ucell.atoms: ~UnitCell() frees it (set_atom_flag = true). - // Do NOT delete[] iat2it/iat2ia either: they are owned by the internal - // Statistics member, whose destructor releases them (see AGENTS.md). + // iat2it/iat2ia are std::vector members owned by UnitCell. } UnitCell ucell; @@ -44,8 +43,8 @@ TEST_F(TestUnitCellIndex, SetIat2itia) TEST_F(TestUnitCellIndex, SetIat2itiaCalledTwice) { - // Cover the delete[] + new[] reallocation path: the second call must - // produce the same mapping as the first. + // Cover the resize-on-repeat path: the second call must produce the + // same mapping as the first. ucell.set_iat2itia(); ucell.set_iat2itia(); EXPECT_EQ(ucell.iat2it[0], 0); diff --git a/source/source_cell/test/unitcell_test.cpp b/source/source_cell/test/unitcell_test.cpp index 630390c5b7d..8f65ba87445 100644 --- a/source/source_cell/test/unitcell_test.cpp +++ b/source/source_cell/test/unitcell_test.cpp @@ -7,6 +7,7 @@ #include "source_base/mathzone.h" #include "prepare_unitcell.h" #include +#include #include #include @@ -52,6 +53,22 @@ class UcellTest : public ::testing::Test std::string output; }; +/// Compile-time guard: UnitCell owns the raw pointer 'atoms' tracked by +/// 'set_atom_flag' and exposes reference aliases into its 'lat' member, so +/// neither copy nor move can be implemented correctly until that ownership +/// model is refactored. Re-enabling either operation must fail this TU. +TEST(UnitCellTypeTraits, NotCopyableOrMovable) +{ + static_assert(!std::is_copy_constructible::value, + "UnitCell must not be copy constructible"); + static_assert(!std::is_copy_assignable::value, + "UnitCell must not be copy assignable"); + static_assert(!std::is_move_constructible::value, + "UnitCell must not be move constructible"); + static_assert(!std::is_move_assignable::value, + "UnitCell must not be move assignable"); +} + using UcellDeathTest = UcellTest; TEST_F(UcellTest, Constructor) @@ -64,10 +81,10 @@ TEST_F(UcellTest, Constructor) EXPECT_EQ(ucell->nat, 0); EXPECT_EQ(ucell->namax, 0); EXPECT_EQ(ucell->nwmax, 0); - EXPECT_EQ(ucell->iat2it, nullptr); - EXPECT_EQ(ucell->iat2ia, nullptr); - EXPECT_EQ(ucell->iwt2iat, nullptr); - EXPECT_EQ(ucell->iwt2iw, nullptr); + EXPECT_TRUE(ucell->iat2it.empty()); + EXPECT_TRUE(ucell->iat2ia.empty()); + EXPECT_TRUE(ucell->iwt2iat.empty()); + EXPECT_TRUE(ucell->iwt2iw.empty()); EXPECT_DOUBLE_EQ(ucell->tpiba, 0.0); EXPECT_DOUBLE_EQ(ucell->tpiba2, 0.0); EXPECT_DOUBLE_EQ(ucell->omega, 0.0); diff --git a/source/source_cell/unitcell.h b/source/source_cell/unitcell.h index 359456737a3..de3c04b8e02 100644 --- a/source/source_cell/unitcell.h +++ b/source/source_cell/unitcell.h @@ -2,6 +2,7 @@ #define UNITCELL_H #include +#include #include "source_base/global_function.h" #include "source_cell/sep_cell.h" #include "source_cell/magnetism.h" @@ -17,6 +18,18 @@ class UnitCell : public BaseCell { UnitCell(); ~UnitCell(); + /// Copy and move are both disabled. The defaulted move constructor + /// shallow-copies the owning pointer 'atoms' and 'set_atom_flag' without + /// resetting the source (double free), and the reference aliases into + /// 'lat' (Coordinate, a1, latvec, ...) make the defaulted move assignment + /// implicitly deleted (-Wdefaulted-function-deleted). Correct move + /// semantics require an owning container for 'atoms' and removing the + /// aliases; until then any accidental move must fail at compile time. + UnitCell(const UnitCell&) = delete; + UnitCell& operator=(const UnitCell&) = delete; + UnitCell(UnitCell&&) = delete; + UnitCell& operator=(UnitCell&&) = delete; + /// @brief Initialize basic cell parameters (latname, ntype, lmaxmax, init_vel) /// from INPUT and parse fixed_axes into lat_axis_free flags. void setup_from_input(const std::string& latname_in, @@ -138,16 +151,15 @@ class UnitCell : public BaseCell { /// @name Statistics /// @{ - Statistics st; - int& ntype = st.ntype; - int& nat = st.nat; - int*& iat2it = st.iat2it; - int*& iat2ia = st.iat2ia; - int*& iwt2iat = st.iwt2iat; - int*& iwt2iw = st.iwt2iw; - ModuleBase::IntArray& itia2iat = st.itia2iat; - int& namax = st.namax; - int& nwmax = st.nwmax; + int ntype = 0; ///< number of atom species in UnitCell + int nat = 0; ///< total number of atoms of all species in unitcell + std::vector iat2it; ///< iat==>it, distinguish a atom belong to which type + std::vector iat2ia; ///< iat==>ia + std::vector iwt2iat; ///< iwt ==> iat + std::vector iwt2iw; ///< iwt ==> iw (Peize Lin add 2018-07-02) + ModuleBase::IntArray itia2iat; ///< (it, ia)==>iat, the index in nat (add 2009-3-2 by mohan) + int namax = 0; ///< the max na among all atom species + int nwmax = 0; ///< the max nw among all atom species /// @} Atom* atoms = nullptr; diff --git a/source/source_cell/unitcell_data.h b/source/source_cell/unitcell_data.h index 5e185773239..5fc98a74540 100644 --- a/source/source_cell/unitcell_data.h +++ b/source/source_cell/unitcell_data.h @@ -34,42 +34,4 @@ struct Lattice ModuleBase::Matrix3 invGGT = ModuleBase::Matrix3(); ///< inverse G }; -/** - * @brief Statistics data and index maps. - * - * Relationships between indices: - * - ntype, it: atom type index - * - nat, iat: total atom index - * - atoms[it].na, ia: atom index within type - * - atoms[it].nw, iw: orbital index within atom - * - * - if know it ==> atoms[it].na; atoms[it].nw - * - if know iat ==> it; ia - * - if know ia, must have known it ==> iat - * - if know iwt, must have known it, ia ==> iwt - */ -struct Statistics -{ - int ntype = 0; ///< number of atom species in UnitCell - int nat = 0; ///< total number of atoms of all species in unitcell - int* iat2it = nullptr; ///< iat==>it, distinguish a atom belong to which type - int* iat2ia = nullptr; ///< iat==>ia - int* iwt2iat = nullptr; ///< iwt ==> iat - int* iwt2iw = nullptr; ///< iwt ==> iw (Peize Lin add 2018-07-02) - ModuleBase::IntArray itia2iat; ///< (it, ia)==>iat, the index in nat (add 2009-3-2 by mohan) - int namax = 0; ///< the max na among all atom species - int nwmax = 0; ///< the max nw among all atom species - - /** - * @brief Destructor. - */ - ~Statistics() - { - delete[] iat2it; - delete[] iat2ia; - delete[] iwt2iat; - delete[] iwt2iw; - } -}; - #endif diff --git a/source/source_cell/unitcell_index.cpp b/source/source_cell/unitcell_index.cpp index edc05240bbf..678416acb80 100644 --- a/source/source_cell/unitcell_index.cpp +++ b/source/source_cell/unitcell_index.cpp @@ -8,10 +8,8 @@ void UnitCell::set_iat2itia() { assert(nat > 0); - delete[] iat2it; - delete[] iat2ia; - this->iat2it = new int[nat]; - this->iat2ia = new int[nat]; + this->iat2it.resize(nat); + this->iat2ia.resize(nat); int iat = 0; for (int it = 0; it < ntype; it++) { for (int ia = 0; ia < atoms[it].na; ia++) { diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index c73f823c226..f6aa0ec8aca 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/source/source_esolver/esolver_fp.cpp @@ -119,7 +119,9 @@ void ESolver_FP::before_all_runners(BaseCell& basecell, const Input_para& inp) if (ModuleSymmetry::Symmetry::symm_flag == 1) { const int cal_symm_repr[2] = {this->inp_->cal_symm_repr[0], this->inp_->cal_symm_repr[1]}; - ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running, + ucell.symm.analy_sys(ucell.lat, ucell.atoms, + ucell.nat, ucell.ntype, ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + GlobalV::ofs_running, this->inp_->symmetry_prec, inp.nspin, this->inp_->calculation, cal_symm_repr); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); } @@ -213,7 +215,9 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) if (update_symmetry) { const int cal_symm_repr[2] = {this->inp_->cal_symm_repr[0], this->inp_->cal_symm_repr[1]}; - ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running, + ucell.symm.analy_sys(ucell.lat, ucell.atoms, + ucell.nat, ucell.ntype, ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + GlobalV::ofs_running, this->inp_->symmetry_prec, this->inp_->nspin, this->inp_->calculation, cal_symm_repr); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); } diff --git a/source/source_esolver/esolver_gets.cpp b/source/source_esolver/esolver_gets.cpp index 82000ddd105..d0a2d727587 100644 --- a/source/source_esolver/esolver_gets.cpp +++ b/source/source_esolver/esolver_gets.cpp @@ -69,8 +69,12 @@ void ESolver_GetS::before_all_runners(BaseCell& basecell, const Input_para& inp) { const int cal_symm_repr[2] = {this->inp_->cal_symm_repr[0], this->inp_->cal_symm_repr[1]}; ucell.symm.analy_sys(ucell.lat, - ucell.st, ucell.atoms, + ucell.nat, + ucell.ntype, + ucell.iat2it, + ucell.iat2ia, + ucell.itia2iat, GlobalV::ofs_running, this->inp_->symmetry_prec, inp.nspin, diff --git a/source/source_esolver/esolver_lr_lcao_tddft.cpp b/source/source_esolver/esolver_lr_lcao_tddft.cpp index e754ab02041..bfb9bf2bb00 100644 --- a/source/source_esolver/esolver_lr_lcao_tddft.cpp +++ b/source/source_esolver/esolver_lr_lcao_tddft.cpp @@ -374,7 +374,9 @@ void ModuleESolver::ESolver_LR::initialize_from_unitcell_(UnitCell& ucell if (ModuleSymmetry::Symmetry::symm_flag == 1) { const int cal_symm_repr[2] = {this->inp_->cal_symm_repr[0], this->inp_->cal_symm_repr[1]}; - ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running, + ucell.symm.analy_sys(ucell.lat, ucell.atoms, + ucell.nat, ucell.ntype, ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + GlobalV::ofs_running, this->inp_->symmetry_prec, this->inp_->nspin, this->inp_->calculation, cal_symm_repr); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); } diff --git a/source/source_esolver/test/esolver_dp_test.cpp b/source/source_esolver/test/esolver_dp_test.cpp index a60b40a0783..4bd3c445015 100644 --- a/source/source_esolver/test/esolver_dp_test.cpp +++ b/source/source_esolver/test/esolver_dp_test.cpp @@ -26,10 +26,10 @@ class ESolverDPTest : public ::testing::Test { // Initialize variables before each test esolver = new ModuleESolver::ESolver_DP("./support/case_1.pb"); - ucell.iat2it = new int[2]; + ucell.iat2it.resize(2); ucell.iat2it[0] = 0; ucell.iat2it[1] = 1; - ucell.iat2ia = new int[2]; + ucell.iat2ia.resize(2); ucell.iat2ia[0] = 0; ucell.iat2ia[1] = 1; ucell.nat = 2; diff --git a/source/source_estate/module_dm/unittests/prepare_unitcell.h b/source/source_estate/module_dm/unittests/prepare_unitcell.h index 43e00e8f765..8faec588075 100644 --- a/source/source_estate/module_dm/unittests/prepare_unitcell.h +++ b/source/source_estate/module_dm/unittests/prepare_unitcell.h @@ -223,14 +223,10 @@ class UcellTestPrepare ucell.namax = namax; // cal_index assert(nlocal > 0); - delete[] ucell.iwt2iat; - delete[] ucell.iwt2iw; - delete[] ucell.iat2it; - delete[] ucell.iat2ia; - ucell.iwt2iat = new int[nlocal]; - ucell.iwt2iw = new int[nlocal]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; // set_iat2itia + ucell.iwt2iat.resize(nlocal); + ucell.iwt2iw.resize(nlocal); + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); // set_iat2itia ucell.itia2iat.create(ucell.ntype, ucell.namax); ucell.set_iat2iwt(1); int iat = 0; diff --git a/source/source_estate/module_dm/unittests/test_cal_dm_r.cpp b/source/source_estate/module_dm/unittests/test_cal_dm_r.cpp index dbad57eb1de..9176463cb3f 100644 --- a/source/source_estate/module_dm/unittests/test_cal_dm_r.cpp +++ b/source/source_estate/module_dm/unittests/test_cal_dm_r.cpp @@ -38,8 +38,8 @@ class DMTest : public testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_estate/module_dm/unittests/test_dm_constructor.cpp b/source/source_estate/module_dm/unittests/test_dm_constructor.cpp index f270439fe6f..a6b0062dd02 100644 --- a/source/source_estate/module_dm/unittests/test_dm_constructor.cpp +++ b/source/source_estate/module_dm/unittests/test_dm_constructor.cpp @@ -37,8 +37,8 @@ class DMTest : public testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_estate/module_dm/unittests/test_dm_from_psi.cpp b/source/source_estate/module_dm/unittests/test_dm_from_psi.cpp index 3988193ff72..49ac90cf4a2 100644 --- a/source/source_estate/module_dm/unittests/test_dm_from_psi.cpp +++ b/source/source_estate/module_dm/unittests/test_dm_from_psi.cpp @@ -38,8 +38,8 @@ class DMTest : public testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau = new ModuleBase::Vector3[ucell.nat]; ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_estate/module_dm/unittests/test_dm_r_init.cpp b/source/source_estate/module_dm/unittests/test_dm_r_init.cpp index 8c1a6b988e5..bce2a29a776 100644 --- a/source/source_estate/module_dm/unittests/test_dm_r_init.cpp +++ b/source/source_estate/module_dm/unittests/test_dm_r_init.cpp @@ -38,8 +38,8 @@ class DMTest : public testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_estate/test/test_occ_mixer.cpp b/source/source_estate/test/test_occ_mixer.cpp index 33f9a36e4e7..2d2f04e0a87 100644 --- a/source/source_estate/test/test_occ_mixer.cpp +++ b/source/source_estate/test/test_occ_mixer.cpp @@ -45,11 +45,9 @@ class OccMatMixerTest : public ::testing::Test cell.ntype = 1; cell.nat = 2; atoms_storage.resize(cell.ntype); - // UnitCell exposes raw pointers; Statistics::~Statistics() owns and - // delete[]s iat2it/iat2ia, so they must come from new[] here. cell.atoms = atoms_storage.data(); - cell.iat2it = new int[cell.nat]; - cell.iat2ia = new int[cell.nat]; + cell.iat2it.resize(cell.nat); + cell.iat2ia.resize(cell.nat); cell.itia2iat.create(cell.ntype, cell.nat); cell.atoms[0].na = 2; cell.atoms[0].nwl = 2; // max angular momentum present diff --git a/source/source_hamilt/module_gint/kernel/gint_gpu_vars.cpp b/source/source_hamilt/module_gint/kernel/gint_gpu_vars.cpp index 1b893b36c40..beba4946e99 100644 --- a/source/source_hamilt/module_gint/kernel/gint_gpu_vars.cpp +++ b/source/source_hamilt/module_gint/kernel/gint_gpu_vars.cpp @@ -96,7 +96,7 @@ GintGpuVars::GintGpuVars(std::shared_ptr biggrid_info, CHECK_CUDA(cudaMemcpy(mgrids_pos_d, mgrids_pos_h.data(), sizeof(double3) * mgrid_num, cudaMemcpyHostToDevice)); CHECK_CUDA(cudaMalloc((void**)&iat2it_d, sizeof(int) * ucell.nat)); - CHECK_CUDA(cudaMemcpy(iat2it_d, ucell.iat2it, sizeof(int) * ucell.nat, cudaMemcpyHostToDevice)); + CHECK_CUDA(cudaMemcpy(iat2it_d, ucell.iat2it.data(), sizeof(int) * ucell.nat, cudaMemcpyHostToDevice)); } GintGpuVars::~GintGpuVars() diff --git a/source/source_hamilt/module_gint/test/test_gint_common.cpp b/source/source_hamilt/module_gint/test/test_gint_common.cpp index 4d6a21607ca..478130d0e30 100644 --- a/source/source_hamilt/module_gint/test/test_gint_common.cpp +++ b/source/source_hamilt/module_gint/test/test_gint_common.cpp @@ -13,8 +13,8 @@ class GintCommonTest : public ::testing::Test ucell_.ntype = 1; ucell_.nat = 2; ucell_.atoms = new Atom[ucell_.ntype]; - ucell_.iat2it = new int[ucell_.nat]; - ucell_.iat2ia = new int[ucell_.nat]; + ucell_.iat2it.resize(ucell_.nat); + ucell_.iat2ia.resize(ucell_.nat); ucell_.atoms[0].na = 2; ucell_.atoms[0].nw = 2; diff --git a/source/source_hamilt/module_hcontainer/test/prepare_unitcell.h b/source/source_hamilt/module_hcontainer/test/prepare_unitcell.h index 834a177bef8..e8123a074dd 100644 --- a/source/source_hamilt/module_hcontainer/test/prepare_unitcell.h +++ b/source/source_hamilt/module_hcontainer/test/prepare_unitcell.h @@ -223,14 +223,10 @@ class UcellTestPrepare ucell.namax = namax; // cal_index assert(nlocal > 0); - delete[] ucell.iwt2iat; - delete[] ucell.iwt2iw; - delete[] ucell.iat2it; - delete[] ucell.iat2ia; - ucell.iwt2iat = new int[nlocal]; - ucell.iwt2iw = new int[nlocal]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; // set_iat2itia + ucell.iwt2iat.resize(nlocal); + ucell.iwt2iw.resize(nlocal); + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); // set_iat2itia ucell.itia2iat.create(ucell.ntype, ucell.namax); ucell.set_iat2iwt(1); int iat = 0; diff --git a/source/source_hamilt/module_hcontainer/test/test_func_folding.cpp b/source/source_hamilt/module_hcontainer/test/test_func_folding.cpp index 12f5c9b987e..3ee04d5e360 100644 --- a/source/source_hamilt/module_hcontainer/test/test_func_folding.cpp +++ b/source/source_hamilt/module_hcontainer/test/test_func_folding.cpp @@ -25,7 +25,7 @@ class FoldingTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) { ucell.iat2it[iat] = 0; diff --git a/source/source_hamilt/module_hcontainer/test/test_hcontainer.cpp b/source/source_hamilt/module_hcontainer/test/test_hcontainer.cpp index 172cc088e5b..362b2a95458 100644 --- a/source/source_hamilt/module_hcontainer/test/test_hcontainer.cpp +++ b/source/source_hamilt/module_hcontainer/test/test_hcontainer.cpp @@ -25,7 +25,7 @@ class HContainerTest : public ::testing::Test ucell.ntype = 1; ucell.nat = 3; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) { ucell.iat2it[iat] = 0; diff --git a/source/source_hamilt/module_hcontainer/test/test_hcontainer_complex.cpp b/source/source_hamilt/module_hcontainer/test/test_hcontainer_complex.cpp index c7920e202dd..d0cf0ed52f0 100644 --- a/source/source_hamilt/module_hcontainer/test/test_hcontainer_complex.cpp +++ b/source/source_hamilt/module_hcontainer/test/test_hcontainer_complex.cpp @@ -25,7 +25,7 @@ class HContainerTest : public ::testing::Test ucell.ntype = 1; ucell.nat = 3; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) { ucell.iat2it[iat] = 0; diff --git a/source/source_hamilt/module_hcontainer/test/test_hcontainer_output.cpp b/source/source_hamilt/module_hcontainer/test/test_hcontainer_output.cpp index 443841ffc54..d37ac7e32c1 100644 --- a/source/source_hamilt/module_hcontainer/test/test_hcontainer_output.cpp +++ b/source/source_hamilt/module_hcontainer/test/test_hcontainer_output.cpp @@ -32,8 +32,8 @@ class OutputHContainerTest : public testing::Test ucell.ntype = 1; ucell.nat = 2; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) { ucell.iat2ia[iat] = iat; @@ -41,8 +41,8 @@ class OutputHContainerTest : public testing::Test } ucell.atoms[0].na = 2; ucell.atoms[0].nw = 2; - ucell.iwt2iat = new int[4]; - ucell.iwt2iw = new int[4]; + ucell.iwt2iat.resize(4); + ucell.iwt2iw.resize(4); ucell.itia2iat.create(ucell.ntype, ucell.nat); ucell.set_iat2iwt(1); ucell.itia2iat(0, 0) = 0; diff --git a/source/source_hamilt/module_hcontainer/test/test_hcontainer_time.cpp b/source/source_hamilt/module_hcontainer/test/test_hcontainer_time.cpp index 0ac84c718ba..6afe0976cb7 100644 --- a/source/source_hamilt/module_hcontainer/test/test_hcontainer_time.cpp +++ b/source/source_hamilt/module_hcontainer/test/test_hcontainer_time.cpp @@ -24,7 +24,7 @@ class HContainerTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) { ucell.iat2it[iat] = 0; diff --git a/source/source_hamilt/module_hcontainer/test/test_transfer.cpp b/source/source_hamilt/module_hcontainer/test/test_transfer.cpp index 385df86b6c1..26500c6d5fd 100644 --- a/source/source_hamilt/module_hcontainer/test/test_transfer.cpp +++ b/source/source_hamilt/module_hcontainer/test/test_transfer.cpp @@ -27,8 +27,8 @@ class TransferTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) { diff --git a/source/source_io/module_dos/unittests/test_write_pdos_text.cpp b/source/source_io/module_dos/unittests/test_write_pdos_text.cpp index 2ac56f31727..2ab6f9009b7 100644 --- a/source/source_io/module_dos/unittests/test_write_pdos_text.cpp +++ b/source/source_io/module_dos/unittests/test_write_pdos_text.cpp @@ -32,8 +32,8 @@ static UnitCell* make_si_ucell() atom.l_nchi = {1, 1}; // iat2it / iat2ia - ucell->iat2it = new int[1]; - ucell->iat2ia = new int[1]; + ucell->iat2it.resize(1); + ucell->iat2ia.resize(1); ucell->iat2it[0] = 0; ucell->iat2ia[0] = 0; @@ -57,8 +57,6 @@ class WritePdosTextTest : public ::testing::Test } void TearDown() override { - // iat2it/iat2ia are owned by UnitCell's internal Statistics member, - // whose destructor releases them; do not delete them here. delete ucell; } diff --git a/source/source_io/module_hs/unittests/test_hsr_writer.cpp b/source/source_io/module_hs/unittests/test_hsr_writer.cpp index b2f0079681f..a5f86e7c74e 100644 --- a/source/source_io/module_hs/unittests/test_hsr_writer.cpp +++ b/source/source_io/module_hs/unittests/test_hsr_writer.cpp @@ -421,7 +421,7 @@ TEST(HsrWriterIo, HContainerBinaryMpiGatherWritesCompleteFiles) ucell.set_atom_flag = true; ucell.atoms[0].na = 1; ucell.atoms[0].nw = 2; - ucell.iat2it = new int[1]; + ucell.iat2it.resize(1); ucell.iat2it[0] = 0; ucell.set_iat2iwt(1); const int* iat2iwt = ucell.get_iat2iwt(); @@ -506,8 +506,8 @@ TEST(HsrWriterIo, WriteHsrTextCsrCarriesPerSpinFermiFromEfermi) init_unitcell(ucell); // The __MPI path in write_hsr calls set_atomic_trace and gatherParallels, // which require a real atom-to-orbital map; nullptr with nat=0 throws on - // that path. iat2it is released by UnitCell's Statistics member. - ucell.iat2it = new int[1]; + // that path. + ucell.iat2it.resize(1); ucell.iat2it[0] = 0; ucell.set_iat2iwt(1); const int* iat2iwt_ptr = ucell.get_iat2iwt(); diff --git a/source/source_io/test/restart_exx_csr_test.cpp b/source/source_io/test/restart_exx_csr_test.cpp index eee5fd96e50..4da9b0cef4e 100644 --- a/source/source_io/test/restart_exx_csr_test.cpp +++ b/source/source_io/test/restart_exx_csr_test.cpp @@ -15,8 +15,8 @@ void init_unitcell_for_ri(UnitCell& ucell) ucell.atoms[0].na = 1; ucell.atoms[0].nw = 2; ucell.atoms[0].stapos_wf = 0; - ucell.iat2it = new int[1]{0}; - ucell.iat2ia = new int[1]{0}; + ucell.iat2it = {0}; + ucell.iat2ia = {0}; } } // namespace diff --git a/source/source_lcao/module_dftu/dftu_nao_occ.cpp b/source/source_lcao/module_dftu/dftu_nao_occ.cpp index e535ed59ec2..47ed453fb20 100644 --- a/source/source_lcao/module_dftu/dftu_nao_occ.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_occ.cpp @@ -362,9 +362,12 @@ void cal_occ_mat_k(const Parallel_Orbitals* pv, if (dftu_spacegroup_symmetry && !dftu_occ_symrot_built) { const std::array period{ kv.nmp[0], kv.nmp[1], kv.nmp[2] }; - dftu_occ_symrot.find_irreducible_sector(ucell.symm, ucell.atoms, ucell.st, + dftu_occ_symrot.find_irred_sector(ucell.symm, ucell.atoms, + ucell.nat, ucell.ntype, ucell.iat2it, ucell.iat2ia, ModuleSymmetry::Symmetry_rotation_k::get_bvk_cells(period), period, ucell.lat); - dftu_occ_symrot.cal_Ms(kv, ucell, *pv, nspin); + dftu_occ_symrot.cal_Ms(kv, ucell.symm, ucell.atoms, + ucell.latvec, ucell.lmax, ucell.nat, ucell.iat2it, ucell.iat2ia, + *pv, nspin); dftu_occ_symrot_built = true; } diff --git a/source/source_lcao/module_dftu/dftu_nao_op.cpp b/source/source_lcao/module_dftu/dftu_nao_op.cpp index 022672f386d..b670df11e87 100644 --- a/source/source_lcao/module_dftu/dftu_nao_op.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_op.cpp @@ -141,9 +141,13 @@ void hamilt::DFTU_onsite>::contributeHR() { const std::array period{ this->kv_->nmp[0], this->kv_->nmp[1], this->kv_->nmp[2] }; // for return_lattice to calculate Ms - this->symrot_.find_irreducible_sector(this->ucell->symm, this->ucell->atoms, this->ucell->st, + this->symrot_.find_irred_sector(this->ucell->symm, this->ucell->atoms, + this->ucell->nat, this->ucell->ntype, this->ucell->iat2it, this->ucell->iat2ia, ModuleSymmetry::Symmetry_rotation_k::get_bvk_cells(period), period, this->ucell->lat); - this->symrot_.cal_Ms(*this->kv_, *this->ucell, *pv, this->nspin); + this->symrot_.cal_Ms(*this->kv_, this->ucell->symm, this->ucell->atoms, + this->ucell->latvec, this->ucell->lmax, + this->ucell->nat, this->ucell->iat2it, this->ucell->iat2ia, + *pv, this->nspin); this->symrot_built_ = true; } const int nspin0 = (this->nspin == 2) ? 2 : 1; diff --git a/source/source_lcao/module_dftu/unittests/test_dftu_nao_fs_accum.cpp b/source/source_lcao/module_dftu/unittests/test_dftu_nao_fs_accum.cpp index f5f01e828a0..ac1a7ba014f 100644 --- a/source/source_lcao/module_dftu/unittests/test_dftu_nao_fs_accum.cpp +++ b/source/source_lcao/module_dftu/unittests/test_dftu_nao_fs_accum.cpp @@ -33,8 +33,7 @@ class DftuFsAccumTest : public ::testing::Test pv.set_serial(3, 3); // iwt2iat: orbital 0,1 -> atom 0; orbital 2 -> atom 1 - // Must use new[] because Statistics destructor will delete[] it. - ucell.iwt2iat = new int[3]{0, 0, 1}; + ucell.iwt2iat = {0, 0, 1}; } Parallel_Orbitals pv; diff --git a/source/source_lcao/module_dftu/unittests/test_dftu_nao_ijr.cpp b/source/source_lcao/module_dftu/unittests/test_dftu_nao_ijr.cpp index d70848daf18..0083c0d7dfb 100644 --- a/source/source_lcao/module_dftu/unittests/test_dftu_nao_ijr.cpp +++ b/source/source_lcao/module_dftu/unittests/test_dftu_nao_ijr.cpp @@ -69,16 +69,14 @@ class AccumulateHrIat0Test : public ::testing::Test ucell.ntype = 1; ucell.nat = nat; ucell.atoms = atoms_buf; - iat2it_buf.assign(nat, 0); - iat2ia_buf.resize(nat); - ucell.iat2it = iat2it_buf.data(); - ucell.iat2ia = iat2ia_buf.data(); + ucell.iat2it.assign(nat, 0); + ucell.iat2ia.resize(nat); ucell.atoms[0].tau.resize(nat); ucell.lat0 = 1.0; ucell.itia2iat.create(ucell.ntype, nat); for (int iat = 0; iat < nat; iat++) { - iat2ia_buf[iat] = iat; + ucell.iat2ia[iat] = iat; ucell.atoms[0].tau[iat] = ModuleBase::Vector3(0.0, 0.0, 0.0); ucell.itia2iat(0, iat) = iat; } @@ -102,8 +100,6 @@ class AccumulateHrIat0Test : public ::testing::Test HR.reset(); paraV.reset(); ucell.atoms = nullptr; - ucell.iat2it = nullptr; - ucell.iat2ia = nullptr; } AdjacentAtomInfo make_adjs() const @@ -142,8 +138,6 @@ class AccumulateHrIat0Test : public ::testing::Test int my_rank = 0; UnitCell ucell; Atom atoms_buf[1]; - std::vector iat2it_buf; - std::vector iat2ia_buf; std::unique_ptr paraV; std::unique_ptr> HR; }; diff --git a/source/source_lcao/module_dftu/unittests/test_dftu_nao_op.cpp b/source/source_lcao/module_dftu/unittests/test_dftu_nao_op.cpp index 0a87a8a1578..132131147c1 100644 --- a/source/source_lcao/module_dftu/unittests/test_dftu_nao_op.cpp +++ b/source/source_lcao/module_dftu/unittests/test_dftu_nao_op.cpp @@ -31,8 +31,8 @@ class DFTUTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.lat0 = 1.0; ucell.itia2iat.create(ucell.ntype, ucell.nat); diff --git a/source/source_lcao/module_lr/ri_benchmark/test/ri_benchmark_test.cpp b/source/source_lcao/module_lr/ri_benchmark/test/ri_benchmark_test.cpp index 2a23cc7b420..407697dd869 100644 --- a/source/source_lcao/module_lr/ri_benchmark/test/ri_benchmark_test.cpp +++ b/source/source_lcao/module_lr/ri_benchmark/test/ri_benchmark_test.cpp @@ -20,7 +20,7 @@ Sep_Cell::Sep_Cell() noexcept {} Sep_Cell::~Sep_Cell() noexcept {} UnitCell::UnitCell() { atoms = new Atom[1]; - iat2it = new int[1]; iat2it[0] = 0; + iat2it = {0}; iat2iwt.resize(1, 0); } UnitCell::~UnitCell() { diff --git a/source/source_lcao/module_operator_lcao/test/test_ekinetic.cpp b/source/source_lcao/module_operator_lcao/test/test_ekinetic.cpp index 9edc265ccbf..d241ca82e03 100644 --- a/source/source_lcao/module_operator_lcao/test/test_ekinetic.cpp +++ b/source/source_lcao/module_operator_lcao/test/test_ekinetic.cpp @@ -33,8 +33,8 @@ class EKineticTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_lcao/module_operator_lcao/test/test_ekinetic_serial.cpp b/source/source_lcao/module_operator_lcao/test/test_ekinetic_serial.cpp index 87ce65d4342..9d035bd4bd0 100644 --- a/source/source_lcao/module_operator_lcao/test/test_ekinetic_serial.cpp +++ b/source/source_lcao/module_operator_lcao/test/test_ekinetic_serial.cpp @@ -33,8 +33,8 @@ class EKineticTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) @@ -191,8 +191,8 @@ TEST_F(EKineticTest, singleAtom) ucell_single.ntype = 1; ucell_single.nat = 1; ucell_single.atoms = new Atom[1]; - ucell_single.iat2it = new int[1]; - ucell_single.iat2ia = new int[1]; + ucell_single.iat2it.resize(1); + ucell_single.iat2ia.resize(1); ucell_single.atoms[0].tau.resize(1); ucell_single.itia2iat.create(1, 1); ucell_single.iat2it[0] = 0; diff --git a/source/source_lcao/module_operator_lcao/test/test_nonlocal.cpp b/source/source_lcao/module_operator_lcao/test/test_nonlocal.cpp index 5eea2ba121a..d604dab5f77 100644 --- a/source/source_lcao/module_operator_lcao/test/test_nonlocal.cpp +++ b/source/source_lcao/module_operator_lcao/test/test_nonlocal.cpp @@ -38,8 +38,8 @@ class NonlocalTest : public ::testing::Test ucell.infoNL.reset(lcao_nl); ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_lcao/module_operator_lcao/test/test_overlap.cpp b/source/source_lcao/module_operator_lcao/test/test_overlap.cpp index e3c017857b2..aa6a14ae5a0 100644 --- a/source/source_lcao/module_operator_lcao/test/test_overlap.cpp +++ b/source/source_lcao/module_operator_lcao/test/test_overlap.cpp @@ -33,8 +33,8 @@ class OverlapTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_lcao/module_operator_lcao/test/test_overlap_cd.cpp b/source/source_lcao/module_operator_lcao/test/test_overlap_cd.cpp index ad2e8daf7e1..eca29b00dda 100644 --- a/source/source_lcao/module_operator_lcao/test/test_overlap_cd.cpp +++ b/source/source_lcao/module_operator_lcao/test/test_overlap_cd.cpp @@ -32,8 +32,8 @@ class OverlapTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_lcao/module_operator_lcao/test/test_overlap_serial.cpp b/source/source_lcao/module_operator_lcao/test/test_overlap_serial.cpp index 98b8e4d7419..9a69b82bb7f 100644 --- a/source/source_lcao/module_operator_lcao/test/test_overlap_serial.cpp +++ b/source/source_lcao/module_operator_lcao/test/test_overlap_serial.cpp @@ -33,8 +33,8 @@ class OverlapTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) @@ -198,8 +198,8 @@ TEST_F(OverlapTest, singleAtom) ucell_single.ntype = 1; ucell_single.nat = 1; ucell_single.atoms = new Atom[1]; - ucell_single.iat2it = new int[1]; - ucell_single.iat2ia = new int[1]; + ucell_single.iat2it.resize(1); + ucell_single.iat2ia.resize(1); ucell_single.atoms[0].tau.resize(1); ucell_single.itia2iat.create(1, 1); ucell_single.iat2it[0] = 0; diff --git a/source/source_lcao/module_operator_lcao/test/test_t_nl_cd.cpp b/source/source_lcao/module_operator_lcao/test/test_t_nl_cd.cpp index 5b9e6c1c572..a83235b2816 100644 --- a/source/source_lcao/module_operator_lcao/test/test_t_nl_cd.cpp +++ b/source/source_lcao/module_operator_lcao/test/test_t_nl_cd.cpp @@ -35,8 +35,8 @@ class TNLTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_lcao/module_rdmft/rdmft.cpp b/source/source_lcao/module_rdmft/rdmft.cpp index 22883e46d6d..5567140ebc2 100644 --- a/source/source_lcao/module_rdmft/rdmft.cpp +++ b/source/source_lcao/module_rdmft/rdmft.cpp @@ -157,9 +157,13 @@ void RDMFT::init(Parallel_Orbitals& ParaV_in, if (exx_spacegroup_symmetry) { const std::array& period = RI_Util::get_Born_vonKarmen_period(*kv); - this->symrot_exx.find_irreducible_sector(ucell->symm, ucell->atoms, ucell->st, + this->symrot_exx.find_irred_sector(ucell->symm, ucell->atoms, + ucell->nat, ucell->ntype, ucell->iat2it, ucell->iat2ia, RI_Util::get_Born_von_Karmen_cells(period), period, ucell->lat, PARAM.globalv.global_out_dir); - this->symrot_exx.cal_Ms(*kv, *ucell, *ParaV, nspin); + this->symrot_exx.cal_Ms(*kv, ucell->symm, ucell->atoms, + ucell->latvec, ucell->lmax, + ucell->nat, ucell->iat2it, ucell->iat2ia, + *ParaV, nspin); } if (this->exx_info_->info_ri.real_number) diff --git a/source/source_lcao/module_ri/exx_lri.hpp b/source/source_lcao/module_ri/exx_lri.hpp index c89336ba5cb..2ba07c692f0 100644 --- a/source/source_lcao/module_ri/exx_lri.hpp +++ b/source/source_lcao/module_ri/exx_lri.hpp @@ -7,6 +7,7 @@ #define EXX_LRI_HPP #include "exx_lri.h" +#include "source_cell/unitcell.h" #include "ri_2d_comm.h" #include "ri_util.h" #include "exx_lri_detail.h" @@ -824,7 +825,7 @@ void Exx_LRI::cal_exx_elec(const std::vector, std::set>> judge = RI_2D_Comm::get_2D_judge(ucell,pv); if(p_symrot) - { this->exx_lri.set_symmetry(true, p_symrot->get_irreducible_sector()); } + { this->exx_lri.set_symmetry(true, p_symrot->get_irred_sector()); } else { this->exx_lri.set_symmetry(false, {}); } @@ -857,7 +858,7 @@ void Exx_LRI::cal_exx_elec(const std::vectorexx_lri.post_2D.set_tensors_map2(this->exx_lri.Hs); // rotate locally without repeat - Hs_a2D = p_symrot->restore_HR(ucell.symm, ucell.atoms, ucell.st, 'H', Hs_a2D); + Hs_a2D = p_symrot->restore_HR(ucell.symm, ucell.atoms, ucell.iat2it, 'H', Hs_a2D); // cal energy using full Hs without repeat this->exx_lri.energy = this->exx_lri.post_2D.cal_energy( this->exx_lri.post_2D.saves["Ds_" + suffix], @@ -899,7 +900,7 @@ void Exx_LRI::cal_exx_elec_soc( // pass 2: spinor-coupled rotation of the 4 channels from the irreducible sector to the full BZ std::array>>, 4> Hs_full = - p_symrot->restore_HR_nspin4(ucell.symm, ucell.atoms, ucell.st, 'H', Hs_irr); + p_symrot->restore_HR_nspin4(ucell.symm, ucell.atoms, ucell.iat2it, 'H', Hs_irr); // pass 3: per-channel energy (full Hs, no repeat), then gather the repeated full Hs for abacus for (int is = 0; is < 4; ++is) diff --git a/source/source_lcao/module_ri/exx_lri_interface.hpp b/source/source_lcao/module_ri/exx_lri_interface.hpp index 13fe79697d5..b08514faad0 100644 --- a/source/source_lcao/module_ri/exx_lri_interface.hpp +++ b/source/source_lcao/module_ri/exx_lri_interface.hpp @@ -1,6 +1,7 @@ #ifndef EXX_LRI_INTERFACE_HPP #define EXX_LRI_INTERFACE_HPP #include "exx_lri_interface.h" +#include "source_cell/unitcell.h" #include "source_base/formatter.h" #include "source_base/parallel_common.h" #include "source_estate/elecstate_lcao.h" @@ -99,11 +100,13 @@ void Exx_LRI_Interface::refresh_symmetry( if (this->exx_spacegroup_symmetry) { const std::array& period = RI_Util::get_Born_vonKarmen_period(kv); - this->symrot_.find_irreducible_sector( - ucell.symm, ucell.atoms, ucell.st, + this->symrot_.find_irred_sector( + ucell.symm, ucell.atoms, + ucell.nat, ucell.ntype, ucell.iat2it, ucell.iat2ia, RI_Util::get_Born_von_Karmen_cells(period), period, ucell.lat, PARAM.globalv.global_out_dir); this->symrot_.set_abfs_Lmax(Exx_Abfs::Construct_Orbs::get_Lmax(this->exx_ptr->abfs)); - this->symrot_.cal_Ms(kv, ucell, pv, PARAM.inp.nspin); + this->symrot_.cal_Ms(kv, ucell.symm, ucell.atoms, ucell.latvec, ucell.lmax, + ucell.nat, ucell.iat2it, ucell.iat2ia, pv, PARAM.inp.nspin); } } diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rot_out.cpp b/source/source_lcao/module_ri/module_exx_symmetry/symm_rot_out.cpp index 0ff3c9be33b..bc77bb97acc 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rot_out.cpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rot_out.cpp @@ -1,6 +1,6 @@ #include "./symm_rotation.h" +#include "source_cell/unitcell.h" #include "source_io/module_parameter/parameter.h" - namespace ModuleSymmetry { std::string mat3_fmt(const ModuleBase::Matrix3& m) @@ -19,7 +19,7 @@ namespace ModuleSymmetry ModuleBase::TITLE("ModuleSymmetry", "print_symrot_info_R"); std::ofstream ofs(PARAM.globalv.global_out_dir + "symrot_R.txt"); // Print the irreducible sector (to be optimized) - ofs << "Number of irreducible sector: " << symrot.get_irreducible_sector().size() << std::endl; + ofs << "Number of irreducible sector: " << symrot.get_irred_sector().size() << std::endl; ofs << "Lmax of AOs: " << lmax_ao << "\n"; ofs << "Lmax of ABFs: " << symrot.abfs_Lmax << "\n"; // print AO rotation matrix T diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp index b8261ad6fe1..d7bc618ceb6 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp @@ -1,13 +1,6 @@ #include "symm_rotation.h" -#include "source_base/constants.h" -#include "source_io/module_parameter/parameter.h" -#include -#include "source_base/parallel_reduce.h" -#include "source_base/module_external/scalapack_connector.h" -#include "source_base/tool_title.h" -#include "source_base/timer.h" -#include "source_base/mathzone.h" -#include "source_lcao/module_ri/ri_util.h" + +#include namespace ModuleSymmetry { @@ -18,44 +11,4 @@ namespace ModuleSymmetry for (auto& abfs_T : abfs_l_nchi) { this->abfs_Lmax_ = std::max(this->abfs_Lmax_, static_cast(abfs_T.size()) - 1); } } - std::vector Symmetry_rotation::get_Rs_from_adjacent_list(const UnitCell& ucell, - const Grid_Driver& gd, - const Parallel_Orbitals& pv) const - { - // find the union set of Rs for all the atom pairs - std::set Rs_set; - for (int iat1 = 0;iat1 < ucell.nat;++iat1) - { - auto tau1 = ucell.get_tau(iat1); - int it1 = ucell.iat2it[iat1], ia1 = ucell.iat2ia[iat1]; - AdjacentAtomInfo adjs; - gd.Find_atom(ucell, tau1, it1, ia1, &adjs); - for (int ad = 0; ad < adjs.adj_num + 1; ++ad) - { - const int it2 = adjs.ntype[ad]; - const int ia2 = adjs.natom[ad]; - int iat2 = ucell.itia2iat(it2, ia2); - if (pv.get_nrow_atom(iat1) && pv.get_ncol_atom(iat2)) - { - const ModuleBase::Vector3& R_index = adjs.box[ad]; - if (ucell.cal_dtau(iat1, iat2, R_index).norm() * ucell.lat0 - < ucell.atoms[it1].Rcut + ucell.atoms[it2].Rcut) { - Rs_set.insert({ R_index.x, R_index.y, R_index.z }); -} - } - } - } - // set to vector - std::vector Rs(Rs_set.size()); - for (auto& R : Rs_set) { Rs.push_back(R); -} - return Rs; - } - - std::vector Symmetry_rotation::get_Rs_from_BvK(const K_Vectors& kv) const - { - const TC& period = RI_Util::get_Born_vonKarmen_period(kv); - return RI_Util::get_Born_von_Karmen_cells(period); - } - } diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h index 16d7ba91928..2adbb5561b2 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h @@ -5,7 +5,8 @@ #include "source_basis/module_ao/parallel_orbitals.h" #include #include "source_hamilt/module_hcontainer/hcontainer.h" -#include "source_cell/module_neighbor/sltk_grid_driver.h" + +class UnitCell; namespace ModuleSymmetry { @@ -29,13 +30,17 @@ namespace ModuleSymmetry //-------------------------------------------------------------------------------- /// The main functions to rotate matrices /// Given H(R) in the irreduceble sector, calculate H(R) for all the atompairs and cells. + /// @param iat2it map from global atom index to its species index template // RI::Tensor type std::map, RI::Tensor>> restore_HR( - const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, + const Symmetry& symm, const Atom* atoms, const std::vector& iat2it, const char mode, const std::map, RI::Tensor>>& HR_irreduceble)const; + /// @param iat2it map from global atom index to its species index + /// @param iat2ia map from global atom index to its index within the species template // HContainer type void restore_HR( - const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, + const Symmetry& symm, const Atom* atoms, + const std::vector& iat2it, const std::vector& iat2ia, const char mode, const hamilt::HContainer& HR_irreduceble, hamilt::HContainer& HR_rotated)const; /// (nspin=4) spinor overload: rotate all 4 spin channels of H(R) together. On top of the /// orbital rotation T1^dagger(.)T2 (mode 'H') / T1^T(.)T2^* (mode 'D') applied to every @@ -47,31 +52,17 @@ namespace ModuleSymmetry /// remap of the 4 channels applied after the SU(2) mixing; see symm_rotation_r.hpp. template // RI::Tensor type std::array, RI::Tensor>>, 4> restore_HR_nspin4( - const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, + const Symmetry& symm, const Atom* atoms, const std::vector& iat2it, const char mode, const std::array, RI::Tensor>>, 4>& HR_irreducible_soc)const; //-------------------------------------------------------------------------------- - /// test functions - /// test H(R) rotation: giver a full H(R), pick out H(R) in the irreducible sector, rotate it, and compare with the original full H(R) - template // RI::Tensor type, using col-major implementation - void test_HR_rotation(const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, - const std::map, RI::Tensor>>& HR_full); - template // test the rotation of RI coefficients - void test_Cs_rotation(const Symmetry& symm, const Atom* atoms, const Statistics& st, - const std::map, RI::Tensor>>& Cs_full)const; - template // HContainer type, using row-major implementation - void test_HR_rotation(const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, - const hamilt::HContainer& HR_full); - template // HContainer type + template // HContainer type void print_HR(const std::map, RI::Tensor>>& HR, const std::string name, const double& threshold = 0.0); //-------------------------------------------------------------------------------- private: //-------------------------------------------------------------------------------- - std::vector get_Rs_from_BvK(const K_Vectors& kv)const; - std::vector get_Rs_from_adjacent_list(const UnitCell& ucell, - const Grid_Driver& gd, - const Parallel_Orbitals& pv) const; + //-------------------------------------------------------------------------------- /// The sub functions to rotate matrices @@ -83,8 +74,9 @@ namespace ModuleSymmetry template // pointer type, blas void rotate_atompair_serial(Tdata* TAT, const Tdata* A, const int& nw1, const int& nw2, const int isym, const Atom& a1, const Atom& a2, const char mode)const; - template // HContainer type, pblas - void rotate_atompair_parallel(const TR* Alocal_in, const int isym, const Atom* atoms, const Statistics& st, + template // HContainer type, pblas + void rotate_atompair_parallel(const TR* Alocal_in, const int isym, const Atom* atoms, + const std::vector& iat2it, const std::vector& iat2ia, const Tap& ap_in, const Tap& ap_out, const char mode, const Parallel_Orbitals& pv, TR* Alocal_out, const bool output = false)const; /// rotate a 3-dim C tensor in RI technique diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r.hpp b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r.hpp index 381aaf5f435..86a20e6364f 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r.hpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r.hpp @@ -77,7 +77,7 @@ namespace ModuleSymmetry template std::map, RI::Tensor>> Symmetry_rotation::restore_HR( - const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, + const Symmetry& symm, const Atom* atoms, const std::vector& iat2it, const char mode, const std::map, RI::Tensor>>& HR_irreducible) const { ModuleBase::TITLE("Symmetry_rotation", "restore_HR"); @@ -102,7 +102,7 @@ namespace ModuleSymmetry const int& ap1 = apR.first.first; const int& ap2 = apR.first.second; const TC& R = apR.second; - HR_full[ap1][{ap2, R}] = rotate_atompair_serial(tmp2.second, isym, atoms[st.iat2it[irap1]], atoms[st.iat2it[irap2]], mode); + HR_full[ap1][{ap2, R}] = rotate_atompair_serial(tmp2.second, isym, atoms[iat2it[irap1]], atoms[iat2it[irap2]], mode); } } else { std::cout << "Warning: not found: irreducible atom pair =(" << irap1 << "," << irap2 << "), irR=(" << irR[0] << "," << irR[1] << "," << irR[2] << ")\n";} @@ -121,7 +121,7 @@ namespace ModuleSymmetry // // H_12(R)=T^\dagger(V)H_1'2'(VR+O_1-O_2)T(V) // if (HR_irreducible.find(irap.first) != HR_irreducible.end() && HR_irreducible.at(irap.first).find({ irap.second, irR }) != HR_irreducible.at(irap.first).end()) // HR_full[ap.first][{ap.second, R}] = rotate_atompair_serial(HR_irreducible.at(irap.first).at({ irap.second, irR }), - // isym, atoms[st.iat2it[irap.first]], atoms[st.iat2it[irap.second]], mode); + // isym, atoms[ucell.iat2it[irap.first]], atoms[ucell.iat2it[irap.second]], mode); // else // std::cout << "not found: current atom pair =(" << ap.first << "," << ap.second << "), R=(" << R[0] << "," << R[1] << "," << R[2] << "), irreducible atom pair =(" << irap.first << "," << irap.second << "), irR=(" << irR[0] << "," << irR[1] << "," << irR[2] << ")\n"; // } @@ -171,7 +171,7 @@ namespace ModuleSymmetry // when it is applied in either direction. template std::array, RI::Tensor>>, 4> Symmetry_rotation::restore_HR_nspin4( - const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, + const Symmetry& symm, const Atom* atoms, const std::vector& iat2it, const char mode, const std::array, RI::Tensor>>, 4>& HR_irreducible_soc)const { ModuleBase::TITLE("Symmetry_rotation", "restore_HR_nspin4"); @@ -203,8 +203,8 @@ namespace ModuleSymmetry std::cout << "Warning: not found: irreducible atom pair =(" << irap1 << "," << irap2 << "), irR=(" << irR[0] << "," << irR[1] << "," << irR[2] << ")\n"; continue; } - const Atom& a1 = atoms[st.iat2it[irap1]]; - const Atom& a2 = atoms[st.iat2it[irap2]]; + const Atom& a1 = atoms[iat2it[irap1]]; + const Atom& a2 = atoms[iat2it[irap2]]; // gather the 4 irreducible spin-channel blocks (zero-filled when absent) std::array, 4> Hir; for (int is = 0;is < 4;++is) @@ -405,77 +405,4 @@ namespace ModuleSymmetry } } - template - void Symmetry_rotation::test_HR_rotation(const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, - const std::map, RI::Tensor>>& HR_full) - { - ModuleBase::TITLE("Symmetry_rotation", "test_HR_rotation"); - - // 1. pick out H(R) in the irreducible sector from full H(R) - std::map, RI::Tensor>> HR_irreducible; - for (auto& irap_Rs : this->irs_.irreducible_sector_) - { - const Tap& irap = irap_Rs.first; - for (auto& irR : irap_Rs.second) - { - const std::pair a2_irR = { irap.second, irR }; - HR_irreducible[irap.first][a2_irR] = (HR_full.at(irap.first).count(a2_irR) != 0) ? - HR_full.at(irap.first).at(a2_irR) - : RI::Tensor(HR_full.at(irap.first).begin()->second.shape); - } - } - // 2. rotate - std::map, RI::Tensor>> HR_rotated = restore_HR(symm, atoms, st, mode, HR_irreducible); - // 3. compare - for (auto& HR_ia1 : HR_rotated) - { - int iat1 = HR_ia1.first; - for (auto& HR_ia12R : HR_ia1.second) - { - int iat2 = HR_ia12R.first.first; - TC R = HR_ia12R.first.second; - const RI::Tensor& HR_rot = HR_ia12R.second; - if (HR_full.at(iat1).count({ iat2, R }) == 0)// rot back but not found - { - std::cout << "R_rot not found in atom pair (" << iat1 << ", " << iat2 << "): R=(" << R[0] << "," << R[1] << "," << R[2] << "):\n"; - continue; - } - const RI::Tensor& HR_ref = HR_full.at(iat1).at({ iat2, R }); - assert(HR_rot.shape[0] == HR_ref.shape[0]); - assert(HR_rot.shape[1] == HR_ref.shape[1]); - // output - std::cout << "atom pair (" << iat1 << ", " << iat2 << "), R=(" << R[0] << "," << R[1] << "," << R[2] << "):\n"; - print_tensor(HR_rot, std::string("R_rot").insert(0, 1, mode)); - print_tensor(HR_ref, std::string("R_ref").insert(0, 1, mode)); - } - } - } - - template - void Symmetry_rotation::test_Cs_rotation(const Symmetry& symm, const Atom* atoms, const Statistics& st, - const std::map, RI::Tensor>>& Cs_full)const - { - for (auto& sector_pair : this->irs_.full_map_to_irreducible_sector_) - { - const TapR& apR = sector_pair.first; - const int& isym = sector_pair.second.first; - const TapR& irapR = sector_pair.second.second; - // if (apR.first != irapR.first || apR.second != irapR.second) - if (apR.first != irapR.first) - { - std::cout << "irapR=(" << irapR.first.first << "," << irapR.first.second << "), (" << irapR.second[0] << "," << irapR.second[1] << "," << irapR.second[2] << "):\n"; - std::cout << "apR=(" << apR.first.first << "," << apR.first.second << "), (" << apR.second[0] << "," << apR.second[1] << "," << apR.second[2] << "):\n"; - const RI::Tensor& Cs_ir = Cs_full.at(irapR.first.first).at({ irapR.first.second,irapR.second }); - const RI::Tensor& Cs_ref = Cs_full.at(apR.first.first).at({ apR.first.second,apR.second }); - const RI::Tensor& Cs_rot = this->rotate_singleC_serial(Cs_ir, isym, - atoms[st.iat2it[irapR.first.first]], atoms[st.iat2it[irapR.first.second]], irapR.first.first); - print_tensor3(Cs_rot, "Cs_rot"); - print_tensor3(Cs_ref, "Cs_ref"); - print_tensor3(Cs_ir, "Cs_irreducible"); - exit(0); - } - } - - } - } \ No newline at end of file diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r_hcontainer.hpp b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r_hcontainer.hpp index 7e96480dabd..ad0f7274952 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r_hcontainer.hpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r_hcontainer.hpp @@ -36,7 +36,8 @@ namespace ModuleSymmetry template // HContainer type void Symmetry_rotation::restore_HR( - const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, + const Symmetry& symm, const Atom* atoms, + const std::vector& iat2it, const std::vector& iat2ia, const char mode, const hamilt::HContainer& HR_irreduceble, hamilt::HContainer& HR_rotated)const { @@ -59,7 +60,7 @@ namespace ModuleSymmetry const int R_hc = ap_hc.find_R(R[0], R[1], R[2]); if (R_hc < 0) continue; TR* ijR_ptr = ap_hc.get_pointer(R_hc); - rotate_atompair_parallel(irijR_ptr, isym, atoms, st, irap, ap, mode, *irap_hc.get_paraV(), ijR_ptr); + rotate_atompair_parallel(irijR_ptr, isym, atoms, iat2it, iat2ia, irap, ap, mode, *irap_hc.get_paraV(), ijR_ptr); } ModuleBase::timer::end("Symmetry_rotation", "restore_HR"); } @@ -77,13 +78,14 @@ namespace ModuleSymmetry obj[starti + i + (startj + j) * nr] = block(j, i); }; template - void Symmetry_rotation::rotate_atompair_parallel(const TR* Alocal_in, const int isym, const Atom* atoms, const Statistics& st, + void Symmetry_rotation::rotate_atompair_parallel(const TR* Alocal_in, const int isym, const Atom* atoms, + const std::vector& iat2it, const std::vector& iat2ia, const Tap& ap_in, const Tap& ap_out, const char mode, const Parallel_Orbitals& pv, TR* Alocal_out, const bool output)const { // all the matrices are row-major (col-contiguous) int iat1 = ap_in.first, iat2 = ap_in.second; - int it1 = st.iat2it[iat1], it2 = st.iat2it[iat2]; - int ia1 = st.iat2ia[iat1], ia2 = st.iat2ia[iat2]; + int it1 = iat2it[iat1], it2 = iat2it[iat2]; + int ia1 = iat2ia[iat1], ia2 = iat2ia[iat2]; // contruct T matrix std::vector T1, T2; auto set_rotation_matrix = [&](const int& it, const int& ia, std::vector& T)->int @@ -149,8 +151,8 @@ namespace ModuleSymmetry // copy back to Alocal_out iat1 = ap_out.first, iat2 = ap_out.second; - it1 = st.iat2it[iat1], it2 = st.iat2it[iat2]; - ia1 = st.iat2ia[iat1], ia2 = st.iat2ia[iat2]; + it1 = iat2it[iat1], it2 = iat2it[iat2]; + ia1 = iat2ia[iat1], ia2 = iat2ia[iat2]; abr1 = pv.atom_begin_row[iat1], abc2 = pv.atom_begin_col[iat2]; iw1start = atoms[it1].stapos_wf + ia1 * atoms[it1].nw; iw2start = atoms[it2].stapos_wf + ia2 * atoms[it2].nw; @@ -162,73 +164,4 @@ namespace ModuleSymmetry = TAT[(pv.local2global_row(i + abr1) - iw1start) * atoms[it2].nw + (pv.local2global_col(j + abc2) - iw2start)]; } } - - template - void Symmetry_rotation::test_HR_rotation(const Symmetry& symm, const Atom* atoms, const Statistics& st, - const char mode, const hamilt::HContainer& HR_full) - { - ModuleBase::TITLE("Symmetry_rotation", "test_HR_rotation"); - auto get_irreducible_ijR_info = [&HR_full, this]() -> std::vector - { - std::vector irreducible_ijR_info; - irreducible_ijR_info.push_back(this->irs_.irreducible_sector_.size()); - for (auto& irap_irR : this->irs_.irreducible_sector_) - { - const int iat1 = irap_irR.first.first, iat2 = irap_irR.first.second; - irreducible_ijR_info.insert(irreducible_ijR_info.end(), { iat1, iat2, 0 }); - int nR = 0; - for (auto& R : irap_irR.second) - if (HR_full.get_atom_pair(iat1, iat2).find_R(R[0], R[1], R[2]) >= 0) - { - irreducible_ijR_info.insert(irreducible_ijR_info.end(), R.begin(), R.end()); - ++nR; - } - irreducible_ijR_info[irreducible_ijR_info.size() - 3 * nR - 1] = nR; - } - return irreducible_ijR_info; - }; - - const Parallel_Orbitals* pv = HR_full.get_atom_pair(0, 0).get_paraV(); - // 1. pick out H(R) in the irreducible sector from full H(R) - const std::vector& irreducible_ijR_info = get_irreducible_ijR_info(); - hamilt::HContainer HR_irreducible(pv, nullptr, &irreducible_ijR_info); - HR_irreducible.set_zero(); - for (auto& irap_irR : this->irs_.irreducible_sector_) - { - const int iat1 = irap_irR.first.first, iat2 = irap_irR.first.second; - const hamilt::AtomPair& irap = HR_irreducible.get_atom_pair(iat1, iat2); - const hamilt::AtomPair& ap_full = HR_full.get_atom_pair(iat1, iat2); - for (auto& R : irap_irR.second) - if (irap.find_R(R[0], R[1], R[2]) >= 0) - { // out-of-range R can be added to irreducible sector to be calculated, but not in this test for lack of reference - TR* irptr = irap.get_HR_values(R[0], R[1], R[2]).get_pointer(); - TR* ptr = ap_full.get_HR_values(R[0], R[1], R[2]).get_pointer(); - std::copy(ptr, ptr + pv->get_nrow_atom(iat1) * pv->get_ncol_atom(iat2), irptr); - } - } - - //2. rotate - hamilt::HContainer HR_rotated(HR_full); - HR_rotated.set_zero(); - this->restore_HR(symm, atoms, st, mode, HR_irreducible, HR_rotated); - //3. compare - for (int iat1 = 0;iat1 < st.nat;++iat1) - { - for (int iat2 = 0;iat2 < st.nat;++iat2) - { - const hamilt::AtomPair& ap_full = HR_full.get_atom_pair(iat1, iat2); - const hamilt::AtomPair& ap_rotated = HR_rotated.get_atom_pair(iat1, iat2); - assert(ap_full.get_R_size() == ap_rotated.get_R_size()); - for (int irR = 0;irR < ap_full.get_R_size();++irR) - { - const TR* full_ptr = ap_full.get_pointer(irR); - int* R = ap_full.get_R_index(irR); - const TR* rotated_ptr = ap_rotated.get_HR_values(R[0], R[1], R[2]).get_pointer(); - GlobalV::ofs_running << "atom pair: (" << iat1 << ", " << iat2 << "), R: (" << R[0] << " " << R[1] << " " << R[2] << ")\n"; - print_atompair_local(*pv, iat1, iat2, rotated_ptr, std::string("R_rot").insert(0, 1, mode)); - print_atompair_local(*pv, iat1, iat2, full_ptr, std::string("R_ref").insert(0, 1, mode)); - } - } - } - } } \ No newline at end of file diff --git a/source/source_lcao/module_ri/module_exx_symmetry/test/symm_rotation_test.cpp b/source/source_lcao/module_ri/module_exx_symmetry/test/symm_rotation_test.cpp index ab93816a077..73ca3b68752 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/test/symm_rotation_test.cpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/test/symm_rotation_test.cpp @@ -1,5 +1,6 @@ #include "mpi.h" #include "../symm_rotation.h" +#include "source_cell/unitcell.h" #include "gtest/gtest.h" #define DOUBLETHRESHOLD 1e-8 diff --git a/source/source_lcao/module_ri/module_exx_symmetry/test/test_symm_rotation.cpp b/source/source_lcao/module_ri/module_exx_symmetry/test/test_symm_rotation.cpp index cee79ba3fb6..131d9d21150 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/test/test_symm_rotation.cpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/test/test_symm_rotation.cpp @@ -174,8 +174,8 @@ TEST(SymmetryDensityRestoration, RebuildAfterCellSymmetryAnalysis) cell.a1 = ModuleBase::Vector3(1, 0, 0); cell.a2 = ModuleBase::Vector3(0, 1, 0); cell.a3 = ModuleBase::Vector3(0, 0, 1); - cell.st.iat2it = new int[3]{0, 0, 1}; - cell.st.iat2ia = new int[3]{0, 1, 0}; + cell.iat2it = {0, 0, 1}; + cell.iat2ia = {0, 1, 0}; atoms[0].label = "A"; atoms[0].na = 2; atoms[0].nw = 1; @@ -193,7 +193,9 @@ TEST(SymmetryDensityRestoration, RebuildAfterCellSymmetryAnalysis) std::ofstream log; const int representation[2] = {0, 0}; const std::string calculation = "cell-relax"; - cell.symm.analy_sys(cell.lat, cell.st, cell.atoms, log, 1e-6, 1, calculation, representation); + cell.symm.analy_sys(cell.lat, cell.atoms, cell.nat, cell.ntype, + cell.iat2it, cell.iat2ia, cell.itia2iat, + log, 1e-6, 1, calculation, representation); const int old_operations = cell.symm.nrotk; K_Vectors kv; kv.set_nks(1); @@ -213,23 +215,31 @@ TEST(SymmetryDensityRestoration, RebuildAfterCellSymmetryAnalysis) const ModuleSymmetry::TC period = {2, 2, 2}; const std::vector cells = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}}; CellSymmetryRotation reused; - reused.find_irreducible_sector(cell.symm, cell.atoms, cell.st, cells, period, cell.lat); - reused.cal_Ms(kv, cell, pv, 1); + reused.find_irred_sector(cell.symm, cell.atoms, cell.nat, cell.ntype, + cell.iat2it, cell.iat2ia, cells, period, cell.lat); + reused.cal_Ms(kv, cell.symm, cell.atoms, cell.latvec, cell.lmax, + cell.nat, cell.iat2it, cell.iat2ia, pv, 1); const auto old_rotations = reused.rotations(); atoms[0].taud = {{0.9, 0.2, 0.3}, {0.2, 0.3, 0.2}}; atoms[1].taud = {{0.05, 0.25, 0.25}}; atoms[0].tau = atoms[0].taud; atoms[1].tau = atoms[1].taud; - cell.symm.analy_sys(cell.lat, cell.st, cell.atoms, log, 1e-6, 1, calculation, representation); + cell.symm.analy_sys(cell.lat, cell.atoms, cell.nat, cell.ntype, + cell.iat2it, cell.iat2ia, cell.itia2iat, + log, 1e-6, 1, calculation, representation); ASSERT_EQ(cell.symm.nrotk, old_operations); reused.reset_symmetry(); - reused.find_irreducible_sector(cell.symm, cell.atoms, cell.st, cells, period, cell.lat); - reused.cal_Ms(kv, cell, pv, 1); + reused.find_irred_sector(cell.symm, cell.atoms, cell.nat, cell.ntype, + cell.iat2it, cell.iat2ia, cells, period, cell.lat); + reused.cal_Ms(kv, cell.symm, cell.atoms, cell.latvec, cell.lmax, + cell.nat, cell.iat2it, cell.iat2ia, pv, 1); CellSymmetryRotation fresh; - fresh.find_irreducible_sector(cell.symm, cell.atoms, cell.st, cells, period, cell.lat); - fresh.cal_Ms(kv, cell, pv, 1); - EXPECT_EQ(reused.get_irreducible_sector(), fresh.get_irreducible_sector()); + fresh.find_irred_sector(cell.symm, cell.atoms, cell.nat, cell.ntype, + cell.iat2it, cell.iat2ia, cells, period, cell.lat); + fresh.cal_Ms(kv, cell.symm, cell.atoms, cell.latvec, cell.lmax, + cell.nat, cell.iat2it, cell.iat2ia, pv, 1); + EXPECT_EQ(reused.get_irred_sector(), fresh.get_irred_sector()); EXPECT_EQ(reused.rotations(), fresh.rotations()); EXPECT_NE(old_rotations, fresh.rotations()); for (int atom = 0; atom < cell.nat; ++atom) @@ -249,7 +259,7 @@ TEST(SymmetryDensityRestoration, ReturnLatticePreservesBoundaryRepresentatives) { ModuleSymmetry::Symmetry symmetry; symmetry.epsilon = 3.2e-5; - ModuleSymmetry::Irreducible_Sector sector; + ModuleSymmetry::Irred_Sector sector; const ModuleBase::Matrix3 reflection(-1, 0, 0, 0, 1, 0, 0, 0, 1); const ModuleBase::Vector3 translation(0.0, 0.0, 0.0); const ModuleBase::Vector3 source(0.99999, 0.25, 0.5); @@ -274,7 +284,7 @@ TEST(SymmetryDensityRestoration, ReturnLatticePreservesTranslationRepresentative { ModuleSymmetry::Symmetry symmetry; symmetry.epsilon = 3.2e-5; - ModuleSymmetry::Irreducible_Sector sector; + ModuleSymmetry::Irred_Sector sector; const ModuleBase::Matrix3 identity(1, 0, 0, 0, 1, 0, 0, 0, 1); const ModuleBase::Vector3 source(0.25, 0.5, 0.75); const ModuleBase::Vector3 translation(0.99999, 0.0, 0.0); diff --git a/source/source_lcao/module_ri/rpa_lri.hpp b/source/source_lcao/module_ri/rpa_lri.hpp index 5f19a66120b..e30641d8f44 100644 --- a/source/source_lcao/module_ri/rpa_lri.hpp +++ b/source/source_lcao/module_ri/rpa_lri.hpp @@ -16,6 +16,7 @@ #include #include #include "source_lcao/module_ri/module_exx_symmetry/symm_rotation.h" +#include "source_cell/unitcell.h" #include "rpa_lri.h" #include "exx_lri.h" @@ -167,13 +168,16 @@ void RPA_LRI::cal_postSCF_exx(const module_dm::DensityMatrix { const std::array period = RI_Util::get_Born_vonKarmen_period(kv); const auto& Rs = RI_Util::get_Born_von_Karmen_cells(period); - symrot.find_irreducible_sector(ucell.symm, ucell.atoms, ucell.st, Rs, period, ucell.lat, PARAM.globalv.global_out_dir); + symrot.find_irred_sector(ucell.symm, ucell.atoms, + ucell.nat, ucell.ntype, ucell.iat2it, ucell.iat2ia, + Rs, period, ucell.lat, PARAM.globalv.global_out_dir); // set Lmax of the rotation matrices to max(l_ao, l_abf), to support rotation under ABF // NOTE: Using Exx_Abfs::Construct_Orbs::get_Lmax() to compute Lmax from the actual ABFs // instead of relying on exx_cut_coulomb->abfs_Lmax() (not yet initialized) or // this->info.abfs_Lmax (defaults to 0). This ensures correct Lmax for symmetry rotation. symrot.set_abfs_Lmax(Exx_Abfs::Construct_Orbs::get_Lmax(abfs_for_lmax)); - symrot.cal_Ms(kv, ucell, parav, PARAM.inp.nspin); + symrot.cal_Ms(kv, ucell.symm, ucell.atoms, ucell.latvec, ucell.lmax, + ucell.nat, ucell.iat2it, ucell.iat2ia, parav, PARAM.inp.nspin); // output Ts (symrot_R.txt) and Ms (symrot_k.txt) ModuleSymmetry::print_symrot_info_R(symrot, ucell.symm, ucell.lmax, Rs); ModuleSymmetry::print_symrot_info_k(symrot, kv, ucell); diff --git a/source/source_lcao/test/test_init_dm_from_file.cpp b/source/source_lcao/test/test_init_dm_from_file.cpp index fb2ab91d67d..99647533ab5 100644 --- a/source/source_lcao/test/test_init_dm_from_file.cpp +++ b/source/source_lcao/test/test_init_dm_from_file.cpp @@ -39,8 +39,8 @@ class InitDMFileTest : public testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.atoms[0].taud.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); diff --git a/source/source_lcao/test/test_output_hcontainer_consistency.cpp b/source/source_lcao/test/test_output_hcontainer_consistency.cpp index 075495d4f1b..bed577dc630 100644 --- a/source/source_lcao/test/test_output_hcontainer_consistency.cpp +++ b/source/source_lcao/test/test_output_hcontainer_consistency.cpp @@ -41,8 +41,8 @@ class OutputHContainerTest : public testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.atoms[0].taud.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); diff --git a/source/source_pw/module_dfpt/test/dfpt_pw_data_test.cpp b/source/source_pw/module_dfpt/test/dfpt_pw_data_test.cpp index dcfc50dfbd7..3a89b83ba29 100644 --- a/source/source_pw/module_dfpt/test/dfpt_pw_data_test.cpp +++ b/source/source_pw/module_dfpt/test/dfpt_pw_data_test.cpp @@ -50,7 +50,9 @@ class DFPT_PW_DataTest : public DFPTStruTestFixture ofs_running.open("tmp_dfpt_qlist"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {2, 2, 2}, true); data.init(&qlist, 1, 2, 3, 0, 1, 1, nullptr); } diff --git a/source/source_pw/module_dfpt/test/dfpt_pw_run_test.cpp b/source/source_pw/module_dfpt/test/dfpt_pw_run_test.cpp index 881eca4d94f..c2b5646666f 100644 --- a/source/source_pw/module_dfpt/test/dfpt_pw_run_test.cpp +++ b/source/source_pw/module_dfpt/test/dfpt_pw_run_test.cpp @@ -46,7 +46,9 @@ class DFPT_PWRunTest : public DFPTStruTestFixture ofs_running.open("tmp_dfpt_run"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); ucell.symm = symm; } diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_q0_serial_test.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_q0_serial_test.cpp index d7f12199275..a02dd2a66c1 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_q0_serial_test.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_q0_serial_test.cpp @@ -556,10 +556,8 @@ TEST_F(DFPTQ0SerialTest, StarRotationCyclicGroup) ucell_.atoms[0].taud[i] = ModuleBase::Vector3(t0[i], t0[(i + 1) % 3], t0[(i + 2) % 3]); ucell_.atoms[0].tau[i] = ucell_.atoms[0].taud[i] * a_; } - delete[] ucell_.iat2it; - delete[] ucell_.iat2ia; - ucell_.iat2it = new int[3]; - ucell_.iat2ia = new int[3]; + ucell_.iat2it.resize(3); + ucell_.iat2ia.resize(3); for (int i = 0; i < 3; ++i) { ucell_.iat2it[i] = 0; diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.cpp index 13d843f72ed..dfbc4f70923 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.cpp @@ -24,8 +24,8 @@ void DFPTSerialBase::SetUpCell() ucell_.tpiba = ModuleBase::TWO_PI / lat0_; ucell_.tpiba2 = ucell_.tpiba * ucell_.tpiba; ucell_.omega = a_ * a_ * a_ * lat0_ * lat0_ * lat0_; - ucell_.iat2it = new int[1]; - ucell_.iat2ia = new int[1]; + ucell_.iat2it.resize(1); + ucell_.iat2ia.resize(1); ucell_.iat2it[0] = 0; ucell_.iat2ia[0] = 0; MakeCoulombAtom(); @@ -61,10 +61,6 @@ void DFPTSerialBase::TearDown() { delete[] ucell_.atoms; ucell_.atoms = nullptr; - delete[] ucell_.iat2it; - ucell_.iat2it = nullptr; - delete[] ucell_.iat2ia; - ucell_.iat2ia = nullptr; } void DFPTSerialBase::MakeCoulombAtom() @@ -143,10 +139,8 @@ void DFPTSerialBase::MakeTwoAtomCell() at.ncpp.nh = 0; at.mass = (it == 0) ? 12.0 : 4.0; } - delete[] ucell_.iat2it; - delete[] ucell_.iat2ia; - ucell_.iat2it = new int[2]; - ucell_.iat2ia = new int[2]; + ucell_.iat2it.resize(2); + ucell_.iat2ia.resize(2); ucell_.iat2it[0] = 0; ucell_.iat2ia[0] = 0; ucell_.iat2it[1] = 1; diff --git a/source/source_pw/module_pwdft/force_pw.cpp b/source/source_pw/module_pwdft/force_pw.cpp index d37e2379954..71dcec97601 100644 --- a/source/source_pw/module_pwdft/force_pw.cpp +++ b/source/source_pw/module_pwdft/force_pw.cpp @@ -367,7 +367,7 @@ void Forces::cal_force_loc(const UnitCell& ucell, resmem_var_op()(forcelc_d, this->nat * 3); resmem_var_op()(vloc_d, vloc.nr * vloc.nc); - syncmem_int_h2d_op()(iat2it_d, ucell.iat2it, this->nat); + syncmem_int_h2d_op()(iat2it_d, ucell.iat2it.data(), this->nat); syncmem_int_h2d_op()(ig2gg_d, rho_basis->ig2igg, rho_basis->npw); syncmem_var_h2d_op()(gcar_d, gcar_h.data(), rho_basis->npw * 3); syncmem_var_h2d_op()(tau_d, tau_h.data(), this->nat * 3); @@ -554,7 +554,7 @@ void Forces::cal_force_ew(const UnitCell& ucell, resmem_complex_op()(aux_d, rho_basis->npw); resmem_var_op()(forceion_d, this->nat * 3); - syncmem_int_h2d_op()(iat2it_d, ucell.iat2it, this->nat); + syncmem_int_h2d_op()(iat2it_d, ucell.iat2it.data(), this->nat); syncmem_var_h2d_op()(gcar_d, gcar_h.data(), rho_basis->npw * 3); syncmem_var_h2d_op()(tau_d, tau_h.data(), this->nat * 3); syncmem_var_h2d_op()(it_fact_d, it_fact_h.data(), ucell.ntype); diff --git a/source/source_pw/module_pwdft/test/dftu_base_test.cpp b/source/source_pw/module_pwdft/test/dftu_base_test.cpp index 2281b7f2186..cb2c3dc3a0d 100644 --- a/source/source_pw/module_pwdft/test/dftu_base_test.cpp +++ b/source/source_pw/module_pwdft/test/dftu_base_test.cpp @@ -32,8 +32,8 @@ class DFTUBaseTest : public testing::Test ucell.ntype = 1; ucell.nat = 1; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.atoms[0].taud.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); diff --git a/source/source_relax/unittests/for_test.h b/source/source_relax/unittests/for_test.h index 043c88ed470..698f57f3d4f 100644 --- a/source/source_relax/unittests/for_test.h +++ b/source/source_relax/unittests/for_test.h @@ -24,11 +24,6 @@ UnitCell::UnitCell() namax = 0; nwmax = 0; - iat2it = nullptr; - iat2ia = nullptr; - iwt2iat = nullptr; - iwt2iw = nullptr; - itia2iat.create(1, 1); latvec = ModuleBase::Matrix3(); diff --git a/source/source_relax/unittests/test_ions_move_bfgs2.cpp b/source/source_relax/unittests/test_ions_move_bfgs2.cpp index 76dad4a9f8b..8d176b780f7 100644 --- a/source/source_relax/unittests/test_ions_move_bfgs2.cpp +++ b/source/source_relax/unittests/test_ions_move_bfgs2.cpp @@ -226,8 +226,8 @@ TEST_F(BFGSTest, UpdateUsesAbsoluteDisplacementThreshold) ucell.atoms = new Atom[ucell.ntype]; ucell.atoms[0].na = 1; ucell.atoms[0].mbl = std::vector>(1, {1, 1, 1}); - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.iat2it[0] = 0; ucell.iat2ia[0] = 0; @@ -266,8 +266,8 @@ TEST_F(BFGSTest, GetPosAndPostaud) ucell.atoms[0].tau[1].x = 1.0; ucell.atoms[0].tau[1].y = 0.0; ucell.atoms[0].tau[1].z = 0.0; // allocate mapping arrays - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); int k = 0; for (int it = 0; it < ucell.ntype; ++it) { for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { @@ -301,8 +301,8 @@ TEST_F(BFGSTest, CalculateLargestGrad) ucell.atoms[0].mbl = std::vector>(2, {1, 1, 1}); // mapping arrays - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); int k = 0; for (int it = 0; it < ucell.ntype; ++it) { for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { @@ -343,8 +343,8 @@ TEST_F(BFGSTest, RelaxStepBasic) ucell.atoms[0].tau = std::vector>(2); ucell.atoms[0].taud = std::vector>(2); ucell.atoms[0].mbl = std::vector>(2, {1, 1, 1}); - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); int k = 0; for (int it = 0; it < ucell.ntype; ++it) { for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { diff --git a/source/source_relax/unittests/test_ions_move_lbfgs.cpp b/source/source_relax/unittests/test_ions_move_lbfgs.cpp index e7f882719b7..93275d92b9e 100644 --- a/source/source_relax/unittests/test_ions_move_lbfgs.cpp +++ b/source/source_relax/unittests/test_ions_move_lbfgs.cpp @@ -38,10 +38,6 @@ UnitCell::UnitCell() ntype = 1; nat = 1; - iat2it = nullptr; - iat2ia = nullptr; - iwt2iat = nullptr; - iwt2iw = nullptr; itia2iat.create(1, 1); atoms = new Atom[ntype]; diff --git a/source/source_relax/unittests/test_relax_nsync.cpp b/source/source_relax/unittests/test_relax_nsync.cpp index fffd6fdd16f..0efd1289544 100644 --- a/source/source_relax/unittests/test_relax_nsync.cpp +++ b/source/source_relax/unittests/test_relax_nsync.cpp @@ -49,10 +49,6 @@ UnitCell::UnitCell() ntype = 1; nat = 1; - iat2it = nullptr; - iat2ia = nullptr; - iwt2iat = nullptr; - iwt2iw = nullptr; itia2iat.create(1, 1); atoms = new Atom[ntype]; @@ -181,13 +177,10 @@ class IonCellOptimizerTest : public ::testing::Test ucell.lat_axis_free[0] = 1; ucell.lat_axis_free[1] = 1; ucell.lat_axis_free[2] = 1; - // The mock UnitCell leaves iat2it/iat2ia null; the L-BFGS update loop - // dereferences them from the second iteration on, so give the single - // atom a valid type/index. Ownership passes to UnitCell's internal - // Statistics member, whose destructor releases them -- do NOT delete - // them here. - ucell.iat2it = new int[natom]; - ucell.iat2ia = new int[natom]; + // The L-BFGS update loop dereferences iat2it/iat2ia from the second + // iteration on, so give the single atom a valid type/index. + ucell.iat2it.resize(natom); + ucell.iat2ia.resize(natom); ucell.iat2it[0] = 0; ucell.iat2ia[0] = 0; } diff --git a/source/source_relax/unittests/test_relax_sync.cpp b/source/source_relax/unittests/test_relax_sync.cpp index ce176d42345..0405bcf59bc 100644 --- a/source/source_relax/unittests/test_relax_sync.cpp +++ b/source/source_relax/unittests/test_relax_sync.cpp @@ -89,8 +89,8 @@ class Test_SETGRAD : public testing::Test ucell.omega = 1.0; ucell.lat0 = 1.0; - ucell.iat2it = new int[nat]; - ucell.iat2ia = new int[nat]; + ucell.iat2it.resize(nat); + ucell.iat2ia.resize(nat); ucell.atoms[0].mbl.resize(nat); ucell.atoms[0].taud.resize(nat); ucell.atoms[0].tau.resize(nat); @@ -296,8 +296,8 @@ class Test_RELAX : public testing::Test ucell.omega = 452.590903143121; ucell.lat0 = 1.8897259886; - ucell.iat2it = new int[nat]; - ucell.iat2ia = new int[nat]; + ucell.iat2it.resize(nat); + ucell.iat2ia.resize(nat); ucell.iat2it[0] = 0; ucell.iat2it[1] = 1; ucell.iat2it[2] = 2; @@ -417,8 +417,8 @@ TEST(RelaxSyncSummary, ConvergedPrintsSummary) ucell.atoms[0].label = "Si"; ucell.omega = 1.0; ucell.lat0 = 1.0; - ucell.iat2it = new int[nat]; - ucell.iat2ia = new int[nat]; + ucell.iat2it.resize(nat); + ucell.iat2ia.resize(nat); ucell.iat2it[0] = 0; ucell.iat2ia[0] = 0; ucell.atoms[0].mbl.resize(nat); @@ -477,8 +477,8 @@ TEST(RelaxSyncSummary, TwoStepConvergence) ucell.atoms[0].label = "Si"; ucell.omega = 1.0; ucell.lat0 = 1.0; - ucell.iat2it = new int[nat]; - ucell.iat2ia = new int[nat]; + ucell.iat2it.resize(nat); + ucell.iat2ia.resize(nat); ucell.iat2it[0] = 0; ucell.iat2ia[0] = 0; ucell.atoms[0].mbl.resize(nat); @@ -547,8 +547,8 @@ TEST(RelaxSyncSummary, CellRelaxConvergedPrintsStressHistory) ucell.atoms[0].label = "Si"; ucell.omega = 1.0; ucell.lat0 = 1.0; - ucell.iat2it = new int[nat]; - ucell.iat2ia = new int[nat]; + ucell.iat2it.resize(nat); + ucell.iat2ia.resize(nat); ucell.iat2it[0] = 0; ucell.iat2ia[0] = 0; ucell.atoms[0].mbl.resize(nat); @@ -622,8 +622,8 @@ void make_two_type_cell(UnitCell& ucell, const int mbl_flat[9]) ucell.atoms[t].vel.resize(na); } // iat: 0 -> (type0, ia0); 1 -> (type1, ia0); 2 -> (type1, ia1) - ucell.iat2it = new int[3]; - ucell.iat2ia = new int[3]; + ucell.iat2it.resize(3); + ucell.iat2ia.resize(3); ucell.iat2it[0] = 0; ucell.iat2it[1] = 1; ucell.iat2it[2] = 1; ucell.iat2ia[0] = 0; ucell.iat2ia[1] = 0; ucell.iat2ia[2] = 1; @@ -648,11 +648,6 @@ void make_two_type_cell(UnitCell& ucell, const int mbl_flat[9]) void free_two_type_cell(UnitCell& ucell) { delete[] ucell.atoms; - // iat2it / iat2ia are owned by the mock; release explicitly. - delete[] ucell.iat2it; - delete[] ucell.iat2ia; - ucell.iat2it = nullptr; - ucell.iat2ia = nullptr; } } // namespace @@ -730,8 +725,8 @@ TEST(RelaxSyncMaxStep, GeometryMovedAndNotConverged) ucell.atoms[0].label = "Si"; ucell.omega = 1.0; ucell.lat0 = 1.0; - ucell.iat2it = new int[nat]; - ucell.iat2ia = new int[nat]; + ucell.iat2it.resize(nat); + ucell.iat2ia.resize(nat); ucell.iat2it[0] = 0; ucell.iat2ia[0] = 0; ucell.atoms[0].mbl.resize(nat);