From 42f33d27db87ce71775d171f73784644b64efe47 Mon Sep 17 00:00:00 2001 From: DanielZhao0432 Date: Sat, 10 Oct 2026 12:02:15 +0800 Subject: [PATCH 1/3] refactor(cell): split unitcell.cpp into index/stats/setup files --- source/source_cell/CMakeLists.txt | 3 + source/source_cell/test/CMakeLists.txt | 6 + source/source_cell/test_pw/CMakeLists.txt | 2 +- source/source_cell/unitcell.cpp | 335 +----------------- source/source_cell/unitcell_index.cpp | 48 +++ source/source_cell/unitcell_setup.cpp | 259 ++++++++++++++ source/source_cell/unitcell_stats.cpp | 39 ++ .../module_deepks/test/CMakeLists.txt | 3 + source/source_md/test/CMakeLists.txt | 3 + .../module_pwdft/test/CMakeLists.txt | 6 + 10 files changed, 371 insertions(+), 333 deletions(-) create mode 100644 source/source_cell/unitcell_index.cpp create mode 100644 source/source_cell/unitcell_setup.cpp create mode 100644 source/source_cell/unitcell_stats.cpp diff --git a/source/source_cell/CMakeLists.txt b/source/source_cell/CMakeLists.txt index b3e574e7144..e70e2aa7654 100644 --- a/source/source_cell/CMakeLists.txt +++ b/source/source_cell/CMakeLists.txt @@ -24,6 +24,9 @@ add_library( mdcell_reader.cpp mdcell.cpp unitcell.cpp + unitcell_index.cpp + unitcell_stats.cpp + unitcell_setup.cpp read_atoms.cpp read_atoms_helper.cpp read_orb.cpp diff --git a/source/source_cell/test/CMakeLists.txt b/source/source_cell/test/CMakeLists.txt index e0db610c5bf..dc436a829d4 100644 --- a/source/source_cell/test/CMakeLists.txt +++ b/source/source_cell/test/CMakeLists.txt @@ -23,6 +23,9 @@ install(FILES bcast_sep_cell_test.sh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) list(APPEND cell_simple_srcs ../unitcell.cpp + ../unitcell_index.cpp + ../unitcell_stats.cpp + ../unitcell_setup.cpp ../update_cell.cpp ../bcast_cell.cpp ../print_cell.cpp @@ -138,6 +141,9 @@ AddTest( ../print_cell.cpp ../atom_spec.cpp ../unitcell.cpp + ../unitcell_index.cpp + ../unitcell_stats.cpp + ../unitcell_setup.cpp ../update_cell.cpp ../bcast_cell.cpp ../atom_pseudo.cpp diff --git a/source/source_cell/test_pw/CMakeLists.txt b/source/source_cell/test_pw/CMakeLists.txt index 2b1bde23d6a..f5051876918 100644 --- a/source/source_cell/test_pw/CMakeLists.txt +++ b/source/source_cell/test_pw/CMakeLists.txt @@ -12,7 +12,7 @@ install(FILES unitcell_test_pw_para.sh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) AddTest( TARGET MODULE_CELL_unitcell_test_pw LIBS parameter base device - SOURCES unitcell_test_pw.cpp ../unitcell.cpp ../read_atoms.cpp ../read_atoms_helper.cpp ../cell_tools.cpp ../atom_spec.cpp ../update_cell.cpp ../bcast_cell.cpp + SOURCES unitcell_test_pw.cpp ../unitcell.cpp ../unitcell_index.cpp ../unitcell_stats.cpp ../unitcell_setup.cpp ../read_atoms.cpp ../read_atoms_helper.cpp ../cell_tools.cpp ../atom_spec.cpp ../update_cell.cpp ../bcast_cell.cpp ../atom_pseudo.cpp ../pseudo.cpp ../read_pp.cpp ../read_pp_complete.cpp ../read_pp_upf201.cpp ../read_pp_upf100.cpp ../read_stru.cpp ../read_atom_species.cpp ../read_pp_vwr.cpp ../read_pp_blps.cpp diff --git a/source/source_cell/unitcell.cpp b/source/source_cell/unitcell.cpp index cc4db426946..f7cc675ecd7 100644 --- a/source/source_cell/unitcell.cpp +++ b/source/source_cell/unitcell.cpp @@ -1,30 +1,10 @@ /** * @file unitcell.cpp - * @brief Implementation of UnitCell class. + * @brief Implementation of UnitCell class: constructor and destructor. + * Indexing helpers live in unitcell_index.cpp, per-species counting + * helpers in unitcell_stats.cpp, and setup routines in unitcell_setup.cpp. */ -#include -#include // Peize Lin fix bug about strcmp 2016-08-02 - -#include "source_base/constants.h" -#include "source_base/global_function.h" -#include "source_base/global_variable.h" #include "unitcell.h" -#include "bcast_cell.h" -#include "source_base/tool_quit.h" -#include "source_base/output.h" - -#include "source_cell/read_stru.h" -#include "source_base/atom_in.h" -#include "source_base/element_elec_config.h" -#include "source_base/global_file.h" -#include "source_base/parallel_common.h" -#include "source_cell/sep_cell.h" - -#ifdef __MPI -#include "mpi.h" -#endif - -#include "update_cell.h" UnitCell::UnitCell() { @@ -38,312 +18,3 @@ UnitCell::~UnitCell() delete[] atoms; } } - - -void UnitCell::set_iat2itia() { - assert(nat > 0); - delete[] iat2it; - delete[] iat2ia; - this->iat2it = new int[nat]; - this->iat2ia = new int[nat]; - int iat = 0; - for (int it = 0; it < ntype; it++) { - for (int ia = 0; ia < atoms[it].na; ia++) { - this->iat2it[iat] = it; - this->iat2ia[iat] = ia; - ++iat; - } - } - return; -} - -std::map UnitCell::get_atom_Counts() const { - std::map atomCounts; - for (int it = 0; it < this->ntype; it++) { - atomCounts.insert(std::pair(it, this->atoms[it].na)); - } - return atomCounts; -} - -std::map UnitCell::get_orbital_Counts() const { - std::map orbitalCounts; - for (int it = 0; it < this->ntype; it++) { - orbitalCounts.insert(std::pair(it, this->atoms[it].nw)); - } - return orbitalCounts; -} - -std::map> UnitCell::get_lnchi_Counts() const { - std::map> lnchiCounts; - for (int it = 0; it < this->ntype; it++) { - for (int L = 0; L < this->atoms[it].nwl + 1; L++) { - // Check if the key 'it' exists in the outer map - if (lnchiCounts.find(it) == lnchiCounts.end()) { - // If it doesn't exist, initialize an empty inner map - lnchiCounts[it] = std::map(); - } - int l_nchi = this->atoms[it].l_nchi[L]; - // Insert the key-value pair into the inner map - lnchiCounts[it].insert(std::pair(L, l_nchi)); - } - } - return lnchiCounts; -} - -//============================================================== -// Calculate various lattice related quantities for given latvec -//============================================================== -void UnitCell::setup_cell(const std::string& fn, std::ofstream& log, const double symmetry_prec, - const int dfthalf_type, const std::string& pseudo_dir, const int nspin, - const std::string& basis_type, const std::string& orbital_dir, const std::string& init_wfc, - const double onsite_radius, const bool deepks_setorb, const bool rpa, - const bool fixed_atoms, const bool noncolin, const std::string& calculation, - const std::string& esolver_type, const int symmetry) -{ - ModuleBase::TITLE("UnitCell", "setup_cell"); - - assert(ntype > 0); - - // (1) init *Atom class array. - this->atoms = new Atom[this->ntype]; // atom species. - this->set_atom_flag = true; - - this->symm.epsilon = symmetry_prec; - this->symm.epsilon_input = symmetry_prec; - - bool ok = true; - bool ok2 = true; - - bool ok3 = true; // for sep potential in DFT-1/2 - - // (3) read in atom information - this->pseudo_fn.resize(ntype); - this->pseudo_type.resize(ntype); - this->orbital_fn.resize(ntype); - - if (GlobalV::MY_RANK == 0) - { - // open "atom_unitcell" file. - std::ifstream ifa(fn.c_str(), std::ios::in); - if (!ifa) - { - GlobalV::ofs_warning << fn; - ok = false; - } - - if (ok) - { - log << "\n\n"; - log << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; - log << " | |" << std::endl; - log << " | #Setup Unitcell# |" << std::endl; - log << " | From the input file and the structure file we know the number of |" << std::endl; - log << " | different elments in this unitcell, then we list the detail |" << std::endl; - log << " | information for each element, especially the zeta and polar atomic |" << std::endl; - log << " | orbital number for each element. The total atom number is counted. |" << std::endl; - log << " | We calculate the nearest atom distance for each atom and show the |" << std::endl; - log << " | Cartesian and Direct coordinates for each atom. We list the file |" << std::endl; - log << " | address for atomic orbitals. The volume and the lattice vectors |" << std::endl; - log << " | in real and reciprocal space is also shown. |" << std::endl; - log << " | |" << std::endl; - log << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; - log << "\n"; - - log << " READING UNITCELL INFORMATION" << std::endl; - //======================== - // call read_atom_species - //======================== - const bool read_atom_species = unitcell::read_atom_species(ifa, log, *this, - basis_type, orbital_dir, init_wfc, onsite_radius, deepks_setorb, rpa); - //======================== - // call read_lattice_constant - //======================== - const bool read_lattice_constant = unitcell::read_lattice_constant(ifa, log ,this->lat); - //========================== - // readl sep potential, currently using the pseudopotential folder (pseudo_dir in INPUT) - //========================== - if (dfthalf_type > 0) { - sep_cell.init(this->ntype); - std::vector atom_labels(this->ntype); - for (int i = 0; i < this->ntype; ++i) - { - atom_labels[i] = this->atoms[i].label; - } - ok3 = sep_cell.read_sep_potentials(ifa, pseudo_dir, GlobalV::ofs_warning, atom_labels); - } - //========================== - // call read_atom_positions - //========================== - ok2 = unitcell::read_atom_positions(*this, ifa, log, GlobalV::ofs_warning, nspin, - basis_type, orbital_dir, init_wfc, onsite_radius, fixed_atoms, noncolin, - calculation, esolver_type, symmetry); - } - } -#ifdef __MPI - Parallel_Common::bcast_bool(ok); - Parallel_Common::bcast_bool(ok2); - Parallel_Common::bcast_bool(ok3); -#endif - if (!ok) { - ModuleBase::WARNING_QUIT( - "UnitCell::setup_cell", - "Can not find the file containing atom positions.!"); - } - if (!ok2) { - ModuleBase::WARNING_QUIT("UnitCell::setup_cell", - "Something wrong during read_atom_positions."); - } - if (!ok3) { - ModuleBase::WARNING_QUIT("UnitCell::setup_cell", "Something wrong during read_sep_potentials"); - } - -#ifdef __MPI - unitcell::bcast_unitcell(*this, nspin); - sep_cell.bcast_sep_cell(); -#endif - - //======================================================== - // Calculate unit cell volume - // the reason to calculate volume here is - // Firstly, latvec must be read in. - //======================================================== - assert(lat0 > 0.0); - this->omega = latvec.Det() * this->lat0 * this->lat0 * this->lat0; - - - if (this->omega < 0) - { - std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - std::cout << " Warning: The lattice vector is left-handed; a right-handed vector is prefered." << std::endl; - std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - GlobalV::ofs_warning << - "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - GlobalV::ofs_warning << - " Warning: The lattice vector is left-handed; a right-handed vector is prefered." << std::endl; - GlobalV::ofs_warning << - "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - this->omega = std::abs(this->omega); - } - else if (this->omega == 0) - { - ModuleBase::WARNING_QUIT("setup_cell", "The volume is zero."); - } - else - { - ModuleBase::GlobalFunc::OUT(log, "Cell volume (Bohr^3)", this->omega); - ModuleBase::GlobalFunc::OUT(log, "Cell volume (A^3)", this->omega * pow(ModuleBase::BOHR_TO_A, 3)); - } - - //========================================================== - // Calculate recip. lattice vectors and dot products - // latvec have the unit of lat0, but G has the unit 2Pi/lat0 - //========================================================== - this->GT = latvec.Inverse(); - this->G = GT.Transpose(); - this->GGT = G * GT; - this->invGGT = GGT.Inverse(); - - log << std::endl; - output::printM3(log, - "Lattice vectors: (Cartesian coordinate: in unit of a_0)", - latvec); - output::printM3( - log, - "Reciprocal vectors: (Cartesian coordinate: in unit of 2 pi/a_0)", - G); - - //=================================== - // set index for iat2it, iat2ia - //=================================== - this->set_iat2itia(); - - sep_cell.set_omega(this->omega, this->tpiba2); - - return; -} - - -void UnitCell::set_iat2iwt(const int& npol_in) -{ -#ifdef __DEBUG - assert(npol_in == 1 || npol_in == 2); - assert(this->nat > 0); - assert(this->ntype > 0); -#endif - this->iat2iwt.resize(this->nat); - this->npol = npol_in; - int iat = 0; - int iwt = 0; - - for (int it = 0; it < this->ntype; it++) - { - for (int ia = 0; ia < atoms[it].na; ia++) - { - this->iat2iwt[iat] = iwt; - iwt += atoms[it].nw * this->npol; - ++iat; - } - } - return; -} - - - -void UnitCell::setup_from_input(const std::string& latname_in, - const int& ntype_in, - const int& lmaxmax_in, - const bool& init_vel_in, - const std::string& fixed_axes_in) { - this->latName = latname_in; - this->ntype = ntype_in; - this->magnet.start_mag.resize(ntype_in, 0.0); - this->lmaxmax = lmaxmax_in; - this->init_vel = init_vel_in; - // pengfei Li add 2018-11-11 - if (fixed_axes_in == "None") { - this->lat_axis_free[0] = 1; - this->lat_axis_free[1] = 1; - this->lat_axis_free[2] = 1; - } else if (fixed_axes_in == "volume") { - this->lat_axis_free[0] = 1; - this->lat_axis_free[1] = 1; - this->lat_axis_free[2] = 1; - } else if (fixed_axes_in == "shape") { - this->lat_axis_free[0] = 1; - this->lat_axis_free[1] = 1; - this->lat_axis_free[2] = 1; - } else if (fixed_axes_in == "a") { - this->lat_axis_free[0] = 0; - this->lat_axis_free[1] = 1; - this->lat_axis_free[2] = 1; - } else if (fixed_axes_in == "b") { - this->lat_axis_free[0] = 1; - this->lat_axis_free[1] = 0; - this->lat_axis_free[2] = 1; - } else if (fixed_axes_in == "c") { - this->lat_axis_free[0] = 1; - this->lat_axis_free[1] = 1; - this->lat_axis_free[2] = 0; - } else if (fixed_axes_in == "ab") { - this->lat_axis_free[0] = 0; - this->lat_axis_free[1] = 0; - this->lat_axis_free[2] = 1; - } else if (fixed_axes_in == "ac") { - this->lat_axis_free[0] = 0; - this->lat_axis_free[1] = 1; - this->lat_axis_free[2] = 0; - } else if (fixed_axes_in == "bc") { - this->lat_axis_free[0] = 1; - this->lat_axis_free[1] = 0; - this->lat_axis_free[2] = 0; - } else if (fixed_axes_in == "abc") { - this->lat_axis_free[0] = 0; - this->lat_axis_free[1] = 0; - this->lat_axis_free[2] = 0; - } else { - ModuleBase::WARNING_QUIT( - "Input", - "fixed_axes should be none, volume, shape, a, b, c, ab, ac, bc or abc!"); - } - return; -} diff --git a/source/source_cell/unitcell_index.cpp b/source/source_cell/unitcell_index.cpp new file mode 100644 index 00000000000..edc05240bbf --- /dev/null +++ b/source/source_cell/unitcell_index.cpp @@ -0,0 +1,48 @@ +/** + * @file unitcell_index.cpp + * @brief Atom/orbital indexing helpers of UnitCell: set_iat2itia, set_iat2iwt. + */ +#include "unitcell.h" + +#include + +void UnitCell::set_iat2itia() { + assert(nat > 0); + delete[] iat2it; + delete[] iat2ia; + this->iat2it = new int[nat]; + this->iat2ia = new int[nat]; + int iat = 0; + for (int it = 0; it < ntype; it++) { + for (int ia = 0; ia < atoms[it].na; ia++) { + this->iat2it[iat] = it; + this->iat2ia[iat] = ia; + ++iat; + } + } + return; +} + +void UnitCell::set_iat2iwt(const int& npol_in) +{ +#ifdef __DEBUG + assert(npol_in == 1 || npol_in == 2); + assert(this->nat > 0); + assert(this->ntype > 0); +#endif + this->iat2iwt.resize(this->nat); + this->npol = npol_in; + int iat = 0; + int iwt = 0; + + for (int it = 0; it < this->ntype; it++) + { + for (int ia = 0; ia < atoms[it].na; ia++) + { + this->iat2iwt[iat] = iwt; + iwt += atoms[it].nw * this->npol; + ++iat; + } + } + return; +} diff --git a/source/source_cell/unitcell_setup.cpp b/source/source_cell/unitcell_setup.cpp new file mode 100644 index 00000000000..11ae80df0c5 --- /dev/null +++ b/source/source_cell/unitcell_setup.cpp @@ -0,0 +1,259 @@ +/** + * @file unitcell_setup.cpp + * @brief Setup routines of UnitCell: setup_cell, setup_from_input. + */ +#include +#include // Peize Lin fix bug about strcmp 2016-08-02 + +#include "source_base/constants.h" +#include "source_base/global_function.h" +#include "source_base/global_variable.h" +#include "unitcell.h" +#include "bcast_cell.h" +#include "source_base/tool_quit.h" +#include "source_base/output.h" + +#include "source_cell/read_stru.h" +#include "source_base/atom_in.h" +#include "source_base/element_elec_config.h" +#include "source_base/global_file.h" +#include "source_base/parallel_common.h" +#include "source_cell/sep_cell.h" + +#ifdef __MPI +#include "mpi.h" +#endif + +#include "update_cell.h" + +//============================================================== +// Calculate various lattice related quantities for given latvec +//============================================================== +void UnitCell::setup_cell(const std::string& fn, std::ofstream& log, const double symmetry_prec, + const int dfthalf_type, const std::string& pseudo_dir, const int nspin, + const std::string& basis_type, const std::string& orbital_dir, const std::string& init_wfc, + const double onsite_radius, const bool deepks_setorb, const bool rpa, + const bool fixed_atoms, const bool noncolin, const std::string& calculation, + const std::string& esolver_type, const int symmetry) +{ + ModuleBase::TITLE("UnitCell", "setup_cell"); + + assert(ntype > 0); + + // (1) init *Atom class array. + this->atoms = new Atom[this->ntype]; // atom species. + this->set_atom_flag = true; + + this->symm.epsilon = symmetry_prec; + this->symm.epsilon_input = symmetry_prec; + + bool ok = true; + bool ok2 = true; + + bool ok3 = true; // for sep potential in DFT-1/2 + + // (3) read in atom information + this->pseudo_fn.resize(ntype); + this->pseudo_type.resize(ntype); + this->orbital_fn.resize(ntype); + + if (GlobalV::MY_RANK == 0) + { + // open "atom_unitcell" file. + std::ifstream ifa(fn.c_str(), std::ios::in); + if (!ifa) + { + GlobalV::ofs_warning << fn; + ok = false; + } + + if (ok) + { + log << "\n\n"; + log << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; + log << " | |" << std::endl; + log << " | #Setup Unitcell# |" << std::endl; + log << " | From the input file and the structure file we know the number of |" << std::endl; + log << " | different elments in this unitcell, then we list the detail |" << std::endl; + log << " | information for each element, especially the zeta and polar atomic |" << std::endl; + log << " | orbital number for each element. The total atom number is counted. |" << std::endl; + log << " | We calculate the nearest atom distance for each atom and show the |" << std::endl; + log << " | Cartesian and Direct coordinates for each atom. We list the file |" << std::endl; + log << " | address for atomic orbitals. The volume and the lattice vectors |" << std::endl; + log << " | in real and reciprocal space is also shown. |" << std::endl; + log << " | |" << std::endl; + log << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; + log << "\n"; + + log << " READING UNITCELL INFORMATION" << std::endl; + //======================== + // call read_atom_species + //======================== + const bool read_atom_species = unitcell::read_atom_species(ifa, log, *this, + basis_type, orbital_dir, init_wfc, onsite_radius, deepks_setorb, rpa); + //======================== + // call read_lattice_constant + //======================== + const bool read_lattice_constant = unitcell::read_lattice_constant(ifa, log ,this->lat); + //========================== + // readl sep potential, currently using the pseudopotential folder (pseudo_dir in INPUT) + //========================== + if (dfthalf_type > 0) { + sep_cell.init(this->ntype); + std::vector atom_labels(this->ntype); + for (int i = 0; i < this->ntype; ++i) + { + atom_labels[i] = this->atoms[i].label; + } + ok3 = sep_cell.read_sep_potentials(ifa, pseudo_dir, GlobalV::ofs_warning, atom_labels); + } + //========================== + // call read_atom_positions + //========================== + ok2 = unitcell::read_atom_positions(*this, ifa, log, GlobalV::ofs_warning, nspin, + basis_type, orbital_dir, init_wfc, onsite_radius, fixed_atoms, noncolin, + calculation, esolver_type, symmetry); + } + } +#ifdef __MPI + Parallel_Common::bcast_bool(ok); + Parallel_Common::bcast_bool(ok2); + Parallel_Common::bcast_bool(ok3); +#endif + if (!ok) { + ModuleBase::WARNING_QUIT( + "UnitCell::setup_cell", + "Can not find the file containing atom positions.!"); + } + if (!ok2) { + ModuleBase::WARNING_QUIT("UnitCell::setup_cell", + "Something wrong during read_atom_positions."); + } + if (!ok3) { + ModuleBase::WARNING_QUIT("UnitCell::setup_cell", "Something wrong during read_sep_potentials"); + } + +#ifdef __MPI + unitcell::bcast_unitcell(*this, nspin); + sep_cell.bcast_sep_cell(); +#endif + + //======================================================== + // Calculate unit cell volume + // the reason to calculate volume here is + // Firstly, latvec must be read in. + //======================================================== + assert(lat0 > 0.0); + this->omega = latvec.Det() * this->lat0 * this->lat0 * this->lat0; + + + if (this->omega < 0) + { + std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + std::cout << " Warning: The lattice vector is left-handed; a right-handed vector is prefered." << std::endl; + std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + GlobalV::ofs_warning << + "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + GlobalV::ofs_warning << + " Warning: The lattice vector is left-handed; a right-handed vector is prefered." << std::endl; + GlobalV::ofs_warning << + "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + this->omega = std::abs(this->omega); + } + else if (this->omega == 0) + { + ModuleBase::WARNING_QUIT("setup_cell", "The volume is zero."); + } + else + { + ModuleBase::GlobalFunc::OUT(log, "Cell volume (Bohr^3)", this->omega); + ModuleBase::GlobalFunc::OUT(log, "Cell volume (A^3)", this->omega * pow(ModuleBase::BOHR_TO_A, 3)); + } + + //========================================================== + // Calculate recip. lattice vectors and dot products + // latvec have the unit of lat0, but G has the unit 2Pi/lat0 + //========================================================== + this->GT = latvec.Inverse(); + this->G = GT.Transpose(); + this->GGT = G * GT; + this->invGGT = GGT.Inverse(); + + log << std::endl; + output::printM3(log, + "Lattice vectors: (Cartesian coordinate: in unit of a_0)", + latvec); + output::printM3( + log, + "Reciprocal vectors: (Cartesian coordinate: in unit of 2 pi/a_0)", + G); + + //=================================== + // set index for iat2it, iat2ia + //=================================== + this->set_iat2itia(); + + sep_cell.set_omega(this->omega, this->tpiba2); + + return; +} + + +void UnitCell::setup_from_input(const std::string& latname_in, + const int& ntype_in, + const int& lmaxmax_in, + const bool& init_vel_in, + const std::string& fixed_axes_in) { + this->latName = latname_in; + this->ntype = ntype_in; + this->magnet.start_mag.resize(ntype_in, 0.0); + this->lmaxmax = lmaxmax_in; + this->init_vel = init_vel_in; + // pengfei Li add 2018-11-11 + if (fixed_axes_in == "None") { + this->lat_axis_free[0] = 1; + this->lat_axis_free[1] = 1; + this->lat_axis_free[2] = 1; + } else if (fixed_axes_in == "volume") { + this->lat_axis_free[0] = 1; + this->lat_axis_free[1] = 1; + this->lat_axis_free[2] = 1; + } else if (fixed_axes_in == "shape") { + this->lat_axis_free[0] = 1; + this->lat_axis_free[1] = 1; + this->lat_axis_free[2] = 1; + } else if (fixed_axes_in == "a") { + this->lat_axis_free[0] = 0; + this->lat_axis_free[1] = 1; + this->lat_axis_free[2] = 1; + } else if (fixed_axes_in == "b") { + this->lat_axis_free[0] = 1; + this->lat_axis_free[1] = 0; + this->lat_axis_free[2] = 1; + } else if (fixed_axes_in == "c") { + this->lat_axis_free[0] = 1; + this->lat_axis_free[1] = 1; + this->lat_axis_free[2] = 0; + } else if (fixed_axes_in == "ab") { + this->lat_axis_free[0] = 0; + this->lat_axis_free[1] = 0; + this->lat_axis_free[2] = 1; + } else if (fixed_axes_in == "ac") { + this->lat_axis_free[0] = 0; + this->lat_axis_free[1] = 1; + this->lat_axis_free[2] = 0; + } else if (fixed_axes_in == "bc") { + this->lat_axis_free[0] = 1; + this->lat_axis_free[1] = 0; + this->lat_axis_free[2] = 0; + } else if (fixed_axes_in == "abc") { + this->lat_axis_free[0] = 0; + this->lat_axis_free[1] = 0; + this->lat_axis_free[2] = 0; + } else { + ModuleBase::WARNING_QUIT( + "Input", + "fixed_axes should be none, volume, shape, a, b, c, ab, ac, bc or abc!"); + } + return; +} diff --git a/source/source_cell/unitcell_stats.cpp b/source/source_cell/unitcell_stats.cpp new file mode 100644 index 00000000000..a550740ba2c --- /dev/null +++ b/source/source_cell/unitcell_stats.cpp @@ -0,0 +1,39 @@ +/** + * @file unitcell_stats.cpp + * @brief Per-species counting helpers of UnitCell: + * get_atom_Counts, get_orbital_Counts, get_lnchi_Counts. + */ +#include "unitcell.h" + +std::map UnitCell::get_atom_Counts() const { + std::map atomCounts; + for (int it = 0; it < this->ntype; it++) { + atomCounts.insert(std::pair(it, this->atoms[it].na)); + } + return atomCounts; +} + +std::map UnitCell::get_orbital_Counts() const { + std::map orbitalCounts; + for (int it = 0; it < this->ntype; it++) { + orbitalCounts.insert(std::pair(it, this->atoms[it].nw)); + } + return orbitalCounts; +} + +std::map> UnitCell::get_lnchi_Counts() const { + std::map> lnchiCounts; + for (int it = 0; it < this->ntype; it++) { + for (int L = 0; L < this->atoms[it].nwl + 1; L++) { + // Check if the key 'it' exists in the outer map + if (lnchiCounts.find(it) == lnchiCounts.end()) { + // If it doesn't exist, initialize an empty inner map + lnchiCounts[it] = std::map(); + } + int l_nchi = this->atoms[it].l_nchi[L]; + // Insert the key-value pair into the inner map + lnchiCounts[it].insert(std::pair(L, l_nchi)); + } + } + return lnchiCounts; +} diff --git a/source/source_lcao/module_deepks/test/CMakeLists.txt b/source/source_lcao/module_deepks/test/CMakeLists.txt index 56748e172c7..b7d0510e0d9 100644 --- a/source/source_lcao/module_deepks/test/CMakeLists.txt +++ b/source/source_lcao/module_deepks/test/CMakeLists.txt @@ -15,6 +15,9 @@ set(DEEPKS_UNIT_COMMON_SOURCES deepks_test_prep.cpp deepks_test.cpp ../../../source_cell/unitcell.cpp + ../../../source_cell/unitcell_index.cpp + ../../../source_cell/unitcell_stats.cpp + ../../../source_cell/unitcell_setup.cpp ../../../source_cell/update_cell.cpp ../../../source_cell/bcast_cell.cpp ../../../source_cell/atom_spec.cpp diff --git a/source/source_md/test/CMakeLists.txt b/source/source_md/test/CMakeLists.txt index 406cd5be5d8..4686c1f9b9a 100644 --- a/source/source_md/test/CMakeLists.txt +++ b/source/source_md/test/CMakeLists.txt @@ -5,6 +5,9 @@ list(APPEND depend_files ../md_func.cpp ../../source_cell/basecell.cpp ../../source_cell/unitcell.cpp + ../../source_cell/unitcell_index.cpp + ../../source_cell/unitcell_stats.cpp + ../../source_cell/unitcell_setup.cpp ../../source_cell/update_cell.cpp ../../source_cell/bcast_cell.cpp ../../source_cell/print_cell.cpp diff --git a/source/source_pw/module_pwdft/test/CMakeLists.txt b/source/source_pw/module_pwdft/test/CMakeLists.txt index 326f27a64e6..93b877b7682 100644 --- a/source/source_pw/module_pwdft/test/CMakeLists.txt +++ b/source/source_pw/module_pwdft/test/CMakeLists.txt @@ -32,6 +32,9 @@ AddTest( ../../../source_estate/occ_matrix.cpp ../../../source_estate/occ_mixer.cpp ../../../source_cell/unitcell.cpp + ../../../source_cell/unitcell_index.cpp + ../../../source_cell/unitcell_stats.cpp + ../../../source_cell/unitcell_setup.cpp ../../../source_cell/atom_spec.cpp ../../../source_cell/atom_pseudo.cpp ../../../source_cell/pseudo.cpp @@ -60,6 +63,9 @@ AddTest( LIBS parameter base device planewave SOURCES structure_factor_test.cpp ../stru_fac.cpp ../../../source_cell/unitcell.cpp + ../../../source_cell/unitcell_index.cpp + ../../../source_cell/unitcell_stats.cpp + ../../../source_cell/unitcell_setup.cpp ../../../source_cell/update_cell.cpp ../../../source_cell/bcast_cell.cpp From 4e7dacaa2b51a934ad55be5200d8a4b47b9d2422 Mon Sep 17 00:00:00 2001 From: DanielZhao0432 Date: Sat, 10 Oct 2026 12:05:00 +0800 Subject: [PATCH 2/3] test(cell): add unit tests for UnitCell::set_iat2itia and set_iat2iwt --- source/source_cell/test/CMakeLists.txt | 6 ++ .../source_cell/test/test_unitcell_index.cpp | 75 +++++++++++++++++++ 2 files changed, 81 insertions(+) create mode 100644 source/source_cell/test/test_unitcell_index.cpp diff --git a/source/source_cell/test/CMakeLists.txt b/source/source_cell/test/CMakeLists.txt index dc436a829d4..f10c87d0f5e 100644 --- a/source/source_cell/test/CMakeLists.txt +++ b/source/source_cell/test/CMakeLists.txt @@ -182,6 +182,12 @@ AddTest( SOURCES unitcell_test.cpp ../cal_ux.cpp ) +AddTest( + TARGET MODULE_CELL_unitcell_index_test + LIBS base device cell_info symmetry + SOURCES test_unitcell_index.cpp ../magnetism.cpp +) + AddTest( TARGET MODULE_CELL_print_cell_test LIBS base device cell_info symmetry diff --git a/source/source_cell/test/test_unitcell_index.cpp b/source/source_cell/test/test_unitcell_index.cpp new file mode 100644 index 00000000000..478fc66b04e --- /dev/null +++ b/source/source_cell/test/test_unitcell_index.cpp @@ -0,0 +1,75 @@ +#include "../unitcell.h" + +#include + +class TestUnitCellIndex : public ::testing::Test +{ + protected: + void SetUp() override + { + // Build a minimal UnitCell: + // 2 species, 2 atoms of the first species, 1 atom of the second + ucell.ntype = 2; + ucell.nat = 3; + ucell.atoms = new Atom[2]; + ucell.atoms[0].na = 2; + ucell.atoms[0].nw = 4; // 4 orbitals + ucell.atoms[1].na = 1; + ucell.atoms[1].nw = 2; // 2 orbitals + ucell.set_atom_flag = true; + } + + 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). + } + + UnitCell ucell; +}; + +TEST_F(TestUnitCellIndex, SetIat2itia) +{ + ucell.set_iat2itia(); + // iat2it: [0, 0, 1] + EXPECT_EQ(ucell.iat2it[0], 0); + EXPECT_EQ(ucell.iat2it[1], 0); + EXPECT_EQ(ucell.iat2it[2], 1); + // iat2ia: [0, 1, 0] + EXPECT_EQ(ucell.iat2ia[0], 0); + EXPECT_EQ(ucell.iat2ia[1], 1); + EXPECT_EQ(ucell.iat2ia[2], 0); +} + +TEST_F(TestUnitCellIndex, SetIat2itiaCalledTwice) +{ + // Cover the delete[] + new[] reallocation path: the second call must + // produce the same mapping as the first. + ucell.set_iat2itia(); + ucell.set_iat2itia(); + EXPECT_EQ(ucell.iat2it[0], 0); + EXPECT_EQ(ucell.iat2it[1], 0); + EXPECT_EQ(ucell.iat2it[2], 1); + EXPECT_EQ(ucell.iat2ia[0], 0); + EXPECT_EQ(ucell.iat2ia[1], 1); + EXPECT_EQ(ucell.iat2ia[2], 0); +} + +TEST_F(TestUnitCellIndex, SetIat2iwtNpol1) +{ + ucell.set_iat2iwt(1); + // iat=0: iwt=0, iat=1: iwt=4, iat=2: iwt=8 + EXPECT_EQ(ucell.get_iat2iwt()[0], 0); + EXPECT_EQ(ucell.get_iat2iwt()[1], 4); + EXPECT_EQ(ucell.get_iat2iwt()[2], 8); +} + +TEST_F(TestUnitCellIndex, SetIat2iwtNpol2) +{ + ucell.set_iat2iwt(2); + // iat=0: iwt=0, iat=1: iwt=8, iat=2: iwt=16 + EXPECT_EQ(ucell.get_iat2iwt()[0], 0); + EXPECT_EQ(ucell.get_iat2iwt()[1], 8); + EXPECT_EQ(ucell.get_iat2iwt()[2], 16); +} From d85a375bfb93f64aad1a0c97d374a89018328218 Mon Sep 17 00:00:00 2001 From: DanielZhao0432 Date: Sat, 10 Oct 2026 14:38:21 +0800 Subject: [PATCH 3/3] fix(build): register unitcell_index/stats/setup objects in Makefile build --- source/Makefile.Objects | 3 +++ 1 file changed, 3 insertions(+) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index f34c7ca601f..150070b799b 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -212,6 +212,9 @@ OBJS_CELL=atom_pseudo.o\ read_pp_blps.o\ read_pp_vwr.o\ unitcell.o\ + unitcell_index.o\ + unitcell_stats.o\ + unitcell_setup.o\ read_atoms.o\ read_atoms_helper.o\ print_cell.o\