diff --git a/source/Makefile.Objects b/source/Makefile.Objects index acbe4e2bb8..6aeaa0c75a 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -198,6 +198,7 @@ OBJS_CELL=atom_pseudo.o\ klist.o\ k_vector_utils.o\ cell_index.o\ + cell_tools.o\ check_atomic_stru.o\ update_cell.o\ bcast_cell.o\ diff --git a/source/source_cell/CMakeLists.txt b/source/source_cell/CMakeLists.txt index 2afc796896..ca64ccf60d 100644 --- a/source/source_cell/CMakeLists.txt +++ b/source/source_cell/CMakeLists.txt @@ -22,6 +22,7 @@ add_library( klist.cpp parallel_kpoints.cpp cell_index.cpp + cell_tools.cpp check_atomic_stru.cpp update_cell.cpp magnetism.cpp diff --git a/source/source_cell/cell_tools.cpp b/source/source_cell/cell_tools.cpp new file mode 100644 index 0000000000..6c380927cb --- /dev/null +++ b/source/source_cell/cell_tools.cpp @@ -0,0 +1,117 @@ +/** + * @file cell_tools.cpp + * @brief Implementation of cell tool free functions. + */ +#include "cell_tools.h" + +namespace unitcell +{ + std::vector get_atomLabels(const Atom* atoms, const int ntype) + { + std::vector atomLabels(ntype); + for (int it = 0; it < ntype; it++) + { + atomLabels[it] = atoms[it].label; + } + return atomLabels; + } + + std::vector get_atomCounts(const Atom* atoms, const int ntype) + { + std::vector atomCounts(ntype); + for (int it = 0; it < ntype; it++) + { + atomCounts[it] = atoms[it].na; + } + return atomCounts; + } + + std::vector> get_lnchiCounts(const Atom* atoms, const int ntype) + { + std::vector> lnchiCounts(ntype); + for (int it = 0; it < ntype; it++) + { + lnchiCounts[it].resize(atoms[it].nwl + 1); + for (int L = 0; L < atoms[it].nwl + 1; L++) + { + lnchiCounts[it][L] = atoms[it].l_nchi[L]; + } + } + return lnchiCounts; + } + + std::vector> get_target_mag(const Atom* atoms, + const int ntype, + const int nat) + { + std::vector> target_mag(nat); + int iat = 0; + for (int it = 0; it < ntype; it++) + { + for (int ia = 0; ia < atoms[it].na; ia++) + { + target_mag[iat] = atoms[it].m_loc_[ia]; + ++iat; + } + } + return target_mag; + } + + std::vector> get_lambda(const Atom* atoms, + const int ntype, + const int nat) + { + std::vector> lambda(nat); + int iat = 0; + for (int it = 0; it < ntype; it++) + { + for (int ia = 0; ia < atoms[it].na; ia++) + { + lambda[iat] = atoms[it].lambda[ia]; + ++iat; + } + } + return lambda; + } + + std::vector> get_constrain(const Atom* atoms, + const int ntype, + const int nat) + { + std::vector> constrain(nat); + int iat = 0; + for (int it = 0; it < ntype; it++) + { + for (int ia = 0; ia < atoms[it].na; ia++) + { + constrain[iat] = atoms[it].constrain[ia]; + ++iat; + } + } + return constrain; + } + + bool if_atoms_can_move(const Atom* atoms, const int ntype) + { + for (int it = 0; it < ntype; it++) + { + for (int ia = 0; ia < atoms[it].na; ia++) + { + if (atoms[it].mbl[ia].x || atoms[it].mbl[ia].y || atoms[it].mbl[ia].z) + { + return true; + } + } + } + return false; + } + + bool if_cell_can_change(const std::vector& lat_axis_free) + { + if (lat_axis_free[0] || lat_axis_free[1] || lat_axis_free[2]) + { + return true; + } + return false; + } +} diff --git a/source/source_cell/cell_tools.h b/source/source_cell/cell_tools.h new file mode 100644 index 0000000000..226416166a --- /dev/null +++ b/source/source_cell/cell_tools.h @@ -0,0 +1,76 @@ +/** + * @file cell_tools.h + * @brief Free function tools for extracting cell/atom information. + */ +#ifndef CELL_TOOLS_H +#define CELL_TOOLS_H + +#include +#include + +#include "source_base/vector3.h" +#include "source_cell/atom_spec.h" + +/** + * @brief Free functions for extracting atom/orbital info from Atom array. + */ +namespace unitcell +{ + /// @brief Get atom labels for each atom type. + /// @param atoms atom pointer [in] + /// @param ntype number of atom types [in] + /// @return vector of atom labels, one per type + std::vector get_atomLabels(const Atom* atoms, const int ntype); + + /// @brief Get atom counts (number of atoms) for each atom type. + /// @param atoms atom pointer [in] + /// @param ntype number of atom types [in] + /// @return vector of atom counts, one per type + std::vector get_atomCounts(const Atom* atoms, const int ntype); + + /// @brief Get lnchi counts (number of chi functions per L) for each atom type. + /// @param atoms atom pointer [in] + /// @param ntype number of atom types [in] + /// @return vector of lnchi counts, one vector per type + std::vector> get_lnchiCounts(const Atom* atoms, const int ntype); + + /// @brief Get target magnetic moment for each atom (used by deltaspin). + /// @param atoms atom pointer [in] + /// @param ntype number of atom types [in] + /// @param nat total number of atoms [in] + /// @return vector of target magnetic moments, one per atom + std::vector> get_target_mag(const Atom* atoms, + const int ntype, + const int nat); + + /// @brief Get Lagrange multiplier for each atom (used by deltaspin). + /// @param atoms atom pointer [in] + /// @param ntype number of atom types [in] + /// @param nat total number of atoms [in] + /// @return vector of Lagrange multipliers, one per atom + std::vector> get_lambda(const Atom* atoms, + const int ntype, + const int nat); + + /// @brief Get constrain flag for each atom (used by deltaspin). + /// @param atoms atom pointer [in] + /// @param ntype number of atom types [in] + /// @param nat total number of atoms [in] + /// @return vector of constrain flags, one per atom + std::vector> get_constrain(const Atom* atoms, + const int ntype, + const int nat); + + /// @brief Judge if any atom can move (any mbl component is non-zero). + /// @param atoms atom pointer [in] + /// @param ntype number of atom types [in] + /// @return true if at least one atom is allowed to move + bool if_atoms_can_move(const Atom* atoms, const int ntype); + + /// @brief Judge if any lattice vector can change. + /// @param lat_axis_free lattice-axis freedom flags (size 3) [in] + /// @return true if at least one lattice axis is free to change + bool if_cell_can_change(const std::vector& lat_axis_free); +} + +#endif // CELL_TOOLS_H diff --git a/source/source_cell/module_neighbor/test/prepare_unitcell.h b/source/source_cell/module_neighbor/test/prepare_unitcell.h index 9cabf07f9e..92a0506b1a 100644 --- a/source/source_cell/module_neighbor/test/prepare_unitcell.h +++ b/source/source_cell/module_neighbor/test/prepare_unitcell.h @@ -71,14 +71,12 @@ class UcellTestPrepare //basic info this->ntype = this->elements.size(); UnitCell* ucell = new UnitCell; - ucell->setup(this->latname, + ucell->setup_from_input(this->latname, this->ntype, this->lmaxmax, this->init_vel, this->fixed_axes); - ucell->atom_label.resize(ucell->ntype); - ucell->atom_mass.resize(ucell->ntype); ucell->pseudo_fn.resize(ucell->ntype); ucell->pseudo_type.resize(ucell->ntype); ucell->orbital_fn.resize(ucell->ntype); @@ -87,8 +85,6 @@ class UcellTestPrepare ucell->magnet.ux_[2] = 0.0; for(int it=0;itntype;++it) { - ucell->atom_label[it] = this->elements[it]; - ucell->atom_mass[it] = this->atomic_mass[it]; ucell->pseudo_fn[it] = this->pp_files[it]; ucell->pseudo_type[it] = this->pp_types[it]; ucell->orbital_fn[it] = this->orb_files[it]; @@ -148,7 +144,7 @@ class UcellTestPrepare ucell->atoms[it].angle2.resize(ucell->atoms[it].na); ucell->atoms[it].m_loc_.resize(ucell->atoms[it].na); ucell->atoms[it].mbl.resize(ucell->atoms[it].na); - ucell->atoms[it].mass = ucell->atom_mass[it]; // mass set here + ucell->atoms[it].mass = this->atomic_mass[it]; for(int ia=0; iaatoms[it].na; ++ia) { diff --git a/source/source_cell/module_symmetry/symm_magnetic.cpp b/source/source_cell/module_symmetry/symm_magnetic.cpp index 168183c4ce..599f612655 100644 --- a/source/source_cell/module_symmetry/symm_magnetic.cpp +++ b/source/source_cell/module_symmetry/symm_magnetic.cpp @@ -2,7 +2,7 @@ using namespace ModuleSymmetry; #include "symmetry_rotation_spin.h" -#include "source_io/module_parameter/parameter.h" +#include "source_base/global_variable.h" #include #include diff --git a/source/source_cell/print_cell.cpp b/source/source_cell/print_cell.cpp index 0c71b8bdbb..ae55fb727b 100644 --- a/source/source_cell/print_cell.cpp +++ b/source/source_cell/print_cell.cpp @@ -5,6 +5,7 @@ #include "source_base/formatter.h" #include "source_base/tool_title.h" #include "source_base/global_variable.h" +#include "source_base/output.h" namespace unitcell { @@ -98,8 +99,8 @@ namespace unitcell for(int it=0; it> ucell.atom_label[i] >> ucell.atom_mass[i]; + ss >> ucell.atoms[i].label >> ucell.atoms[i].mass; ucell.pseudo_fn[i] = "auto"; ucell.pseudo_type[i] = "auto"; @@ -73,8 +73,8 @@ bool read_atom_species(std::ifstream& ifa, // Peize Lin test for bsse 2021.04.07 const std::string bsse_label = "empty"; ucell.atoms[i].flag_empty_element = - (search( ucell.atom_label[i].begin(), ucell.atom_label[i].end(), - bsse_label.begin(), bsse_label.end() ) != ucell.atom_label[i].end()) + (search( ucell.atoms[i].label.begin(), ucell.atoms[i].label.end(), + bsse_label.begin(), bsse_label.end() ) != ucell.atoms[i].label.end()) ? true : false; } } diff --git a/source/source_cell/read_atoms.cpp b/source/source_cell/read_atoms.cpp index 359f8b6202..d3991d662e 100644 --- a/source/source_cell/read_atoms.cpp +++ b/source/source_cell/read_atoms.cpp @@ -72,7 +72,7 @@ bool unitcell::read_atom_positions(UnitCell& ucell, if (na > 0) { - unitcell::allocate_atom_properties(ucell.atoms[it], na, ucell.atom_mass[it]); + unitcell::allocate_atom_properties(ucell.atoms[it], na); for (int ia = 0;ia < na; ia++) { // modify the reading of frozen ions and velocities -- Yuanbo Li 2021/8/20 diff --git a/source/source_cell/read_atoms_helper.cpp b/source/source_cell/read_atoms_helper.cpp index 896311ae70..8e5e75e6c9 100644 --- a/source/source_cell/read_atoms_helper.cpp +++ b/source/source_cell/read_atoms_helper.cpp @@ -5,6 +5,7 @@ #include "read_stru.h" #include "print_cell.h" #include "read_orb.h" +#include "cell_tools.h" #include #include #include @@ -47,7 +48,7 @@ bool validate_coordinate_system(const std::string& Coordinate, return true; } -void allocate_atom_properties(Atom& atom, int na, double mass) +void allocate_atom_properties(Atom& atom, int na) { atom.tau.resize(na, ModuleBase::Vector3(0,0,0)); atom.dis.resize(na, ModuleBase::Vector3(0,0,0)); @@ -61,7 +62,6 @@ void allocate_atom_properties(Atom& atom, int na, double mass) atom.m_loc_.resize(na, ModuleBase::Vector3(0,0,0)); atom.lambda.resize(na, ModuleBase::Vector3(0,0,0)); atom.constrain.resize(na, ModuleBase::Vector3(0,0,0)); - atom.mass = mass; } void set_atom_movement_flags(Atom& atom, int ia, @@ -138,7 +138,7 @@ bool finalize_atom_positions(UnitCell& ucell, const std::string& esolver_type) { // Check if any atom can move in MD - if(!ucell.if_atoms_can_move() && calculation=="md" && esolver_type!="tddft") + if(!unitcell::if_atoms_can_move(ucell.atoms, ucell.ntype) && calculation=="md" && esolver_type!="tddft") { ModuleBase::WARNING("read_atoms", "no atoms can move in MD simulations!"); return false; @@ -502,13 +502,14 @@ bool read_atom_type_header(int it, UnitCell& ucell, // (1) read in atom label // start magnetization //======================================= + const std::string label_from_species = ucell.atoms[it].label; ModuleBase::GlobalFunc::READ_VALUE(ifpos, ucell.atoms[it].label); - if(ucell.atoms[it].label != ucell.atom_label[it]) + if(ucell.atoms[it].label != label_from_species) { ofs_warning << " Label orders in ATOMIC_POSITIONS and ATOMIC_SPECIES sections do not match!" << std::endl; ofs_warning << " Label read from ATOMIC_POSITIONS is " << ucell.atoms[it].label << std::endl; - ofs_warning << " Label from ATOMIC_SPECIES is " << ucell.atom_label[it] << std::endl; + ofs_warning << " Label from ATOMIC_SPECIES is " << label_from_species << std::endl; return false; } ModuleBase::GlobalFunc::OUT(ofs_running, "Atom label", ucell.atoms[it].label); diff --git a/source/source_cell/read_atoms_helper.h b/source/source_cell/read_atoms_helper.h index 505e5d1236..05d6091cd8 100644 --- a/source/source_cell/read_atoms_helper.h +++ b/source/source_cell/read_atoms_helper.h @@ -24,7 +24,7 @@ bool validate_coordinate_system(const std::string& Coordinate, * @param na Number of atoms * @param mass Atomic mass */ -void allocate_atom_properties(Atom& atom, int na, double mass); +void allocate_atom_properties(Atom& atom, int na); /** * @brief Set atom movement constraints based on fixed_atoms parameter diff --git a/source/source_cell/read_pseudo.cpp b/source/source_cell/read_pseudo.cpp index 4fcad45d13..a2210509b8 100644 --- a/source/source_cell/read_pseudo.cpp +++ b/source/source_cell/read_pseudo.cpp @@ -7,6 +7,8 @@ #include "cal_atoms_info.h" #include "read_pp.h" #include "bcast_cell.h" +#include "print_cell.h" +#include "source_base/atom_in.h" #include "source_base/element_elec_config.h" #include "source_base/parallel_common.h" @@ -60,7 +62,7 @@ AtomsInfoResult read_pseudo(std::ofstream& ofs, UnitCell& ucell, Atom* atom = &ucell.atoms[it]; if (!(atom->label_orb.empty())) { - ucell.compare_atom_labels(atom->label_orb, atom->ncpp.psd); + unitcell::compare_atom_labels(atom->label_orb, atom->ncpp.psd); } } @@ -393,7 +395,7 @@ void print_unitcell_pseudo(const std::string& fn, UnitCell& ucell) ModuleBase::TITLE("unitcell", "print_unitcell_pseudo"); std::ofstream ofs(fn.c_str()); - ucell.print_cell(ofs); + unitcell::print_cell(ucell, ofs); for (int i = 0; i < ucell.ntype; i++) { ucell.atoms[i].print_Atom(ofs); @@ -403,4 +405,77 @@ void print_unitcell_pseudo(const std::string& fn, UnitCell& ucell) return; } +void compare_atom_labels(const std::string& label1, const std::string& label2) +{ + if (label1!= label2) //'!( "Ag" == "Ag" || "47" == "47" || "Silver" == Silver" )' + { + atom_in ai; + if (!(std::to_string(ai.atom_Z[label1]) == label2 + || // '!( "Ag" == "47" )' + ai.atom_symbol[label1] == label2 || // '!( "Ag" == "Silver" )' + label1 == std::to_string(ai.atom_Z[label2]) + || // '!( "47" == "Ag" )' + label1 == std::to_string(ai.symbol_Z[label2]) + || // '!( "47" == "Silver" )' + label1 == ai.atom_symbol[label2] || // '!( "Silver" == "Ag" )' + std::to_string(ai.symbol_Z[label1]) + == label2)) // '!( "Silver" == "47" )' + { + std::string stru_label = ""; + std::string psuedo_label = ""; + for (int ip = 0; ip < label1.length(); ip++) + { + if (!(isdigit(label1[ip]) || label1[ip] == '_')) + { + stru_label += label1[ip]; + } + else + { + break; + } + } + stru_label[0] = toupper(stru_label[0]); + + for (int ip = 0; ip < label2.length(); ip++) + { + if (!(isdigit(label2[ip]) || label2[ip] == '_')) + { + psuedo_label += label2[ip]; + } + else + { + break; + } + } + psuedo_label[0] = toupper(psuedo_label[0]); + + if (!(stru_label == psuedo_label + || //' !("Ag1" == "ag_locpsp" || "47" == "47" || "Silver" == + //Silver" )' + std::to_string(ai.atom_Z[stru_label]) == psuedo_label + || // ' !("Ag1" == "47" )' + ai.atom_symbol[stru_label] == psuedo_label + || // ' !("Ag1" == "Silver")' + stru_label == std::to_string(ai.atom_Z[psuedo_label]) + || // ' !("47" == "Ag1" )' + stru_label == std::to_string(ai.symbol_Z[psuedo_label]) + || // ' !("47" == "Silver1" )' + stru_label == ai.atom_symbol[psuedo_label] + || // ' !("Silver1" == "Ag" )' + std::to_string(ai.symbol_Z[stru_label]) + == psuedo_label)) // ' !("Silver1" == "47" )' + + { + std::string atom_label_in_orbtial + = "atom label in orbital file "; + std::string mismatch_with_pseudo + = " mismatch with pseudo file of "; + ModuleBase::WARNING_QUIT("UnitCell::read_pseudo", + atom_label_in_orbtial + label1 + + mismatch_with_pseudo + label2); + } + } + } +} + } diff --git a/source/source_cell/read_pseudo.h b/source/source_cell/read_pseudo.h index 40b94ef945..f4512ade70 100644 --- a/source/source_cell/read_pseudo.h +++ b/source/source_cell/read_pseudo.h @@ -129,6 +129,15 @@ namespace unitcell { */ void cal_natomwfc(std::ofstream& log,int& natomwfc,const int ntype,const Atom* atoms,const int nspin); + /** + * @brief Check consistency between two atom labels from STRU and pseudo or + * orb file. + * + * @param label1 atom label from STRU [in] + * @param label2 atom label from pseudo or orbital file [in] + */ + void compare_atom_labels(const std::string& label1, const std::string& label2); + } #endif \ No newline at end of file diff --git a/source/source_cell/test/CMakeLists.txt b/source/source_cell/test/CMakeLists.txt index a67e86a40d..9ce431e3aa 100644 --- a/source/source_cell/test/CMakeLists.txt +++ b/source/source_cell/test/CMakeLists.txt @@ -46,6 +46,7 @@ list(APPEND cell_simple_srcs ../read_orb.cpp ../sep.cpp ../sep_cell.cpp + ../cell_tools.cpp ) add_library(cell_info OBJECT ${cell_simple_srcs}) @@ -109,6 +110,7 @@ AddTest( SOURCES read_atoms_helper_test.cpp ../read_atoms.cpp ../read_atoms_helper.cpp + ../cell_tools.cpp ../read_orb.cpp ../read_stru.cpp ../print_cell.cpp diff --git a/source/source_cell/test/prepare_unitcell.h b/source/source_cell/test/prepare_unitcell.h index 00e2c08019..b8118fd992 100644 --- a/source/source_cell/test/prepare_unitcell.h +++ b/source/source_cell/test/prepare_unitcell.h @@ -70,14 +70,12 @@ class UcellTestPrepare //basic info this->ntype = this->elements.size(); std::unique_ptr ucell(new UnitCell); - ucell->setup(this->latname, + ucell->setup_from_input(this->latname, this->ntype, this->lmaxmax, this->init_vel, this->fixed_axes); - ucell->atom_label.resize(ucell->ntype); - ucell->atom_mass.resize(ucell->ntype); ucell->pseudo_fn.resize(ucell->ntype); ucell->pseudo_type.resize(ucell->ntype); ucell->orbital_fn.resize(ucell->ntype); @@ -86,8 +84,6 @@ class UcellTestPrepare ucell->magnet.ux_[2] = 0.0; for(int it=0;itntype;++it) { - ucell->atom_label[it] = this->elements[it]; - ucell->atom_mass[it] = this->atomic_mass[it]; ucell->pseudo_fn[it] = this->pp_files[it]; ucell->pseudo_type[it] = this->pp_types[it]; ucell->orbital_fn[it] = this->orb_files[it]; @@ -149,7 +145,7 @@ class UcellTestPrepare ucell->atoms[it].mbl.resize(ucell->atoms[it].na); ucell->atoms[it].lambda.resize(ucell->atoms[it].na); ucell->atoms[it].constrain.resize(ucell->atoms[it].na); - ucell->atoms[it].mass = ucell->atom_mass[it]; // mass set here + ucell->atoms[it].mass = this->atomic_mass[it]; for(int ia=0; iaatoms[it].na; ++ia) { if (ucell->Coordinate == "Direct") diff --git a/source/source_cell/test/read_atoms_helper_test.cpp b/source/source_cell/test/read_atoms_helper_test.cpp index 7c232a4da4..77e6c49dc1 100644 --- a/source/source_cell/test/read_atoms_helper_test.cpp +++ b/source/source_cell/test/read_atoms_helper_test.cpp @@ -151,9 +151,9 @@ TEST_F(ReadAtomsHelperTest, AllocateAtomProperties) { Atom atom; int na = 5; - double mass = 12.0; + atom.mass = 12.0; - unitcell::allocate_atom_properties(atom, na, mass); + unitcell::allocate_atom_properties(atom, na); EXPECT_EQ(atom.tau.size(), na); EXPECT_EQ(atom.dis.size(), na); @@ -167,7 +167,7 @@ TEST_F(ReadAtomsHelperTest, AllocateAtomProperties) EXPECT_EQ(atom.m_loc_.size(), na); EXPECT_EQ(atom.lambda.size(), na); EXPECT_EQ(atom.constrain.size(), na); - EXPECT_DOUBLE_EQ(atom.mass, mass); + EXPECT_DOUBLE_EQ(atom.mass, 12.0); } // Test transform_atom_coordinates for Direct coordinates diff --git a/source/source_cell/test/sepcell_test.cpp b/source/source_cell/test/sepcell_test.cpp index c58933ac51..811c1277a2 100644 --- a/source/source_cell/test/sepcell_test.cpp +++ b/source/source_cell/test/sepcell_test.cpp @@ -63,9 +63,6 @@ class SepCellTest : public ::testing::Test // Initialize UnitCell for tests that need it. // This setup is common for many read_sep_potentials tests. ucell.ntype = 2; - ucell.atom_label.resize(ucell.ntype); - ucell.atom_label[0] = "Li"; - ucell.atom_label[1] = "F"; ucell.atoms = new Atom[ucell.ntype]; ucell.atoms[0].label = "Li"; ucell.atoms[0].na = 1; // Number of atoms of this type @@ -121,7 +118,8 @@ TEST_F(SepCellTest, ReadSepPotentialsSuccess) sep_cell.init(ucell.ntype); std::ofstream ofs_running_dummy("dummy_ofs_running.tmp"); - int result = sep_cell.read_sep_potentials(ifs, pp_dir, ofs_running_dummy, ucell.atom_label); + std::vector atom_labels = {ucell.atoms[0].label, ucell.atoms[1].label}; + int result = sep_cell.read_sep_potentials(ifs, pp_dir, ofs_running_dummy, atom_labels); ifs.close(); std::remove("dummy_ofs_running.tmp"); @@ -166,7 +164,8 @@ TEST_F(SepCellTest, ReadSepPotentialsNoSepFilesSection) std::ofstream ofs_running_dummy("dummy_ofs_running.tmp"); sep_cell.init(ucell.ntype); - int result = sep_cell.read_sep_potentials(ifs, pp_dir, ofs_running_dummy, ucell.atom_label); + std::vector atom_labels = {ucell.atoms[0].label, ucell.atoms[1].label}; + int result = sep_cell.read_sep_potentials(ifs, pp_dir, ofs_running_dummy, atom_labels); ifs.close(); std::remove("dummy_ofs_running.tmp"); @@ -189,7 +188,8 @@ TEST_F(SepCellTest, BcastSepCell) sep_cell.init(ucell.ntype); std::ofstream ofs_running_dummy("dummy_ofs_running.tmp"); - int result = sep_cell.read_sep_potentials(ifs, pp_dir, ofs_running_dummy, ucell.atom_label); + std::vector atom_labels = {ucell.atoms[0].label, ucell.atoms[1].label}; + int result = sep_cell.read_sep_potentials(ifs, pp_dir, ofs_running_dummy, atom_labels); ifs.close(); std::remove("dummy_ofs_running.tmp"); diff --git a/source/source_cell/test/support/mock_unitcell.cpp b/source/source_cell/test/support/mock_unitcell.cpp index 33dced94e4..0b761a34a7 100644 --- a/source/source_cell/test/support/mock_unitcell.cpp +++ b/source/source_cell/test/support/mock_unitcell.cpp @@ -20,8 +20,6 @@ SepPot::~SepPot(){} Sep_Cell::Sep_Cell() noexcept {} Sep_Cell::~Sep_Cell() noexcept {} -void UnitCell::print_cell(std::ofstream& ofs) const {} - void UnitCell::set_iat2itia() {} 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, @@ -30,11 +28,7 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log, const doubl const bool fixed_atoms, const bool noncolin, const std::string& calculation, const std::string& esolver_type, const int symmetry) {} -bool UnitCell::if_atoms_can_move() const { return true; } - -bool UnitCell::if_cell_can_change() const { return true; } - -void UnitCell::setup(const std::string& latname_in, +void UnitCell::setup_from_input(const std::string& latname_in, const int& ntype_in, const int& lmaxmax_in, const bool& init_vel_in, @@ -43,5 +37,3 @@ void UnitCell::setup(const std::string& latname_in, namespace unitcell { void cal_nelec(const Atom* atoms, const int& ntype, double& nelec, const double nelec_delta) {} } - -void UnitCell::compare_atom_labels(const std::string &label1, const std::string &label2) const {} diff --git a/source/source_cell/test/unitcell_test.cpp b/source/source_cell/test/unitcell_test.cpp index 210cb606ac..beb8d8fc69 100644 --- a/source/source_cell/test/unitcell_test.cpp +++ b/source/source_cell/test/unitcell_test.cpp @@ -5,6 +5,7 @@ #include "source_cell/read_orb.h" #include "source_cell/read_pseudo.h" #include "source_cell/read_stru.h" +#include "source_cell/cell_tools.h" #include "source_cell/print_cell.h" #include "memory" #include "source_cell/read_stru.h" @@ -36,7 +37,7 @@ Magnetism::~Magnetism() * - Constructor: * - UnitCell() and ~UnitCell() * - Setup: - * - setup(): to set latname, ntype, lmaxmax, init_vel, and lc + * - setup_from_input(): to set latname, ntype, lmaxmax, init_vel, and lc * - if_cell_can_change(): judge if any lattice vector can change * - RemakeCell * - remake_cell(): rebuild cell according to its latName @@ -49,10 +50,6 @@ Magnetism::~Magnetism() * - iat2iait(): depends on the above function, but can find both ia & it from iat * - ijat2iaitjajt(): find ia, it, ja, jt from ijat (ijat_max = nat*nat) * which collapses it, ia, jt, ja loop into a single loop - * - step_ia(): periodically set ia to 0 when ia reaches atom[it].na - 1 - * - step_it(): periodically set it to 0 when it reaches ntype -1 - * - step_iait(): return true only the above two conditions are true - * - step_jajtiait(): return ture only two of the above function (for i and j) are true * - GetAtomCounts * - get_atomCounts(): get atomCounts, which is a map from atom type to atom number * - GetOrbitalCounts @@ -73,7 +70,7 @@ Magnetism::~Magnetism() * - PrintTauCartesian * - print_tau(): print atomic coordinates, magmom and initial velocities * - PrintUnitcellPseudo - * - Actually an integrated function to call UnitCell::print_cell and Atom::print_Atom + * - Actually an integrated function to call unitcell::print_cell and Atom::print_Atom * - UpdateVel * - update_vel(const ModuleBase::Vector3* vel_in) * - CalUx @@ -167,7 +164,7 @@ TEST_F(UcellTest, Setup) std::vector fixed_axes_in = {"None", "volume", "shape", "a", "b", "c", "ab", "ac", "bc", "abc"}; for (int i = 0; i < fixed_axes_in.size(); ++i) { - ucell->setup(latname_in, ntype_in, lmaxmax_in, init_vel_in, fixed_axes_in[i]); + ucell->setup_from_input(latname_in, ntype_in, lmaxmax_in, init_vel_in, fixed_axes_in[i]); EXPECT_EQ(ucell->latName, latname_in); EXPECT_EQ(ucell->ntype, ntype_in); EXPECT_EQ(ucell->lmaxmax, lmaxmax_in); @@ -177,56 +174,56 @@ TEST_F(UcellTest, Setup) EXPECT_EQ(ucell->lat_axis_free[0], 1); EXPECT_EQ(ucell->lat_axis_free[1], 1); EXPECT_EQ(ucell->lat_axis_free[2], 1); - EXPECT_TRUE(ucell->if_cell_can_change()); + EXPECT_TRUE(unitcell::if_cell_can_change(ucell->lat_axis_free)); } else if (fixed_axes_in[i] == "a") { EXPECT_EQ(ucell->lat_axis_free[0], 0); EXPECT_EQ(ucell->lat_axis_free[1], 1); EXPECT_EQ(ucell->lat_axis_free[2], 1); - EXPECT_TRUE(ucell->if_cell_can_change()); + EXPECT_TRUE(unitcell::if_cell_can_change(ucell->lat_axis_free)); } else if (fixed_axes_in[i] == "b") { EXPECT_EQ(ucell->lat_axis_free[0], 1); EXPECT_EQ(ucell->lat_axis_free[1], 0); EXPECT_EQ(ucell->lat_axis_free[2], 1); - EXPECT_TRUE(ucell->if_cell_can_change()); + EXPECT_TRUE(unitcell::if_cell_can_change(ucell->lat_axis_free)); } else if (fixed_axes_in[i] == "c") { EXPECT_EQ(ucell->lat_axis_free[0], 1); EXPECT_EQ(ucell->lat_axis_free[1], 1); EXPECT_EQ(ucell->lat_axis_free[2], 0); - EXPECT_TRUE(ucell->if_cell_can_change()); + EXPECT_TRUE(unitcell::if_cell_can_change(ucell->lat_axis_free)); } else if (fixed_axes_in[i] == "ab") { EXPECT_EQ(ucell->lat_axis_free[0], 0); EXPECT_EQ(ucell->lat_axis_free[1], 0); EXPECT_EQ(ucell->lat_axis_free[2], 1); - EXPECT_TRUE(ucell->if_cell_can_change()); + EXPECT_TRUE(unitcell::if_cell_can_change(ucell->lat_axis_free)); } else if (fixed_axes_in[i] == "ac") { EXPECT_EQ(ucell->lat_axis_free[0], 0); EXPECT_EQ(ucell->lat_axis_free[1], 1); EXPECT_EQ(ucell->lat_axis_free[2], 0); - EXPECT_TRUE(ucell->if_cell_can_change()); + EXPECT_TRUE(unitcell::if_cell_can_change(ucell->lat_axis_free)); } else if (fixed_axes_in[i] == "bc") { EXPECT_EQ(ucell->lat_axis_free[0], 1); EXPECT_EQ(ucell->lat_axis_free[1], 0); EXPECT_EQ(ucell->lat_axis_free[2], 0); - EXPECT_TRUE(ucell->if_cell_can_change()); + EXPECT_TRUE(unitcell::if_cell_can_change(ucell->lat_axis_free)); } else if (fixed_axes_in[i] == "abc") { EXPECT_EQ(ucell->lat_axis_free[0], 0); EXPECT_EQ(ucell->lat_axis_free[1], 0); EXPECT_EQ(ucell->lat_axis_free[2], 0); - EXPECT_FALSE(ucell->if_cell_can_change()); + EXPECT_FALSE(unitcell::if_cell_can_change(ucell->lat_axis_free)); } } } @@ -239,14 +236,14 @@ TEST_F(UcellDeathTest, CompareAatomLabel) = {"Ag", "47", "Silver", "Ag", "47", "Silver", "Ag", "47", "Silver", "Ag1", "ag", "ag_locpsp", "Ag"}; for (int it = 0; it < 12; it++) { - ucell->compare_atom_labels(stru_label[it], pseudo_label[it]); + unitcell::compare_atom_labels(stru_label[it], pseudo_label[it]); } stru_label[0] = "Fe"; pseudo_label[0] = "O"; std::string atom_label_in_orbtial = "atom label in orbital file "; std::string mismatch_with_pseudo = " mismatch with pseudo file of "; testing::internal::CaptureStdout(); - EXPECT_EXIT(ucell->compare_atom_labels(stru_label[0], pseudo_label[0]), ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(unitcell::compare_atom_labels(stru_label[0], pseudo_label[0]), ::testing::ExitedWithCode(1), ""); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr(atom_label_in_orbtial + stru_label[0] + mismatch_with_pseudo + pseudo_label[0])); @@ -575,15 +572,11 @@ TEST_F(UcellTest, Index) int it_beg2; long long iat2 = ucell->nat + 1; EXPECT_FALSE(ucell->iat2iait(iat2, &ia_beg2, &it_beg2)); - // test ijat2iaitjajt, step_jajtiait, step_iat, step_ia, step_it + // test ijat2iaitjajt int ia_test; int it_test; int ja_test; int jt_test; - int ia_test2 = 0; - int it_test2 = 0; - int ja_test2 = 0; - int jt_test2 = 0; long long ijat = 0; for (int it = 0; it < utp.natom.size(); ++it) { @@ -599,15 +592,6 @@ TEST_F(UcellTest, Index) EXPECT_EQ(ja_test, ja); EXPECT_EQ(jt_test, jt); ++ijat; - if (it_test == utp.natom.size() - 1 && ia_test == utp.natom[it] - 1 - && jt_test == utp.natom.size() - 1 && ja_test == utp.natom[jt] - 1) - { - EXPECT_TRUE(ucell->step_jajtiait(&ja_test, &jt_test, &ia_test, &it_test)); - } - else - { - EXPECT_FALSE(ucell->step_jajtiait(&ja_test, &jt_test, &ia_test, &it_test)); - } } } } @@ -624,7 +608,7 @@ TEST_F(UcellTest, GetAtomCounts) EXPECT_EQ(atomCounts[0], 1); EXPECT_EQ(atomCounts[1], 2); /// atomCounts as vector - std::vector atomCounts2 = ucell->get_atomCounts(); + std::vector atomCounts2 = unitcell::get_atomCounts(ucell->atoms, ucell->ntype); EXPECT_EQ(atomCounts2[0], 1); EXPECT_EQ(atomCounts2[1], 2); } @@ -654,7 +638,7 @@ TEST_F(UcellTest, GetLnchiCounts) EXPECT_EQ(LnchiCounts[1][1], 1); EXPECT_EQ(LnchiCounts[1][2], 1); /// LnchiCounts as vector - std::vector> LnchiCounts2 = ucell->get_lnchiCounts(); + std::vector> LnchiCounts2 = unitcell::get_lnchiCounts(ucell->atoms, ucell->ntype); EXPECT_EQ(LnchiCounts2[0][0], 1); EXPECT_EQ(LnchiCounts2[0][1], 1); EXPECT_EQ(LnchiCounts2[0][2], 1); @@ -727,7 +711,7 @@ TEST_F(UcellTest, SelectiveDynamics) { UcellTestPrepare utp = UcellTestLib["C1H2-SD"]; ucell = utp.SetUcellInfo(); - EXPECT_TRUE(ucell->if_atoms_can_move()); + EXPECT_TRUE(unitcell::if_atoms_can_move(ucell->atoms, ucell->ntype)); } @@ -760,7 +744,7 @@ TEST_F(UcellTest, PrintCell) ucell = utp.SetUcellInfo(); std::ofstream ofs; ofs.open("printcell.log"); - ucell->print_cell(ofs); + unitcell::print_cell(*ucell, ofs); ofs.close(); std::ifstream ifs; ifs.open("printcell.log"); @@ -1049,8 +1033,6 @@ class UcellTestReadStru : public ::testing::Test void SetUp() override { ucell->ntype = 2; - ucell->atom_mass.resize(ucell->ntype); - ucell->atom_label.resize(ucell->ntype); ucell->pseudo_fn.resize(ucell->ntype); ucell->pseudo_type.resize(ucell->ntype); ucell->orbital_fn.resize(ucell->ntype); diff --git a/source/source_cell/test/unitcell_test_para.cpp b/source/source_cell/test/unitcell_test_para.cpp index 9a40ceb853..ebe3e9dbe5 100644 --- a/source/source_cell/test/unitcell_test_para.cpp +++ b/source/source_cell/test/unitcell_test_para.cpp @@ -7,6 +7,7 @@ #include "source_base/global_variable.h" #include "source_base/mathzone.h" #include "source_cell/unitcell.h" +#include "source_cell/cell_tools.h" #include "source_cell/read_pseudo.h" #include #include @@ -77,7 +78,7 @@ TEST_F(UcellTest, BcastUnitcell) EXPECT_EQ(ucell->atoms[0].na, 1); EXPECT_EQ(ucell->atoms[1].na, 2); /// this is to ensure all processes have the atom label info - auto atom_labels = ucell->get_atomLabels(); + auto atom_labels = unitcell::get_atomLabels(ucell->atoms, ucell->ntype); std::string atom_type1_expected = "C"; std::string atom_type2_expected = "H"; EXPECT_EQ(atom_labels[0], atom_type1_expected); @@ -94,7 +95,7 @@ TEST_F(UcellTest, BcastLattice) EXPECT_EQ(ucell->atoms[0].na, 1); EXPECT_EQ(ucell->atoms[1].na, 2); /// this is to ensure all processes have the atom label info - auto atom_labels = ucell->get_atomLabels(); + auto atom_labels = unitcell::get_atomLabels(ucell->atoms, ucell->ntype); std::string atom_type1_expected = "C"; std::string atom_type2_expected = "H"; EXPECT_EQ(atom_labels[0], atom_type1_expected); diff --git a/source/source_cell/test/unitcell_test_setupcell.cpp b/source/source_cell/test/unitcell_test_setupcell.cpp index 5e37bc096e..509a079c2a 100644 --- a/source/source_cell/test/unitcell_test_setupcell.cpp +++ b/source/source_cell/test/unitcell_test_setupcell.cpp @@ -65,8 +65,6 @@ class UcellTest : public ::testing::Test { ucell->lmaxmax = 2; ucell->ntype = 2; - ucell->atom_mass.resize(ucell->ntype); - ucell->atom_label.resize(ucell->ntype); ucell->pseudo_fn.resize(ucell->ntype); ucell->pseudo_type.resize(ucell->ntype); ucell->orbital_fn.resize(ucell->ntype); diff --git a/source/source_cell/test_pw/CMakeLists.txt b/source/source_cell/test_pw/CMakeLists.txt index b852cf2da3..d1cdbb5190 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 ../atom_spec.cpp ../update_cell.cpp ../bcast_cell.cpp + 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 ../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/test_pw/unitcell_test_pw.cpp b/source/source_cell/test_pw/unitcell_test_pw.cpp index d9c980d7e4..7a9c8d6862 100644 --- a/source/source_cell/test_pw/unitcell_test_pw.cpp +++ b/source/source_cell/test_pw/unitcell_test_pw.cpp @@ -57,8 +57,6 @@ class UcellTest : public ::testing::Test { ucell->lmaxmax = 2; ucell->ntype = 2; - ucell->atom_mass.resize(ucell->ntype); - ucell->atom_label.resize(ucell->ntype); ucell->pseudo_fn.resize(ucell->ntype); ucell->pseudo_type.resize(ucell->ntype); ucell->orbital_fn.resize(ucell->ntype); diff --git a/source/source_cell/unitcell.cpp b/source/source_cell/unitcell.cpp index e05b859dc2..710f1db31c 100644 --- a/source/source_cell/unitcell.cpp +++ b/source/source_cell/unitcell.cpp @@ -40,28 +40,6 @@ UnitCell::~UnitCell() } -void UnitCell::print_cell(std::ofstream& ofs) const { - - ModuleBase::GlobalFunc::OUT(ofs, "print_unitcell()"); - - ModuleBase::GlobalFunc::OUT(ofs, "latName", latName); - ModuleBase::GlobalFunc::OUT(ofs, "ntype", ntype); - ModuleBase::GlobalFunc::OUT(ofs, "nat", nat); - ModuleBase::GlobalFunc::OUT(ofs, "lat0", lat0); - ModuleBase::GlobalFunc::OUT(ofs, "lat0_angstrom", lat0_angstrom); - ModuleBase::GlobalFunc::OUT(ofs, "tpiba", tpiba); - ModuleBase::GlobalFunc::OUT(ofs, "omega", omega); - - output::printM3(ofs, "Lattices Vector (R) : ", latvec); - output::printM3(ofs, "Supercell lattice vector : ", latvec_supercell); - output::printM3(ofs, "Reciprocal lattice Vector (G): ", G); - output::printM3(ofs, "GGT : ", GGT); - - ofs << std::endl; - return; -} - - void UnitCell::set_iat2itia() { assert(nat > 0); delete[] iat2it; @@ -112,83 +90,15 @@ std::map> UnitCell::get_lnchi_Counts() const { return lnchiCounts; } -std::vector UnitCell::get_atomLabels() const { - std::vector atomLabels(this->ntype); - for (int it = 0; it < this->ntype; it++) { - atomLabels[it] = this->atoms[it].label; - } - return atomLabels; -} - -std::vector UnitCell::get_atomCounts() const { - std::vector atomCounts(this->ntype); - for (int it = 0; it < this->ntype; it++) { - atomCounts[it] = this->atoms[it].na; - } - return atomCounts; -} - -std::vector> UnitCell::get_lnchiCounts() const { - std::vector> lnchiCounts(this->ntype); - for (int it = 0; it < this->ntype; it++) { - lnchiCounts[it].resize(this->atoms[it].nwl + 1); - for (int L = 0; L < this->atoms[it].nwl + 1; L++) { - lnchiCounts[it][L] = this->atoms[it].l_nchi[L]; - } - } - return lnchiCounts; -} - -std::vector> UnitCell::get_target_mag() const -{ - std::vector> target_mag(this->nat); - for (int it = 0; it < this->ntype; it++) - { - for (int ia = 0; ia < this->atoms[it].na; ia++) - { - int iat = itia2iat(it, ia); - target_mag[iat] = this->atoms[it].m_loc_[ia]; - } - } - return target_mag; -} - -std::vector> UnitCell::get_lambda() const -{ - std::vector> lambda(this->nat); - for (int it = 0; it < this->ntype; it++) - { - for (int ia = 0; ia < this->atoms[it].na; ia++) - { - int iat = itia2iat(it, ia); - lambda[iat] = this->atoms[it].lambda[ia]; - } - } - return lambda; -} - -std::vector> UnitCell::get_constrain() const -{ - std::vector> constrain(this->nat); - for (int it = 0; it < this->ntype; it++) - { - for (int ia = 0; ia < this->atoms[it].na; ia++) - { - int iat = itia2iat(it, ia); - constrain[iat] = this->atoms[it].constrain[ia]; - } - } - return constrain; -} - //============================================================== // 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, +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) + const bool fixed_atoms, const bool noncolin, const std::string& calculation, + const std::string& esolver_type, const int symmetry) { ModuleBase::TITLE("UnitCell", "setup_cell"); @@ -207,8 +117,6 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log, const doubl bool ok3 = true; // for sep potential in DFT-1/2 // (3) read in atom information - this->atom_mass.resize(ntype); - this->atom_label.resize(ntype); this->pseudo_fn.resize(ntype); this->pseudo_type.resize(ntype); this->orbital_fn.resize(ntype); @@ -256,7 +164,12 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log, const doubl //========================== if (dfthalf_type > 0) { sep_cell.init(this->ntype); - ok3 = sep_cell.read_sep_potentials(ifa, pseudo_dir, GlobalV::ofs_warning, this->atom_label); + 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 @@ -330,12 +243,6 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log, const doubl this->GGT = G * GT; this->invGGT = GGT.Inverse(); - // LiuXh add 20180515 - this->GT0 = latvec.Inverse(); - this->G0 = GT.Transpose(); - this->GGT0 = G * GT; - this->invGGT0 = GGT.Inverse(); - log << std::endl; output::printM3(log, "Lattice vectors: (Cartesian coordinate: in unit of a_0)", @@ -382,35 +289,7 @@ void UnitCell::set_iat2iwt(const int& npol_in) -// check if any atom can be moved -bool UnitCell::if_atoms_can_move() const -{ - for (int it = 0; it < this->ntype; it++) - { - Atom* atom = &atoms[it]; - for (int ia = 0; ia < atom->na; ia++) - { - if (atom->mbl[ia].x || atom->mbl[ia].y || atom->mbl[ia].z) - { - return true; - } - } - } - return false; -} - -// check if lattice vector can be changed -bool UnitCell::if_cell_can_change() const -{ - // need to be fixed next - if (this->lat_axis_free[0] || this->lat_axis_free[1] || this->lat_axis_free[2]) - { - return true; - } - return false; -} - -void UnitCell::setup(const std::string& latname_in, +void UnitCell::setup_from_input(const std::string& latname_in, const int& ntype_in, const int& lmaxmax_in, const bool& init_vel_in, @@ -468,77 +347,3 @@ void UnitCell::setup(const std::string& latname_in, } return; } - - -void UnitCell::compare_atom_labels(const std::string& label1, const std::string& label2) const -{ - if (label1!= label2) //'!( "Ag" == "Ag" || "47" == "47" || "Silver" == Silver" )' - { - atom_in ai; - if (!(std::to_string(ai.atom_Z[label1]) == label2 - || // '!( "Ag" == "47" )' - ai.atom_symbol[label1] == label2 || // '!( "Ag" == "Silver" )' - label1 == std::to_string(ai.atom_Z[label2]) - || // '!( "47" == "Ag" )' - label1 == std::to_string(ai.symbol_Z[label2]) - || // '!( "47" == "Silver" )' - label1 == ai.atom_symbol[label2] || // '!( "Silver" == "Ag" )' - std::to_string(ai.symbol_Z[label1]) - == label2)) // '!( "Silver" == "47" )' - { - std::string stru_label = ""; - std::string psuedo_label = ""; - for (int ip = 0; ip < label1.length(); ip++) - { - if (!(isdigit(label1[ip]) || label1[ip] == '_')) - { - stru_label += label1[ip]; - } - else - { - break; - } - } - stru_label[0] = toupper(stru_label[0]); - - for (int ip = 0; ip < label2.length(); ip++) - { - if (!(isdigit(label2[ip]) || label2[ip] == '_')) - { - psuedo_label += label2[ip]; - } - else - { - break; - } - } - psuedo_label[0] = toupper(psuedo_label[0]); - - if (!(stru_label == psuedo_label - || //' !("Ag1" == "ag_locpsp" || "47" == "47" || "Silver" == - //Silver" )' - std::to_string(ai.atom_Z[stru_label]) == psuedo_label - || // ' !("Ag1" == "47" )' - ai.atom_symbol[stru_label] == psuedo_label - || // ' !("Ag1" == "Silver")' - stru_label == std::to_string(ai.atom_Z[psuedo_label]) - || // ' !("47" == "Ag1" )' - stru_label == std::to_string(ai.symbol_Z[psuedo_label]) - || // ' !("47" == "Silver1" )' - stru_label == ai.atom_symbol[psuedo_label] - || // ' !("Silver1" == "Ag" )' - std::to_string(ai.symbol_Z[stru_label]) - == psuedo_label)) // ' !("Silver1" == "47" )' - - { - std::string atom_label_in_orbtial - = "atom label in orbital file "; - std::string mismatch_with_pseudo - = " mismatch with pseudo file of "; - ModuleBase::WARNING_QUIT("UnitCell::read_pseudo", - atom_label_in_orbtial + label1 - + mismatch_with_pseudo + label2); - } - } - } -} diff --git a/source/source_cell/unitcell.h b/source/source_cell/unitcell.h index 5113dd4cd7..2be50a9877 100644 --- a/source/source_cell/unitcell.h +++ b/source/source_cell/unitcell.h @@ -15,6 +15,11 @@ */ class UnitCell : public AtomProvider, public BaseCell { public: + UnitCell(); + ~UnitCell(); + + /// @name BaseCell / AtomProvider interface overrides + /// @{ double get_lat0() const override { return lat0; } @@ -42,72 +47,43 @@ class UnitCell : public AtomProvider, public BaseCell { ModuleBase::Vector3 get_tau(int i, int j) const override { return atoms[i].tau[j]; } + /// @} - Atom* atoms = nullptr; - Sep_Cell sep_cell; - - bool set_atom_flag = false; ///< added on 2009-3-8 by mohan - Magnetism magnet; ///< magnetism Yu Liu 2021-07-03 - std::vector> atom_mulliken; ///< [nat][nspin] - int n_mag_at = 0; - - Lattice lat; - std::string& Coordinate = lat.Coordinate; - std::string& latName = lat.latName; - double& lat0 = lat.lat0; - double& lat0_angstrom = lat.lat0_angstrom; - double& tpiba = lat.tpiba; - double& tpiba2 = lat.tpiba2; - double& omega = lat.omega; - std::vector& lat_axis_free = lat.lat_axis_free; + /// @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, + const int& ntype_in, + const int& lmaxmax_in, + const bool& init_vel_in, + const std::string& fixed_axes_in); - ModuleBase::Matrix3& latvec = lat.latvec; - ModuleBase::Vector3&a1 = lat.a1, &a2 = lat.a2, &a3 = lat.a3; - ModuleBase::Vector3& latcenter = lat.latcenter; - ModuleBase::Matrix3& latvec_supercell = lat.latvec_supercell; - ModuleBase::Matrix3& G = lat.G; - ModuleBase::Matrix3& GT = lat.GT; - ModuleBase::Matrix3& GGT = lat.GGT; - ModuleBase::Matrix3& invGGT = lat.invGGT; + void 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); - 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; + void set_iat2itia(); - ModuleSymmetry::Symmetry symm; + void set_iat2iwt(const int& npol_in); /// iat2iwt is the atom index iat to the first global index for orbital of /// this atom the size of iat2iwt is nat, the value should be /// sum_{i=0}^{iat-1} atoms[it].nw * npol where the npol is the number of /// polarizations, 1 for non-magnetic(NSPIN=1 or 2), 2 for magnetic(only /// NSPIN=4) this part only used for Atomic Orbital based calculation - public: /// @brief Indexing tool for find orbital global index from it,ia,iw template inline Tiait itiaiw2iwt(const Tiait& it, const Tiait& ia, const Tiait& iw) const { return Tiait(this->iat2iwt[this->itia2iat(it, ia)] + iw); } - /// @brief Initialize iat2iwt - void set_iat2iwt(const int& npol_in); /// @brief Get iat2iwt inline const int* get_iat2iwt() const { return iat2iwt.data(); } /// @brief Get npol inline const int& get_npol() const { return npol; } - private: - std::vector iat2iwt; ///< iat ==> iwt, the first global index for orbital of this atom - int npol = 1; ///< number of spin polarizations, initialized in set_iat2iwt - /// ----------------- END of iat2iwt part ----------------- - - public: /// @brief Indexing tools for ia and it /// @return true if the last out is reset template @@ -135,41 +111,6 @@ class UnitCell : public AtomProvider, public BaseCell { return true; } - template - inline bool step_it(Tiait* it) const { - if (++(*it) >= ntype) { - *it = 0; - return true; - } - return false; - } - - template - inline bool step_ia(const Tiait it, Tiait* ia) const { - if (++(*ia) >= atoms[it].na) { - *ia = 0; - return true; - } - return false; - } - - template - inline bool step_iait(Tiait* ia, Tiait* it) const { - if (step_ia(*it, ia)) { - return step_it(it); - } - return false; - } - - template - inline bool - step_jajtiait(Tiait* ja, Tiait* jt, Tiait* ia, Tiait* it) const { - if (step_iait(ja, jt)) { - return step_iait(ia, it); - } - return false; - } - /// @brief Get tau for atom iat inline const ModuleBase::Vector3& get_tau(const int& iat) const { return atoms[iat2it[iat]].tau[iat2ia[iat]]; @@ -184,26 +125,78 @@ class UnitCell : public AtomProvider, public BaseCell { + double(R.z) * a3 - get_tau(iat1); } - /// @brief LiuXh add 20180515 - ModuleBase::Matrix3 G0; - ModuleBase::Matrix3 GT0; - ModuleBase::Matrix3 GGT0; - ModuleBase::Matrix3 invGGT0; + /// @brief get atomCounts, which is a map from element type to atom number + std::map get_atom_Counts() const; + /// @brief get orbitalCounts, which is a map from element type to orbital + /// number + std::map get_orbital_Counts() const; + /// @brief get lnchiCounts, which is a map from element type to the l:nchi + /// map + std::map> get_lnchi_Counts() const; + /// these are newly added functions, the three above functions are + /// deprecated and will be removed in the future - /// @todo Encapsulate ionic_position_updated and cell_parameter_updated with - /// setters that enforce state invariants; currently exposed as mutable - /// flags that can be toggled from anywhere. - bool ionic_position_updated - = false; ///< whether the ionic position has been updated - bool cell_parameter_updated - = false; ///< whether the cell parameters are updated + public: + + /// @name Lattice + /// @{ + Lattice lat; + std::string& Coordinate = lat.Coordinate; + std::string& latName = lat.latName; + double& lat0 = lat.lat0; + double& lat0_angstrom = lat.lat0_angstrom; + double& tpiba = lat.tpiba; + double& tpiba2 = lat.tpiba2; + double& omega = lat.omega; + std::vector& lat_axis_free = lat.lat_axis_free; + + ModuleBase::Matrix3& latvec = lat.latvec; + ModuleBase::Vector3&a1 = lat.a1, &a2 = lat.a2, &a3 = lat.a3; + ModuleBase::Vector3& latcenter = lat.latcenter; + ModuleBase::Matrix3& latvec_supercell = lat.latvec_supercell; + ModuleBase::Matrix3& G = lat.G; + ModuleBase::Matrix3& GT = lat.GT; + ModuleBase::Matrix3& GGT = lat.GGT; + ModuleBase::Matrix3& invGGT = lat.invGGT; + /// @} + + /// @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; + /// @} + + Atom* atoms = nullptr; + Sep_Cell sep_cell; + + /// @name Magnetism + /// @{ + Magnetism magnet; ///< magnetism Yu Liu 2021-07-03 + std::vector> atom_mulliken; ///< [nat][nspin] + int n_mag_at = 0; + /// @} + + /// @name Symmetry + /// @{ + ModuleSymmetry::Symmetry symm; + /// @} + /// @name Pseudopotential / Orbital parameters /// @brief meshx : max number of mesh point in pseudopotential file /// @brief natomwfc : number of starting wavefunctions /// @brief lmax : Max L used for localized orbital /// @brief nmax : Max N used for localized orbital /// @brief lmax_ppwf : Max L of pseudo wave functions /// @brief lmaxmax : revert from INPUT + /// @{ int meshx = 0; int natomwfc = 0; int lmax = 0; @@ -213,17 +206,10 @@ class UnitCell : public AtomProvider, public BaseCell { int lmaxmax = 0; ///< liuyu 2021-07-04 bool init_vel = false; ///< liuyu 2021-07-15 // double nelec; + /// @} - private: - ModuleBase::Matrix3 stress; ///< calculate stress on the cell - - public: - UnitCell(); - ~UnitCell(); - void print_cell(std::ofstream& ofs) const; - - std::vector atom_mass; - std::vector atom_label; + /// @name File lists + /// @{ std::vector pseudo_fn; std::vector pseudo_type; @@ -231,15 +217,22 @@ class UnitCell : public AtomProvider, public BaseCell { std::string descriptor_file; ///< filenames of descriptor_file, liuyu add 2023-04-06 std::vector abfs_orbital_files; ///< ABFS orbital filenames read from STRU "ABFS_ORBITAL" (used by LCAO EXX) std::vector jle_orbital_files; ///< JLE orbital filenames read from STRU "ABFS_JLES_ORBITAL" (used by LCAO EXX) + /// @} - void set_iat2itia(); - - void 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); + /// @name State flags + /// @todo Encapsulate ionic_position_updated and cell_parameter_updated with + /// setters that enforce state invariants; currently exposed as mutable + /// flags that can be toggled from anywhere. + /// @{ + bool set_atom_flag = false; ///< added on 2009-3-8 by mohan + bool ionic_position_updated + = false; ///< whether the ionic position has been updated + bool cell_parameter_updated + = false; ///< whether the cell parameters are updated + /// @} + /// @name Nonlocal info + /// @{ /** * @brief Pointer to non-local pseudopotential information. * @@ -247,49 +240,19 @@ class UnitCell : public AtomProvider, public BaseCell { * to non-local projector data. It is null for non-LCAO calculations. */ std::unique_ptr infoNL; + /// @} - /// @brief For constrained vc-relaxation where type of lattice is fixed, adjust the lattice vectors + private: + // --------------------- Private Data --------------------- - /// @brief cal_natomwfc : calculate total number of atomic wavefunctions - /// @brief cal_nwfc : calculate total number of local basis and lmax - /// @brief cal_meshx : calculate max number of mesh points in pp file - bool if_atoms_can_move() const; - bool if_cell_can_change() const; - void setup(const std::string& latname_in, - const int& ntype_in, - const int& lmaxmax_in, - const bool& init_vel_in, - const std::string& fixed_axes_in); + std::vector iat2iwt; ///< iat ==> iwt, the first global index for orbital of this atom + int npol = 1; ///< number of spin polarizations, initialized in set_iat2iwt + /// ----------------- END of iat2iwt part ----------------- - /// @brief check consistency between two atom labels from STRU and pseudo or - /// orb file - void compare_atom_labels(const std::string& label1, const std::string& label2) const; - /// @brief get atomCounts, which is a map from element type to atom number - std::map get_atom_Counts() const; - /// @brief get orbitalCounts, which is a map from element type to orbital - /// number - std::map get_orbital_Counts() const; - /// @brief get lnchiCounts, which is a map from element type to the l:nchi - /// map - std::map> get_lnchi_Counts() const; - /// these are newly added functions, the three above functions are - /// deprecated and will be removed in the future - /// @brief get atom labels - std::vector get_atomLabels() const; - /// @brief get atomCounts, which is a vector of element type with atom - /// number - std::vector get_atomCounts() const; - /// @brief get lnchiCounts, which is a vector of element type with the - /// l:nchi vector - std::vector> get_lnchiCounts() const; - /// @brief get target magnetic moment for deltaspin - std::vector> get_target_mag() const; - /// @brief get lagrange multiplier for deltaspin - std::vector> get_lambda() const; - /// @brief get constrain for deltaspin - std::vector> get_constrain() const; + ModuleBase::Matrix3 stress; ///< calculate stress on the cell - private: + /// @name BaseCell private overrides + /// @{ Kind get_kind() const override { return Kind::unit_cell; @@ -304,6 +267,7 @@ class UnitCell : public AtomProvider, public BaseCell { { return GT; } + /// @} }; #endif // unitcell class diff --git a/source/source_esolver/test/esolver_dp_test.cpp b/source/source_esolver/test/esolver_dp_test.cpp index f25b89ab75..c4daaef119 100644 --- a/source/source_esolver/test/esolver_dp_test.cpp +++ b/source/source_esolver/test/esolver_dp_test.cpp @@ -43,9 +43,8 @@ class ESolverDPTest : public ::testing::Test ucell.atoms[0].taud.resize(1, ModuleBase::Vector3(0.0, 0.0, 0.0)); ucell.atoms[1].taud.resize(1, ModuleBase::Vector3(0.0, 0.0, 0.0)); - ucell.atom_label.resize(ucell.ntype); - ucell.atom_label[0] = "Cu"; - ucell.atom_label[1] = "Al"; + ucell.atoms[0].label = "Cu"; + ucell.atoms[1].label = "Al"; esolver->before_all_runners(ucell, inp); } diff --git a/source/source_esolver/test/for_test.h b/source/source_esolver/test/for_test.h index ba8c9bb415..dd7448b956 100644 --- a/source/source_esolver/test/for_test.h +++ b/source/source_esolver/test/for_test.h @@ -22,12 +22,10 @@ UnitCell::UnitCell() ntype = 2; nat = 2; - atom_label.resize(ntype); - atom_label[0] = "Al"; - atom_label[1] = "Cu"; - atoms = new Atom[ntype]; set_atom_flag = true; + atoms[0].label = "Al"; + atoms[1].label = "Cu"; for (int it = 0; it < ntype; it++) { diff --git a/source/source_estate/module_dm/test/prepare_unitcell.h b/source/source_estate/module_dm/test/prepare_unitcell.h index ae8a582f00..43e00e8f76 100644 --- a/source/source_estate/module_dm/test/prepare_unitcell.h +++ b/source/source_estate/module_dm/test/prepare_unitcell.h @@ -72,11 +72,9 @@ class UcellTestPrepare // basic info this->ntype = this->elements.size(); static UnitCell ucell; - ucell.setup(this->latname, this->ntype, this->lmaxmax, this->init_vel, this->fixed_axes); + ucell.setup_from_input(this->latname, this->ntype, this->lmaxmax, this->init_vel, this->fixed_axes); - ucell.atom_label.resize(ucell.ntype); - ucell.atom_mass.resize(ucell.ntype); - ucell.pseudo_fn.resize(ucell.ntype); + ucell.pseudo_fn.resize(ucell.ntype); ucell.pseudo_type.resize(ucell.ntype); ucell.orbital_fn.resize(ucell.ntype); ucell.magnet.ux_[0] = 0.0; // ux_ set here @@ -84,8 +82,6 @@ class UcellTestPrepare ucell.magnet.ux_[2] = 0.0; for (int it = 0; it < ucell.ntype; ++it) { - ucell.atom_label[it] = this->elements[it]; - ucell.atom_mass[it] = this->atomic_mass[it]; ucell.pseudo_fn[it] = this->pp_files[it]; ucell.pseudo_type[it] = this->pp_types[it]; ucell.orbital_fn[it] = this->orb_files[it]; @@ -148,7 +144,7 @@ class UcellTestPrepare ucell.atoms[it].angle2.resize(ucell.atoms[it].na); ucell.atoms[it].m_loc_.resize(ucell.atoms[it].na); ucell.atoms[it].mbl.resize(ucell.atoms[it].na); - ucell.atoms[it].mass = ucell.atom_mass[it]; // mass set here + ucell.atoms[it].mass = this->atomic_mass[it]; for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { if (ucell.Coordinate == "Direct") diff --git a/source/source_estate/test/prepare_unitcell.h b/source/source_estate/test/prepare_unitcell.h index 5d0d19e0d0..9ec33230a5 100644 --- a/source/source_estate/test/prepare_unitcell.h +++ b/source/source_estate/test/prepare_unitcell.h @@ -55,13 +55,11 @@ class UcellTestPrepare //basic info this->ntype = this->elements.size(); std::unique_ptr ucell(new UnitCell); - ucell->setup(this->latname, + ucell->setup_from_input(this->latname, this->ntype, this->lmaxmax, this->init_vel, this->fixed_axes); - ucell->atom_label.resize(ucell->ntype); - ucell->atom_mass.resize(ucell->ntype); ucell->pseudo_fn.resize(ucell->ntype); ucell->pseudo_type.resize(ucell->ntype); ucell->orbital_fn.resize(ucell->ntype); @@ -70,8 +68,6 @@ class UcellTestPrepare ucell->magnet.ux_[2] = 0.0; for(int it=0;itntype;++it) { - ucell->atom_label[it] = this->elements[it]; - ucell->atom_mass[it] = this->atomic_mass[it]; ucell->pseudo_fn[it] = this->pp_files[it]; ucell->pseudo_type[it] = this->pp_types[it]; ucell->orbital_fn[it] = this->orb_files[it]; @@ -131,7 +127,7 @@ class UcellTestPrepare ucell->atoms[it].angle2.resize(ucell->atoms[it].na); ucell->atoms[it].m_loc_.resize(ucell->atoms[it].na); ucell->atoms[it].mbl.resize(ucell->atoms[it].na); - ucell->atoms[it].mass = ucell->atom_mass[it]; // mass set here + ucell->atoms[it].mass = this->atomic_mass[it]; for(int ia=0; iaatoms[it].na; ++ia) { if (ucell->Coordinate == "Direct") diff --git a/source/source_hamilt/module_hcontainer/test/prepare_unitcell.h b/source/source_hamilt/module_hcontainer/test/prepare_unitcell.h index 8708c6a41e..834a177bef 100644 --- a/source/source_hamilt/module_hcontainer/test/prepare_unitcell.h +++ b/source/source_hamilt/module_hcontainer/test/prepare_unitcell.h @@ -71,10 +71,8 @@ class UcellTestPrepare //basic info this->ntype = this->elements.size(); static UnitCell ucell; - ucell.setup(this->latname, this->ntype, this->lmaxmax, this->init_vel, this->fixed_axes); - ucell.atom_label.resize(ucell.ntype); - ucell.atom_mass.resize(ucell.ntype); - ucell.pseudo_fn.resize(ucell.ntype); + ucell.setup_from_input(this->latname, this->ntype, this->lmaxmax, this->init_vel, this->fixed_axes); + ucell.pseudo_fn.resize(ucell.ntype); ucell.pseudo_type.resize(ucell.ntype); ucell.orbital_fn.resize(ucell.ntype); @@ -83,8 +81,6 @@ class UcellTestPrepare ucell.magnet.ux_[2] = 0.0; for (int it = 0; it < ucell.ntype; ++it) { - ucell.atom_label[it] = this->elements[it]; - ucell.atom_mass[it] = this->atomic_mass[it]; ucell.pseudo_fn[it] = this->pp_files[it]; ucell.pseudo_type[it] = this->pp_types[it]; ucell.orbital_fn[it] = this->orb_files[it]; @@ -147,7 +143,7 @@ class UcellTestPrepare ucell.atoms[it].angle2.resize(ucell.atoms[it].na); ucell.atoms[it].m_loc_.resize(ucell.atoms[it].na); ucell.atoms[it].mbl.resize(ucell.atoms[it].na); - ucell.atoms[it].mass = ucell.atom_mass[it]; // mass set here + ucell.atoms[it].mass = this->atomic_mass[it]; for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { if (ucell.Coordinate == "Direct") diff --git a/source/source_io/module_json/init_info.cpp b/source/source_io/module_json/init_info.cpp index 83ff731c1e..9eafceb4bf 100644 --- a/source/source_io/module_json/init_info.cpp +++ b/source/source_io/module_json/init_info.cpp @@ -99,7 +99,6 @@ void gen_stru(UnitCell* ucell) // atom coordinate, mag and label const double lat0_angstrom = ucell->lat0_angstrom; - std::string* label = ucell->atom_label.data(); for (int i = 0; i < ntype; i++) { ModuleBase::Vector3* tau = ucell->atoms[i].tau.data(); @@ -117,7 +116,7 @@ void gen_stru(UnitCell* ucell) // Json::AbacusJson::add_Json(ucell->atoms[i].mag[j],true,"init","mag"); - std::string str = label[i]; + std::string str = ucell->atoms[i].label; Json::AbacusJson::add_json({"init", "label"}, str, true); // Json::AbacusJson::add_Json(str,true,"init","label"); } diff --git a/source/source_io/module_json/test/para_json_test.cpp b/source/source_io/module_json/test/para_json_test.cpp index ca98e61159..c3e4390973 100644 --- a/source/source_io/module_json/test/para_json_test.cpp +++ b/source/source_io/module_json/test/para_json_test.cpp @@ -325,7 +325,6 @@ TEST(AbacusJsonTest, Init_stru_test) ucell.pseudo_fn.resize(1); ucell.orbital_fn.resize(1); ucell.atoms = atomlist; - ucell.atom_label.resize(1); ucell.lat0 = lat0; ucell.lat0_angstrom = lat0 * ModuleBase::BOHR_TO_A; @@ -337,9 +336,8 @@ TEST(AbacusJsonTest, Init_stru_test) // fill ucell for (int i = 0; i < 1; i++) { - ucell.atom_label[i] = "Si"; + ucell.atoms[i].label = "Si"; atomlist[i].na = 2; - atomlist[i].label = "Fe"; ucell.pseudo_fn[i] = "si.ufp"; ucell.atoms[i].tau.resize(2); atomlist[i].mag.resize(2); @@ -358,9 +356,9 @@ TEST(AbacusJsonTest, Init_stru_test) ASSERT_EQ(Json::AbacusJson::doc["init"]["mag"][0].GetDouble(), 0); ASSERT_EQ(Json::AbacusJson::doc["init"]["mag"][1].GetDouble(), 131.0); - ASSERT_STREQ(Json::AbacusJson::doc["init"]["pp"]["Fe"].GetString(), "si.ufp"); + ASSERT_STREQ(Json::AbacusJson::doc["init"]["pp"]["Si"].GetString(), "si.ufp"); ASSERT_STREQ(Json::AbacusJson::doc["init"]["label"][0].GetString(), "Si"); - ASSERT_STREQ(Json::AbacusJson::doc["init"]["element"]["Fe"].GetString(), ""); + ASSERT_STREQ(Json::AbacusJson::doc["init"]["element"]["Si"].GetString(), ""); ASSERT_EQ(Json::AbacusJson::doc["init"]["coordinate"][0][0].GetDouble(), 0); ASSERT_EQ(Json::AbacusJson::doc["init"]["coordinate"][0][1].GetDouble(), 0); diff --git a/source/source_io/module_mulliken/cal_mag.h b/source/source_io/module_mulliken/cal_mag.h index b0ec41974a..63377231c7 100644 --- a/source/source_io/module_mulliken/cal_mag.h +++ b/source/source_io/module_mulliken/cal_mag.h @@ -5,6 +5,7 @@ #include "source_basis/module_ao/ORB_read.h" #include "source_basis/module_nao/two_center_bundle.h" #include "source_cell/cell_index.h" +#include "source_cell/cell_tools.h" #include "source_cell/klist.h" #include "source_cell/module_neighbor/sltk_grid_driver.h" #include "source_cell/unitcell.h" @@ -36,8 +37,9 @@ void cal_mag(Parallel_Orbitals* pv, if (PARAM.inp.out_mul) { auto cell_index - = CellIndex(ucell.get_atomLabels(), - ucell.get_atomCounts(), ucell.get_lnchiCounts(), PARAM.inp.nspin); + = CellIndex(unitcell::get_atomLabels(ucell.atoms, ucell.ntype), + unitcell::get_atomCounts(ucell.atoms, ucell.ntype), + unitcell::get_lnchiCounts(ucell.atoms, ucell.ntype), PARAM.inp.nspin); auto out_s_k = ModuleIO::Output_Sk(p_ham, pv, PARAM.inp.nspin, kv.get_nks()); auto out_dm_k = ModuleIO::Output_DMK(dm, pv, PARAM.inp.nspin, kv.get_nks()); @@ -61,7 +63,7 @@ void cal_mag(Parallel_Orbitals* pv, std::vector mag_x(ucell.nat, 0.0); std::vector mag_y(ucell.nat, 0.0); std::vector mag_z(ucell.nat, 0.0); - auto atomLabels = ucell.get_atomLabels(); + auto atomLabels = unitcell::get_atomLabels(ucell.atoms, ucell.ntype); if(PARAM.inp.nspin == 2) { diff --git a/source/source_io/module_output/cif_io.cpp b/source/source_io/module_output/cif_io.cpp index 4b9a0e0fea..540ea5c36e 100644 --- a/source/source_io/module_output/cif_io.cpp +++ b/source/source_io/module_output/cif_io.cpp @@ -307,7 +307,6 @@ void ModuleIO::CifParser::_unpack_ucell(const UnitCell& ucell, for (int i = 0; i < natom; ++i) { atom_site_labels[i] = ucell.atoms[ucell.iat2it[i]].ncpp.psd; // the most standard label - atom_site_labels[i] = atom_site_labels[i].empty() ? ucell.atom_label[ucell.iat2it[i]]: atom_site_labels[i]; atom_site_labels[i] = atom_site_labels[i].empty() ? ucell.atoms[ucell.iat2it[i]].label: atom_site_labels[i]; assert(!atom_site_labels[i].empty()); // ensure the label is not empty atom_site_fract_coords[3 * i] = ucell.atoms[ucell.iat2it[i]].taud[ucell.iat2ia[i]].x; diff --git a/source/source_io/test/for_testing_input_conv.h b/source/source_io/test/for_testing_input_conv.h index be67c4af01..5f34903207 100644 --- a/source/source_io/test/for_testing_input_conv.h +++ b/source/source_io/test/for_testing_input_conv.h @@ -149,9 +149,9 @@ void Occupy::decision(const std::string& name, const double& smearing_sigma) { return; } -// void UnitCell::setup(const std::string&,const int&,const int&,const +// void UnitCell::setup_from_input(const std::string&,const int&,const int&,const // bool&,const std::string&){return;} -void UnitCell::setup(const std::string& latname_in, +void UnitCell::setup_from_input(const std::string& latname_in, const int& ntype_in, const int& lmaxmax_in, const bool& init_vel_in, diff --git a/source/source_io/test/prepare_unitcell.h b/source/source_io/test/prepare_unitcell.h index 152a018a02..7164f317c2 100644 --- a/source/source_io/test/prepare_unitcell.h +++ b/source/source_io/test/prepare_unitcell.h @@ -71,13 +71,11 @@ class UcellTestPrepare //basic info this->ntype = this->elements.size(); UnitCell* ucell = new UnitCell; - ucell->setup(this->latname, + ucell->setup_from_input(this->latname, this->ntype, this->lmaxmax, this->init_vel, this->fixed_axes); - ucell->atom_label.resize(ucell->ntype); - ucell->atom_mass.resize(ucell->ntype); ucell->pseudo_fn.resize(ucell->ntype); ucell->pseudo_type.resize(ucell->ntype); @@ -87,8 +85,6 @@ class UcellTestPrepare ucell->magnet.ux_[2] = 0.0; for(int it=0;itntype;++it) { - ucell->atom_label[it] = this->elements[it]; - ucell->atom_mass[it] = this->atomic_mass[it]; ucell->pseudo_fn[it] = this->pp_files[it]; ucell->pseudo_type[it] = this->pp_types[it]; ucell->orbital_fn[it] = this->orb_files[it]; @@ -148,7 +144,7 @@ class UcellTestPrepare ucell->atoms[it].angle2.resize(ucell->atoms[it].na); ucell->atoms[it].m_loc_.resize(ucell->atoms[it].na); ucell->atoms[it].mbl.resize(ucell->atoms[it].na); - ucell->atoms[it].mass = ucell->atom_mass[it]; // mass set here + ucell->atoms[it].mass = this->atomic_mass[it]; for(int ia=0; iaatoms[it].na; ++ia) { if (ucell->Coordinate == "Direct") diff --git a/source/source_io/test_serial/prepare_unitcell.h b/source/source_io/test_serial/prepare_unitcell.h index 476fb8af7f..d07e358f2f 100644 --- a/source/source_io/test_serial/prepare_unitcell.h +++ b/source/source_io/test_serial/prepare_unitcell.h @@ -71,13 +71,11 @@ class UcellTestPrepare //basic info this->ntype = this->elements.size(); UnitCell* ucell = new UnitCell; - ucell->setup(this->latname, + ucell->setup_from_input(this->latname, this->ntype, this->lmaxmax, this->init_vel, this->fixed_axes); - ucell->atom_label.resize(ucell->ntype); - ucell->atom_mass.resize(ucell->ntype); ucell->pseudo_fn.resize(ucell->ntype); ucell->pseudo_type.resize(ucell->ntype); @@ -87,8 +85,6 @@ class UcellTestPrepare ucell->magnet.ux_[2] = 0.0; for(int it=0;itntype;++it) { - ucell->atom_label[it] = this->elements[it]; - ucell->atom_mass[it] = this->atomic_mass[it]; ucell->pseudo_fn[it] = this->pp_files[it]; ucell->pseudo_type[it] = this->pp_types[it]; ucell->orbital_fn[it] = this->orb_files[it]; @@ -148,7 +144,7 @@ class UcellTestPrepare ucell->atoms[it].angle2.resize(ucell->atoms[it].na); ucell->atoms[it].m_loc_.resize(ucell->atoms[it].na); ucell->atoms[it].mbl.resize(ucell->atoms[it].na); - ucell->atoms[it].mass = ucell->atom_mass[it]; // mass set here + ucell->atoms[it].mass = this->atomic_mass[it]; for(int ia=0; iaatoms[it].na; ++ia) { if (ucell->Coordinate == "Direct") diff --git a/source/source_lcao/module_deepks/deepks_basic.cpp b/source/source_lcao/module_deepks/deepks_basic.cpp index b9e6d63f4a..7eeb6c56ef 100644 --- a/source/source_lcao/module_deepks/deepks_basic.cpp +++ b/source/source_lcao/module_deepks/deepks_basic.cpp @@ -361,7 +361,7 @@ void DeePKS_domain::prepare_atom(const UnitCell& ucell, torch::Tensor& atom_out) { for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { - atom_out[index][0] = AtomInfo.atom_Z[ucell.atom_label[it]]; + atom_out[index][0] = AtomInfo.atom_Z[ucell.atoms[it].label]; // use bohr as unit atom_out[index][1] = ucell.atoms[it].tau[ia].x * ucell.lat0; diff --git a/source/source_lcao/module_deepks/test/CMakeLists.txt b/source/source_lcao/module_deepks/test/CMakeLists.txt index d35286d71f..beefdfecc5 100644 --- a/source/source_lcao/module_deepks/test/CMakeLists.txt +++ b/source/source_lcao/module_deepks/test/CMakeLists.txt @@ -21,6 +21,7 @@ set(DEEPKS_UNIT_COMMON_SOURCES ../../../source_cell/atom_pseudo.cpp ../../../source_cell/read_atoms.cpp ../../../source_cell/read_atoms_helper.cpp + ../../../source_cell/cell_tools.cpp ../../../source_cell/read_stru.cpp ../../../source_cell/print_cell.cpp ../../../source_cell/read_atom_species.cpp diff --git a/source/source_lcao/module_deltaspin/init_sc.cpp b/source/source_lcao/module_deltaspin/init_sc.cpp index 73da9388ec..e5636c1b76 100644 --- a/source/source_lcao/module_deltaspin/init_sc.cpp +++ b/source/source_lcao/module_deltaspin/init_sc.cpp @@ -1,4 +1,5 @@ #include "spin_constrain.h" +#include "source_cell/cell_tools.h" /** * @file init_sc.cpp @@ -68,9 +69,9 @@ void spinconstrain::SpinConstrain::init_sc(double sc_thr_in, // Step 4: Load target magnetic moments and initial lambda from UnitCell // These are parsed from the STRU file's "sc_mag" and "lambda" keywords - this->set_target_mag(ucell.get_target_mag()); - this->lambda_ = ucell.get_lambda(); - this->constrain_ = ucell.get_constrain(); + this->set_target_mag(unitcell::get_target_mag(ucell.atoms, ucell.ntype, ucell.nat)); + this->lambda_ = unitcell::get_lambda(ucell.atoms, ucell.ntype, ucell.nat); + this->constrain_ = unitcell::get_constrain(ucell.atoms, ucell.ntype, ucell.nat); // Step 5: CRITICAL FIX for collinear spin (nspin=2) // In collinear mode, spins are constrained along the z-axis only. @@ -89,7 +90,7 @@ void spinconstrain::SpinConstrain::init_sc(double sc_thr_in, } // Step 6: Set auxiliary parameters - this->atomLabels_ = ucell.get_atomLabels(); // "Fe_0", "Fe_1", etc. + this->atomLabels_ = unitcell::get_atomLabels(ucell.atoms, ucell.ntype); // "Fe_0", "Fe_1", etc. this->direction_only_ = direction_only_in; // Only optimize spin direction this->tpiba = ucell.tpiba; // 2*pi/a lattice scaling this->pw_wfc_ = pw_wfc_in; // PW basis (PW mode only) diff --git a/source/source_lcao/module_deltaspin/test/prepare_unitcell.h b/source/source_lcao/module_deltaspin/test/prepare_unitcell.h index 7abe206dbf..b0a2d0c85d 100644 --- a/source/source_lcao/module_deltaspin/test/prepare_unitcell.h +++ b/source/source_lcao/module_deltaspin/test/prepare_unitcell.h @@ -71,11 +71,9 @@ class UcellTestPrepare //basic info this->ntype = this->elements.size(); static UnitCell ucell; - ucell.setup(this->latname, this->ntype, this->lmaxmax, this->init_vel, this->fixed_axes); + ucell.setup_from_input(this->latname, this->ntype, this->lmaxmax, this->init_vel, this->fixed_axes); delete[] ucell.orbital_fn; delete[] ucell.magnet.start_magnetization; // mag set here - ucell->atom_label.resize(ucell->ntype); - ucell->atom_mass.resize(ucell->ntype); ucell->pseudo_fn.resize(ucell->ntype); ucell->pseudo_type.resize(ucell->ntype); @@ -86,8 +84,6 @@ class UcellTestPrepare ucell.magnet.ux_[2] = 0.0; for (int it = 0; it < ucell.ntype; ++it) { - ucell.atom_label[it] = this->elements[it]; - ucell.atom_mass[it] = this->atomic_mass[it]; ucell.pseudo_fn[it] = this->pp_files[it]; ucell.pseudo_type[it] = this->pp_types[it]; ucell.orbital_fn[it] = this->orb_files[it]; @@ -158,7 +154,7 @@ class UcellTestPrepare ucell.atoms[it].angle2 = new double[ucell.atoms[it].na]; ucell.atoms[it].m_loc_ = new ModuleBase::Vector3[ucell.atoms[it].na]; ucell.atoms[it].mbl = new ModuleBase::Vector3[ucell.atoms[it].na]; - ucell.atoms[it].mass = ucell.atom_mass[it]; // mass set here + ucell.atoms[it].mass = this->atomic_mass[it]; for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { if (ucell.Coordinate == "Cartesian") diff --git a/source/source_main/driver_run.cpp b/source/source_main/driver_run.cpp index c9faecf562..d1d945fc68 100644 --- a/source/source_main/driver_run.cpp +++ b/source/source_main/driver_run.cpp @@ -49,7 +49,7 @@ void Driver::driver_run() // the life of ucell should begin here, mohan 2024-05-12 UnitCell ucell; - ucell.setup(PARAM.inp.latname, + ucell.setup_from_input(PARAM.inp.latname, PARAM.inp.ntype, PARAM.inp.lmaxmax, PARAM.inp.init_vel, diff --git a/source/source_md/md_func.cpp b/source/source_md/md_func.cpp index 6bd1b60dd5..43e5590722 100644 --- a/source/source_md/md_func.cpp +++ b/source/source_md/md_func.cpp @@ -387,7 +387,7 @@ void dump_info(const int& step, { for (int ia = 0; ia < unit_in.atoms[it].na; ++ia) { - ofs << " " << index << " " << unit_in.atom_label[it] << " " << unit_in.atoms[it].tau[ia].x * unit_pos + ofs << " " << index << " " << unit_in.atoms[it].label << " " << unit_in.atoms[it].tau[ia].x * unit_pos << " " << unit_in.atoms[it].tau[ia].y * unit_pos << " " << unit_in.atoms[it].tau[ia].z * unit_pos; if (param_in.mdp.dump_force) diff --git a/source/source_md/test/CMakeLists.txt b/source/source_md/test/CMakeLists.txt index 31e456d2ea..9f74a519a7 100644 --- a/source/source_md/test/CMakeLists.txt +++ b/source/source_md/test/CMakeLists.txt @@ -14,6 +14,7 @@ list(APPEND depend_files ../../source_cell/atom_pseudo.cpp ../../source_cell/read_atoms.cpp ../../source_cell/read_atoms_helper.cpp + ../../source_cell/cell_tools.cpp ../../source_cell/pseudo.cpp ../../source_cell/read_pp.cpp ../../source_cell/read_pp_complete.cpp diff --git a/source/source_md/test/setcell.h b/source/source_md/test/setcell.h index b36151ee76..572e417ed7 100644 --- a/source/source_md/test/setcell.h +++ b/source/source_md/test/setcell.h @@ -28,10 +28,8 @@ class Setcell ucell.atoms = new Atom[ucell.ntype]; ucell.set_atom_flag = true; - ucell.atom_mass.resize(ucell.ntype); - ucell.atom_label.resize(ucell.ntype); - ucell.atom_mass[0] = 39.948; - ucell.atom_label[0] = "Ar"; + ucell.atoms[0].mass = 39.948; + ucell.atoms[0].label = "Ar"; ucell.lat0 = 1; ucell.lat0_angstrom = ucell.lat0 * ModuleBase::BOHR_TO_A; @@ -63,7 +61,6 @@ class Setcell ucell.atoms[0].taud.resize(4); ucell.atoms[0].vel.resize(4); ucell.atoms[0].mbl.resize(4); - ucell.atoms[0].mass = ucell.atom_mass[0]; ucell.atoms[0].angle1.resize(4); ucell.atoms[0].angle2.resize(4); @@ -90,11 +87,6 @@ class Setcell ucell.GGT = ucell.G * ucell.GT; ucell.invGGT = ucell.GGT.Inverse(); - ucell.GT0 = ucell.latvec.Inverse(); - ucell.G0 = ucell.GT.Transpose(); - ucell.GGT0 = ucell.G * ucell.GT; - ucell.invGGT0 = ucell.GGT.Inverse(); - ucell.set_iat2itia(); }; diff --git a/source/source_psi/test/psi_initializer_unit_test.cpp b/source/source_psi/test/psi_initializer_unit_test.cpp index 342e04a5a6..d2301597ca 100644 --- a/source/source_psi/test/psi_initializer_unit_test.cpp +++ b/source/source_psi/test/psi_initializer_unit_test.cpp @@ -138,9 +138,6 @@ class PsiIntializerUnitTest : public ::testing::Test { this->p_ucell->tpiba = 2.0 * M_PI / this->p_ucell->lat0; this->p_ucell->tpiba2 = this->p_ucell->tpiba * this->p_ucell->tpiba; // atom - this->p_ucell->atom_label.shrink_to_fit(); - this->p_ucell->atom_label.resize(1); - this->p_ucell->atom_label[0] = "Si"; // atom properties this->p_ucell->nat = 1; this->p_ucell->ntype = 1; diff --git a/source/source_pw/module_pwdft/onsite_proj.cpp b/source/source_pw/module_pwdft/onsite_proj.cpp index d4af4a7c2f..67458d7e7b 100644 --- a/source/source_pw/module_pwdft/onsite_proj.cpp +++ b/source/source_pw/module_pwdft/onsite_proj.cpp @@ -8,6 +8,7 @@ #include "source_pw/module_pwdft/onsite_proj_print.h" #include "source_lcao/module_dftu/dftu.h" #include "source_lcao/module_deltaspin/spin_constrain.h" +#include "source_cell/cell_tools.h" #include "source_io/module_parameter/parameter.h" #include "source_base/projgen.h" @@ -580,7 +581,7 @@ void projectors::OnsiteProjector::cal_occupations( Parallel_Reduce::reduce_double_allpool(npool, GlobalV::NPROC_IN_POOL, (double*)(&(occs[0])), occs.size()*2); // occ has been reduced and calculate mag // Print orbital charge analysis - auto atom_labels = this->ucell->get_atomLabels(); + auto atom_labels = unitcell::get_atomLabels(this->ucell->atoms, this->ucell->ntype); print::print_orb_chg(this->ucell, occs, this->iat_nh, atom_labels); // print charge diff --git a/source/source_pw/module_pwdft/test/CMakeLists.txt b/source/source_pw/module_pwdft/test/CMakeLists.txt index 08a5c9dd8e..123e0ccd61 100644 --- a/source/source_pw/module_pwdft/test/CMakeLists.txt +++ b/source/source_pw/module_pwdft/test/CMakeLists.txt @@ -44,6 +44,7 @@ AddTest( ../../../source_cell/read_atom_species.cpp ../../../source_cell/read_atoms.cpp ../../../source_cell/read_atoms_helper.cpp + ../../../source_cell/cell_tools.cpp ../../../source_cell/read_pp.cpp ../../../source_cell/read_pp_complete.cpp ../../../source_cell/read_pp_upf100.cpp diff --git a/source/source_relax/relax_nsync.cpp b/source/source_relax/relax_nsync.cpp index 71e342ef66..3e5f4bf114 100644 --- a/source/source_relax/relax_nsync.cpp +++ b/source/source_relax/relax_nsync.cpp @@ -2,6 +2,7 @@ #include "source_base/global_function.h" #include "source_base/global_variable.h" #include "source_io/module_parameter/parameter.h" +#include "source_cell/cell_tools.h" #include "source_cell/update_cell.h" /** @@ -81,8 +82,8 @@ bool IonCellOptimizer::relax_step(const int& istep, } // Determine what relaxation steps are needed - const bool need_atom_relax = (is_relax || is_cell_relax) && ucell.if_atoms_can_move(); - const bool need_cell_relax = is_cell_relax && ucell.if_cell_can_change(); + const bool need_atom_relax = (is_relax || is_cell_relax) && unitcell::if_atoms_can_move(ucell.atoms, ucell.ntype); + const bool need_cell_relax = is_cell_relax && unitcell::if_cell_can_change(ucell.lat_axis_free); // Atomic relaxation branch if (need_atom_relax) @@ -138,7 +139,7 @@ bool IonCellOptimizer::relax_step(const int& istep, return converged; } - else if (is_cell_relax && !ucell.if_cell_can_change()) + else if (is_cell_relax && !unitcell::if_cell_can_change(ucell.lat_axis_free)) { ModuleBase::WARNING("IonCellOptimizer", "Lattice vectors are not allowed to change!"); return true; diff --git a/source/source_relax/test/for_test.h b/source/source_relax/test/for_test.h index 3cdd6af573..ed16ca04cc 100644 --- a/source/source_relax/test/for_test.h +++ b/source/source_relax/test/for_test.h @@ -42,8 +42,6 @@ UnitCell::UnitCell() tpiba2 = 0.0; omega = 0.0; - atom_mass.shrink_to_fit(); - atom_label.resize(1); pseudo_fn.resize(1); pseudo_type.resize(1); orbital_fn.resize(1);