diff --git a/docs/CONTRIBUTING.md b/docs/CONTRIBUTING.md index 0d0dda27e7..caa29c6f24 100644 --- a/docs/CONTRIBUTING.md +++ b/docs/CONTRIBUTING.md @@ -301,6 +301,49 @@ To add a unit test: in the `abacus-develop` directory. +## Adding an integrate test +The integrate test is a test suite for testing the whole ABACUS package. The examples are located in the `tests/integrate` directory. Before adding a new test, please firstly read `README.md` in `tests/integrate` to understand the structure of the integrate test. To add an integrate test: +1. Add a new directory under `tests/integrate` for the new test. +2. Prepare the input files for the new test. + - The input files should be placed in the new directory. Pseudopotential files and orbital files should be placed in `tests/PP_ORB`. You should define the correct `pseudo_dir` and `orb_dir`(if need orbital files) in INPUT with the relative path to the `tests/PP_ORB` directory, and be sure the new test can be run successfully. + - The running time of the new test should not exceed 20 seconds. You can try to reduce the time by below methods (on the premise of ensuring the effectiveness of the test): + - Reduce the number of atoms in the unit cell (1~2 atoms). + - Reduce the number of k-points (`1 1 1` or `2 2 2`). + - Reduce ecutwfc (20~50 Ry). + - Reduce the number of steps for relax or md job (2~3 steps). + - Reduce the basis set for LCAO calculations (DZP orbital and 6 a.u. cutoff). + - For PW calculations, should set `pw_seed 1` in INPUT file to ensure the reproducibility of the test. +3. Generate the reference results for the new test. + - Run the new test with GNU compiler and 4 MPI processes 2 OpenMP threads. The command is `OMP_NUM_THREADS=2 mpirun -np 4 abacus > log.txt`. + - Execute tests/integrate/tools/catch_properties.sh script to generate the reference results. At the new test directory, run `bash ../tools/catch_properties.sh result.ref`. + A `result.ref` file may be like: + ```text + etotref -3439.007931317310 + etotperatomref -3439.0079313173 + totaltimeref 2.78 + ``` + - If you want to test the correctness of some output files, you need to do extra below steps: + 1. add the corresponding comparison method in `catch_properties.sh`. For example, to verify whether the output of the SPIN1_CHG.cube file is correct, you need to add the following code in `catch_properties.sh`: + ```bash + has_band=$(awk '$1=="out_band" {a=$2} END{print a}' INPUT) # check if the BAND is outputed + if ! test -z "$has_band" && [ $has_band == 1 ]; then # if band is outputed, then check if the band is correct + bandref=refBANDS_1.dat # this file should be prepared in new test directory + bandcal=OUT.autotest/BANDS_1.dat # this file is generated by each run of test + python3 ../tools/CompareFile.py $bandref $bandcal 8 # compare the new band file with the reference file + echo "CompareBand_pass $?" >>$1 # record the comparison result, $? is the return value of last command + fi + ``` + `CompareFile.py` is used to determine if two files are identical. It accepts three arguments: the first two are the files to be compared, and the third specifies the precision for comparing numerical values. The comparison fails if the difference between any two corresponding numerical values exceeds 1e-{precision} (such as: 1e-8 in previous case). If the files are identical, the script returns 0; otherwise, it returns 1. + + 2. Add the reference file (such as: `refBANDS_1.dat` in previous case) to the new test directory. + + 3. Add the reference comparison result to the `result.ref` file. For example, `CompareBand_pass 0` means the comparison of the band file is passed. (This statement should be added before the `totaltimeref` line) +4. Add a `jd` file in the new test directory, which is one setence to describe the new test. +5. Add the new test to `tests/integrate/CASES_CPU.txt` file (or `tests/integrate/CASES_GPU.txt` file if it is for GPU verion). +6. Enter directory tests/integrate and run `bash Autotest.sh -r ` to check if the new test can be run successfully. + + + ## Debugging the codes For the unexpected results when developing ABACUS, [GDB](https://www.sourceware.org/gdb/) will come in handy. diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 07c6c6811f..c58a651601 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -172,6 +172,7 @@ OBJS_CELL=atom_pseudo.o\ setup_nonlocal.o\ klist.o\ cell_index.o\ + check_atomic_stru.o\ OBJS_DEEPKS=LCAO_deepks.o\ LCAO_deepks_fdelta.o\ @@ -232,7 +233,6 @@ OBJS_GINT=gint.o\ gint_tau.o\ gint_vl.o\ gint_k_env.o\ - gint_k_sparse.o\ gint_k_sparse1.o\ gint_k_pvpr.o\ gint_k_pvdpr.o\ diff --git a/source/driver_run.cpp b/source/driver_run.cpp index 7a64fbe487..e379673e63 100644 --- a/source/driver_run.cpp +++ b/source/driver_run.cpp @@ -1,48 +1,53 @@ #include "driver.h" +#include "module_cell/check_atomic_stru.h" #include "module_cell/module_neighbor/sltk_atom_arrange.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/input.h" +#include "module_io/para_json.h" #include "module_io/print_info.h" #include "module_io/winput.h" #include "module_md/run_md.h" -#include "module_io/para_json.h" /** - * @brief This is the driver function which defines the workflow of ABACUS calculations. - * It relies on the class Esolver, which is a class that organizes workflows of single point calculations. - * - * For calculations involving change of configuration (lattice parameter & ionic motion), - * this driver calls Esolver::Run and the configuration-changing subroutine in a alternating manner. - * + * @brief This is the driver function which defines the workflow of ABACUS + * calculations. It relies on the class Esolver, which is a class that organizes + * workflows of single point calculations. + * + * For calculations involving change of configuration (lattice parameter & ionic + * motion), this driver calls Esolver::Run and the configuration-changing + * subroutine in a alternating manner. + * * Information is passed between the two subroutines by class UnitCell - * - * Esolver::Run takes in a configuration and provides force and stress, - * the configuration-changing subroutine takes force and stress and updates the configuration + * + * Esolver::Run takes in a configuration and provides force and stress, + * the configuration-changing subroutine takes force and stress and updates the + * configuration */ -void Driver::driver_run(void) -{ +void Driver::driver_run() { ModuleBase::TITLE("Driver", "driver_line"); ModuleBase::timer::tick("Driver", "driver_line"); - //! 1: initialize the ESolver - ModuleESolver::ESolver *p_esolver = nullptr; + //! 1: initialize the ESolver + ModuleESolver::ESolver* p_esolver = nullptr; ModuleESolver::init_esolver(p_esolver); //! 2: setup cell and atom information // this warning should not be here, mohan 2024-05-22 #ifndef __LCAO - if(GlobalV::BASIS_TYPE == "lcao_in_pw" || GlobalV::BASIS_TYPE == "lcao") - { - ModuleBase::WARNING_QUIT("driver","to use LCAO basis, compile with __LCAO"); + if (GlobalV::BASIS_TYPE == "lcao_in_pw" || GlobalV::BASIS_TYPE == "lcao") { + ModuleBase::WARNING_QUIT("driver", + "to use LCAO basis, compile with __LCAO"); } #endif // the life of ucell should begin here, mohan 2024-05-12 // delete ucell as a GlobalC in near future GlobalC::ucell.setup_cell(GlobalV::stru_file, GlobalV::ofs_running); + Check_Atomic_Stru::check_atomic_stru(GlobalC::ucell, + GlobalV::MIN_DIST_COEF); - //! 3: initialize Esolver and fill json-structure + //! 3: initialize Esolver and fill json-structure p_esolver->before_all_runners(INPUT, GlobalC::ucell); // this Json part should be moved to before_all_runners, mohan 2024-05-12 @@ -52,32 +57,25 @@ void Driver::driver_run(void) const std::string cal_type = GlobalV::CALCULATION; - //! 4: different types of calculations - if(cal_type == "md") - { + //! 4: different types of calculations + if (cal_type == "md") { Run_MD::md_line(GlobalC::ucell, p_esolver, INPUT.mdp); - } - else if(cal_type == "scf" - || cal_type == "relax" - || cal_type == "cell-relax") - { + } else if (cal_type == "scf" || cal_type == "relax" + || cal_type == "cell-relax") { Relax_Driver rl_driver; rl_driver.relax_driver(p_esolver); - } - else - { + } else { //! supported "other" functions: - //! nscf(PW,LCAO), - //! get_pchg(LCAO), - //! test_memory(PW,LCAO), + //! nscf(PW,LCAO), + //! get_pchg(LCAO), + //! test_memory(PW,LCAO), //! test_neighbour(LCAO), - //! get_S(LCAO), + //! get_S(LCAO), //! gen_bessel(PW), et al. const int istep = 0; p_esolver->others(istep); } - //! 5: clean up esolver p_esolver->after_all_runners(); ModuleESolver::clean_esolver(p_esolver); diff --git a/source/module_cell/CMakeLists.txt b/source/module_cell/CMakeLists.txt index e3d24bd05e..527177d8c5 100644 --- a/source/module_cell/CMakeLists.txt +++ b/source/module_cell/CMakeLists.txt @@ -22,6 +22,7 @@ add_library( klist.cpp parallel_kpoints.cpp cell_index.cpp + check_atomic_stru.cpp ) if(ENABLE_COVERAGE) diff --git a/source/module_cell/check_atomic_stru.cpp b/source/module_cell/check_atomic_stru.cpp new file mode 100644 index 0000000000..caa662097f --- /dev/null +++ b/source/module_cell/check_atomic_stru.cpp @@ -0,0 +1,159 @@ +#include "check_atomic_stru.h" + +#include "module_base/element_covalent_radius.h" + +void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, double& factor) { + // First we calculate all bond length in the structure, + // and compare with the covalent_bond_length, + // if there has bond length is shorter than covalent_bond_length * factor, + // we think this structure is unreasonable. + const double warning_coef = 0.6; + assert(ucell.ntype > 0); + std::stringstream errorlog; + bool all_pass = true; + bool no_warning = true; + for (int it1 = 0; it1 < ucell.ntype; it1++) { + std::string symbol1 = ""; + for (char ch: ucell.atoms[it1].label) { + if (std::isalpha(ch)) { + symbol1.push_back(ch); + } + } + // std::string symbol1 = ucell.atoms[it1].label; + double symbol1_covalent_radius; + if (ModuleBase::CovalentRadius.find(symbol1) + != ModuleBase::CovalentRadius.end()) { + symbol1_covalent_radius = ModuleBase::CovalentRadius.at(symbol1); + } else { + std::stringstream mess; + mess << "Notice: symbol '" << symbol1 + << "' is not an element symbol!!!! "; + mess << "set the covalent radius to be 0." << std::endl; + GlobalV::ofs_running << mess.str(); + std::cout << mess.str(); + symbol1_covalent_radius = 0.0; + } + + for (int ia1 = 0; ia1 < ucell.atoms[it1].na; ia1++) { + double x1 = ucell.atoms[it1].taud[ia1].x; + double y1 = ucell.atoms[it1].taud[ia1].y; + double z1 = ucell.atoms[it1].taud[ia1].z; + + for (int it2 = 0; it2 < ucell.ntype; it2++) { + std::string symbol2 = ucell.atoms[it2].label; + double symbol2_covalent_radius; + if (ModuleBase::CovalentRadius.find(symbol2) + != ModuleBase::CovalentRadius.end()) { + symbol2_covalent_radius + = ModuleBase::CovalentRadius.at(symbol2); + } else { + symbol2_covalent_radius = 0.0; + } + + double covalent_length + = (symbol1_covalent_radius + symbol2_covalent_radius) + / ModuleBase::BOHR_TO_A; + + for (int ia2 = 0; ia2 < ucell.atoms[it2].na; ia2++) { + for (int a = -1; a < 2; a++) { + for (int b = -1; b < 2; b++) { + for (int c = -1; c < 2; c++) { + if (it1 > it2) + continue; + else if (it1 == it2 && ia1 > ia2) + continue; + else if (it1 == it2 && ia1 == ia2 && a == 0 + && b == 0 && c == 0) + continue; + + double x2 = ucell.atoms[it2].taud[ia2].x + a; + double y2 = ucell.atoms[it2].taud[ia2].y + b; + double z2 = ucell.atoms[it2].taud[ia2].z + c; + + double bond_length + = sqrt(pow((x2 - x1) * ucell.a1.x + + (y2 - y1) * ucell.a2.x + + (z2 - z1) * ucell.a3.x, + 2) + + pow((x2 - x1) * ucell.a1.y + + (y2 - y1) * ucell.a2.y + + (z2 - z1) * ucell.a3.y, + 2) + + pow((x2 - x1) * ucell.a1.z + + (y2 - y1) * ucell.a2.z + + (z2 - z1) * ucell.a3.z, + 2)) + * ucell.lat0; + + if (bond_length < covalent_length * factor + || bond_length + < covalent_length * warning_coef) { + errorlog.setf(std::ios_base::fixed, + std::ios_base::floatfield); + errorlog << std::setw(3) << ia1 + 1 + << "-th " << std::setw(3) + << ucell.atoms[it1].label << ", "; + errorlog << std::setw(3) << ia2 + 1 + << "-th " << std::setw(3) + << ucell.atoms[it2].label; + errorlog << " (cell:" << std::setw(2) << a + << " " << std::setw(2) << b << " " + << std::setw(2) << c << ")"; + errorlog << ", distance= " + << std::setprecision(3) + << bond_length << " Bohr ("; + errorlog + << bond_length * ModuleBase::BOHR_TO_A + << " Angstrom)" << std::endl; + + if (bond_length + < covalent_length * factor) { + all_pass = false; + } else { + no_warning = false; + } + } + } // c + } // b + } // a + } // ia2 + } // it2 + } // ia1 + } // it1 + + if (!all_pass || !no_warning) { + std::stringstream mess; + mess << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" + << std::endl; + mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" + << std::endl; + mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" + << std::endl; + mess << "!!! WARNING: Some atoms are too close!!!" << std::endl; + mess << "!!! Please check the nearest-neighbor list in log file." + << std::endl; + mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" + << std::endl; + mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" + << std::endl; + mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" + << std::endl; + + GlobalV::ofs_running << mess.str() << mess.str() << mess.str() + << errorlog.str(); + std::cout << mess.str() << mess.str() << mess.str() << std::endl; + + if (!all_pass) { + mess.clear(); + mess.str(""); + mess << "If this structure is what you want, you can set " + "'min_dist_coef'" + << std::endl; + mess << "as a smaller value (the current value is " << factor + << ") in INPUT file." << std::endl; + GlobalV::ofs_running << mess.str(); + std::cout << mess.str(); + ModuleBase::WARNING_QUIT("Input", "The structure is unreasonable!"); + } + } +} diff --git a/source/module_cell/check_atomic_stru.h b/source/module_cell/check_atomic_stru.h new file mode 100644 index 0000000000..832b5d54bc --- /dev/null +++ b/source/module_cell/check_atomic_stru.h @@ -0,0 +1,11 @@ +#ifndef CHECK_ATOMIC_STRU_H +#define CHECK_ATOMIC_STRU_H + +#include "unitcell.h" + +class Check_Atomic_Stru { + public: + static void check_atomic_stru(UnitCell& ucell, double& factor); +}; + +#endif \ No newline at end of file diff --git a/source/module_cell/test/CMakeLists.txt b/source/module_cell/test/CMakeLists.txt index 0a0be0f2e9..b84b473862 100644 --- a/source/module_cell/test/CMakeLists.txt +++ b/source/module_cell/test/CMakeLists.txt @@ -24,6 +24,7 @@ list(APPEND cell_simple_srcs ../read_pp_upf100.cpp ../read_pp_vwr.cpp ../read_pp_blps.cpp + ../check_atomic_stru.cpp ) add_library(cell_info OBJECT ${cell_simple_srcs}) diff --git a/source/module_cell/test/support/mock_unitcell.cpp b/source/module_cell/test/support/mock_unitcell.cpp index dd91c3cefb..1da3c1aad4 100644 --- a/source/module_cell/test/support/mock_unitcell.cpp +++ b/source/module_cell/test/support/mock_unitcell.cpp @@ -1,14 +1,16 @@ #include "module_cell/unitcell.h" /* README: - This file supports idea like "I dont need any functions of UnitCell, I want to avoid using UnitCell functions because there is GLobalC, which will bring - endless compile troubles like undefined behavior" + This file supports idea like "I dont need any functions of UnitCell, I want + to avoid using UnitCell functions because there is GLobalC, which will bring + endless compile troubles like undefined behavior" */ void UnitCell::cal_ux() {} -bool UnitCell::judge_parallel(double a[3], ModuleBase::Vector3 b) {return true;} +bool UnitCell::judge_parallel(double a[3], ModuleBase::Vector3 b) { + return true; +} void UnitCell::set_iat2iwt(const int& npol_in) {} -UnitCell::UnitCell() -{ +UnitCell::UnitCell() { if (GlobalV::test_unitcell) ModuleBase::TITLE("unitcell", "Constructor"); Coordinate = "Direct"; @@ -48,8 +50,7 @@ UnitCell::UnitCell() set_atom_flag = false; } -UnitCell::~UnitCell() -{ +UnitCell::~UnitCell() { delete[] atom_label; delete[] atom_mass; delete[] pseudo_fn; @@ -60,37 +61,47 @@ UnitCell::~UnitCell() delete[] iwt2iat; delete[] iwt2iw; delete[] lc; - if(set_atom_flag) - { - delete[] atoms; - } + if (set_atom_flag) { + delete[] atoms; + } } void UnitCell::print_cell(std::ofstream& ofs) const {} void UnitCell::print_cell_xyz(const std::string& fn) const {} void UnitCell::print_cell_cif(const std::string& fn) const {} -int UnitCell::read_atom_species(std::ifstream &ifa, std::ofstream &ofs_running) {return 0;} -void UnitCell::read_cell_pseudopots(const std::string &pp_dir, std::ofstream &log) {} -bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_running, std::ofstream &ofs_warning) {return true;} +int UnitCell::read_atom_species(std::ifstream& ifa, + std::ofstream& ofs_running) { + return 0; +} +void UnitCell::read_cell_pseudopots(const std::string& pp_dir, + std::ofstream& log) {} +bool UnitCell::read_atom_positions(std::ifstream& ifpos, + std::ofstream& ofs_running, + std::ofstream& ofs_warning) { + return true; +} void UnitCell::update_pos_tau(const double* pos) {} void UnitCell::update_pos_taud(double* posd_in) {} void UnitCell::update_pos_taud(const ModuleBase::Vector3* posd_in) {} void UnitCell::update_vel(const ModuleBase::Vector3* vel_in) {} void UnitCell::periodic_boundary_adjustment() {} void UnitCell::bcast_atoms_tau() {} -bool UnitCell::judge_big_cell() const {return true;} -void UnitCell::update_stress(ModuleBase::matrix &scs) {} -void UnitCell::update_force(ModuleBase::matrix &fcs) {} +bool UnitCell::judge_big_cell() const { return true; } +void UnitCell::update_stress(ModuleBase::matrix& scs) {} +void UnitCell::update_force(ModuleBase::matrix& fcs) {} #ifdef __MPI void UnitCell::bcast_unitcell() {} void UnitCell::bcast_unitcell2() {} #endif -void UnitCell::set_iat2itia(void) {} -void UnitCell::setup_cell(const std::string &fn, std::ofstream &log) {} -void UnitCell::read_orb_file(int it, std::string &orb_file, std::ofstream &ofs_running, Atom *atom) {} -void UnitCell::read_pseudo(std::ofstream &ofs) {} -int UnitCell::find_type(const std::string &label) {return 0;} -void UnitCell::print_tau(void) const {} -void UnitCell::print_stru_file(const std::string& fn, +void UnitCell::set_iat2itia() {} +void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) {} +void UnitCell::read_orb_file(int it, + std::string& orb_file, + std::ofstream& ofs_running, + Atom* atom) {} +void UnitCell::read_pseudo(std::ofstream& ofs) {} +int UnitCell::find_type(const std::string& label) { return 0; } +void UnitCell::print_tau() const {} +void UnitCell::print_stru_file(const std::string& fn, const int& nspin, const bool& direct, const bool& vel, @@ -98,17 +109,20 @@ void UnitCell::print_stru_file(const std::string& fn, const bool& orb, const bool& dpks_desc, const int& iproc) const {} -void UnitCell::check_dtau(void) {} -void UnitCell::setup_cell_after_vc(std::ofstream &log) {} +void UnitCell::check_dtau() {} +void UnitCell::setup_cell_after_vc(std::ofstream& log) {} void UnitCell::remake_cell() {} -void UnitCell::cal_nwfc(std::ofstream &log) {} +void UnitCell::cal_nwfc(std::ofstream& log) {} void UnitCell::cal_meshx() {} -void UnitCell::cal_natomwfc(std::ofstream &log) {} -void UnitCell::print_unitcell_pseudo(const std::string &fn) {} -bool UnitCell::check_tau(void)const {return true;} -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, const int &ntype_in, const int &lmaxmax_in, const bool &init_vel_in, const std::string &fixed_axes_in) {} -void UnitCell::check_structure(double factor) {} +void UnitCell::cal_natomwfc(std::ofstream& log) {} +void UnitCell::print_unitcell_pseudo(const std::string& fn) {} +bool UnitCell::check_tau() const { return true; } +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, + const int& ntype_in, + const int& lmaxmax_in, + const bool& init_vel_in, + const std::string& fixed_axes_in) {} void UnitCell::cal_nelec(double& nelec) {} void UnitCell::compare_atom_labels(std::string label1, std::string label2) {} \ No newline at end of file diff --git a/source/module_cell/test/unitcell_test_readpp.cpp b/source/module_cell/test/unitcell_test_readpp.cpp index 8ece29e72e..167fb8e5b7 100644 --- a/source/module_cell/test/unitcell_test_readpp.cpp +++ b/source/module_cell/test/unitcell_test_readpp.cpp @@ -1,30 +1,28 @@ -#include "gtest/gtest.h" -#include "gmock/gmock.h" #include "memory" -#include "module_base/mathzone.h" #include "module_base/global_variable.h" +#include "module_base/mathzone.h" +#include "module_cell/check_atomic_stru.h" #include "module_cell/unitcell.h" -#include + +#include "gmock/gmock.h" +#include "gtest/gtest.h" #include +#include #ifdef __MPI #include "mpi.h" #endif #include "prepare_unitcell.h" #ifdef __LCAO -InfoNonlocal::InfoNonlocal(){} -InfoNonlocal::~InfoNonlocal(){} +InfoNonlocal::InfoNonlocal() {} +InfoNonlocal::~InfoNonlocal() {} #endif -Magnetism::Magnetism() -{ - this->tot_magnetization = 0.0; - this->abs_magnetization = 0.0; - this->start_magnetization = nullptr; -} -Magnetism::~Magnetism() -{ - delete[] this->start_magnetization; +Magnetism::Magnetism() { + this->tot_magnetization = 0.0; + this->abs_magnetization = 0.0; + this->start_magnetization = nullptr; } +Magnetism::~Magnetism() { delete[] this->start_magnetization; } /************************************************ * unit test of class UnitCell @@ -33,14 +31,16 @@ Magnetism::~Magnetism() /** * - Tested Functions: * - ReadCellPPWarning1 - * - read_cell_pseudopots(): error when average the pseudopotential: error_ap + * - read_cell_pseudopots(): error when average the pseudopotential: + * error_ap * - ReadCellPPWarning2 * - read_cell_pseudopots(): Couldn't find pseudopotential file:: error == 1 * - ReadCellPPWarning3 * - read_cell_pseudopots(): Pseudopotential data do not match: error ==2 * - error==3 is currently difficult to reach in read_pseudo_vwr * - ReadCellPPWarning4 - * - read_cell_pseudopots(): dft_functional from INPUT does not match that in pseudopot file + * - read_cell_pseudopots(): dft_functional from INPUT does not match that + * in pseudopot file * - ReadCellPPWarning5 * - read_cell_pseudopots(): Unknown pseudopotential type * - ReadCellPP @@ -48,19 +48,24 @@ Magnetism::~Magnetism() * - CalMeshx * - cal_meshx(): calculate max mesh info from atomic pseudo potential file * - CalNatomwfc1 - * - cal_natomwfc(): calculate total number of atomic orbitals in pseudo potential file + * - cal_natomwfc(): calculate total number of atomic orbitals in pseudo + * potential file * - NSPIN != 4 - * - this corresponds to number_of_wfc, PP_CHI in pp file, and atoms[it].ncpp.lchi[ncpp.nchi] + * - this corresponds to number_of_wfc, PP_CHI in pp file, and + * atoms[it].ncpp.lchi[ncpp.nchi] * - setup the total number of PAOs: pseudopotential atomic orbitals * - CalNatomwfc2 - * - cal_natomwfc(): calculate total number of atomic orbitals in pseudo potential file + * - cal_natomwfc(): calculate total number of atomic orbitals in pseudo + * potential file * - NSPIN ==4, has_so = false * - CalNatomwfc3 - * - cal_natomwfc(): calculate total number of atomic orbitals in pseudo potential file + * - cal_natomwfc(): calculate total number of atomic orbitals in pseudo + * potential file * - NSPIN ==4, has_so = true * - CalNwfc1 * - cal_nwfc(): calcuate the total number of local basis: NSPIN != 4 - * - this corresponds to number_of_proj, PP_BETA in pp file, and atoms[it].l_nchi[nw], nw from orb file + * - this corresponds to number_of_proj, PP_BETA in pp file, and + * atoms[it].l_nchi[nw], nw from orb file * - setup GlobalV::NLOCAL * - interfaces initialed in this function: * - itia2iat @@ -71,327 +76,339 @@ Magnetism::~Magnetism() * - CalNwfc2 * - cal_nwfc(): calcuate the total number of local basis: NSPIN == 4 * - CheckStructure - * - check_structure(): check if too atoms are two close + * - check_atomic_stru(): check if too atoms are two close * - ReadPseudoWarning1 * - read_pseudo(): All DFT functional must consistent. * - ReadPseudoWarning2 - * - read_pseudo(): number valence electrons > corresponding minimum possible of an element + * - read_pseudo(): number valence electrons > corresponding minimum + * possible of an element * - CalNelec: UnitCell::cal_nelec * - calculate the total number of valence electrons from psp files */ -//mock function +// mock function #ifdef __LCAO -void LCAO_Orbitals::bcast_files( - const int &ntype_in, - const int &my_rank) -{ - return; +void LCAO_Orbitals::bcast_files(const int& ntype_in, const int& my_rank) { + return; } #endif -class UcellTest : public ::testing::Test -{ -protected: - UcellTestPrepare utp = UcellTestLib["C1H2-Read"]; - std::unique_ptr ucell; - std::ofstream ofs; - std::string pp_dir; - std::string output; - void SetUp() - { - ofs.open("running.log"); - GlobalV::relax_new = utp.relax_new; - GlobalV::global_out_dir = "./"; - ucell = utp.SetUcellInfo(); - GlobalV::LSPINORB = false; - pp_dir = "./support/"; - GlobalV::PSEUDORCUT = 15.0; - GlobalV::DFT_FUNCTIONAL = "default"; - GlobalV::test_unitcell = 1; - GlobalV::test_pseudo_cell = 1; - GlobalV::NSPIN = 1; - GlobalV::BASIS_TYPE = "pw"; - } - void TearDown() - { - ofs.close(); - } +class UcellTest : public ::testing::Test { + protected: + UcellTestPrepare utp = UcellTestLib["C1H2-Read"]; + std::unique_ptr ucell; + std::ofstream ofs; + std::string pp_dir; + std::string output; + void SetUp() { + ofs.open("running.log"); + GlobalV::relax_new = utp.relax_new; + GlobalV::global_out_dir = "./"; + ucell = utp.SetUcellInfo(); + GlobalV::LSPINORB = false; + pp_dir = "./support/"; + GlobalV::PSEUDORCUT = 15.0; + GlobalV::DFT_FUNCTIONAL = "default"; + GlobalV::test_unitcell = 1; + GlobalV::test_pseudo_cell = 1; + GlobalV::NSPIN = 1; + GlobalV::BASIS_TYPE = "pw"; + } + void TearDown() { ofs.close(); } }; using UcellDeathTest = UcellTest; -TEST_F(UcellDeathTest,ReadCellPPWarning1) -{ - GlobalV::LSPINORB = true; - ucell->pseudo_fn[1] = "H_sr.upf"; - testing::internal::CaptureStdout(); - EXPECT_EXIT(ucell->read_cell_pseudopots(pp_dir,ofs), - ::testing::ExitedWithCode(0),""); - output = testing::internal::GetCapturedStdout(); - EXPECT_THAT(output,testing::HasSubstr("error when average the pseudopotential.")); +TEST_F(UcellDeathTest, ReadCellPPWarning1) { + GlobalV::LSPINORB = true; + ucell->pseudo_fn[1] = "H_sr.upf"; + testing::internal::CaptureStdout(); + EXPECT_EXIT(ucell->read_cell_pseudopots(pp_dir, ofs), + ::testing::ExitedWithCode(0), + ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, + testing::HasSubstr("error when average the pseudopotential.")); } -TEST_F(UcellDeathTest,ReadCellPPWarning2) -{ - pp_dir = "./arbitrary/"; - testing::internal::CaptureStdout(); - EXPECT_EXIT(ucell->read_cell_pseudopots(pp_dir,ofs), - ::testing::ExitedWithCode(0),""); - output = testing::internal::GetCapturedStdout(); - EXPECT_THAT(output,testing::HasSubstr("Couldn't find pseudopotential file")); +TEST_F(UcellDeathTest, ReadCellPPWarning2) { + pp_dir = "./arbitrary/"; + testing::internal::CaptureStdout(); + EXPECT_EXIT(ucell->read_cell_pseudopots(pp_dir, ofs), + ::testing::ExitedWithCode(0), + ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, + testing::HasSubstr("Couldn't find pseudopotential file")); } -TEST_F(UcellDeathTest,ReadCellPPWarning3) -{ - ucell->pseudo_type[0] = "upf"; - testing::internal::CaptureStdout(); - EXPECT_EXIT(ucell->read_cell_pseudopots(pp_dir,ofs), - ::testing::ExitedWithCode(0),""); - output = testing::internal::GetCapturedStdout(); - EXPECT_THAT(output,testing::HasSubstr("Pseudopotential data do not match.")); +TEST_F(UcellDeathTest, ReadCellPPWarning3) { + ucell->pseudo_type[0] = "upf"; + testing::internal::CaptureStdout(); + EXPECT_EXIT(ucell->read_cell_pseudopots(pp_dir, ofs), + ::testing::ExitedWithCode(0), + ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, + testing::HasSubstr("Pseudopotential data do not match.")); } -TEST_F(UcellDeathTest,ReadCellPPWarning4) -{ - GlobalV::DFT_FUNCTIONAL = "LDA"; - testing::internal::CaptureStdout(); - EXPECT_NO_THROW(ucell->read_cell_pseudopots(pp_dir,ofs)); - output = testing::internal::GetCapturedStdout(); +TEST_F(UcellDeathTest, ReadCellPPWarning4) { + GlobalV::DFT_FUNCTIONAL = "LDA"; + testing::internal::CaptureStdout(); + EXPECT_NO_THROW(ucell->read_cell_pseudopots(pp_dir, ofs)); + output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("dft_functional readin is: LDA")); - EXPECT_THAT(output, testing::HasSubstr("dft_functional in pseudopot file is: PBE")); - EXPECT_THAT(output, testing::HasSubstr("Please make sure this is what you need")); + EXPECT_THAT(output, + testing::HasSubstr("dft_functional in pseudopot file is: PBE")); + EXPECT_THAT(output, + testing::HasSubstr("Please make sure this is what you need")); } -TEST_F(UcellDeathTest, ReadCellPPWarning5) -{ +TEST_F(UcellDeathTest, ReadCellPPWarning5) { ucell->pseudo_type[0] = "upf0000"; testing::internal::CaptureStdout(); - EXPECT_EXIT(ucell->read_cell_pseudopots(pp_dir, ofs), ::testing::ExitedWithCode(0), ""); + EXPECT_EXIT(ucell->read_cell_pseudopots(pp_dir, ofs), + ::testing::ExitedWithCode(0), + ""); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("Unknown pseudopotential type.")); } -TEST_F(UcellTest,ReadCellPP) -{ - ucell->atoms[1].flag_empty_element = true; - ucell->read_cell_pseudopots(pp_dir,ofs); - EXPECT_EQ(ucell->atoms[0].ncpp.pp_type,"NC"); - EXPECT_FALSE(ucell->atoms[0].ncpp.has_so); //becomes false in average_p - EXPECT_FALSE(ucell->atoms[1].ncpp.has_so); - EXPECT_EQ(ucell->atoms[0].ncpp.nchi,2); //3=>2 in average_p - EXPECT_EQ(ucell->atoms[1].ncpp.nchi,1); - ofs.close(); - std::ifstream ifs; - ifs.open("running.log"); - std::string str((std::istreambuf_iterator(ifs)),std::istreambuf_iterator()); - EXPECT_THAT(str, testing::HasSubstr("Read in pseudopotential file is C.upf")); - EXPECT_THAT(str, testing::HasSubstr("pseudopotential type = NC")); - EXPECT_THAT(str, testing::HasSubstr("exchange-correlation functional = PBE")); - EXPECT_THAT(str, testing::HasSubstr("valence electrons = 4")); - EXPECT_THAT(str, testing::HasSubstr("Read in pseudopotential file is H.upf")); - EXPECT_THAT(str, testing::HasSubstr("valence electrons = 0")); +TEST_F(UcellTest, ReadCellPP) { + ucell->atoms[1].flag_empty_element = true; + ucell->read_cell_pseudopots(pp_dir, ofs); + EXPECT_EQ(ucell->atoms[0].ncpp.pp_type, "NC"); + EXPECT_FALSE(ucell->atoms[0].ncpp.has_so); // becomes false in average_p + EXPECT_FALSE(ucell->atoms[1].ncpp.has_so); + EXPECT_EQ(ucell->atoms[0].ncpp.nchi, 2); // 3=>2 in average_p + EXPECT_EQ(ucell->atoms[1].ncpp.nchi, 1); + ofs.close(); + std::ifstream ifs; + ifs.open("running.log"); + std::string str((std::istreambuf_iterator(ifs)), + std::istreambuf_iterator()); + EXPECT_THAT(str, + testing::HasSubstr("Read in pseudopotential file is C.upf")); + EXPECT_THAT(str, testing::HasSubstr("pseudopotential type = NC")); + EXPECT_THAT(str, + testing::HasSubstr("exchange-correlation functional = PBE")); + EXPECT_THAT(str, testing::HasSubstr("valence electrons = 4")); + EXPECT_THAT(str, + testing::HasSubstr("Read in pseudopotential file is H.upf")); + EXPECT_THAT(str, testing::HasSubstr("valence electrons = 0")); } -TEST_F(UcellTest,CalMeshx) -{ - ucell->read_cell_pseudopots(pp_dir,ofs); - ucell->cal_meshx(); - EXPECT_EQ(ucell->atoms[0].ncpp.msh,1247); - EXPECT_EQ(ucell->atoms[1].ncpp.msh,1165); - EXPECT_EQ(ucell->meshx,1247); +TEST_F(UcellTest, CalMeshx) { + ucell->read_cell_pseudopots(pp_dir, ofs); + ucell->cal_meshx(); + EXPECT_EQ(ucell->atoms[0].ncpp.msh, 1247); + EXPECT_EQ(ucell->atoms[1].ncpp.msh, 1165); + EXPECT_EQ(ucell->meshx, 1247); } -TEST_F(UcellTest,CalNatomwfc1) -{ - ucell->read_cell_pseudopots(pp_dir,ofs); - EXPECT_FALSE(ucell->atoms[0].ncpp.has_so); - EXPECT_FALSE(ucell->atoms[1].ncpp.has_so); - ucell->cal_natomwfc(ofs); - EXPECT_EQ(ucell->atoms[0].ncpp.nchi,2); - EXPECT_EQ(ucell->atoms[1].ncpp.nchi,1); - EXPECT_EQ(ucell->atoms[0].na,1); - EXPECT_EQ(ucell->atoms[1].na,2); - EXPECT_EQ(ucell->natomwfc,(1+3)*1+1*2); +TEST_F(UcellTest, CalNatomwfc1) { + ucell->read_cell_pseudopots(pp_dir, ofs); + EXPECT_FALSE(ucell->atoms[0].ncpp.has_so); + EXPECT_FALSE(ucell->atoms[1].ncpp.has_so); + ucell->cal_natomwfc(ofs); + EXPECT_EQ(ucell->atoms[0].ncpp.nchi, 2); + EXPECT_EQ(ucell->atoms[1].ncpp.nchi, 1); + EXPECT_EQ(ucell->atoms[0].na, 1); + EXPECT_EQ(ucell->atoms[1].na, 2); + EXPECT_EQ(ucell->natomwfc, (1 + 3) * 1 + 1 * 2); } -TEST_F(UcellTest,CalNatomwfc2) -{ - GlobalV::LSPINORB = false; - GlobalV::NSPIN = 4; - ucell->read_cell_pseudopots(pp_dir,ofs); - EXPECT_FALSE(ucell->atoms[0].ncpp.has_so); - EXPECT_FALSE(ucell->atoms[1].ncpp.has_so); - ucell->cal_natomwfc(ofs); - EXPECT_EQ(ucell->atoms[0].ncpp.nchi,2); - EXPECT_EQ(ucell->atoms[1].ncpp.nchi,1); - EXPECT_EQ(ucell->atoms[0].na,1); - EXPECT_EQ(ucell->atoms[1].na,2); - EXPECT_EQ(ucell->natomwfc,((1+3)*1+1*2)*2); +TEST_F(UcellTest, CalNatomwfc2) { + GlobalV::LSPINORB = false; + GlobalV::NSPIN = 4; + ucell->read_cell_pseudopots(pp_dir, ofs); + EXPECT_FALSE(ucell->atoms[0].ncpp.has_so); + EXPECT_FALSE(ucell->atoms[1].ncpp.has_so); + ucell->cal_natomwfc(ofs); + EXPECT_EQ(ucell->atoms[0].ncpp.nchi, 2); + EXPECT_EQ(ucell->atoms[1].ncpp.nchi, 1); + EXPECT_EQ(ucell->atoms[0].na, 1); + EXPECT_EQ(ucell->atoms[1].na, 2); + EXPECT_EQ(ucell->natomwfc, ((1 + 3) * 1 + 1 * 2) * 2); } -TEST_F(UcellTest,CalNatomwfc3) -{ - GlobalV::LSPINORB = true; - GlobalV::NSPIN = 4; - ucell->read_cell_pseudopots(pp_dir,ofs); - EXPECT_TRUE(ucell->atoms[0].ncpp.has_so); - EXPECT_TRUE(ucell->atoms[1].ncpp.has_so); - ucell->cal_natomwfc(ofs); - EXPECT_EQ(ucell->atoms[0].ncpp.nchi,3); - EXPECT_EQ(ucell->atoms[1].ncpp.nchi,1); - EXPECT_EQ(ucell->atoms[0].na,1); - EXPECT_EQ(ucell->atoms[1].na,2); - EXPECT_EQ(ucell->natomwfc,((2*0+2)+(2*1+2)+(2*1))*1+(2*0+2)*2); +TEST_F(UcellTest, CalNatomwfc3) { + GlobalV::LSPINORB = true; + GlobalV::NSPIN = 4; + ucell->read_cell_pseudopots(pp_dir, ofs); + EXPECT_TRUE(ucell->atoms[0].ncpp.has_so); + EXPECT_TRUE(ucell->atoms[1].ncpp.has_so); + ucell->cal_natomwfc(ofs); + EXPECT_EQ(ucell->atoms[0].ncpp.nchi, 3); + EXPECT_EQ(ucell->atoms[1].ncpp.nchi, 1); + EXPECT_EQ(ucell->atoms[0].na, 1); + EXPECT_EQ(ucell->atoms[1].na, 2); + EXPECT_EQ(ucell->natomwfc, + ((2 * 0 + 2) + (2 * 1 + 2) + (2 * 1)) * 1 + (2 * 0 + 2) * 2); } -TEST_F(UcellTest,CalNwfc1) -{ - ucell->read_cell_pseudopots(pp_dir,ofs); - EXPECT_FALSE(ucell->atoms[0].ncpp.has_so); - EXPECT_FALSE(ucell->atoms[1].ncpp.has_so); - ucell->cal_nwfc(ofs); - EXPECT_EQ(ucell->atoms[0].iw2l[8],2); - EXPECT_EQ(ucell->atoms[0].iw2n[8],0); - EXPECT_EQ(ucell->atoms[0].iw2m[8],4); - EXPECT_EQ(ucell->atoms[1].iw2l[8],2); - EXPECT_EQ(ucell->atoms[1].iw2n[8],0); - EXPECT_EQ(ucell->atoms[1].iw2m[8],4); - EXPECT_EQ(ucell->atoms[1].iw2_ylm[8],8); - //here is the default table for pw basis calculation - // nw = 1*1 + 3*1 + 5*1 = 9 - // L N m L*L+m - // 0 0 0 0 0 - // 1 1 0 0 0 - // 2 1 0 1 2 - // 3 1 0 2 2 - // 4 2 0 0 4 - // 5 2 0 1 5 - // 6 2 0 2 6 - // 7 2 0 3 7 - // 8 2 0 4 8 - EXPECT_EQ(ucell->atoms[0].na,1); - EXPECT_EQ(ucell->atoms[1].na,2); - EXPECT_EQ(ucell->namax,2); - EXPECT_EQ(ucell->atoms[0].nw,9); - EXPECT_EQ(ucell->atoms[1].nw,9); - EXPECT_EQ(ucell->nwmax,9); - EXPECT_EQ(GlobalV::NLOCAL,3*9); - //check itia2iat - EXPECT_EQ(ucell->itia2iat.getSize(), 4); - EXPECT_EQ(ucell->itia2iat(0,0), 0); - EXPECT_EQ(ucell->itia2iat(0,1), 0); - EXPECT_EQ(ucell->itia2iat(1,0), 1); - EXPECT_EQ(ucell->itia2iat(1,1), 2); - //check iat2iwt - EXPECT_EQ(ucell->get_npol(), 1); - EXPECT_EQ(ucell->get_iat2iwt()[0], 0); - EXPECT_EQ(ucell->get_iat2iwt()[1], 9); - EXPECT_EQ(ucell->get_iat2iwt()[2], 18); - //check itiaiw2iwt - EXPECT_EQ(ucell->itiaiw2iwt(0, 0, 0), 0); - EXPECT_EQ(ucell->itiaiw2iwt(0, 0, 1), 1); - EXPECT_EQ(ucell->itiaiw2iwt(0, 0, 8), 8); - EXPECT_EQ(ucell->itiaiw2iwt(1, 0, 0), 9); - EXPECT_EQ(ucell->itiaiw2iwt(1, 1, 0), 18); - //check itia2iat - EXPECT_EQ(ucell->itia2iat.getSize(), 4); - EXPECT_EQ(ucell->itia2iat(0, 0), 0); - EXPECT_EQ(ucell->itia2iat(0, 1), 0); - EXPECT_EQ(ucell->itia2iat(1, 0), 1); - EXPECT_EQ(ucell->itia2iat(1, 1), 2); - //check iwt2iat - EXPECT_EQ(ucell->iwt2iat[0], 0); - EXPECT_EQ(ucell->iwt2iat[10], 1); - EXPECT_EQ(ucell->iwt2iat[20], 2); - //check iwt2iw - EXPECT_EQ(ucell->iwt2iw[0], 0); - EXPECT_EQ(ucell->iwt2iw[10], 1); - EXPECT_EQ(ucell->iwt2iw[20], 2); +TEST_F(UcellTest, CalNwfc1) { + ucell->read_cell_pseudopots(pp_dir, ofs); + EXPECT_FALSE(ucell->atoms[0].ncpp.has_so); + EXPECT_FALSE(ucell->atoms[1].ncpp.has_so); + ucell->cal_nwfc(ofs); + EXPECT_EQ(ucell->atoms[0].iw2l[8], 2); + EXPECT_EQ(ucell->atoms[0].iw2n[8], 0); + EXPECT_EQ(ucell->atoms[0].iw2m[8], 4); + EXPECT_EQ(ucell->atoms[1].iw2l[8], 2); + EXPECT_EQ(ucell->atoms[1].iw2n[8], 0); + EXPECT_EQ(ucell->atoms[1].iw2m[8], 4); + EXPECT_EQ(ucell->atoms[1].iw2_ylm[8], 8); + // here is the default table for pw basis calculation + // nw = 1*1 + 3*1 + 5*1 = 9 + // L N m L*L+m + // 0 0 0 0 0 + // 1 1 0 0 0 + // 2 1 0 1 2 + // 3 1 0 2 2 + // 4 2 0 0 4 + // 5 2 0 1 5 + // 6 2 0 2 6 + // 7 2 0 3 7 + // 8 2 0 4 8 + EXPECT_EQ(ucell->atoms[0].na, 1); + EXPECT_EQ(ucell->atoms[1].na, 2); + EXPECT_EQ(ucell->namax, 2); + EXPECT_EQ(ucell->atoms[0].nw, 9); + EXPECT_EQ(ucell->atoms[1].nw, 9); + EXPECT_EQ(ucell->nwmax, 9); + EXPECT_EQ(GlobalV::NLOCAL, 3 * 9); + // check itia2iat + EXPECT_EQ(ucell->itia2iat.getSize(), 4); + EXPECT_EQ(ucell->itia2iat(0, 0), 0); + EXPECT_EQ(ucell->itia2iat(0, 1), 0); + EXPECT_EQ(ucell->itia2iat(1, 0), 1); + EXPECT_EQ(ucell->itia2iat(1, 1), 2); + // check iat2iwt + EXPECT_EQ(ucell->get_npol(), 1); + EXPECT_EQ(ucell->get_iat2iwt()[0], 0); + EXPECT_EQ(ucell->get_iat2iwt()[1], 9); + EXPECT_EQ(ucell->get_iat2iwt()[2], 18); + // check itiaiw2iwt + EXPECT_EQ(ucell->itiaiw2iwt(0, 0, 0), 0); + EXPECT_EQ(ucell->itiaiw2iwt(0, 0, 1), 1); + EXPECT_EQ(ucell->itiaiw2iwt(0, 0, 8), 8); + EXPECT_EQ(ucell->itiaiw2iwt(1, 0, 0), 9); + EXPECT_EQ(ucell->itiaiw2iwt(1, 1, 0), 18); + // check itia2iat + EXPECT_EQ(ucell->itia2iat.getSize(), 4); + EXPECT_EQ(ucell->itia2iat(0, 0), 0); + EXPECT_EQ(ucell->itia2iat(0, 1), 0); + EXPECT_EQ(ucell->itia2iat(1, 0), 1); + EXPECT_EQ(ucell->itia2iat(1, 1), 2); + // check iwt2iat + EXPECT_EQ(ucell->iwt2iat[0], 0); + EXPECT_EQ(ucell->iwt2iat[10], 1); + EXPECT_EQ(ucell->iwt2iat[20], 2); + // check iwt2iw + EXPECT_EQ(ucell->iwt2iw[0], 0); + EXPECT_EQ(ucell->iwt2iw[10], 1); + EXPECT_EQ(ucell->iwt2iw[20], 2); } -TEST_F(UcellTest,CalNwfc2) -{ - GlobalV::NSPIN = 4; - GlobalV::BASIS_TYPE = "lcao"; - ucell->read_cell_pseudopots(pp_dir,ofs); - EXPECT_FALSE(ucell->atoms[0].ncpp.has_so); - EXPECT_FALSE(ucell->atoms[1].ncpp.has_so); - ucell->cal_nwfc(ofs); - EXPECT_EQ(GlobalV::NLOCAL,3*9*2); +TEST_F(UcellTest, CalNwfc2) { + GlobalV::NSPIN = 4; + GlobalV::BASIS_TYPE = "lcao"; + ucell->read_cell_pseudopots(pp_dir, ofs); + EXPECT_FALSE(ucell->atoms[0].ncpp.has_so); + EXPECT_FALSE(ucell->atoms[1].ncpp.has_so); + ucell->cal_nwfc(ofs); + EXPECT_EQ(GlobalV::NLOCAL, 3 * 9 * 2); } -TEST_F(UcellDeathTest,CheckStructure) -{ - ucell->read_cell_pseudopots(pp_dir,ofs); - EXPECT_FALSE(ucell->atoms[0].ncpp.has_so); - EXPECT_FALSE(ucell->atoms[1].ncpp.has_so); - //trial 1 - testing::internal::CaptureStdout(); - EXPECT_NO_THROW(ucell->check_structure(0.2)); - output = testing::internal::GetCapturedStdout(); - EXPECT_THAT(output,testing::HasSubstr("WARNING: Some atoms are too close!!!")); - //trial 2 - testing::internal::CaptureStdout(); - EXPECT_EXIT(ucell->check_structure(0.4),::testing::ExitedWithCode(0),""); - output = testing::internal::GetCapturedStdout(); - EXPECT_THAT(output,testing::HasSubstr("The structure is unreasonable!")); - //trial 3 - ucell->atoms[0].ncpp.psd = "arbitrary"; - testing::internal::CaptureStdout(); - EXPECT_NO_THROW(ucell->check_structure(0.2)); - output = testing::internal::GetCapturedStdout(); - EXPECT_THAT(output,testing::HasSubstr("Notice: symbol 'arbitrary' is not an element symbol!!!! set the covalent radius to be 0.")); +TEST_F(UcellDeathTest, CheckStructure) { + ucell->read_cell_pseudopots(pp_dir, ofs); + EXPECT_FALSE(ucell->atoms[0].ncpp.has_so); + EXPECT_FALSE(ucell->atoms[1].ncpp.has_so); + // trial 1 + testing::internal::CaptureStdout(); + double factor = 0.2; + EXPECT_NO_THROW(Check_Atomic_Stru::check_atomic_stru(*ucell, factor)); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, + testing::HasSubstr("WARNING: Some atoms are too close!!!")); + // trial 2 + testing::internal::CaptureStdout(); + factor = 0.4; + EXPECT_EXIT(Check_Atomic_Stru::check_atomic_stru(*ucell, factor), + ::testing::ExitedWithCode(0), + ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("The structure is unreasonable!")); + // trial 3 + ucell->atoms[0].label = "arbitrary"; + testing::internal::CaptureStdout(); + factor = 0.2; + EXPECT_NO_THROW(Check_Atomic_Stru::check_atomic_stru(*ucell, factor)); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT( + output, + testing::HasSubstr("Notice: symbol 'arbitrary' is not an element " + "symbol!!!! set the covalent radius to be 0.")); + // trial 4 + ucell->atoms[0].label = "Fe1"; + testing::internal::CaptureStdout(); + factor = 0.2; + EXPECT_NO_THROW(Check_Atomic_Stru::check_atomic_stru(*ucell, factor)); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, + testing::HasSubstr("WARNING: Some atoms are too close!!!")); } -TEST_F(UcellDeathTest,ReadPseudoWarning1) -{ - GlobalV::global_pseudo_dir = pp_dir; - GlobalV::out_element_info = 1; - GlobalV::MIN_DIST_COEF = 0.2; - ucell->pseudo_fn[1] = "H_sr_lda.upf"; - testing::internal::CaptureStdout(); - EXPECT_EXIT(ucell->read_pseudo(ofs),::testing::ExitedWithCode(0),""); - output = testing::internal::GetCapturedStdout(); - EXPECT_THAT(output,testing::HasSubstr("All DFT functional must consistent.")); +TEST_F(UcellDeathTest, ReadPseudoWarning1) { + GlobalV::global_pseudo_dir = pp_dir; + GlobalV::out_element_info = 1; + GlobalV::MIN_DIST_COEF = 0.2; + ucell->pseudo_fn[1] = "H_sr_lda.upf"; + testing::internal::CaptureStdout(); + EXPECT_EXIT(ucell->read_pseudo(ofs), ::testing::ExitedWithCode(0), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, + testing::HasSubstr("All DFT functional must consistent.")); } -TEST_F(UcellDeathTest,ReadPseudoWarning2) -{ - GlobalV::global_pseudo_dir = pp_dir; - GlobalV::out_element_info = 1; - GlobalV::MIN_DIST_COEF = 0.2; - ucell->pseudo_fn[0] = "Al_ONCV_PBE-1.0.upf"; - testing::internal::CaptureStdout(); - EXPECT_NO_THROW(ucell->read_pseudo(ofs)); - output = testing::internal::GetCapturedStdout(); - EXPECT_THAT(output,testing::HasSubstr("Warning: the number of valence electrons in pseudopotential > 3 for Al: [Ne] 3s2 3p1")); +TEST_F(UcellDeathTest, ReadPseudoWarning2) { + GlobalV::global_pseudo_dir = pp_dir; + GlobalV::out_element_info = 1; + GlobalV::MIN_DIST_COEF = 0.2; + ucell->pseudo_fn[0] = "Al_ONCV_PBE-1.0.upf"; + testing::internal::CaptureStdout(); + EXPECT_NO_THROW(ucell->read_pseudo(ofs)); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT( + output, + testing::HasSubstr("Warning: the number of valence electrons in " + "pseudopotential > 3 for Al: [Ne] 3s2 3p1")); } -TEST_F(UcellTest, CalNelec) -{ +TEST_F(UcellTest, CalNelec) { ucell->read_cell_pseudopots(pp_dir, ofs); - EXPECT_EQ(4,ucell->atoms[0].ncpp.zv); - EXPECT_EQ(1,ucell->atoms[1].ncpp.zv); - EXPECT_EQ(1,ucell->atoms[0].na); - EXPECT_EQ(2,ucell->atoms[1].na); + EXPECT_EQ(4, ucell->atoms[0].ncpp.zv); + EXPECT_EQ(1, ucell->atoms[1].ncpp.zv); + EXPECT_EQ(1, ucell->atoms[0].na); + EXPECT_EQ(2, ucell->atoms[1].na); double nelec = 0; ucell->cal_nelec(nelec); - EXPECT_DOUBLE_EQ(6,nelec); + EXPECT_DOUBLE_EQ(6, nelec); } #ifdef __MPI #include "mpi.h" -int main(int argc, char **argv) -{ - MPI_Init(&argc, &argv); - testing::InitGoogleTest(&argc, argv); +int main(int argc, char** argv) { + MPI_Init(&argc, &argv); + testing::InitGoogleTest(&argc, argv); - MPI_Comm_size(MPI_COMM_WORLD,&GlobalV::NPROC); - MPI_Comm_rank(MPI_COMM_WORLD,&GlobalV::MY_RANK); + MPI_Comm_size(MPI_COMM_WORLD, &GlobalV::NPROC); + MPI_Comm_rank(MPI_COMM_WORLD, &GlobalV::MY_RANK); - int result = RUN_ALL_TESTS(); - MPI_Finalize(); - return result; + int result = RUN_ALL_TESTS(); + MPI_Finalize(); + return result; } #endif diff --git a/source/module_cell/unitcell.cpp b/source/module_cell/unitcell.cpp index 836f3c825e..dcd564f345 100644 --- a/source/module_cell/unitcell.cpp +++ b/source/module_cell/unitcell.cpp @@ -12,7 +12,6 @@ #include "../module_basis/module_ao/ORB_read.h" // to use 'ORB' -- mohan 2021-01-30 #endif #include "module_base/atom_in.h" -#include "module_base/element_covalent_radius.h" #include "module_base/element_elec_config.h" #include "module_base/global_file.h" #include "module_base/parallel_common.h" @@ -27,8 +26,7 @@ #include "module_ri/serialization_cereal.h" #endif -UnitCell::UnitCell() -{ +UnitCell::UnitCell() { if (GlobalV::test_unitcell) ModuleBase::TITLE("unitcell", "Constructor"); Coordinate = "Direct"; @@ -69,8 +67,7 @@ UnitCell::UnitCell() set_atom_flag = false; } -UnitCell::~UnitCell() -{ +UnitCell::~UnitCell() { delete[] atom_label; delete[] atom_mass; delete[] pseudo_fn; @@ -81,16 +78,14 @@ UnitCell::~UnitCell() delete[] iwt2iat; delete[] iwt2iw; delete[] lc; - if (set_atom_flag) - { + if (set_atom_flag) { delete[] atoms; } } #include "module_base/parallel_common.h" #ifdef __MPI -void UnitCell::bcast_unitcell() -{ +void UnitCell::bcast_unitcell() { if (GlobalV::test_unitcell) ModuleBase::TITLE("UnitCell", "bcast_unitcell"); Parallel_Common::bcast_string(Coordinate); @@ -144,36 +139,33 @@ void UnitCell::bcast_unitcell() Parallel_Common::bcast_double(latvec_supercell.e33); Parallel_Common::bcast_double(magnet.start_magnetization, ntype); - if (GlobalV::NSPIN == 4) - { + if (GlobalV::NSPIN == 4) { Parallel_Common::bcast_double(magnet.ux_[0]); Parallel_Common::bcast_double(magnet.ux_[1]); Parallel_Common::bcast_double(magnet.ux_[2]); } - for (int i = 0; i < ntype; i++) - { + for (int i = 0; i < ntype; i++) { atoms[i].bcast_atom(); // init tau and mbl array } #ifdef __EXX - ModuleBase::bcast_data_cereal(GlobalC::exx_info.info_ri.files_abfs, MPI_COMM_WORLD, 0); + ModuleBase::bcast_data_cereal(GlobalC::exx_info.info_ri.files_abfs, + MPI_COMM_WORLD, + 0); #endif return; } -void UnitCell::bcast_unitcell2() -{ - for (int i = 0; i < ntype; i++) - { +void UnitCell::bcast_unitcell2() { + for (int i = 0; i < ntype; i++) { atoms[i].bcast_atom2(); } return; } #endif -void UnitCell::print_cell(std::ofstream& ofs) const -{ +void UnitCell::print_cell(std::ofstream& ofs) const { if (GlobalV::test_unitcell) ModuleBase::TITLE("UnitCell", "print_cell"); @@ -196,8 +188,7 @@ void UnitCell::print_cell(std::ofstream& ofs) const return; } -void UnitCell::print_cell_cif(const std::string& fn) const -{ +void UnitCell::print_cell_cif(const std::string& fn) const { if (GlobalV::test_unitcell) ModuleBase::TITLE("UnitCell", "print_cell_cif"); @@ -224,9 +215,12 @@ void UnitCell::print_cell_cif(const std::string& fn) const // ofs << "_cell_angle_beta " << "90.00" << std::endl; // ofs << "_cell_angle_gamma " << "90.00" << std::endl; // xiaohui modify alpha, beta and gamma 2015-09-29 - double angle_alpha = acos((a2 * a3) / (a2.norm() * a3.norm())) / ModuleBase::PI * 180.0; - double angle_beta = acos((a1 * a3) / (a1.norm() * a3.norm())) / ModuleBase::PI * 180.0; - double angle_gamma = acos((a1 * a2) / (a1.norm() * a2.norm())) / ModuleBase::PI * 180.0; + double angle_alpha + = acos((a2 * a3) / (a2.norm() * a3.norm())) / ModuleBase::PI * 180.0; + double angle_beta + = acos((a1 * a3) / (a1.norm() * a3.norm())) / ModuleBase::PI * 180.0; + double angle_gamma + = acos((a1 * a2) / (a1.norm() * a2.norm())) / ModuleBase::PI * 180.0; ofs << "_cell_angle_alpha " << angle_alpha << std::endl; ofs << "_cell_angle_beta " << angle_beta << std::endl; ofs << "_cell_angle_gamma " << angle_gamma << std::endl; @@ -241,12 +235,11 @@ void UnitCell::print_cell_cif(const std::string& fn) const ofs << "_atom_site_fract_x" << std::endl; ofs << "_atom_site_fract_y" << std::endl; ofs << "_atom_site_fract_z" << std::endl; - for (int it = 0; it < ntype; it++) - { - for (int ia = 0; ia < atoms[it].na; ia++) - { - ofs << atoms[it].label << " " << atoms[it].taud[ia].x << " " << atoms[it].taud[ia].y << " " - << atoms[it].taud[ia].z << std::endl; + for (int it = 0; it < ntype; it++) { + for (int ia = 0; ia < atoms[it].na; ia++) { + ofs << atoms[it].label << " " << atoms[it].taud[ia].x << " " + << atoms[it].taud[ia].y << " " << atoms[it].taud[ia].z + << std::endl; } } ofs.close(); @@ -272,8 +265,10 @@ void UnitCell::print_cell_xyz(const std::string& fn) const { for (int ia = 0; ia < atoms[it].na; ia++) { - ofs << atoms[it].label << " " << atoms[it].tau[ia].x * lat0 * 0.529177 << " " - << atoms[it].tau[ia].y * lat0 * 0.529177 << " " << atoms[it].tau[ia].z * lat0 * 0.529177 << std::endl; + ofs << atoms[it].label << " " << atoms[it].tau[ia].x * lat0 * +0.529177 << " " + << atoms[it].tau[ia].y * lat0 * 0.529177 << " " << +atoms[it].tau[ia].z * lat0 * 0.529177 << std::endl; } } @@ -282,18 +277,15 @@ void UnitCell::print_cell_xyz(const std::string& fn) const } */ -void UnitCell::set_iat2itia() -{ +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++) - { + 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; @@ -302,36 +294,28 @@ void UnitCell::set_iat2itia() return; } -std::map UnitCell::get_atom_Counts() const -{ +std::map UnitCell::get_atom_Counts() const { std::map atomCounts; - for (int it = 0; it < this->ntype; it++) - { + 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 UnitCell::get_orbital_Counts() const { std::map orbitalCounts; - for (int it = 0; it < this->ntype; it++) - { + 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> 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++) - { + 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 (lnchiCounts.find(it) == lnchiCounts.end()) { // If it doesn't exist, initialize an empty inner map lnchiCounts[it] = std::map(); } @@ -343,53 +327,42 @@ std::map> UnitCell::get_lnchi_Counts() const return lnchiCounts; } -std::vector UnitCell::get_atomLabels() const -{ +std::vector UnitCell::get_atomLabels() const { std::vector atomLabels(this->ntype); - for (int it = 0; it < this->ntype; it++) - { + for (int it = 0; it < this->ntype; it++) { atomLabels[it] = this->atoms[it].label; } return atomLabels; } -std::vector UnitCell::get_atomCounts() const -{ +std::vector UnitCell::get_atomCounts() const { std::vector atomCounts(this->ntype); - for (int it = 0; it < this->ntype; it++) - { + for (int it = 0; it < this->ntype; it++) { atomCounts[it] = this->atoms[it].na; } return atomCounts; } -std::vector> UnitCell::get_lnchiCounts() const -{ +std::vector> UnitCell::get_lnchiCounts() const { std::vector> lnchiCounts(this->ntype); - for (int it = 0; it < this->ntype; it++) - { + 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++) - { + for (int L = 0; L < this->atoms[it].nwl + 1; L++) { lnchiCounts[it][L] = this->atoms[it].l_nchi[L]; } } return lnchiCounts; } -void UnitCell::update_pos_tau(const double* pos) -{ +void UnitCell::update_pos_tau(const double* pos) { int iat = 0; - for (int it = 0; it < this->ntype; it++) - { + for (int it = 0; it < this->ntype; it++) { Atom* atom = &this->atoms[it]; - for (int ia = 0; ia < atom->na; ia++) - { - for (int ik = 0; ik < 3; ++ik) - { - if (atom->mbl[ia][ik]) - { - atom->dis[ia][ik] = pos[3 * iat + ik] / this->lat0 - atom->tau[ia][ik]; + for (int ia = 0; ia < atom->na; ia++) { + for (int ik = 0; ik < 3; ++ik) { + if (atom->mbl[ia][ik]) { + atom->dis[ia][ik] + = pos[3 * iat + ik] / this->lat0 - atom->tau[ia][ik]; atom->tau[ia][ik] = pos[3 * iat + ik] / this->lat0; } } @@ -405,16 +378,12 @@ void UnitCell::update_pos_tau(const double* pos) this->bcast_atoms_tau(); } -void UnitCell::update_pos_taud(double* posd_in) -{ +void UnitCell::update_pos_taud(double* posd_in) { int iat = 0; - for (int it = 0; it < this->ntype; it++) - { + for (int it = 0; it < this->ntype; it++) { Atom* atom = &this->atoms[it]; - for (int ia = 0; ia < atom->na; ia++) - { - for (int ik = 0; ik < 3; ++ik) - { + for (int ia = 0; ia < atom->na; ia++) { + for (int ik = 0; ik < 3; ++ik) { atom->taud[ia][ik] += posd_in[3 * iat + ik]; atom->dis[ia][ik] = posd_in[3 * iat + ik]; } @@ -427,16 +396,12 @@ void UnitCell::update_pos_taud(double* posd_in) } // posd_in is atomic displacements here liuyu 2023-03-22 -void UnitCell::update_pos_taud(const ModuleBase::Vector3* posd_in) -{ +void UnitCell::update_pos_taud(const ModuleBase::Vector3* posd_in) { int iat = 0; - for (int it = 0; it < this->ntype; it++) - { + for (int it = 0; it < this->ntype; it++) { Atom* atom = &this->atoms[it]; - for (int ia = 0; ia < atom->na; ia++) - { - for (int ik = 0; ik < 3; ++ik) - { + for (int ia = 0; ia < atom->na; ia++) { + for (int ik = 0; ik < 3; ++ik) { atom->taud[ia][ik] += posd_in[iat][ik]; atom->dis[ia][ik] = posd_in[iat][ik]; } @@ -448,14 +413,11 @@ void UnitCell::update_pos_taud(const ModuleBase::Vector3* posd_in) this->bcast_atoms_tau(); } -void UnitCell::update_vel(const ModuleBase::Vector3* vel_in) -{ +void UnitCell::update_vel(const ModuleBase::Vector3* vel_in) { int iat = 0; - for (int it = 0; it < this->ntype; ++it) - { + for (int it = 0; it < this->ntype; ++it) { Atom* atom = &this->atoms[it]; - for (int ia = 0; ia < atom->na; ++ia) - { + for (int ia = 0; ia < atom->na; ++ia) { this->atoms[it].vel[ia] = vel_in[iat]; ++iat; } @@ -463,19 +425,16 @@ void UnitCell::update_vel(const ModuleBase::Vector3* vel_in) assert(iat == this->nat); } -void UnitCell::periodic_boundary_adjustment() -{ +void UnitCell::periodic_boundary_adjustment() { //---------------------------------------------- // because of the periodic boundary condition // we need to adjust the atom positions, // first adjust direct coordinates, // then update them into cartesian coordinates, //---------------------------------------------- - for (int it = 0; it < this->ntype; it++) - { + for (int it = 0; it < this->ntype; it++) { Atom* atom = &this->atoms[it]; - for (int ia = 0; ia < atom->na; ia++) - { + for (int ia = 0; ia < atom->na; ia++) { // mohan update 2011-03-21 if (atom->taud[ia].x < 0) atom->taud[ia].x += 1.0; @@ -490,14 +449,17 @@ void UnitCell::periodic_boundary_adjustment() if (atom->taud[ia].z >= 1.0) atom->taud[ia].z -= 1.0; - if (atom->taud[ia].x < 0 || atom->taud[ia].y < 0 || atom->taud[ia].z < 0 || atom->taud[ia].x >= 1.0 - || atom->taud[ia].y >= 1.0 || atom->taud[ia].z >= 1.0) - { - GlobalV::ofs_warning << " it=" << it + 1 << " ia=" << ia + 1 << std::endl; - GlobalV::ofs_warning << "d=" << atom->taud[ia].x << " " << atom->taud[ia].y << " " << atom->taud[ia].z + if (atom->taud[ia].x < 0 || atom->taud[ia].y < 0 + || atom->taud[ia].z < 0 || atom->taud[ia].x >= 1.0 + || atom->taud[ia].y >= 1.0 || atom->taud[ia].z >= 1.0) { + GlobalV::ofs_warning << " it=" << it + 1 << " ia=" << ia + 1 << std::endl; - ModuleBase::WARNING_QUIT("Ions_Move_Basic::move_ions", - "the movement of atom is larger than the length of cell."); + GlobalV::ofs_warning << "d=" << atom->taud[ia].x << " " + << atom->taud[ia].y << " " + << atom->taud[ia].z << std::endl; + ModuleBase::WARNING_QUIT( + "Ions_Move_Basic::move_ions", + "the movement of atom is larger than the length of cell."); } atom->tau[ia] = atom->taud[ia] * this->latvec; @@ -506,19 +468,16 @@ void UnitCell::periodic_boundary_adjustment() return; } -void UnitCell::bcast_atoms_tau() -{ +void UnitCell::bcast_atoms_tau() { #ifdef __MPI MPI_Barrier(MPI_COMM_WORLD); - for (int i = 0; i < ntype; i++) - { + for (int i = 0; i < ntype; i++) { atoms[i].bcast_atom(); // bcast tau array } #endif } -void UnitCell::cal_ux() -{ +void UnitCell::cal_ux() { double amag, uxmod; int starting_it = 0; int starting_ia = 0; @@ -527,13 +486,12 @@ void UnitCell::cal_ux() magnet.lsign_ = false; ModuleBase::GlobalFunc::ZEROS(magnet.ux_, 3); - for (int it = 0; it < ntype; it++) - { - for (int ia = 0; ia < atoms[it].na; ia++) - { - amag = pow(atoms[it].m_loc_[ia].x, 2) + pow(atoms[it].m_loc_[ia].y, 2) + pow(atoms[it].m_loc_[ia].z, 2); - if (amag > 1e-6) - { + for (int it = 0; it < ntype; it++) { + for (int ia = 0; ia < atoms[it].na; ia++) { + amag = pow(atoms[it].m_loc_[ia].x, 2) + + pow(atoms[it].m_loc_[ia].y, 2) + + pow(atoms[it].m_loc_[ia].z, 2); + if (amag > 1e-6) { magnet.ux_[0] = atoms[it].m_loc_[ia].x; magnet.ux_[1] = atoms[it].m_loc_[ia].y; magnet.ux_[2] = atoms[it].m_loc_[ia].z; @@ -547,38 +505,37 @@ void UnitCell::cal_ux() break; } // whether the initial magnetizations is parallel - for (int it = starting_it; it < ntype; it++) - { - for (int ia = 0; ia < atoms[it].na; ia++) - { - if (it > starting_it || ia > starting_ia) - { - magnet.lsign_ = magnet.lsign_ && judge_parallel(magnet.ux_, atoms[it].m_loc_[ia]); + for (int it = starting_it; it < ntype; it++) { + for (int ia = 0; ia < atoms[it].na; ia++) { + if (it > starting_it || ia > starting_ia) { + magnet.lsign_ + = magnet.lsign_ + && judge_parallel(magnet.ux_, atoms[it].m_loc_[ia]); } } } - if (magnet.lsign_) - { - uxmod = pow(magnet.ux_[0], 2) + pow(magnet.ux_[1], 2) + pow(magnet.ux_[2], 2); - if (uxmod < 1e-6) - { + if (magnet.lsign_) { + uxmod = pow(magnet.ux_[0], 2) + pow(magnet.ux_[1], 2) + + pow(magnet.ux_[2], 2); + if (uxmod < 1e-6) { ModuleBase::WARNING_QUIT("cal_ux", "wrong uxmod"); } - for (int i = 0; i < 3; i++) - { + for (int i = 0; i < 3; i++) { magnet.ux_[i] *= 1 / sqrt(uxmod); } // std::cout<<" Fixed quantization axis for GGA: " - //< b) -{ +bool UnitCell::judge_parallel(double a[3], ModuleBase::Vector3 b) { bool jp = false; double cross; - cross = pow((a[1] * b.z - a[2] * b.y), 2) + pow((a[2] * b.x - a[0] * b.z), 2) + pow((a[0] * b.y - a[1] * b.x), 2); + cross = pow((a[1] * b.z - a[2] * b.y), 2) + + pow((a[2] * b.x - a[0] * b.z), 2) + + pow((a[0] * b.y - a[1] * b.x), 2); jp = (fabs(cross) < 1e-6); return jp; } @@ -586,8 +543,7 @@ bool UnitCell::judge_parallel(double a[3], ModuleBase::Vector3 b) //============================================================== // Calculate various lattice related quantities for given latvec //============================================================== -void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) -{ +void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) { ModuleBase::TITLE("UnitCell", "setup_cell"); // (1) init mag assert(ntype > 0); @@ -602,33 +558,56 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) bool ok2 = true; // (3) read in atom information - if (GlobalV::MY_RANK == 0) - { + if (GlobalV::MY_RANK == 0) { // open "atom_unitcell" file. std::ifstream ifa(fn.c_str(), std::ios::in); - if (!ifa) - { + if (!ifa) { GlobalV::ofs_warning << fn; ok = false; } - if (ok) - { + if (ok) { log << "\n\n\n\n"; - log << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; - log << " | |" << std::endl; - log << " | Reading atom information in 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 << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" + ">>>>>>>>>>>>" + << std::endl; + log << " | " + " |" + << std::endl; + log << " | Reading atom information in 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\n\n\n"; log << " READING UNITCELL INFORMATION" << std::endl; @@ -646,19 +625,19 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) #ifdef __MPI Parallel_Common::bcast_bool(ok); Parallel_Common::bcast_bool(ok2); - if (GlobalV::NSPIN == 4) - { + if (GlobalV::NSPIN == 4) { Parallel_Common::bcast_bool(GlobalV::DOMAG); Parallel_Common::bcast_bool(GlobalV::DOMAG_Z); } #endif - if (!ok) - { - ModuleBase::WARNING_QUIT("UnitCell::setup_cell", "Can not find the file containing atom positions.!"); + 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 (!ok2) { + ModuleBase::WARNING_QUIT("UnitCell::setup_cell", + "Something wrong during read_atom_positions."); } #ifdef __MPI @@ -666,16 +645,14 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) #endif // after read STRU, calculate initial total magnetization when NSPIN=2 - if (GlobalV::NSPIN == 2 && !GlobalV::TWO_EFERMI) - { - for (int it = 0; it < this->ntype; it++) - { - for (int ia = 0; ia < this->atoms[it].na; ia++) - { + if (GlobalV::NSPIN == 2 && !GlobalV::TWO_EFERMI) { + for (int it = 0; it < this->ntype; it++) { + for (int ia = 0; ia < this->atoms[it].na; ia++) { GlobalV::nupdown += this->atoms[it].mag[ia]; } } - GlobalV::ofs_running << " The readin total magnetization is " << GlobalV::nupdown << std::endl; + GlobalV::ofs_running << " The readin total magnetization is " + << GlobalV::nupdown << std::endl; } //======================================================== @@ -685,16 +662,16 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) //======================================================== assert(lat0 > 0.0); this->omega = std::abs(latvec.Det()) * this->lat0 * lat0 * lat0; - if (this->omega <= 0) - { + if (this->omega <= 0) { std::cout << "The volume is negative: " << this->omega << std::endl; ModuleBase::WARNING_QUIT("setup_cell", "omega <= 0 ."); - } - else - { + } else { log << std::endl; ModuleBase::GlobalFunc::OUT(log, "Volume (Bohr^3)", this->omega); - ModuleBase::GlobalFunc::OUT(log, "Volume (A^3)", this->omega * pow(ModuleBase::BOHR_TO_A, 3)); + ModuleBase::GlobalFunc::OUT(log, + "Volume (A^3)", + this->omega + * pow(ModuleBase::BOHR_TO_A, 3)); } //========================================================== @@ -713,8 +690,13 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) this->invGGT0 = 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); + 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); // OUT(log,"lattice center x",latcenter.x); // OUT(log,"lattice center y",latcenter.y); // OUT(log,"lattice center z",latcenter.z); @@ -725,8 +707,7 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) this->set_iat2itia(); #ifdef USE_PAW - if (GlobalV::use_paw) - { + if (GlobalV::use_paw) { GlobalC::paw_cell.set_libpaw_cell(latvec, lat0); int* typat; @@ -736,10 +717,8 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) xred = new double[nat * 3]; int iat = 0; - for (int it = 0; it < ntype; it++) - { - for (int ia = 0; ia < atoms[it].na; ia++) - { + for (int it = 0; it < ntype; it++) { + for (int ia = 0; ia < atoms[it].na; ia++) { typat[iat] = it + 1; // Fortran index starts from 1 !!!! xred[iat * 3 + 0] = atoms[it].taud[ia].x; xred[iat * 3 + 1] = atoms[it].taud[ia].y; @@ -761,50 +740,72 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) return; } -void UnitCell::read_pseudo(std::ofstream& ofs) -{ +void UnitCell::read_pseudo(std::ofstream& ofs) { // read in non-local pseudopotential and ouput the projectors. ofs << "\n\n\n\n"; - ofs << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; - ofs << " | |" << std::endl; - ofs << " | Reading pseudopotentials files: |" << std::endl; - ofs << " | The pseudopotential file is in UPF format. The 'NC' indicates that |" << std::endl; - ofs << " | the type of pseudopotential is 'norm conserving'. Functional of |" << std::endl; - ofs << " | exchange and correlation is decided by 4 given parameters in UPF |" << std::endl; - ofs << " | file. We also read in the 'core correction' if there exists. |" << std::endl; - ofs << " | Also we can read the valence electrons number and the maximal |" << std::endl; - ofs << " | angular momentum used in this pseudopotential. We also read in the |" << std::endl; - ofs << " | trail wave function, trail atomic density and local-pseudopotential|" << std::endl; - ofs << " | on logrithmic grid. The non-local pseudopotential projector is also|" << std::endl; - ofs << " | read in if there is any. |" << std::endl; - ofs << " | |" << std::endl; - ofs << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; + ofs << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" + ">>>>" + << std::endl; + ofs << " | " + " |" + << std::endl; + ofs << " | Reading pseudopotentials files: " + " |" + << std::endl; + ofs << " | The pseudopotential file is in UPF format. The 'NC' indicates " + "that |" + << std::endl; + ofs << " | the type of pseudopotential is 'norm conserving'. Functional of " + " |" + << std::endl; + ofs << " | exchange and correlation is decided by 4 given parameters in " + "UPF |" + << std::endl; + ofs << " | file. We also read in the 'core correction' if there exists. " + " |" + << std::endl; + ofs << " | Also we can read the valence electrons number and the maximal " + " |" + << std::endl; + ofs << " | angular momentum used in this pseudopotential. We also read in " + "the |" + << std::endl; + ofs << " | trail wave function, trail atomic density and " + "local-pseudopotential|" + << std::endl; + ofs << " | on logrithmic grid. The non-local pseudopotential projector is " + "also|" + << std::endl; + ofs << " | read in if there is any. " + " |" + << std::endl; + ofs << " | " + " |" + << std::endl; + ofs << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" + "<<<<" + << std::endl; ofs << "\n\n\n\n"; read_cell_pseudopots(GlobalV::global_pseudo_dir, ofs); - if (GlobalV::MY_RANK == 0) - { - for (int it = 0; it < this->ntype; it++) - { + if (GlobalV::MY_RANK == 0) { + for (int it = 0; it < this->ntype; it++) { Atom* atom = &atoms[it]; - if (!(atom->label_orb.empty())) - { + if (!(atom->label_orb.empty())) { compare_atom_labels(atom->label_orb, atom->ncpp.psd); } } - if (GlobalV::out_element_info) - { - for (int i = 0; i < this->ntype; i++) - { + if (GlobalV::out_element_info) { + for (int i = 0; i < this->ntype; i++) { ModuleBase::Global_File::make_dir_atom(this->atoms[i].label); } - for (int it = 0; it < ntype; it++) - { + for (int it = 0; it < ntype; it++) { Atom* atom = &atoms[it]; std::stringstream ss; - ss << GlobalV::global_out_dir << atom->label << "/" << atom->label << ".NONLOCAL"; + ss << GlobalV::global_out_dir << atom->label << "/" + << atom->label << ".NONLOCAL"; std::ofstream ofs(ss.str().c_str()); ofs << "
" << std::endl; @@ -819,18 +820,16 @@ void UnitCell::read_pseudo(std::ofstream& ofs) ofs << "\n" << std::endl; ofs << std::setw(10) << atom->ncpp.nbeta << "\t" << "nummber of projectors." << std::endl; - for (int ib = 0; ib < atom->ncpp.nbeta; ib++) - { - for (int ib2 = 0; ib2 < atom->ncpp.nbeta; ib2++) - { - ofs << std::setw(10) << atom->ncpp.lll[ib] << " " << atom->ncpp.lll[ib2] << " " + for (int ib = 0; ib < atom->ncpp.nbeta; ib++) { + for (int ib2 = 0; ib2 < atom->ncpp.nbeta; ib2++) { + ofs << std::setw(10) << atom->ncpp.lll[ib] << " " + << atom->ncpp.lll[ib2] << " " << atom->ncpp.dion(ib, ib2) << std::endl; } } ofs << "" << std::endl; - for (int i = 0; i < atom->ncpp.nbeta; i++) - { + for (int i = 0; i < atom->ncpp.nbeta; i++) { ofs << "" << std::endl; ofs << std::setw(10) << i << "\t" << "the index of projectors." << std::endl; @@ -840,10 +839,8 @@ void UnitCell::read_pseudo(std::ofstream& ofs) // mohan add // only keep the nonzero part. int cut_mesh = atom->ncpp.mesh; - for (int j = atom->ncpp.mesh - 1; j >= 0; --j) - { - if (std::abs(atom->ncpp.betar(i, j)) > 1.0e-10) - { + for (int j = atom->ncpp.mesh - 1; j >= 0; --j) { + if (std::abs(atom->ncpp.betar(i, j)) > 1.0e-10) { cut_mesh = j; break; } @@ -854,10 +851,10 @@ void UnitCell::read_pseudo(std::ofstream& ofs) ofs << std::setw(10) << cut_mesh << "\t" << "the number of mesh points." << std::endl; - for (int j = 0; j < cut_mesh; ++j) - { - ofs << std::setw(15) << atom->ncpp.r[j] << std::setw(15) << atom->ncpp.betar(i, j) - << std::setw(15) << atom->ncpp.rab[j] << std::endl; + for (int j = 0; j < cut_mesh; ++j) { + ofs << std::setw(15) << atom->ncpp.r[j] << std::setw(15) + << atom->ncpp.betar(i, j) << std::setw(15) + << atom->ncpp.rab[j] << std::endl; } ofs << "" << std::endl; } @@ -871,25 +868,23 @@ void UnitCell::read_pseudo(std::ofstream& ofs) bcast_unitcell2(); #endif - for (int it = 0; it < ntype; it++) - { - if (atoms[0].ncpp.xc_func != atoms[it].ncpp.xc_func) - { - GlobalV::ofs_warning << "\n type " << atoms[0].label << " functional is " << atoms[0].ncpp.xc_func; + for (int it = 0; it < ntype; it++) { + if (atoms[0].ncpp.xc_func != atoms[it].ncpp.xc_func) { + GlobalV::ofs_warning << "\n type " << atoms[0].label + << " functional is " << atoms[0].ncpp.xc_func; - GlobalV::ofs_warning << "\n type " << atoms[it].label << " functional is " << atoms[it].ncpp.xc_func + GlobalV::ofs_warning << "\n type " << atoms[it].label + << " functional is " << atoms[it].ncpp.xc_func << std::endl; - ModuleBase::WARNING_QUIT("setup_cell", "All DFT functional must consistent."); + ModuleBase::WARNING_QUIT("setup_cell", + "All DFT functional must consistent."); } - if (atoms[it].ncpp.tvanp) - { + if (atoms[it].ncpp.tvanp) { GlobalV::use_uspp = true; } } - check_structure(GlobalV::MIN_DIST_COEF); - // setup the total number of PAOs cal_natomwfc(ofs); @@ -897,54 +892,63 @@ void UnitCell::read_pseudo(std::ofstream& ofs) cal_nwfc(ofs); // Check whether the number of valence is minimum - if (GlobalV::MY_RANK == 0) - { + if (GlobalV::MY_RANK == 0) { int abtype = 0; - for (int it = 0; it < ntype; it++) - { - if (ModuleBase::MinZval.find(atoms[it].ncpp.psd) != ModuleBase::MinZval.end()) - { - if (atoms[it].ncpp.zv > ModuleBase::MinZval.at(atoms[it].ncpp.psd)) - { + for (int it = 0; it < ntype; it++) { + if (ModuleBase::MinZval.find(atoms[it].ncpp.psd) + != ModuleBase::MinZval.end()) { + if (atoms[it].ncpp.zv + > ModuleBase::MinZval.at(atoms[it].ncpp.psd)) { abtype += 1; - if (abtype == 1) - { - std::cout << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" + if (abtype == 1) { + std::cout << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" + "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" "%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - ofs << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" + ofs << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" + "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" "%%%%%%%%%%%%%%%%%%%%%" << std::endl; } - std::cout << " Warning: the number of valence electrons in pseudopotential > " + std::cout << " Warning: the number of valence electrons in " + "pseudopotential > " << ModuleBase::MinZval.at(atoms[it].ncpp.psd); - std::cout << " for " << atoms[it].ncpp.psd << ": " << ModuleBase::EleConfig.at(atoms[it].ncpp.psd) + std::cout << " for " << atoms[it].ncpp.psd << ": " + << ModuleBase::EleConfig.at(atoms[it].ncpp.psd) << std::endl; - ofs << " Warning: the number of valence electrons in pseudopotential > " + ofs << " Warning: the number of valence electrons in " + "pseudopotential > " << ModuleBase::MinZval.at(atoms[it].ncpp.psd); - ofs << " for " << atoms[it].ncpp.psd << ": " << ModuleBase::EleConfig.at(atoms[it].ncpp.psd) + ofs << " for " << atoms[it].ncpp.psd << ": " + << ModuleBase::EleConfig.at(atoms[it].ncpp.psd) << std::endl; } } } - if (abtype > 0) - { - std::cout << " Pseudopotentials with additional electrons can yield (more) accurate outcomes, but may be " + if (abtype > 0) { + std::cout << " Pseudopotentials with additional electrons can " + "yield (more) accurate outcomes, but may be " "less efficient." << std::endl; - std::cout << " If you're confident that your chosen pseudopotential is appropriate, you can safely ignore " - "this warning." - << std::endl; - std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" + std::cout + << " If you're confident that your chosen pseudopotential is " + "appropriate, you can safely ignore " + "this warning." + << std::endl; + std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" + "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" "%%%%%%%%%%%%\n" << std::endl; - ofs << " Pseudopotentials with additional electrons can yield (more) accurate outcomes, but may be less " + ofs << " Pseudopotentials with additional electrons can yield " + "(more) accurate outcomes, but may be less " "efficient." << std::endl; - ofs << " If you're confident that your chosen pseudopotential is appropriate, you can safely ignore this " + ofs << " If you're confident that your chosen pseudopotential is " + "appropriate, you can safely ignore this " "warning." << std::endl; - ofs << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" + ofs << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" + "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" "%%%%%%%"; ModuleBase::GlobalFunc::OUT(ofs, ""); } @@ -966,8 +970,7 @@ void UnitCell::read_pseudo(std::ofstream& ofs) // atoms[].stapos_wf // GlobalV::NBANDS //=========================================== -void UnitCell::cal_nwfc(std::ofstream& log) -{ +void UnitCell::cal_nwfc(std::ofstream& log) { ModuleBase::TITLE("UnitCell", "cal_nwfc"); assert(ntype > 0); assert(nat > 0); @@ -975,8 +978,7 @@ void UnitCell::cal_nwfc(std::ofstream& log) //=========================== // (1) set iw2l, iw2n, iw2m //=========================== - for (int it = 0; it < ntype; it++) - { + for (int it = 0; it < ntype; it++) { this->atoms[it].set_index(); } @@ -985,8 +987,7 @@ void UnitCell::cal_nwfc(std::ofstream& log) //=========================== this->namax = 0; this->nwmax = 0; - for (int it = 0; it < ntype; it++) - { + for (int it = 0; it < ntype; it++) { this->namax = std::max(atoms[it].na, namax); this->nwmax = std::max(atoms[it].nw, nwmax); } @@ -999,22 +1000,19 @@ void UnitCell::cal_nwfc(std::ofstream& log) // (3) set nwfc and stapos_wf //=========================== GlobalV::NLOCAL = 0; - for (int it = 0; it < ntype; it++) - { + for (int it = 0; it < ntype; it++) { atoms[it].stapos_wf = GlobalV::NLOCAL; const int nlocal_it = atoms[it].nw * atoms[it].na; - if (GlobalV::NSPIN != 4) - { + if (GlobalV::NSPIN != 4) { GlobalV::NLOCAL += nlocal_it; - } - else - { + } else { GlobalV::NLOCAL += nlocal_it * 2; // zhengdy-soc } // for tests // OUT(GlobalV::ofs_running,ss1.str(),nlocal_it); - // OUT(GlobalV::ofs_running,"start position of local orbitals",atoms[it].stapos_wf); + // OUT(GlobalV::ofs_running,"start position of local + //orbitals",atoms[it].stapos_wf); } // OUT(GlobalV::ofs_running,"NLOCAL",GlobalV::NLOCAL); @@ -1036,14 +1034,11 @@ void UnitCell::cal_nwfc(std::ofstream& log) this->set_iat2iwt(GlobalV::NPOL); int iat = 0; int iwt = 0; - for (int it = 0; it < ntype; it++) - { - for (int ia = 0; ia < atoms[it].na; ia++) - { + for (int it = 0; it < ntype; it++) { + for (int ia = 0; ia < atoms[it].na; ia++) { this->itia2iat(it, ia) = iat; // this->iat2ia[iat] = ia; - for (int iw = 0; iw < atoms[it].nw * GlobalV::NPOL; iw++) - { + for (int iw = 0; iw < atoms[it].nw * GlobalV::NPOL; iw++) { // this->itiaiw2iwt(it, ia, iw) = iwt; this->iwt2iat[iwt] = iat; this->iwt2iw[iwt] = iw; @@ -1058,17 +1053,14 @@ void UnitCell::cal_nwfc(std::ofstream& log) //======================== this->lmax = 0; this->nmax = 0; - for (int it = 0; it < ntype; it++) - { + for (int it = 0; it < ntype; it++) { lmax = std::max(lmax, atoms[it].nwl); - for (int l = 0; l < atoms[it].nwl + 1; l++) - { + for (int l = 0; l < atoms[it].nwl + 1; l++) { nmax = std::max(nmax, atoms[it].l_nchi[l]); } int nchi = 0; - for (int l = 0; l < atoms[it].nwl + 1; l++) - { + for (int l = 0; l < atoms[it].nwl + 1; l++) { nchi += atoms[it].l_nchi[l]; } this->nmax_total = std::max(nmax_total, nchi); @@ -1078,12 +1070,9 @@ void UnitCell::cal_nwfc(std::ofstream& log) // (6) set lmax_ppwf //======================= this->lmax_ppwf = 0; - for (int it = 0; it < ntype; it++) - { - for (int ic = 0; ic < atoms[it].ncpp.nchi; ic++) - { - if (lmax_ppwf < atoms[it].ncpp.lchi[ic]) - { + for (int it = 0; it < ntype; it++) { + for (int ic = 0; ic < atoms[it].ncpp.nchi; ic++) { + if (lmax_ppwf < atoms[it].ncpp.lchi[ic]) { this->lmax_ppwf = atoms[it].ncpp.lchi[ic]; } } @@ -1092,8 +1081,8 @@ void UnitCell::cal_nwfc(std::ofstream& log) /* for(int it=0; it< ntype; it++) { - std::cout << " label=" << it << " nbeta=" << atoms[it].nbeta << std::endl; - for(int ic=0; ic GlobalV::NLOCAL !"); + // ModuleBase::WARNING_QUIT("cal_nwfc","NBANDS must > GlobalV::NLOCAL + //!"); // } // } } @@ -1125,8 +1115,7 @@ void UnitCell::cal_nwfc(std::ofstream& log) return; } -void UnitCell::set_iat2iwt(const int& npol_in) -{ +void UnitCell::set_iat2iwt(const int& npol_in) { #ifdef __DEBUG assert(npol_in == 1 || npol_in == 2); assert(this->nat > 0); @@ -1136,10 +1125,8 @@ void UnitCell::set_iat2iwt(const int& npol_in) 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++) - { + 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; @@ -1152,16 +1139,13 @@ void UnitCell::set_iat2iwt(const int& npol_in) // Target : meshx // Demand : atoms[].msh //====================== -void UnitCell::cal_meshx() -{ +void UnitCell::cal_meshx() { if (GlobalV::test_pseudo_cell) ModuleBase::TITLE("UnitCell", "cal_meshx"); this->meshx = 0; - for (int it = 0; it < this->ntype; it++) - { + for (int it = 0; it < this->ntype; it++) { const int mesh = this->atoms[it].ncpp.msh; - if (mesh > this->meshx) - { + if (mesh > this->meshx) { this->meshx = mesh; } } @@ -1175,61 +1159,55 @@ void UnitCell::cal_meshx() // atoms[].oc // atoms[].na //========================= -void UnitCell::cal_natomwfc(std::ofstream& log) -{ +void UnitCell::cal_natomwfc(std::ofstream& log) { if (GlobalV::test_pseudo_cell) ModuleBase::TITLE("UnitCell", "cal_natomwfc"); this->natomwfc = 0; - for (int it = 0; it < ntype; it++) - { + for (int it = 0; it < ntype; it++) { //============================ // Use pseudo-atomic orbitals //============================ int tmp = 0; - for (int l = 0; l < atoms[it].ncpp.nchi; l++) - { - if (atoms[it].ncpp.oc[l] >= 0) - { - if (GlobalV::NSPIN == 4) - { - if (atoms[it].ncpp.has_so) - { + for (int l = 0; l < atoms[it].ncpp.nchi; l++) { + if (atoms[it].ncpp.oc[l] >= 0) { + if (GlobalV::NSPIN == 4) { + if (atoms[it].ncpp.has_so) { tmp += 2 * atoms[it].ncpp.lchi[l]; - if (fabs(atoms[it].ncpp.jchi[l] - atoms[it].ncpp.lchi[l] - 0.5) < 1e-6) + if (fabs(atoms[it].ncpp.jchi[l] - atoms[it].ncpp.lchi[l] + - 0.5) + < 1e-6) tmp += 2; - } - else - { + } else { tmp += 2 * (2 * atoms[it].ncpp.lchi[l] + 1); } - } - else + } else tmp += 2 * atoms[it].ncpp.lchi[l] + 1; } } natomwfc += tmp * atoms[it].na; } - ModuleBase::GlobalFunc::OUT(log, "initial pseudo atomic orbital number", natomwfc); + ModuleBase::GlobalFunc::OUT(log, + "initial pseudo atomic orbital number", + natomwfc); return; } // LiuXh add a new function here, // 20180515 -void UnitCell::setup_cell_after_vc(std::ofstream& log) -{ +void UnitCell::setup_cell_after_vc(std::ofstream& log) { ModuleBase::TITLE("UnitCell", "setup_cell_after_vc"); assert(lat0 > 0.0); this->omega = std::abs(latvec.Det()) * this->lat0 * lat0 * lat0; - if (this->omega <= 0) - { + if (this->omega <= 0) { ModuleBase::WARNING_QUIT("setup_cell_after_vc", "omega <= 0 ."); - } - else - { + } else { log << std::endl; ModuleBase::GlobalFunc::OUT(log, "Volume (Bohr^3)", this->omega); - ModuleBase::GlobalFunc::OUT(log, "Volume (A^3)", this->omega * pow(ModuleBase::BOHR_TO_A, 3)); + ModuleBase::GlobalFunc::OUT(log, + "Volume (A^3)", + this->omega + * pow(ModuleBase::BOHR_TO_A, 3)); } lat0_angstrom = lat0 * 0.529177; @@ -1258,11 +1236,9 @@ void UnitCell::setup_cell_after_vc(std::ofstream& log) this->GGT = G * GT; this->invGGT = GGT.Inverse(); - for (int it = 0; it < ntype; it++) - { + for (int it = 0; it < ntype; it++) { Atom* atom = &atoms[it]; - for (int ia = 0; ia < atom->na; ia++) - { + for (int ia = 0; ia < atom->na; ia++) { atom->tau[ia] = atom->taud[ia] * latvec; } } @@ -1272,20 +1248,22 @@ void UnitCell::setup_cell_after_vc(std::ofstream& log) #endif 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); + 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); return; } // check if any atom can be moved -bool UnitCell::if_atoms_can_move() const -{ - for (int it = 0; it < this->ntype; it++) - { +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++) - { + for (int ia = 0; ia < atom->na; ia++) { if (atom->mbl[ia].x || atom->mbl[ia].y || atom->mbl[ia].z) return true; } @@ -1294,11 +1272,9 @@ bool UnitCell::if_atoms_can_move() const } // check if lattice vector can be changed -bool UnitCell::if_cell_can_change() const -{ +bool UnitCell::if_cell_can_change() const { // need to be fixed next - if (this->lc[0] || this->lc[1] || this->lc[2]) - { + if (this->lc[0] || this->lc[1] || this->lc[2]) { return true; } return false; @@ -1308,250 +1284,94 @@ void UnitCell::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) -{ + const std::string& fixed_axes_in) { this->latName = latname_in; this->ntype = ntype_in; this->lmaxmax = lmaxmax_in; this->init_vel = init_vel_in; // pengfei Li add 2018-11-11 - if (fixed_axes_in == "None") - { + if (fixed_axes_in == "None") { this->lc[0] = 1; this->lc[1] = 1; this->lc[2] = 1; - } - else if (fixed_axes_in == "volume") - { + } else if (fixed_axes_in == "volume") { this->lc[0] = 1; this->lc[1] = 1; this->lc[2] = 1; - if (!GlobalV::relax_new) - { + if (!GlobalV::relax_new) { ModuleBase::WARNING_QUIT( "Input", - "there are bugs in the old implementation; set relax_new to be 1 for fixed_volume relaxation"); + "there are bugs in the old implementation; set relax_new to be " + "1 for fixed_volume relaxation"); } - } - else if (fixed_axes_in == "shape") - { - if (!GlobalV::relax_new) - { - ModuleBase::WARNING_QUIT("Input", "set relax_new to be 1 for fixed_shape relaxation"); + } else if (fixed_axes_in == "shape") { + if (!GlobalV::relax_new) { + ModuleBase::WARNING_QUIT( + "Input", + "set relax_new to be 1 for fixed_shape relaxation"); } this->lc[0] = 1; this->lc[1] = 1; this->lc[2] = 1; - } - else if (fixed_axes_in == "a") - { + } else if (fixed_axes_in == "a") { this->lc[0] = 0; this->lc[1] = 1; this->lc[2] = 1; - } - else if (fixed_axes_in == "b") - { + } else if (fixed_axes_in == "b") { this->lc[0] = 1; this->lc[1] = 0; this->lc[2] = 1; - } - else if (fixed_axes_in == "c") - { + } else if (fixed_axes_in == "c") { this->lc[0] = 1; this->lc[1] = 1; this->lc[2] = 0; - } - else if (fixed_axes_in == "ab") - { + } else if (fixed_axes_in == "ab") { this->lc[0] = 0; this->lc[1] = 0; this->lc[2] = 1; - } - else if (fixed_axes_in == "ac") - { + } else if (fixed_axes_in == "ac") { this->lc[0] = 0; this->lc[1] = 1; this->lc[2] = 0; - } - else if (fixed_axes_in == "bc") - { + } else if (fixed_axes_in == "bc") { this->lc[0] = 1; this->lc[1] = 0; this->lc[2] = 0; - } - else if (fixed_axes_in == "abc") - { + } else if (fixed_axes_in == "abc") { this->lc[0] = 0; this->lc[1] = 0; this->lc[2] = 0; - } - else - { - ModuleBase::WARNING_QUIT("Input", "fixed_axes should be None,volume,shape,a,b,c,ab,ac,bc or abc!"); + } else { + ModuleBase::WARNING_QUIT( + "Input", + "fixed_axes should be None,volume,shape,a,b,c,ab,ac,bc or abc!"); } return; } -void UnitCell::check_structure(double factor) -{ - // First we calculate all bond length in the structure, - // and compare with the covalent_bond_length, - // if there has bond length is shorter than covalent_bond_length * factor, - // we think this structure is unreasonable. - const double warning_coef = 0.6; - assert(ntype > 0); - std::stringstream errorlog; - bool all_pass = true; - bool no_warning = true; - for (int it1 = 0; it1 < ntype; it1++) - { - std::string symbol1 = this->atoms[it1].ncpp.psd; - double symbol1_covalent_radius; - if (ModuleBase::CovalentRadius.find(symbol1) != ModuleBase::CovalentRadius.end()) - { - symbol1_covalent_radius = ModuleBase::CovalentRadius.at(symbol1); - } - else - { - std::stringstream mess; - mess << "Notice: symbol '" << symbol1 << "' is not an element symbol!!!! "; - mess << "set the covalent radius to be 0." << std::endl; - GlobalV::ofs_running << mess.str(); - std::cout << mess.str(); - symbol1_covalent_radius = 0.0; - } - - for (int ia1 = 0; ia1 < this->atoms[it1].na; ia1++) - { - double x1 = this->atoms[it1].taud[ia1].x; - double y1 = this->atoms[it1].taud[ia1].y; - double z1 = this->atoms[it1].taud[ia1].z; - - for (int it2 = 0; it2 < ntype; it2++) - { - std::string symbol2 = this->atoms[it2].ncpp.psd; - double symbol2_covalent_radius; - if (ModuleBase::CovalentRadius.find(symbol2) != ModuleBase::CovalentRadius.end()) - { - symbol2_covalent_radius = ModuleBase::CovalentRadius.at(symbol2); - } - else - { - symbol2_covalent_radius = 0.0; - } - - double covalent_length = (symbol1_covalent_radius + symbol2_covalent_radius) / ModuleBase::BOHR_TO_A; - - for (int ia2 = 0; ia2 < this->atoms[it2].na; ia2++) - { - for (int a = -1; a < 2; a++) - { - for (int b = -1; b < 2; b++) - { - for (int c = -1; c < 2; c++) - { - if (it1 > it2) - continue; - else if (it1 == it2 && ia1 > ia2) - continue; - else if (it1 == it2 && ia1 == ia2 && a == 0 && b == 0 && c == 0) - continue; - - double x2 = this->atoms[it2].taud[ia2].x + a; - double y2 = this->atoms[it2].taud[ia2].y + b; - double z2 = this->atoms[it2].taud[ia2].z + c; - - double bond_length - = sqrt(pow((x2 - x1) * this->a1.x + (y2 - y1) * this->a2.x + (z2 - z1) * this->a3.x, - 2) - + pow((x2 - x1) * this->a1.y + (y2 - y1) * this->a2.y - + (z2 - z1) * this->a3.y, - 2) - + pow((x2 - x1) * this->a1.z + (y2 - y1) * this->a2.z - + (z2 - z1) * this->a3.z, - 2)) - * this->lat0; - - if (bond_length < covalent_length * factor - || bond_length < covalent_length * warning_coef) - { - errorlog.setf(std::ios_base::fixed, std::ios_base::floatfield); - errorlog << std::setw(3) << ia1 + 1 << "-th " << std::setw(3) - << this->atoms[it1].label << ", "; - errorlog << std::setw(3) << ia2 + 1 << "-th " << std::setw(3) - << this->atoms[it2].label; - errorlog << " (cell:" << std::setw(2) << a << " " << std::setw(2) << b << " " - << std::setw(2) << c << ")"; - errorlog << ", distance= " << std::setprecision(3) << bond_length << " Bohr ("; - errorlog << bond_length * ModuleBase::BOHR_TO_A << " Angstrom)" << std::endl; - - if (bond_length < covalent_length * factor) - { - all_pass = false; - } - else - { - no_warning = false; - } - } - } // c - } // b - } // a - } // ia2 - } // it2 - } // ia1 - } // it1 - - if (!all_pass || !no_warning) - { - std::stringstream mess; - mess << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl; - mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - mess << "!!! WARNING: Some atoms are too close!!!" << std::endl; - mess << "!!! Please check the nearest-neighbor list in log file." << std::endl; - mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl; - mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - - GlobalV::ofs_running << mess.str() << mess.str() << mess.str() << errorlog.str(); - std::cout << mess.str() << mess.str() << mess.str() << std::endl; - - if (!all_pass) - { - mess.clear(); - mess.str(""); - mess << "If this structure is what you want, you can set 'min_dist_coef'" << std::endl; - mess << "as a smaller value (the current value is " << factor << ") in INPUT file." << std::endl; - GlobalV::ofs_running << mess.str(); - std::cout << mess.str(); - ModuleBase::WARNING_QUIT("Input", "The structure is unreasonable!"); - } - } -} - -void UnitCell::remake_cell() -{ +void UnitCell::remake_cell() { ModuleBase::TITLE("UnitCell", "rmake_cell"); // The idea is as follows: for each type of lattice, first calculate // from current latvec the lattice parameters, then use the parameters // to reconstruct latvec - if (latName == "none") - { - ModuleBase::WARNING_QUIT("UnitCell", "to use fixed_ibrav, latname must be provided"); - } - else if (latName == "sc") // ibrav = 1 + if (latName == "none") { + ModuleBase::WARNING_QUIT( + "UnitCell", + "to use fixed_ibrav, latname must be provided"); + } else if (latName == "sc") // ibrav = 1 { - double celldm = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + pow(latvec.e13, 2)); + double celldm = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + + pow(latvec.e13, 2)); latvec.Zero(); latvec.e11 = latvec.e22 = latvec.e33 = celldm; - } - else if (latName == "fcc") // ibrav = 2 + } else if (latName == "fcc") // ibrav = 2 { - double celldm = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + pow(latvec.e13, 2)) / std::sqrt(2.0); + double celldm = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + + pow(latvec.e13, 2)) + / std::sqrt(2.0); latvec.e11 = -celldm; latvec.e12 = 0.0; @@ -1562,10 +1382,11 @@ void UnitCell::remake_cell() latvec.e31 = -celldm; latvec.e32 = celldm; latvec.e33 = 0.0; - } - else if (latName == "bcc") // ibrav = 3 + } else if (latName == "bcc") // ibrav = 3 { - double celldm = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + pow(latvec.e13, 2)) / std::sqrt(3.0); + double celldm = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + + pow(latvec.e13, 2)) + / std::sqrt(3.0); latvec.e11 = celldm; latvec.e12 = celldm; @@ -1576,11 +1397,12 @@ void UnitCell::remake_cell() latvec.e31 = -celldm; latvec.e32 = -celldm; latvec.e33 = celldm; - } - else if (latName == "hexagonal") // ibrav = 4 + } else if (latName == "hexagonal") // ibrav = 4 { - double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + pow(latvec.e13, 2)); - double celldm3 = std::sqrt(pow(latvec.e31, 2) + pow(latvec.e32, 2) + pow(latvec.e33, 2)); + double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + + pow(latvec.e13, 2)); + double celldm3 = std::sqrt(pow(latvec.e31, 2) + pow(latvec.e32, 2) + + pow(latvec.e33, 2)); double e22 = sqrt(3.0) / 2.0; latvec.e11 = celldm1; @@ -1592,16 +1414,17 @@ void UnitCell::remake_cell() latvec.e31 = 0.0; latvec.e32 = 0.0; latvec.e33 = celldm3; - } - else if (latName == "trigonal") // ibrav = 5 - { - double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + pow(latvec.e13, 2)); - double celldm2 = std::sqrt(pow(latvec.e21, 2) + pow(latvec.e22, 2) + pow(latvec.e23, 2)); - double celldm12 = (latvec.e11 * latvec.e21 + latvec.e12 * latvec.e22 + latvec.e13 * latvec.e23); + } else if (latName == "trigonal") // ibrav = 5 + { + double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + + pow(latvec.e13, 2)); + double celldm2 = std::sqrt(pow(latvec.e21, 2) + pow(latvec.e22, 2) + + pow(latvec.e23, 2)); + double celldm12 = (latvec.e11 * latvec.e21 + latvec.e12 * latvec.e22 + + latvec.e13 * latvec.e23); double cos12 = celldm12 / celldm1 / celldm2; - if (cos12 <= -0.5 || cos12 >= 1.0) - { + if (cos12 <= -0.5 || cos12 >= 1.0) { ModuleBase::WARNING_QUIT("unitcell", "wrong cos12!"); } double t1 = sqrt(1.0 + 2.0 * cos12); @@ -1621,11 +1444,12 @@ void UnitCell::remake_cell() latvec.e31 = -e11; latvec.e32 = e12; latvec.e33 = e13; - } - else if (latName == "st") // ibrav = 6 + } else if (latName == "st") // ibrav = 6 { - double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + pow(latvec.e13, 2)); - double celldm3 = std::sqrt(pow(latvec.e31, 2) + pow(latvec.e32, 2) + pow(latvec.e33, 2)); + double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + + pow(latvec.e13, 2)); + double celldm3 = std::sqrt(pow(latvec.e31, 2) + pow(latvec.e32, 2) + + pow(latvec.e33, 2)); latvec.e11 = celldm1; latvec.e12 = 0.0; latvec.e13 = 0.0; @@ -1635,8 +1459,7 @@ void UnitCell::remake_cell() latvec.e31 = 0.0; latvec.e32 = 0.0; latvec.e33 = celldm3; - } - else if (latName == "bct") // ibrav = 7 + } else if (latName == "bct") // ibrav = 7 { double celldm1 = std::abs(latvec.e11); double celldm2 = std::abs(latvec.e13); @@ -1650,12 +1473,14 @@ void UnitCell::remake_cell() latvec.e31 = -celldm1; latvec.e32 = -celldm1; latvec.e33 = celldm2; - } - else if (latName == "so") // ibrav = 8 + } else if (latName == "so") // ibrav = 8 { - double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + pow(latvec.e13, 2)); - double celldm2 = std::sqrt(pow(latvec.e21, 2) + pow(latvec.e22, 2) + pow(latvec.e23, 2)); - double celldm3 = std::sqrt(pow(latvec.e31, 2) + pow(latvec.e32, 2) + pow(latvec.e33, 2)); + double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + + pow(latvec.e13, 2)); + double celldm2 = std::sqrt(pow(latvec.e21, 2) + pow(latvec.e22, 2) + + pow(latvec.e23, 2)); + double celldm3 = std::sqrt(pow(latvec.e31, 2) + pow(latvec.e32, 2) + + pow(latvec.e33, 2)); latvec.e11 = celldm1; latvec.e12 = 0.0; @@ -1666,8 +1491,7 @@ void UnitCell::remake_cell() latvec.e31 = 0.0; latvec.e32 = 0.0; latvec.e33 = celldm3; - } - else if (latName == "baco") // ibrav = 9 + } else if (latName == "baco") // ibrav = 9 { double celldm1 = std::abs(latvec.e11); double celldm2 = std::abs(latvec.e22); @@ -1682,8 +1506,7 @@ void UnitCell::remake_cell() latvec.e31 = 0.0; latvec.e32 = 0.0; latvec.e33 = celldm3; - } - else if (latName == "fco") // ibrav = 10 + } else if (latName == "fco") // ibrav = 10 { double celldm1 = std::abs(latvec.e11); double celldm2 = std::abs(latvec.e22); @@ -1698,8 +1521,7 @@ void UnitCell::remake_cell() latvec.e31 = 0.0; latvec.e32 = celldm2; latvec.e33 = celldm3; - } - else if (latName == "bco") // ibrav = 11 + } else if (latName == "bco") // ibrav = 11 { double celldm1 = std::abs(latvec.e11); double celldm2 = std::abs(latvec.e12); @@ -1714,13 +1536,16 @@ void UnitCell::remake_cell() latvec.e31 = -celldm1; latvec.e32 = -celldm2; latvec.e33 = celldm3; - } - else if (latName == "sm") // ibrav = 12 - { - double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + pow(latvec.e13, 2)); - double celldm2 = std::sqrt(pow(latvec.e21, 2) + pow(latvec.e22, 2) + pow(latvec.e23, 2)); - double celldm3 = std::sqrt(pow(latvec.e31, 2) + pow(latvec.e32, 2) + pow(latvec.e33, 2)); - double celldm12 = (latvec.e11 * latvec.e21 + latvec.e12 * latvec.e22 + latvec.e13 * latvec.e23); + } else if (latName == "sm") // ibrav = 12 + { + double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + + pow(latvec.e13, 2)); + double celldm2 = std::sqrt(pow(latvec.e21, 2) + pow(latvec.e22, 2) + + pow(latvec.e23, 2)); + double celldm3 = std::sqrt(pow(latvec.e31, 2) + pow(latvec.e32, 2) + + pow(latvec.e33, 2)); + double celldm12 = (latvec.e11 * latvec.e21 + latvec.e12 * latvec.e22 + + latvec.e13 * latvec.e23); double cos12 = celldm12 / celldm1 / celldm2; double e21 = celldm2 * cos12; @@ -1735,16 +1560,15 @@ void UnitCell::remake_cell() latvec.e31 = 0.0; latvec.e32 = 0.0; latvec.e33 = celldm3; - } - else if (latName == "bacm") // ibrav = 13 + } else if (latName == "bacm") // ibrav = 13 { double celldm1 = std::abs(latvec.e11); - double celldm2 = std::sqrt(pow(latvec.e21, 2) + pow(latvec.e22, 2) + pow(latvec.e23, 2)); + double celldm2 = std::sqrt(pow(latvec.e21, 2) + pow(latvec.e22, 2) + + pow(latvec.e23, 2)); double celldm3 = std::abs(latvec.e13); double cos12 = latvec.e21 / celldm2; - if (cos12 >= 1.0) - { + if (cos12 >= 1.0) { ModuleBase::WARNING_QUIT("unitcell", "wrong cos12!"); } @@ -1760,22 +1584,26 @@ void UnitCell::remake_cell() latvec.e31 = celldm1; latvec.e32 = 0.0; latvec.e33 = celldm3; - } - else if (latName == "triclinic") // ibrav = 14 - { - double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + pow(latvec.e13, 2)); - double celldm2 = std::sqrt(pow(latvec.e21, 2) + pow(latvec.e22, 2) + pow(latvec.e23, 2)); - double celldm3 = std::sqrt(pow(latvec.e31, 2) + pow(latvec.e32, 2) + pow(latvec.e33, 2)); - double celldm12 = (latvec.e11 * latvec.e21 + latvec.e12 * latvec.e22 + latvec.e13 * latvec.e23); + } else if (latName == "triclinic") // ibrav = 14 + { + double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2) + + pow(latvec.e13, 2)); + double celldm2 = std::sqrt(pow(latvec.e21, 2) + pow(latvec.e22, 2) + + pow(latvec.e23, 2)); + double celldm3 = std::sqrt(pow(latvec.e31, 2) + pow(latvec.e32, 2) + + pow(latvec.e33, 2)); + double celldm12 = (latvec.e11 * latvec.e21 + latvec.e12 * latvec.e22 + + latvec.e13 * latvec.e23); double cos12 = celldm12 / celldm1 / celldm2; - double celldm13 = (latvec.e11 * latvec.e31 + latvec.e12 * latvec.e32 + latvec.e13 * latvec.e33); + double celldm13 = (latvec.e11 * latvec.e31 + latvec.e12 * latvec.e32 + + latvec.e13 * latvec.e33); double cos13 = celldm13 / celldm1 / celldm3; - double celldm23 = (latvec.e21 * latvec.e31 + latvec.e22 * latvec.e32 + latvec.e23 * latvec.e33); + double celldm23 = (latvec.e21 * latvec.e31 + latvec.e22 * latvec.e32 + + latvec.e23 * latvec.e33); double cos23 = celldm23 / celldm2 / celldm3; double sin12 = std::sqrt(1.0 - cos12 * cos12); - if (cos12 >= 1.0) - { + if (cos12 >= 1.0) { ModuleBase::WARNING_QUIT("unitcell", "wrong cos12!"); } @@ -1787,61 +1615,69 @@ void UnitCell::remake_cell() latvec.e23 = 0.0; latvec.e31 = celldm3 * cos13; latvec.e32 = celldm3 * (cos23 - cos13 * cos12) / sin12; - double term = 1.0 + 2.0 * cos12 * cos13 * cos23 - cos12 * cos12 - cos13 * cos13 - cos23 * cos23; + double term = 1.0 + 2.0 * cos12 * cos13 * cos23 - cos12 * cos12 + - cos13 * cos13 - cos23 * cos23; term = sqrt(term) / sin12; latvec.e33 = celldm3 * term; - } - else - { + } else { std::cout << "latname is : " << latName << std::endl; - ModuleBase::WARNING_QUIT("UnitCell::read_atom_species", "latname not supported!"); + ModuleBase::WARNING_QUIT("UnitCell::read_atom_species", + "latname not supported!"); } } -void UnitCell::cal_nelec(double& nelec) -{ +void UnitCell::cal_nelec(double& nelec) { ModuleBase::TITLE("UnitCell", "cal_nelec"); GlobalV::ofs_running << "\n SETUP THE ELECTRONS NUMBER" << std::endl; - if (nelec == 0) - { - if (GlobalV::use_paw) - { + if (nelec == 0) { + if (GlobalV::use_paw) { #ifdef USE_PAW - for (int it = 0; it < this->ntype; it++) - { + for (int it = 0; it < this->ntype; it++) { std::stringstream ss1, ss2; - ss1 << " electron number of element " << GlobalC::paw_cell.get_zat(it) << std::endl; - const int nelec_it = GlobalC::paw_cell.get_val(it) * this->atoms[it].na; + ss1 << " electron number of element " + << GlobalC::paw_cell.get_zat(it) << std::endl; + const int nelec_it + = GlobalC::paw_cell.get_val(it) * this->atoms[it].na; nelec += nelec_it; - ss2 << "total electron number of element " << GlobalC::paw_cell.get_zat(it); - - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, ss1.str(), GlobalC::paw_cell.get_val(it)); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, ss2.str(), nelec_it); + ss2 << "total electron number of element " + << GlobalC::paw_cell.get_zat(it); + + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, + ss1.str(), + GlobalC::paw_cell.get_val(it)); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, + ss2.str(), + nelec_it); } #endif - } - else - { - for (int it = 0; it < this->ntype; it++) - { + } else { + for (int it = 0; it < this->ntype; it++) { std::stringstream ss1, ss2; ss1 << "electron number of element " << this->atoms[it].label; - const int nelec_it = this->atoms[it].ncpp.zv * this->atoms[it].na; + const int nelec_it + = this->atoms[it].ncpp.zv * this->atoms[it].na; nelec += nelec_it; - ss2 << "total electron number of element " << this->atoms[it].label; - - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, ss1.str(), this->atoms[it].ncpp.zv); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, ss2.str(), nelec_it); + ss2 << "total electron number of element " + << this->atoms[it].label; + + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, + ss1.str(), + this->atoms[it].ncpp.zv); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, + ss2.str(), + nelec_it); } - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "AUTOSET number of electrons: ", nelec); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, + "AUTOSET number of electrons: ", + nelec); } } - if (GlobalV::nelec_delta != 0) - { + if (GlobalV::nelec_delta != 0) { ModuleBase::GlobalFunc::OUT( GlobalV::ofs_running, - "nelec_delta is NOT zero, please make sure you know what you are doing! nelec_delta: ", + "nelec_delta is NOT zero, please make sure you know what you are " + "doing! nelec_delta: ", GlobalV::nelec_delta); nelec += GlobalV::nelec_delta; ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nelec now: ", nelec); @@ -1849,59 +1685,66 @@ void UnitCell::cal_nelec(double& nelec) return; } -void UnitCell::compare_atom_labels(std::string label1, std::string label2) -{ - if (label1 != label2) //'!( "Ag" == "Ag" || "47" == "47" || "Silver" == Silver" )' +void UnitCell::compare_atom_labels(std::string label1, 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" )' + 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] == '_')) - { + for (int ip = 0; ip < label1.length(); ip++) { + if (!(isdigit(label1[ip]) || label1[ip] == '_')) { stru_label += label1[ip]; - } - else - { + } else { break; } } stru_label[0] = toupper(stru_label[0]); - for (int ip = 0; ip < label2.length(); ip++) - { - if (!(isdigit(label2[ip]) || label2[ip] == '_')) - { + for (int ip = 0; ip < label2.length(); ip++) { + if (!(isdigit(label2[ip]) || label2[ip] == '_')) { psuedo_label += label2[ip]; - } - else - { + } 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" )' + 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 "; + 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); + atom_label_in_orbtial + label1 + + mismatch_with_pseudo + label2); } } } diff --git a/source/module_cell/unitcell.h b/source/module_cell/unitcell.h index 9abadc83cf..ff25373930 100644 --- a/source/module_cell/unitcell.h +++ b/source/module_cell/unitcell.h @@ -13,8 +13,7 @@ #endif // provide the basic information about unitcell. -class UnitCell -{ +class UnitCell { public: Atom* atoms; @@ -58,35 +57,31 @@ class UnitCell ModuleSymmetry::Symmetry symm; // ======================================================== - // 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 + // 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: // 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 - { + inline Tiait + itiaiw2iwt(const Tiait& it, const Tiait& ia, const Tiait& iw) const { return Tiait(this->iat2iwt[this->itia2iat(it, ia)] + iw); } // initialize iat2iwt void set_iat2iwt(const int& npol_in); // get iat2iwt - inline const int* get_iat2iwt() const - { - return iat2iwt.data(); - } + inline const int* get_iat2iwt() const { return iat2iwt.data(); } // get npol - inline const int& get_npol() const - { - return 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 ----------------- + 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: //======================================================== @@ -94,10 +89,8 @@ class UnitCell // return true if the last out is reset //======================================================== template - inline bool iat2iait(const Tiat iat, Tiait* ia, Tiait* it) const - { - if (iat >= nat) - { + inline bool iat2iait(const Tiat iat, Tiait* ia, Tiait* it) const { + if (iat >= nat) { *ia = 0; *it = ntype; return false; @@ -108,8 +101,11 @@ class UnitCell } template - inline bool ijat2iaitjajt(const Tiat ijat, Tiait* ia, Tiait* it, Tiait* ja, Tiait* jt) const - { + inline bool ijat2iaitjajt(const Tiat ijat, + Tiait* ia, + Tiait* it, + Tiait* ja, + Tiait* jt) const { Tiat iat = ijat / nat; Tiat jat = ijat % nat; iat2iait(iat, ia, it); @@ -118,10 +114,8 @@ class UnitCell } template - inline bool step_it(Tiait* it) const - { - if (++(*it) >= ntype) - { + inline bool step_it(Tiait* it) const { + if (++(*it) >= ntype) { *it = 0; return true; } @@ -129,10 +123,8 @@ class UnitCell } template - inline bool step_ia(const Tiait it, Tiait* ia) const - { - if (++(*ia) >= atoms[it].na) - { + inline bool step_ia(const Tiait it, Tiait* ia) const { + if (++(*ia) >= atoms[it].na) { *ia = 0; return true; } @@ -140,37 +132,34 @@ class UnitCell } template - inline bool step_iait(Tiait* ia, Tiait* it) const - { - if (step_ia(*it, ia)) - { + 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)) - { + 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; } // get tau for atom iat - inline const ModuleBase::Vector3& get_tau(const int& iat) const - { + inline const ModuleBase::Vector3& get_tau(const int& iat) const { return atoms[iat2it[iat]].tau[iat2ia[iat]]; } // calculate vector between two atoms with R cell - inline const ModuleBase::Vector3 cal_dtau(const int& iat1, - const int& iat2, - const ModuleBase::Vector3& R) const - { - return get_tau(iat2) + double(R.x) * a1 + double(R.y) * a2 + double(R.z) * a3 - get_tau(iat1); + inline const ModuleBase::Vector3 + cal_dtau(const int& iat1, + const int& iat2, + const ModuleBase::Vector3& R) const { + return get_tau(iat2) + double(R.x) * a1 + double(R.y) * a2 + + double(R.z) * a3 - get_tau(iat1); } // LiuXh add 20180515 @@ -180,8 +169,10 @@ class UnitCell ModuleBase::Matrix3 invGGT0; // I'm doing a bad thing here! Will change later - bool ionic_position_updated = false; // whether the ionic position has been updated - bool cell_parameter_updated = false; // whether the cell parameters are updated + bool ionic_position_updated + = false; // whether the ionic position has been updated + bool cell_parameter_updated + = false; // whether the cell parameters are updated //============================================================ // meshx : max number of mesh point in pseudopotential file @@ -218,7 +209,7 @@ class UnitCell void update_vel(const ModuleBase::Vector3* vel_in); void periodic_boundary_adjustment(); void bcast_atoms_tau(); - bool judge_big_cell(void) const; + bool judge_big_cell() const; void update_stress(ModuleBase::matrix& scs); // updates stress void update_force(ModuleBase::matrix& fcs); // updates force in Atom @@ -226,21 +217,24 @@ class UnitCell double* atom_mass; std::string* atom_label; std::string* pseudo_fn; - std::string* pseudo_type; // pseudopotential types for each elements, sunliang added 2022-09-15. - std::string* orbital_fn; // filenames of orbitals, liuyu add 2022-10-19 - std::string descriptor_file; // filenames of descriptor_file, liuyu add 2023-04-06 + std::string* pseudo_type; // pseudopotential types for each elements, + // sunliang added 2022-09-15. + std::string* orbital_fn; // filenames of orbitals, liuyu add 2022-10-19 + std::string + descriptor_file; // filenames of descriptor_file, liuyu add 2023-04-06 #ifdef __MPI - void bcast_unitcell(void); - void bcast_unitcell2(void); + void bcast_unitcell(); + void bcast_unitcell2(); #endif - void set_iat2itia(void); + void set_iat2itia(); void setup_cell(const std::string& fn, std::ofstream& log); #ifdef __LCAO - InfoNonlocal infoNL; // store nonlocal information of lcao, added by zhengdy 2021-09-07 + InfoNonlocal infoNL; // store nonlocal information of lcao, added by zhengdy + // 2021-09-07 #endif /// @brief read number of numerical orbitals for each angular momentum @@ -248,25 +242,32 @@ class UnitCell /// @param orb_file orbital filename /// @param ofs_running ofstream /// @param atom Atom instance stored in UnitCell - void read_orb_file(int it, std::string& orb_file, std::ofstream& ofs_running, Atom* atom); - int read_atom_species(std::ifstream& ifa, - std::ofstream& ofs_running); // read in the atom information for each type of atom - bool read_atom_positions(std::ifstream& ifpos, - std::ofstream& ofs_running, - std::ofstream& ofs_warning); // read in atomic positions + void read_orb_file(int it, + std::string& orb_file, + std::ofstream& ofs_running, + Atom* atom); + int read_atom_species( + std::ifstream& ifa, + std::ofstream& + ofs_running); // read in the atom information for each type of atom + bool read_atom_positions( + std::ifstream& ifpos, + std::ofstream& ofs_running, + std::ofstream& ofs_warning); // read in atomic positions void read_pseudo(std::ofstream& ofs); int find_type(const std::string& label); - void print_tau(void) const; + void print_tau() const; /** - * @brief UnitCell class is too heavy, this function would be moved elsewhere. Print STRU file respect to given - * setting + * @brief UnitCell class is too heavy, this function would be moved + * elsewhere. Print STRU file respect to given setting * * @param fn STRU file name * @param nspin GlobalV::NSPIN feed in * @param direct true for direct coords, false for cartesian coords * @param vol true for printing velocities - * @param magmom true for printing Mulliken population analysis produced magmom + * @param magmom true for printing Mulliken population analysis produced + * magmom * @param orb true for printing NUMERICAL_ORBITAL section * @param dpks_desc true for printing NUMERICAL_DESCRIPTOR section * @param iproc GlobalV::MY_RANK feed in @@ -279,7 +280,7 @@ class UnitCell const bool& orb = false, const bool& dpks_desc = false, const int& iproc = 0) const; - void check_dtau(void); + void check_dtau(); void setup_cell_after_vc(std::ofstream& log); // LiuXh add 20180515 // for constrained vc-relaxation where type of lattice @@ -298,7 +299,7 @@ class UnitCell void cal_meshx(); void cal_natomwfc(std::ofstream& log); void print_unitcell_pseudo(const std::string& fn); - bool check_tau(void) const; // mohan add 2011-03-03 + bool check_tau() const; // mohan add 2011-03-03 bool if_atoms_can_move() const; bool if_cell_can_change() const; void setup(const std::string& latname_in, @@ -307,26 +308,30 @@ class UnitCell const bool& init_vel_in, const std::string& fixed_axes_in); - void check_structure(double factor); - - /// @brief calculate the total number of electrons in system (GlobalV::nelec) + /// @brief calculate the total number of electrons in system + /// (GlobalV::nelec) void cal_nelec(double& nelec); - /// @brief check consistency between two atom labels from STRU and pseudo or orb file + /// @brief check consistency between two atom labels from STRU and pseudo or + /// orb file void compare_atom_labels(std::string label1, std::string label2); /// @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 + /// @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 + /// @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 + /// 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 + /// @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 + /// @brief get lnchiCounts, which is a vector of element type with the + /// l:nchi vector std::vector> get_lnchiCounts() const; }; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_HS_arrays.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_HS_arrays.hpp index 58f9c57dd9..c325fbf6e1 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_HS_arrays.hpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_HS_arrays.hpp @@ -6,8 +6,7 @@ #include #include -class LCAO_HS_Arrays -{ +class LCAO_HS_Arrays { public: LCAO_HS_Arrays(){}; ~LCAO_HS_Arrays(){}; @@ -23,17 +22,45 @@ class LCAO_HS_Arrays //------------------------------ std::vector Hloc_fixedR; - std::map, std::map>> dHRx_sparse[2]; - std::map, std::map>> dHRy_sparse[2]; - std::map, std::map>> dHRz_sparse[2]; + // For HR_sparse[2], when nspin=1, only 0 is valid, when nspin=2, 0 means + // spin up, 1 means spin down + std::map, + std::map>> + HR_sparse[2]; + std::map, + std::map>> + SR_sparse; + std::map, + std::map>> + TR_sparse; + + std::map, + std::map>> + dHRx_sparse[2]; + std::map, + std::map>> + dHRy_sparse[2]; + std::map, + std::map>> + dHRz_sparse[2]; // For nspin = 4 - std::map, std::map>>> HR_soc_sparse; - std::map, std::map>>> SR_soc_sparse; + std::map, + std::map>>> + HR_soc_sparse; + std::map, + std::map>>> + SR_soc_sparse; - std::map, std::map>>> dHRx_soc_sparse; - std::map, std::map>>> dHRy_soc_sparse; - std::map, std::map>>> dHRz_soc_sparse; + std::map, + std::map>>> + dHRx_soc_sparse; + std::map, + std::map>>> + dHRy_soc_sparse; + std::map, + std::map>>> + dHRz_soc_sparse; }; #endif diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.hpp index a1a79f1d5f..741cfbbc0c 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.hpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.hpp @@ -31,102 +31,99 @@ void sparse_format::cal_HR_exx( const int& current_spin, const double& sparse_threshold, const int (&nmp)[3], - const std::vector>, RI::Tensor>>>& Hexxs) -{ + const std::vector>, + RI::Tensor>>>& Hexxs) { ModuleBase::TITLE("sparse_format", "cal_HR_exx"); ModuleBase::timer::tick("sparse_format", "cal_HR_exx"); const Tdata frac = GlobalC::exx_info.info_global.hybrid_alpha; std::map> atoms_pos; - for (int iat = 0; iat < GlobalC::ucell.nat; ++iat) - { + for (int iat = 0; iat < GlobalC::ucell.nat; ++iat) { atoms_pos[iat] = RI_Util::Vector3_to_array3( - GlobalC::ucell.atoms[GlobalC::ucell.iat2it[iat]].tau[GlobalC::ucell.iat2ia[iat]]); + GlobalC::ucell.atoms[GlobalC::ucell.iat2it[iat]] + .tau[GlobalC::ucell.iat2ia[iat]]); } - const std::array, 3> latvec = {RI_Util::Vector3_to_array3(GlobalC::ucell.a1), - RI_Util::Vector3_to_array3(GlobalC::ucell.a2), - RI_Util::Vector3_to_array3(GlobalC::ucell.a3)}; + const std::array, 3> latvec + = {RI_Util::Vector3_to_array3(GlobalC::ucell.a1), + RI_Util::Vector3_to_array3(GlobalC::ucell.a2), + RI_Util::Vector3_to_array3(GlobalC::ucell.a3)}; const std::array Rs_period = {nmp[0], nmp[1], nmp[2]}; RI::Cell_Nearest cell_nearest; cell_nearest.init(atoms_pos, latvec, Rs_period); - const std::vector is_list - = (GlobalV::NSPIN != 4) ? std::vector{current_spin} : std::vector{0, 1, 2, 3}; + const std::vector is_list = (GlobalV::NSPIN != 4) + ? std::vector{current_spin} + : std::vector{0, 1, 2, 3}; - for (const int is: is_list) - { + for (const int is: is_list) { int is0_b = 0; int is1_b = 0; std::tie(is0_b, is1_b) = RI_2D_Comm::split_is_block(is); - if (Hexxs.empty()) - { + if (Hexxs.empty()) { break; } - for (const auto& HexxA: Hexxs[is]) - { + for (const auto& HexxA: Hexxs[is]) { const int iat0 = HexxA.first; - for (const auto& HexxB: HexxA.second) - { + for (const auto& HexxB: HexxA.second) { const int iat1 = HexxB.first.first; const Abfs::Vector3_Order R = RI_Util::array3_to_Vector3( - cell_nearest.get_cell_nearest_discrete(iat0, iat1, HexxB.first.second)); + cell_nearest.get_cell_nearest_discrete(iat0, + iat1, + HexxB.first.second)); lm.all_R_coor.insert(R); const RI::Tensor& Hexx = HexxB.second; - for (size_t iw0 = 0; iw0 < Hexx.shape[0]; ++iw0) - { + for (size_t iw0 = 0; iw0 < Hexx.shape[0]; ++iw0) { const int iwt0 = RI_2D_Comm::get_iwt(iat0, iw0, is0_b); const int iwt0_local = lm.ParaV->global2local_row(iwt0); - if (iwt0_local < 0) - { + if (iwt0_local < 0) { continue; } - for (size_t iw1 = 0; iw1 < Hexx.shape[1]; ++iw1) - { + for (size_t iw1 = 0; iw1 < Hexx.shape[1]; ++iw1) { const int iwt1 = RI_2D_Comm::get_iwt(iat1, iw1, is1_b); const int iwt1_local = lm.ParaV->global2local_col(iwt1); - if (iwt1_local < 0) - { + if (iwt1_local < 0) { continue; } - if (std::abs(Hexx(iw0, iw1)) > sparse_threshold) - { - if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 2) - { - auto& HR_sparse_ptr = lm.HR_sparse[current_spin][R][iwt0]; + if (std::abs(Hexx(iw0, iw1)) > sparse_threshold) { + if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 2) { + auto& HR_sparse_ptr + = HS_Arrays + .HR_sparse[current_spin][R][iwt0]; double& HR_sparse = HR_sparse_ptr[iwt1]; - HR_sparse += RI::Global_Func::convert(frac * Hexx(iw0, iw1)); - if (std::abs(HR_sparse) <= sparse_threshold) - { + HR_sparse += RI::Global_Func::convert( + frac * Hexx(iw0, iw1)); + if (std::abs(HR_sparse) <= sparse_threshold) { HR_sparse_ptr.erase(iwt1); } - } - else if (GlobalV::NSPIN == 4) - { - auto& HR_sparse_ptr = HS_Arrays.HR_soc_sparse[R][iwt0]; - std::complex& HR_sparse = HR_sparse_ptr[iwt1]; - HR_sparse += RI::Global_Func::convert>(frac * Hexx(iw0, iw1)); - if (std::abs(HR_sparse) <= sparse_threshold) - { + } else if (GlobalV::NSPIN == 4) { + auto& HR_sparse_ptr + = HS_Arrays.HR_soc_sparse[R][iwt0]; + std::complex& HR_sparse + = HR_sparse_ptr[iwt1]; + HR_sparse += RI::Global_Func::convert< + std::complex>(frac + * Hexx(iw0, iw1)); + if (std::abs(HR_sparse) <= sparse_threshold) { HR_sparse_ptr.erase(iwt1); } - } - else - { - throw std::invalid_argument(std::string(__FILE__) + " line " - + std::to_string(__LINE__)); + } else { + throw std::invalid_argument( + std::string(__FILE__) + " line " + + std::to_string(__LINE__)); } } } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.cpp index 8ad1fc8405..d19ce2c1a3 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.cpp @@ -7,54 +7,46 @@ #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" #endif -LCAO_Matrix::LCAO_Matrix() -{ -} +LCAO_Matrix::LCAO_Matrix() {} -LCAO_Matrix::~LCAO_Matrix() -{ -} +LCAO_Matrix::~LCAO_Matrix() {} -void LCAO_Matrix::divide_HS_in_frag(const bool isGamma, Parallel_Orbitals& pv, const int& nks) -{ +void LCAO_Matrix::divide_HS_in_frag(const bool isGamma, + Parallel_Orbitals& pv, + const int& nks) { ModuleBase::TITLE("LCAO_Matrix", "divide_HS_in_frag"); //(1), (2): set up matrix division have been moved into ORB_control // just pass `ParaV` as pointer is enough this->ParaV = &pv; // (3) allocate for S, H_fixed, H, and S_diag - if (isGamma) - { + if (isGamma) { allocate_HS_gamma(this->ParaV->nloc); - } - else - { + } else { allocate_HS_k(this->ParaV->nloc); } #ifdef __DEEPKS // wenfei 2021-12-19 // preparation for DeePKS - if (GlobalV::deepks_out_labels || GlobalV::deepks_scf) - { + if (GlobalV::deepks_out_labels || GlobalV::deepks_scf) { // allocate relevant data structures for calculating descriptors std::vector na; na.resize(GlobalC::ucell.ntype); - for (int it = 0; it < GlobalC::ucell.ntype; it++) - { + for (int it = 0; it < GlobalC::ucell.ntype; it++) { na[it] = GlobalC::ucell.atoms[it].na; } - GlobalC::ld.init(GlobalC::ORB, GlobalC::ucell.nat, GlobalC::ucell.ntype, pv, na); + GlobalC::ld.init(GlobalC::ORB, + GlobalC::ucell.nat, + GlobalC::ucell.ntype, + pv, + na); - if (GlobalV::deepks_scf) - { - if (isGamma) - { + if (GlobalV::deepks_scf) { + if (isGamma) { GlobalC::ld.allocate_V_delta(GlobalC::ucell.nat); - } - else - { + } else { GlobalC::ld.allocate_V_delta(GlobalC::ucell.nat, nks); } } @@ -63,14 +55,12 @@ void LCAO_Matrix::divide_HS_in_frag(const bool isGamma, Parallel_Orbitals& pv, c return; } -void LCAO_Matrix::allocate_HS_gamma(const long& nloc) -{ +void LCAO_Matrix::allocate_HS_gamma(const long& nloc) { ModuleBase::TITLE("LCAO_Matrix", "allocate_HS_gamma"); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nloc", nloc); - if (nloc == 0) - { + if (nloc == 0) { return; } @@ -89,14 +79,12 @@ void LCAO_Matrix::allocate_HS_gamma(const long& nloc) return; } -void LCAO_Matrix::allocate_HS_k(const long& nloc) -{ +void LCAO_Matrix::allocate_HS_k(const long& nloc) { ModuleBase::TITLE("LCAO_Matrix", "allocate_HS_k"); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nloc", nloc); - if (nloc == 0) - { + if (nloc == 0) { return; // mohan fix bug 2012-05-25 } @@ -114,18 +102,18 @@ void LCAO_Matrix::allocate_HS_k(const long& nloc) return; } -void LCAO_Matrix::set_HSgamma(const int& iw1_all, const int& iw2_all, const double& v, double* HSloc) -{ +void LCAO_Matrix::set_HSgamma(const int& iw1_all, + const int& iw2_all, + const double& v, + double* HSloc) { LCAO_Matrix::set_mat2d(iw1_all, iw2_all, v, *this->ParaV, HSloc); return; } -void LCAO_Matrix::zeros_HSgamma(const char& mtype) -{ +void LCAO_Matrix::zeros_HSgamma(const char& mtype) { auto zeros_HSgamma_ker = [&](int num_threads, int thread_id) { long long beg, len; - if (mtype == 'S') - { + if (mtype == 'S') { ModuleBase::BLOCK_TASK_DIST_1D(num_threads, thread_id, (long long)this->Sloc.size(), @@ -134,9 +122,7 @@ void LCAO_Matrix::zeros_HSgamma(const char& mtype) len); ModuleBase::GlobalFunc::ZEROS(this->Sloc.data() + beg, len); - } - else if (mtype == 'T') - { + } else if (mtype == 'T') { ModuleBase::BLOCK_TASK_DIST_1D(num_threads, thread_id, (long long)this->Hloc_fixed.size(), @@ -145,9 +131,7 @@ void LCAO_Matrix::zeros_HSgamma(const char& mtype) len); ModuleBase::GlobalFunc::ZEROS(this->Hloc_fixed.data() + beg, len); - } - else if (mtype == 'H') - { + } else if (mtype == 'H') { ModuleBase::BLOCK_TASK_DIST_1D(num_threads, thread_id, (long long)this->Hloc.size(), @@ -162,12 +146,10 @@ void LCAO_Matrix::zeros_HSgamma(const char& mtype) return; } -void LCAO_Matrix::zeros_HSk(const char& mtype) -{ +void LCAO_Matrix::zeros_HSk(const char& mtype) { auto zeros_HSk_ker = [&](int num_threads, int thread_id) { long long beg, len; - if (mtype == 'S') - { + if (mtype == 'S') { ModuleBase::BLOCK_TASK_DIST_1D(num_threads, thread_id, (long long)this->Sloc2.size(), @@ -175,9 +157,7 @@ void LCAO_Matrix::zeros_HSk(const char& mtype) beg, len); ModuleBase::GlobalFunc::ZEROS(this->Sloc2.data() + beg, len); - } - else if (mtype == 'T') - { + } else if (mtype == 'T') { ModuleBase::BLOCK_TASK_DIST_1D(num_threads, thread_id, (long long)this->Hloc_fixed2.size(), @@ -185,9 +165,7 @@ void LCAO_Matrix::zeros_HSk(const char& mtype) beg, len); ModuleBase::GlobalFunc::ZEROS(this->Hloc_fixed2.data() + beg, len); - } - else if (mtype == 'H') - { + } else if (mtype == 'H') { ModuleBase::BLOCK_TASK_DIST_1D(num_threads, thread_id, (long long)this->Hloc2.size(), @@ -202,31 +180,26 @@ void LCAO_Matrix::zeros_HSk(const char& mtype) } // becareful! Update Hloc, we add new members to it. -void LCAO_Matrix::update_Hloc() -{ +void LCAO_Matrix::update_Hloc() { ModuleBase::TITLE("LCAO_Matrix", "update_Hloc"); #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for (long i = 0; i < this->ParaV->nloc; i++) - { + for (long i = 0; i < this->ParaV->nloc; i++) { Hloc[i] += Hloc_fixed[i]; } return; } -void LCAO_Matrix::update_Hloc2(const int& ik) -{ +void LCAO_Matrix::update_Hloc2(const int& ik) { ModuleBase::TITLE("LCAO_Matrix", "update_Hloc2"); #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for (long i = 0; i < this->ParaV->nloc; i++) - { + for (long i = 0; i < this->ParaV->nloc; i++) { Hloc2[i] += Hloc_fixed2[i]; #ifdef __DEEPKS - if (GlobalV::deepks_scf) - { + if (GlobalV::deepks_scf) { Hloc2[i] += GlobalC::ld.H_V_delta_k[ik][i]; } #endif @@ -235,24 +208,24 @@ void LCAO_Matrix::update_Hloc2(const int& ik) return; } -void LCAO_Matrix::output_HSk(const char& mtype, std::string& fn) -{ +void LCAO_Matrix::output_HSk(const char& mtype, std::string& fn) { ModuleBase::TITLE("LCAO_Matrix", "output_HSk"); std::stringstream ss; ss << GlobalV::global_out_dir << fn; std::ofstream ofs(ss.str().c_str()); ofs << GlobalV::NLOCAL << std::endl; - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - for (int j = 0; j < GlobalV::NLOCAL; j++) - { + for (int i = 0; i < GlobalV::NLOCAL; i++) { + for (int j = 0; j < GlobalV::NLOCAL; j++) { const int index = i * GlobalV::NLOCAL + j; if (mtype == 'S') - ofs << Sloc2[index].real() << " " << Sloc2[index].imag() << std::endl; + ofs << Sloc2[index].real() << " " << Sloc2[index].imag() + << std::endl; else if (mtype == 'T') - ofs << Hloc_fixed2[index].real() << " " << Hloc_fixed2[index].imag() << std::endl; + ofs << Hloc_fixed2[index].real() << " " + << Hloc_fixed2[index].imag() << std::endl; else if (mtype == 'H') - ofs << Hloc2[index].real() << " " << Hloc2[index].imag() << std::endl; + ofs << Hloc2[index].real() << " " << Hloc2[index].imag() + << std::endl; } } ofs.close(); @@ -264,21 +237,17 @@ void LCAO_Matrix::set_HR_tr(const int& Rx, const int& Rz, const int& iw1_all, const int& iw2_all, - const double& v) -{ + const double& v) { const int ir = this->ParaV->global2local_row(iw1_all); const int ic = this->ParaV->global2local_col(iw2_all); // std::cout<<"ir: "<ParaV->nrow + ir; // std::cout<<"index: "<ParaV->ncol + ic; // std::cout<<"index: "<& v) -{ + const std::complex& v) { const int ir = this->ParaV->global2local_row(iw1_all); const int ic = this->ParaV->global2local_col(iw2_all); long index = 0; - if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER()) - { + if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER()) { index = ic * this->ParaV->nrow + ir; - } - else - { + } else { index = ir * this->ParaV->ncol + ic; } assert(index < this->ParaV->nloc); HR_tr_soc[Rx][Ry][Rz][index] = Hloc_fixedR_tr_soc[Rx][Ry][Rz][index] + v; - return; -} - -void LCAO_Matrix::destroy_HS_R_sparse(LCAO_HS_Arrays& HS_Arrays) -{ - ModuleBase::TITLE("LCAO_Matrix", "destroy_HS_R_sparse"); - - if (GlobalV::NSPIN != 4) - { - std::map, std::map>> empty_HR_sparse_up; - std::map, std::map>> empty_HR_sparse_down; - std::map, std::map>> empty_SR_sparse; - HR_sparse[0].swap(empty_HR_sparse_up); - HR_sparse[1].swap(empty_HR_sparse_down); - SR_sparse.swap(empty_SR_sparse); - } - else - { - std::map, std::map>>> - empty_HR_soc_sparse; - std::map, std::map>>> - empty_SR_soc_sparse; - HS_Arrays.HR_soc_sparse.swap(empty_HR_soc_sparse); - HS_Arrays.SR_soc_sparse.swap(empty_SR_soc_sparse); - } - - // 'all_R_coor' has a small memory requirement and does not need to be deleted. - // std::set> empty_all_R_coor; - // all_R_coor.swap(empty_all_R_coor); - - return; -} - -void LCAO_Matrix::destroy_T_R_sparse() -{ - ModuleBase::TITLE("LCAO_Matrix", "destroy_T_R_sparse"); - - if (GlobalV::NSPIN != 4) - { - std::map, std::map>> empty_TR_sparse; - TR_sparse.swap(empty_TR_sparse); - } return; } \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h index 9ff897ab1a..4a95c03984 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h @@ -14,24 +14,28 @@ #include #endif -class LCAO_Matrix -{ +class LCAO_Matrix { public: LCAO_Matrix(); ~LCAO_Matrix(); - void divide_HS_in_frag(const bool isGamma, Parallel_Orbitals& pv, const int& nks); + void divide_HS_in_frag(const bool isGamma, + Parallel_Orbitals& pv, + const int& nks); // folding the fixed Hamiltonian (T+Vnl) if // k-point algorithm is used. - void folding_fixedH(const int& ik, const std::vector>& kvec_d, bool cal_syns = false); + void folding_fixedH(const int& ik, + const std::vector>& kvec_d, + bool cal_syns = false); Parallel_Orbitals* ParaV; #ifdef __EXX using TAC = std::pair>; std::vector>>>* Hexxd; - std::vector>>>>* Hexxc; + std::vector>>>>* + Hexxc; /// @brief Hexxk for all k-points, only for the 1st scf loop ofrestart load std::vector> Hexxd_k_load; std::vector>> Hexxc_k_load; @@ -74,25 +78,25 @@ class LCAO_Matrix std::complex**** SlocR_tr_soc; std::complex**** HR_tr_soc; - // jingan add 2021-6-4, modify 2021-12-2 - // Sparse form of HR and SR, the format is [R_direct_coor][orbit_row][orbit_col] - - // For HR_sparse[2], when nspin=1, only 0 is valid, when nspin=2, 0 means spin up, 1 means spin down - std::map, std::map>> HR_sparse[2]; - std::map, std::map>> SR_sparse; - std::map, std::map>> TR_sparse; - - // Record all R direct coordinate information, even if HR or SR is a zero matrix + // Record all R direct coordinate information, even if HR or SR is a zero + // matrix std::set> all_R_coor; - // Records the R direct coordinates of HR and SR output, This variable will be filled with data when HR and SR files - // are output. + // Records the R direct coordinates of HR and SR output, This variable will + // be filled with data when HR and SR files are output. std::set> output_R_coor; template - static void set_mat2d(const int& global_ir, const int& global_ic, const T& v, const Parallel_Orbitals& pv, T* mat); + static void set_mat2d(const int& global_ir, + const int& global_ic, + const T& v, + const Parallel_Orbitals& pv, + T* mat); - void set_HSgamma(const int& iw1_all, const int& iw2_all, const double& v, double* HSloc); + void set_HSgamma(const int& iw1_all, + const int& iw2_all, + const double& v, + double* HSloc); void set_HR_tr(const int& Rx, const int& Ry, @@ -112,16 +116,11 @@ class LCAO_Matrix void zeros_HSk(const char& mtype); - void update_Hloc(void); + void update_Hloc(); void update_Hloc2(const int& ik); void output_HSk(const char& mtype, std::string& fn); - - // jingan add 2021-6-4, modify 2021-12-2 - void destroy_HS_R_sparse(LCAO_HS_Arrays& HS_Arrays); - - void destroy_T_R_sparse(void); }; #include "LCAO_matrix.hpp" diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp index 17cab1f47f..34c11a3c41 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp @@ -13,8 +13,7 @@ void sparse_format::cal_HSR(const Parallel_Orbitals& pv, const int& current_spin, const double& sparse_thr, const int (&nmp)[3], - hamilt::Hamilt>* p_ham) -{ + hamilt::Hamilt>* p_ham) { ModuleBase::TITLE("sparse_format", "cal_HSR"); sparse_format::set_R_range(lm.all_R_coor, grid); @@ -22,62 +21,69 @@ void sparse_format::cal_HSR(const Parallel_Orbitals& pv, const int nspin = GlobalV::NSPIN; // cal_STN_R_sparse(current_spin, sparse_thr); - if (nspin == 1 || nspin == 2) - { + if (nspin == 1 || nspin == 2) { hamilt::HamiltLCAO, double>* p_ham_lcao - = dynamic_cast, double>*>(p_ham); + = dynamic_cast, double>*>( + p_ham); - if (TD_Velocity::tddft_velocity) - { - sparse_format::cal_HContainer_td(pv, - current_spin, - sparse_thr, - *(p_ham_lcao->getHR()), - TD_Velocity::td_vel_op->HR_sparse_td_vel[current_spin]); - } - else - { + if (TD_Velocity::tddft_velocity) { + sparse_format::cal_HContainer_td( + pv, + current_spin, + sparse_thr, + *(p_ham_lcao->getHR()), + TD_Velocity::td_vel_op->HR_sparse_td_vel[current_spin]); + } else { sparse_format::cal_HContainer_d(pv, current_spin, sparse_thr, *(p_ham_lcao->getHR()), - lm.HR_sparse[current_spin]); + HS_Arrays.HR_sparse[current_spin]); } - sparse_format::cal_HContainer_d(pv, current_spin, sparse_thr, *(p_ham_lcao->getSR()), lm.SR_sparse); - } - else if (nspin == 4) - { - hamilt::HamiltLCAO, std::complex>* p_ham_lcao - = dynamic_cast, std::complex>*>(p_ham); + sparse_format::cal_HContainer_d(pv, + current_spin, + sparse_thr, + *(p_ham_lcao->getSR()), + HS_Arrays.SR_sparse); + } else if (nspin == 4) { + hamilt::HamiltLCAO, std::complex>* + p_ham_lcao + = dynamic_cast, + std::complex>*>(p_ham); - sparse_format::cal_HContainer_cd(pv, current_spin, sparse_thr, *(p_ham_lcao->getHR()), HS_Arrays.HR_soc_sparse); + sparse_format::cal_HContainer_cd(pv, + current_spin, + sparse_thr, + *(p_ham_lcao->getHR()), + HS_Arrays.HR_soc_sparse); - sparse_format::cal_HContainer_cd(pv, current_spin, sparse_thr, *(p_ham_lcao->getSR()), HS_Arrays.SR_soc_sparse); - } - else - { + sparse_format::cal_HContainer_cd(pv, + current_spin, + sparse_thr, + *(p_ham_lcao->getSR()), + HS_Arrays.SR_soc_sparse); + } else { ModuleBase::WARNING_QUIT("cal_HSR", "check the value of nspin."); } // only old DFT+U method need to cal extra contribution to HR - if (GlobalV::dft_plus_u == 2) - { - if (nspin == 1 || nspin == 2) - { - cal_HR_dftu(pv, lm.all_R_coor, lm.SR_sparse, lm.HR_sparse, current_spin, sparse_thr); - } - else if (nspin == 4) - { + if (GlobalV::dft_plus_u == 2) { + if (nspin == 1 || nspin == 2) { + cal_HR_dftu(pv, + lm.all_R_coor, + HS_Arrays.SR_sparse, + HS_Arrays.HR_sparse, + current_spin, + sparse_thr); + } else if (nspin == 4) { cal_HR_dftu_soc(pv, lm.all_R_coor, HS_Arrays.SR_soc_sparse, HS_Arrays.HR_soc_sparse, current_spin, sparse_thr); - } - else - { + } else { ModuleBase::WARNING_QUIT("cal_HSR", "check the value of nspin."); } } @@ -85,15 +91,21 @@ void sparse_format::cal_HSR(const Parallel_Orbitals& pv, #ifdef __EXX #ifdef __MPI // if EXX is considered - if (GlobalC::exx_info.info_global.cal_exx) - { - if (GlobalC::exx_info.info_ri.real_number) - { - sparse_format::cal_HR_exx(lm, HS_Arrays, current_spin, sparse_thr, nmp, *lm.Hexxd); - } - else - { - sparse_format::cal_HR_exx(lm, HS_Arrays, current_spin, sparse_thr, nmp, *lm.Hexxc); + if (GlobalC::exx_info.info_global.cal_exx) { + if (GlobalC::exx_info.info_ri.real_number) { + sparse_format::cal_HR_exx(lm, + HS_Arrays, + current_spin, + sparse_thr, + nmp, + *lm.Hexxd); + } else { + sparse_format::cal_HR_exx(lm, + HS_Arrays, + current_spin, + sparse_thr, + nmp, + *lm.Hexxc); } } #endif // __MPI @@ -109,34 +121,30 @@ void sparse_format::cal_HContainer_d( const int& current_spin, const double& sparse_thr, const hamilt::HContainer& hR, - std::map, std::map>>& target) -{ + std::map, + std::map>>& target) { ModuleBase::TITLE("sparse_format", "cal_HContainer_d"); auto row_indexes = pv.get_indexes_row(); auto col_indexes = pv.get_indexes_col(); - for (int iap = 0; iap < hR.size_atom_pairs(); ++iap) - { + for (int iap = 0; iap < hR.size_atom_pairs(); ++iap) { int atom_i = hR.get_atom_pair(iap).get_atom_i(); int atom_j = hR.get_atom_pair(iap).get_atom_j(); int start_i = pv.atom_begin_row[atom_i]; int start_j = pv.atom_begin_col[atom_j]; int row_size = pv.get_row_size(atom_i); int col_size = pv.get_col_size(atom_j); - for (int iR = 0; iR < hR.get_atom_pair(iap).get_R_size(); ++iR) - { + for (int iR = 0; iR < hR.get_atom_pair(iap).get_R_size(); ++iR) { auto& matrix = hR.get_atom_pair(iap).get_HR_values(iR); - const ModuleBase::Vector3 r_index = hR.get_atom_pair(iap).get_R_index(iR); + const ModuleBase::Vector3 r_index + = hR.get_atom_pair(iap).get_R_index(iR); Abfs::Vector3_Order dR(r_index.x, r_index.y, r_index.z); - for (int i = 0; i < row_size; ++i) - { + for (int i = 0; i < row_size; ++i) { int mu = row_indexes[start_i + i]; - for (int j = 0; j < col_size; ++j) - { + for (int j = 0; j < col_size; ++j) { int nu = col_indexes[start_j + j]; const auto& value_tmp = matrix.get_value(i, j); - if (std::abs(value_tmp) > sparse_thr) - { + if (std::abs(value_tmp) > sparse_thr) { target[dR][mu][nu] = value_tmp; } } @@ -152,34 +160,31 @@ void sparse_format::cal_HContainer_cd( const int& current_spin, const double& sparse_thr, const hamilt::HContainer>& hR, - std::map, std::map>>>& target) -{ + std::map, + std::map>>>& + target) { ModuleBase::TITLE("sparse_format", "cal_HContainer_cd"); auto row_indexes = pv.get_indexes_row(); auto col_indexes = pv.get_indexes_col(); - for (int iap = 0; iap < hR.size_atom_pairs(); ++iap) - { + for (int iap = 0; iap < hR.size_atom_pairs(); ++iap) { int atom_i = hR.get_atom_pair(iap).get_atom_i(); int atom_j = hR.get_atom_pair(iap).get_atom_j(); int start_i = pv.atom_begin_row[atom_i]; int start_j = pv.atom_begin_col[atom_j]; int row_size = pv.get_row_size(atom_i); int col_size = pv.get_col_size(atom_j); - for (int iR = 0; iR < hR.get_atom_pair(iap).get_R_size(); ++iR) - { + for (int iR = 0; iR < hR.get_atom_pair(iap).get_R_size(); ++iR) { auto& matrix = hR.get_atom_pair(iap).get_HR_values(iR); - const ModuleBase::Vector3 r_index = hR.get_atom_pair(iap).get_R_index(iR); + const ModuleBase::Vector3 r_index + = hR.get_atom_pair(iap).get_R_index(iR); Abfs::Vector3_Order dR(r_index.x, r_index.y, r_index.z); - for (int i = 0; i < row_size; ++i) - { + for (int i = 0; i < row_size; ++i) { int mu = row_indexes[start_i + i]; - for (int j = 0; j < col_size; ++j) - { + for (int j = 0; j < col_size; ++j) { int nu = col_indexes[start_j + j]; const auto& value_tmp = matrix.get_value(i, j); - if (std::abs(value_tmp) > sparse_thr) - { + if (std::abs(value_tmp) > sparse_thr) { target[dR][mu][nu] = value_tmp; } } @@ -195,34 +200,32 @@ void sparse_format::cal_HContainer_td( const int& current_spin, const double& sparse_thr, const hamilt::HContainer& hR, - std::map, std::map>>>& target) -{ + std::map, + std::map>>>& + target) { ModuleBase::TITLE("sparse_format", "cal_HContainer_td"); auto row_indexes = pv.get_indexes_row(); auto col_indexes = pv.get_indexes_col(); - for (int iap = 0; iap < hR.size_atom_pairs(); ++iap) - { + for (int iap = 0; iap < hR.size_atom_pairs(); ++iap) { int atom_i = hR.get_atom_pair(iap).get_atom_i(); int atom_j = hR.get_atom_pair(iap).get_atom_j(); int start_i = pv.atom_begin_row[atom_i]; int start_j = pv.atom_begin_col[atom_j]; int row_size = pv.get_row_size(atom_i); int col_size = pv.get_col_size(atom_j); - for (int iR = 0; iR < hR.get_atom_pair(iap).get_R_size(); ++iR) - { + for (int iR = 0; iR < hR.get_atom_pair(iap).get_R_size(); ++iR) { auto& matrix = hR.get_atom_pair(iap).get_HR_values(iR); - const ModuleBase::Vector3 r_index = hR.get_atom_pair(iap).get_R_index(iR); + const ModuleBase::Vector3 r_index + = hR.get_atom_pair(iap).get_R_index(iR); Abfs::Vector3_Order dR(r_index.x, r_index.y, r_index.z); - for (int i = 0; i < row_size; ++i) - { + for (int i = 0; i < row_size; ++i) { int mu = row_indexes[start_i + i]; - for (int j = 0; j < col_size; ++j) - { + for (int j = 0; j < col_size; ++j) { int nu = col_indexes[start_j + j]; - const auto& value_tmp = std::complex(matrix.get_value(i, j), 0.0); - if (std::abs(value_tmp) > sparse_thr) - { + const auto& value_tmp + = std::complex(matrix.get_value(i, j), 0.0); + if (std::abs(value_tmp) > sparse_thr) { target[dR][mu][nu] += value_tmp; } } @@ -237,47 +240,33 @@ void sparse_format::cal_HContainer_td( void sparse_format::clear_zero_elements(LCAO_Matrix& lm, LCAO_HS_Arrays& HS_Arrays, const int& current_spin, - const double& sparse_thr) -{ + const double& sparse_thr) { ModuleBase::TITLE("sparse_format", "clear_zero_elements"); - if (GlobalV::NSPIN != 4) - { - for (auto& R_loop: lm.HR_sparse[current_spin]) - { - for (auto& row_loop: R_loop.second) - { + if (GlobalV::NSPIN != 4) { + for (auto& R_loop: HS_Arrays.HR_sparse[current_spin]) { + for (auto& row_loop: R_loop.second) { auto& col_map = row_loop.second; auto iter = col_map.begin(); - while (iter != col_map.end()) - { - if (std::abs(iter->second) <= sparse_thr) - { + while (iter != col_map.end()) { + if (std::abs(iter->second) <= sparse_thr) { col_map.erase(iter++); - } - else - { + } else { iter++; } } } } - if (TD_Velocity::tddft_velocity) - { - for (auto& R_loop: TD_Velocity::td_vel_op->HR_sparse_td_vel[current_spin]) - { - for (auto& row_loop: R_loop.second) - { + if (TD_Velocity::tddft_velocity) { + for (auto& R_loop: + TD_Velocity::td_vel_op->HR_sparse_td_vel[current_spin]) { + for (auto& row_loop: R_loop.second) { auto& col_map = row_loop.second; auto iter = col_map.begin(); - while (iter != col_map.end()) - { - if (std::abs(iter->second) <= sparse_thr) - { + while (iter != col_map.end()) { + if (std::abs(iter->second) <= sparse_thr) { col_map.erase(iter++); - } - else - { + } else { iter++; } } @@ -285,62 +274,42 @@ void sparse_format::clear_zero_elements(LCAO_Matrix& lm, } } - for (auto& R_loop: lm.SR_sparse) - { - for (auto& row_loop: R_loop.second) - { + for (auto& R_loop: HS_Arrays.SR_sparse) { + for (auto& row_loop: R_loop.second) { auto& col_map = row_loop.second; auto iter = col_map.begin(); - while (iter != col_map.end()) - { - if (std::abs(iter->second) <= sparse_thr) - { + while (iter != col_map.end()) { + if (std::abs(iter->second) <= sparse_thr) { col_map.erase(iter++); - } - else - { + } else { iter++; } } } } - } - else - { - for (auto& R_loop: HS_Arrays.HR_soc_sparse) - { - for (auto& row_loop: R_loop.second) - { + } else { + for (auto& R_loop: HS_Arrays.HR_soc_sparse) { + for (auto& row_loop: R_loop.second) { auto& col_map = row_loop.second; auto iter = col_map.begin(); - while (iter != col_map.end()) - { - if (std::abs(iter->second) <= sparse_thr) - { + while (iter != col_map.end()) { + if (std::abs(iter->second) <= sparse_thr) { col_map.erase(iter++); - } - else - { + } else { iter++; } } // end while iter } // end row loop } // end R loop - for (auto& R_loop: HS_Arrays.SR_soc_sparse) - { - for (auto& row_loop: R_loop.second) - { + for (auto& R_loop: HS_Arrays.SR_soc_sparse) { + for (auto& row_loop: R_loop.second) { auto& col_map = row_loop.second; auto iter = col_map.begin(); - while (iter != col_map.end()) - { - if (std::abs(iter->second) <= sparse_thr) - { + while (iter != col_map.end()) { + if (std::abs(iter->second) <= sparse_thr) { col_map.erase(iter++); - } - else - { + } else { iter++; } } // end while iter @@ -348,5 +317,39 @@ void sparse_format::clear_zero_elements(LCAO_Matrix& lm, } // end R_loop } + return; +} + +void sparse_format::destroy_HS_R_sparse(LCAO_HS_Arrays& HS_Arrays) { + ModuleBase::TITLE("sparse_format", "destroy_HS_R_sparse"); + + if (GlobalV::NSPIN != 4) { + std::map, + std::map>> + empty_HR_sparse_up; + std::map, + std::map>> + empty_HR_sparse_down; + std::map, + std::map>> + empty_SR_sparse; + HS_Arrays.HR_sparse[0].swap(empty_HR_sparse_up); + HS_Arrays.HR_sparse[1].swap(empty_HR_sparse_down); + HS_Arrays.SR_sparse.swap(empty_SR_sparse); + } else { + std::map, + std::map>>> + empty_HR_soc_sparse; + std::map, + std::map>>> + empty_SR_soc_sparse; + HS_Arrays.HR_soc_sparse.swap(empty_HR_soc_sparse); + HS_Arrays.SR_soc_sparse.swap(empty_SR_soc_sparse); + } + + // 'all_R_coor' has a small memory requirement and does not need to be + // deleted. std::set> empty_all_R_coor; + // all_R_coor.swap(empty_all_R_coor); + return; } \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h index 766a641889..13dd667019 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h @@ -4,8 +4,7 @@ #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h" #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" -namespace sparse_format -{ +namespace sparse_format { void cal_HSR(const Parallel_Orbitals& pv, LCAO_Matrix& lm, @@ -16,27 +15,36 @@ void cal_HSR(const Parallel_Orbitals& pv, const int (&nmp)[3], hamilt::Hamilt>* p_ham); -void cal_HContainer_d(const Parallel_Orbitals& pv, - const int& current_spin, - const double& sparse_threshold, - const hamilt::HContainer& hR, - std::map, std::map>>& target); +void cal_HContainer_d( + const Parallel_Orbitals& pv, + const int& current_spin, + const double& sparse_threshold, + const hamilt::HContainer& hR, + std::map, + std::map>>& target); void cal_HContainer_cd( const Parallel_Orbitals& pv, const int& current_spin, const double& sparse_threshold, const hamilt::HContainer>& hR, - std::map, std::map>>>& target); + std::map, + std::map>>>& target); void cal_HContainer_td( const Parallel_Orbitals& pv, const int& current_spin, const double& sparse_threshold, const hamilt::HContainer& hR, - std::map, std::map>>>& target); + std::map, + std::map>>>& target); + +void clear_zero_elements(LCAO_Matrix& lm, + LCAO_HS_Arrays& HS_Arrays, + const int& current_spin, + const double& sparse_thr); -void clear_zero_elements(LCAO_Matrix& lm, LCAO_HS_Arrays& HS_Arrays, const int& current_spin, const double& sparse_thr); +void destroy_HS_R_sparse(LCAO_HS_Arrays& HS_Arrays); } // namespace sparse_format diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_st.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/spar_st.cpp index 109e73221c..5b4bd48a3e 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_st.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_st.cpp @@ -11,12 +11,14 @@ void sparse_format::cal_SR( const Parallel_Orbitals& pv, std::set>& all_R_coor, - std::map, std::map>>& SR_sparse, - std::map, std::map>>>& SR_soc_sparse, + std::map, + std::map>>& SR_sparse, + std::map, + std::map>>>& + SR_soc_sparse, Grid_Driver& grid, const double& sparse_thr, - hamilt::Hamilt>* p_ham) -{ + hamilt::Hamilt>* p_ham) { ModuleBase::TITLE("sparse_format", "cal_SR"); sparse_format::set_R_range(all_R_coor, grid); @@ -24,19 +26,27 @@ void sparse_format::cal_SR( const int nspin = GlobalV::NSPIN; // cal_STN_R_sparse(current_spin, sparse_thr); - if (nspin == 1 || nspin == 2) - { + if (nspin == 1 || nspin == 2) { hamilt::HamiltLCAO, double>* p_ham_lcao - = dynamic_cast, double>*>(p_ham); + = dynamic_cast, double>*>( + p_ham); const int cspin = 0; - sparse_format::cal_HContainer_d(pv, cspin, sparse_thr, *(p_ham_lcao->getSR()), SR_sparse); - } - else if (nspin == 4) - { - hamilt::HamiltLCAO, std::complex>* p_ham_lcao - = dynamic_cast, std::complex>*>(p_ham); + sparse_format::cal_HContainer_d(pv, + cspin, + sparse_thr, + *(p_ham_lcao->getSR()), + SR_sparse); + } else if (nspin == 4) { + hamilt::HamiltLCAO, std::complex>* + p_ham_lcao + = dynamic_cast, + std::complex>*>(p_ham); const int cspin = 0; - sparse_format::cal_HContainer_cd(pv, cspin, sparse_thr, *(p_ham_lcao->getSR()), SR_soc_sparse); + sparse_format::cal_HContainer_cd(pv, + cspin, + sparse_thr, + *(p_ham_lcao->getSR()), + SR_soc_sparse); } return; @@ -48,8 +58,7 @@ void sparse_format::cal_TR(const UnitCell& ucell, LCAO_HS_Arrays& HS_arrays, Grid_Driver& grid, const TwoCenterBundle& two_center_bundle, - const double& sparse_thr) -{ + const double& sparse_thr) { ModuleBase::TITLE("sparse_format", "cal_TR"); // need to rebuild T(R) @@ -84,8 +93,7 @@ void sparse_format::cal_STN_R_for_T(const UnitCell& ucell, LCAO_Matrix& lm, LCAO_HS_Arrays& HS_arrays, Grid_Driver& grid, - const double& sparse_thr) -{ + const double& sparse_thr) { ModuleBase::TITLE("sparse_format", "cal_STN_R_for_T"); const int nspin = GlobalV::NSPIN; @@ -97,18 +105,15 @@ void sparse_format::cal_STN_R_for_T(const UnitCell& ucell, double tmp = 0.0; std::complex tmpc = std::complex(0.0, 0.0); - for (int T1 = 0; T1 < ucell.ntype; ++T1) - { + for (int T1 = 0; T1 < ucell.ntype; ++T1) { Atom* atom1 = &ucell.atoms[T1]; - for (int I1 = 0; I1 < atom1->na; ++I1) - { + for (int I1 = 0; I1 < atom1->na; ++I1) { tau1 = atom1->tau[I1]; grid.Find_atom(ucell, tau1, T1, I1); Atom* atom1 = &ucell.atoms[T1]; const int start = ucell.itiaiw2iwt(T1, I1, 0); - for (int ad = 0; ad < grid.getAdjacentNum() + 1; ++ad) - { + for (int ad = 0; ad < grid.getAdjacentNum() + 1; ++ad) { const int T2 = grid.getType(ad); const int I2 = grid.getNatom(ad); Atom* atom2 = &ucell.atoms[T2]; @@ -116,19 +121,17 @@ void sparse_format::cal_STN_R_for_T(const UnitCell& ucell, tau2 = grid.getAdjacentTau(ad); dtau = tau2 - tau1; double distance = dtau.norm() * ucell.lat0; - double rcut = GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ORB.Phi[T2].getRcut(); + double rcut = GlobalC::ORB.Phi[T1].getRcut() + + GlobalC::ORB.Phi[T2].getRcut(); bool adj = false; - if (distance < rcut) - { + if (distance < rcut) { adj = true; } - else if (distance >= rcut) - { - for (int ad0 = 0; ad0 < grid.getAdjacentNum() + 1; ++ad0) - { + else if (distance >= rcut) { + for (int ad0 = 0; ad0 < grid.getAdjacentNum() + 1; ++ad0) { const int T0 = grid.getType(ad0); tau0 = grid.getAdjacentTau(ad0); @@ -138,49 +141,46 @@ void sparse_format::cal_STN_R_for_T(const UnitCell& ucell, double distance1 = dtau1.norm() * ucell.lat0; double distance2 = dtau2.norm() * ucell.lat0; - double rcut1 = GlobalC::ORB.Phi[T1].getRcut() + ucell.infoNL.Beta[T0].get_rcut_max(); - double rcut2 = GlobalC::ORB.Phi[T2].getRcut() + ucell.infoNL.Beta[T0].get_rcut_max(); + double rcut1 = GlobalC::ORB.Phi[T1].getRcut() + + ucell.infoNL.Beta[T0].get_rcut_max(); + double rcut2 = GlobalC::ORB.Phi[T2].getRcut() + + ucell.infoNL.Beta[T0].get_rcut_max(); - if (distance1 < rcut1 && distance2 < rcut2) - { + if (distance1 < rcut1 && distance2 < rcut2) { adj = true; break; } } } - if (adj) - { + if (adj) { const int start2 = ucell.itiaiw2iwt(T2, I2, 0); - Abfs::Vector3_Order dR(grid.getBox(ad).x, grid.getBox(ad).y, grid.getBox(ad).z); + Abfs::Vector3_Order dR(grid.getBox(ad).x, + grid.getBox(ad).y, + grid.getBox(ad).z); - for (int ii = 0; ii < atom1->nw * GlobalV::NPOL; ii++) - { + for (int ii = 0; ii < atom1->nw * GlobalV::NPOL; ii++) { const int iw1_all = start + ii; const int mu = pv.global2local_row(iw1_all); - if (mu < 0) - { + if (mu < 0) { continue; } - for (int jj = 0; jj < atom2->nw * GlobalV::NPOL; jj++) - { + for (int jj = 0; jj < atom2->nw * GlobalV::NPOL; jj++) { int iw2_all = start2 + jj; const int nu = pv.global2local_col(iw2_all); - if (nu < 0) - { + if (nu < 0) { continue; } - if (nspin == 1 || nspin == 2) - { + if (nspin == 1 || nspin == 2) { tmp = HS_arrays.Hloc_fixedR[index]; - if (std::abs(tmp) > sparse_thr) - { - lm.TR_sparse[dR][iw1_all][iw2_all] = tmp; + if (std::abs(tmp) > sparse_thr) { + HS_arrays.TR_sparse[dR][iw1_all][iw2_all] + = tmp; } } @@ -194,3 +194,15 @@ void sparse_format::cal_STN_R_for_T(const UnitCell& ucell, return; } + +void sparse_format::destroy_T_R_sparse(LCAO_HS_Arrays& HS_Arrays) { + ModuleBase::TITLE("sparse_format", "destroy_T_R_sparse"); + + if (GlobalV::NSPIN != 4) { + std::map, + std::map>> + empty_TR_sparse; + HS_Arrays.TR_sparse.swap(empty_TR_sparse); + } + return; +} diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_st.h b/source/module_hamilt_lcao/hamilt_lcaodft/spar_st.h index 39b5cce533..cad9e0d708 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_st.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_st.h @@ -5,13 +5,15 @@ #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h" #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" -namespace sparse_format -{ +namespace sparse_format { //! calculate overlap matrix with lattice vector R void cal_SR(const Parallel_Orbitals& pv, std::set>& all_R_coor, - std::map, std::map>>& SR_sparse, - std::map, std::map>>>& SR_soc_sparse, + std::map, + std::map>>& SR_sparse, + std::map, + std::map>>>& + SR_soc_sparse, Grid_Driver& grid, const double& sparse_thr, hamilt::Hamilt>* p_ham); @@ -32,6 +34,8 @@ void cal_STN_R_for_T(const UnitCell& ucell, LCAO_HS_Arrays& HS_arrays, Grid_Driver& grid, const double& sparse_thr); + +void destroy_T_R_sparse(LCAO_HS_Arrays& HS_Arrays); } // namespace sparse_format #endif diff --git a/source/module_hamilt_lcao/module_deepks/test/Makefile.Objects b/source/module_hamilt_lcao/module_deepks/test/Makefile.Objects index d4b56b496f..71d2599b8f 100644 --- a/source/module_hamilt_lcao/module_deepks/test/Makefile.Objects +++ b/source/module_hamilt_lcao/module_deepks/test/Makefile.Objects @@ -55,6 +55,7 @@ read_pp_upf201.o\ read_pp_vwr.o\ read_pp_blps.o\ unitcell.o\ +check_atomic_stru.o\ read_atoms.o\ read_cell_pseudopots.o\ setup_nonlocal.o diff --git a/source/module_hamilt_lcao/module_gint/CMakeLists.txt b/source/module_hamilt_lcao/module_gint/CMakeLists.txt index 9243693bce..7fca9dbca9 100644 --- a/source/module_hamilt_lcao/module_gint/CMakeLists.txt +++ b/source/module_hamilt_lcao/module_gint/CMakeLists.txt @@ -9,7 +9,6 @@ list(APPEND objects gint_tau.cpp gint_vl.cpp gint_k_env.cpp - gint_k_sparse.cpp gint_k_sparse1.cpp gint_k_pvpr.cpp gint_k_pvdpr.cpp diff --git a/source/module_hamilt_lcao/module_gint/gint_k.h b/source/module_hamilt_lcao/module_gint/gint_k.h index ac8736be0c..470a108ffa 100644 --- a/source/module_hamilt_lcao/module_gint/gint_k.h +++ b/source/module_hamilt_lcao/module_gint/gint_k.h @@ -11,13 +11,9 @@ // add by jingan for map<> in 2021-12-2, will be deleted in the future #include "module_base/abfs-vector3_order.h" -class Gint_k : public Gint -{ +class Gint_k : public Gint { public: - ~Gint_k() - { - destroy_pvpR(); - } + ~Gint_k() { destroy_pvpR(); } //------------------------------------------------------ // in gint_k_pvpr.cpp //------------------------------------------------------ @@ -25,25 +21,17 @@ class Gint_k : public Gint // for calculating hamiltonian // reset the spin. - void reset_spin(const int& spin_now_in) - { - this->spin_now = spin_now_in; - }; + void reset_spin(const int& spin_now_in) { this->spin_now = spin_now_in; }; // get the spin. - int get_spin() const - { - return spin_now; - } + int get_spin() const { return spin_now; } // renew gint index for new iteration - void renew(const bool& soft = false) - { - if (soft && this->spin_now == 0) - { // in this case, gint will not be recalculated + void renew(const bool& soft = false) { + if (soft + && this->spin_now + == 0) { // in this case, gint will not be recalculated return; - } - else if (this->spin_now != -1) - { + } else if (this->spin_now != -1) { int start_spin = -1; this->reset_spin(start_spin); this->destroy_pvpR(); @@ -66,8 +54,12 @@ class Gint_k : public Gint * @brief transfer pvpR to this->hRGint * then pass this->hRGint to Veff::hR */ - void transfer_pvpR(hamilt::HContainer* hR, const UnitCell* ucell_in, Grid_Driver* gd); - void transfer_pvpR(hamilt::HContainer>* hR, const UnitCell* ucell_in, Grid_Driver* gd); + void transfer_pvpR(hamilt::HContainer* hR, + const UnitCell* ucell_in, + Grid_Driver* gd); + void transfer_pvpR(hamilt::HContainer>* hR, + const UnitCell* ucell_in, + Grid_Driver* gd); //------------------------------------------------------ // in gint_k_env.cpp @@ -80,18 +72,6 @@ class Gint_k : public Gint const std::vector>& kvec_d, UnitCell& ucell); - //------------------------------------------------------ - // in gint_k_sparse.cpp - //------------------------------------------------------ - // related to sparse matrix - // jingan add 2021-6-4, modify 2021-12-2 - void distribute_pvpR_sparseMatrix( - const int current_spin, - const double& sparse_threshold, - const std::map, std::map>>& pvpR_sparseMatrix, - LCAO_Matrix* LM, - Parallel_Orbitals* pv); - //------------------------------------------------------ // in gint_k_sparse1.cpp //------------------------------------------------------ @@ -100,7 +80,9 @@ class Gint_k : public Gint const int current_spin, const int dim, const double& sparse_threshold, - const std::map, std::map>>& pvdpR_sparseMatrix, + const std::map, + std::map>>& + pvdpR_sparseMatrix, LCAO_Matrix* LM, LCAO_HS_Arrays& HS_Arrays, Parallel_Orbitals* pv); @@ -108,7 +90,9 @@ class Gint_k : public Gint void distribute_pvdpR_soc_sparseMatrix( const int dim, const double& sparse_threshold, - const std::map, std::map>>>& + const std::map< + Abfs::Vector3_Order, + std::map>>>& pvdpR_soc_sparseMatrix, LCAO_Matrix* LM, LCAO_HS_Arrays& HS_Arrays, diff --git a/source/module_hamilt_lcao/module_gint/gint_k_sparse.cpp b/source/module_hamilt_lcao/module_gint/gint_k_sparse.cpp deleted file mode 100644 index b729d41912..0000000000 --- a/source/module_hamilt_lcao/module_gint/gint_k_sparse.cpp +++ /dev/null @@ -1,143 +0,0 @@ -#include "gint_k.h" -#include "grid_technique.h" -#include "module_base/global_function.h" -#include "module_base/global_variable.h" -#include "module_base/memory.h" -#include "module_base/parallel_reduce.h" -#include "module_base/timer.h" -#include "module_base/ylm.h" -#include "module_basis/module_ao/ORB_read.h" -#include "module_cell/module_neighbor/sltk_grid_driver.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" - -void Gint_k::distribute_pvpR_sparseMatrix( - const int current_spin, - const double& sparse_threshold, - const std::map, std::map>>& pvpR_sparseMatrix, - LCAO_Matrix* LM, - Parallel_Orbitals* pv) -{ - ModuleBase::TITLE("Gint_k", "distribute_pvpR_sparseMatrix"); - - int total_R_num = LM->all_R_coor.size(); - int* nonzero_num = new int[total_R_num]; - int* minus_nonzero_num = new int[total_R_num]; - ModuleBase::GlobalFunc::ZEROS(nonzero_num, total_R_num); - ModuleBase::GlobalFunc::ZEROS(minus_nonzero_num, total_R_num); - int count = 0; - for (auto& R_coor: LM->all_R_coor) - { - auto iter = pvpR_sparseMatrix.find(R_coor); - if (iter != pvpR_sparseMatrix.end()) - { - for (auto& row_loop: iter->second) - { - nonzero_num[count] += row_loop.second.size(); - } - } - - auto minus_R_coor = -1 * R_coor; - - iter = pvpR_sparseMatrix.find(minus_R_coor); - if (iter != pvpR_sparseMatrix.end()) - { - for (auto& row_loop: iter->second) - { - minus_nonzero_num[count] += row_loop.second.size(); - } - } - - count++; - } - - Parallel_Reduce::reduce_all(nonzero_num, total_R_num); - Parallel_Reduce::reduce_all(minus_nonzero_num, total_R_num); - // Parallel_Reduce::reduce_pool(nonzero_num, total_R_num); - // Parallel_Reduce::reduce_pool(minus_nonzero_num, total_R_num); - - double* tmp = nullptr; - tmp = new double[GlobalV::NLOCAL]; - - count = 0; - for (auto& R_coor: LM->all_R_coor) - { - if (nonzero_num[count] != 0 || minus_nonzero_num[count] != 0) - { - auto minus_R_coor = -1 * R_coor; - - for (int row = 0; row < GlobalV::NLOCAL; ++row) - { - ModuleBase::GlobalFunc::ZEROS(tmp, GlobalV::NLOCAL); - - auto iter = pvpR_sparseMatrix.find(R_coor); - if (iter != pvpR_sparseMatrix.end()) - { - - if (this->gridt->trace_lo[row] >= 0) - { - auto row_iter = iter->second.find(row); - if (row_iter != iter->second.end()) - { - for (auto& value: row_iter->second) - { - tmp[value.first] = value.second; - } - } - } - } - - auto minus_R_iter = pvpR_sparseMatrix.find(minus_R_coor); - if (minus_R_iter != pvpR_sparseMatrix.end()) - { - for (int col = 0; col < row; ++col) - { - if (this->gridt->trace_lo[col] >= 0) - { - auto row_iter = minus_R_iter->second.find(col); - if (row_iter != minus_R_iter->second.end()) - { - auto col_iter = row_iter->second.find(row); - if (col_iter != row_iter->second.end()) - { - tmp[col] = col_iter->second; - } - } - } - } - } - - Parallel_Reduce::reduce_pool(tmp, GlobalV::NLOCAL); - - if (pv->global2local_row(row) >= 0) - { - for (int col = 0; col < GlobalV::NLOCAL; ++col) - { - if (pv->global2local_col(col) >= 0) - { - if (std::abs(tmp[col]) > sparse_threshold) - { - double& value = LM->HR_sparse[current_spin][R_coor][row][col]; - value += tmp[col]; - if (std::abs(value) <= sparse_threshold) - { - LM->HR_sparse[current_spin][R_coor][row].erase(col); - } - } - } - } - } - } - } - - count++; - } - - delete[] nonzero_num; - delete[] minus_nonzero_num; - delete[] tmp; - nonzero_num = nullptr; - minus_nonzero_num = nullptr; - tmp = nullptr; - - return; -} \ No newline at end of file diff --git a/source/module_io/write_HS_R.cpp b/source/module_io/write_HS_R.cpp index 6bbd860c83..fb75d6e28f 100644 --- a/source/module_io/write_HS_R.cpp +++ b/source/module_io/write_HS_R.cpp @@ -9,7 +9,8 @@ // if 'binary=true', output binary file. // The 'sparse_thr' is the accuracy of the sparse matrix. -// If the absolute value of the matrix element is less than or equal to the 'sparse_thr', it will be ignored. +// If the absolute value of the matrix element is less than or equal to the +// 'sparse_thr', it will be ignored. void ModuleIO::output_HSR(const int& istep, const ModuleBase::matrix& v_eff, const Parallel_Orbitals& pv, @@ -21,8 +22,7 @@ void ModuleIO::output_HSR(const int& istep, const std::string& HR_filename_up, const std::string HR_filename_down, const bool& binary, - const double& sparse_thr) -{ + const double& sparse_thr) { ModuleBase::TITLE("ModuleIO", "output_HSR"); ModuleBase::timer::tick("ModuleIO", "output_HSR"); @@ -30,34 +30,58 @@ void ModuleIO::output_HSR(const int& istep, LCAO_HS_Arrays HS_Arrays; - if (nspin == 1 || nspin == 4) - { + if (nspin == 1 || nspin == 4) { const int spin_now = 0; // jingan add 2021-6-4, modify 2021-12-2 - sparse_format::cal_HSR(pv, lm, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham); - } - else if (nspin == 2) - { + sparse_format::cal_HSR(pv, + lm, + HS_Arrays, + grid, + spin_now, + sparse_thr, + kv.nmp, + p_ham); + } else if (nspin == 2) { int spin_now = 1; // save HR of spin down first (the current spin always be down) - sparse_format::cal_HSR(pv, lm, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham); + sparse_format::cal_HSR(pv, + lm, + HS_Arrays, + grid, + spin_now, + sparse_thr, + kv.nmp, + p_ham); // cal HR of the spin up - if (GlobalV::VL_IN_H) - { + if (GlobalV::VL_IN_H) { const int ik = 0; p_ham->refresh(); p_ham->updateHk(ik); spin_now = 0; } - sparse_format::cal_HSR(pv, lm, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham); + sparse_format::cal_HSR(pv, + lm, + HS_Arrays, + grid, + spin_now, + sparse_thr, + kv.nmp, + p_ham); } - ModuleIO::save_HSR_sparse(istep, lm, HS_Arrays, sparse_thr, binary, SR_filename, HR_filename_up, HR_filename_down); + ModuleIO::save_HSR_sparse(istep, + lm, + HS_Arrays, + sparse_thr, + binary, + SR_filename, + HR_filename_up, + HR_filename_down); - lm.destroy_HS_R_sparse(HS_Arrays); + sparse_format::destroy_HS_R_sparse(HS_Arrays); ModuleBase::timer::tick("ModuleIO", "output_HSR"); return; @@ -71,8 +95,7 @@ void ModuleIO::output_dHR(const int& istep, const TwoCenterBundle& two_center_bundle, const K_Vectors& kv, const bool& binary, - const double& sparse_thr) -{ + const double& sparse_thr) { ModuleBase::TITLE("ModuleIO", "output_dHR"); ModuleBase::timer::tick("ModuleIO", "output_dHR"); @@ -82,31 +105,40 @@ void ModuleIO::output_dHR(const int& istep, const int nspin = GlobalV::NSPIN; - if (nspin == 1 || nspin == 4) - { + if (nspin == 1 || nspin == 4) { // mohan add 2024-04-01 const int cspin = 0; - sparse_format::cal_dH(lm, HS_Arrays, grid, two_center_bundle, cspin, sparse_thr, gint_k); - } - else if (nspin == 2) - { - for (int cspin = 0; cspin < 2; cspin++) - { - // note: some MPI process will not have grids when MPI cores are too many, - // v_eff in these processes are empty - const double* vr_eff1 = v_eff.nc * v_eff.nr > 0 ? &(v_eff(cspin, 0)) : nullptr; - - if (!GlobalV::GAMMA_ONLY_LOCAL) - { - if (GlobalV::VL_IN_H) - { - Gint_inout inout(vr_eff1, cspin, Gint_Tools::job_type::dvlocal); + sparse_format::cal_dH(lm, + HS_Arrays, + grid, + two_center_bundle, + cspin, + sparse_thr, + gint_k); + } else if (nspin == 2) { + for (int cspin = 0; cspin < 2; cspin++) { + // note: some MPI process will not have grids when MPI cores are too + // many, v_eff in these processes are empty + const double* vr_eff1 + = v_eff.nc * v_eff.nr > 0 ? &(v_eff(cspin, 0)) : nullptr; + + if (!GlobalV::GAMMA_ONLY_LOCAL) { + if (GlobalV::VL_IN_H) { + Gint_inout inout(vr_eff1, + cspin, + Gint_Tools::job_type::dvlocal); gint_k.cal_gint(&inout); } } - sparse_format::cal_dH(lm, HS_Arrays, grid, two_center_bundle, cspin, sparse_thr, gint_k); + sparse_format::cal_dH(lm, + HS_Arrays, + grid, + two_center_bundle, + cspin, + sparse_thr, + gint_k); } } // mohan update 2024-04-01 @@ -126,20 +158,32 @@ void ModuleIO::output_SR(Parallel_Orbitals& pv, hamilt::Hamilt>* p_ham, const std::string& SR_filename, const bool& binary, - const double& sparse_thr) -{ + const double& sparse_thr) { ModuleBase::TITLE("ModuleIO", "output_SR"); ModuleBase::timer::tick("ModuleIO", "output_SR"); LCAO_HS_Arrays HS_Arrays; - sparse_format::cal_SR(pv, lm.all_R_coor, lm.SR_sparse, HS_Arrays.SR_soc_sparse, grid, sparse_thr, p_ham); + sparse_format::cal_SR(pv, + lm.all_R_coor, + HS_Arrays.SR_sparse, + HS_Arrays.SR_soc_sparse, + grid, + sparse_thr, + p_ham); const int istep = 0; - ModuleIO::save_sparse(lm.SR_sparse, lm.all_R_coor, sparse_thr, binary, SR_filename, *lm.ParaV, "S", istep); + ModuleIO::save_sparse(HS_Arrays.SR_sparse, + lm.all_R_coor, + sparse_thr, + binary, + SR_filename, + *lm.ParaV, + "S", + istep); - lm.destroy_HS_R_sparse(HS_Arrays); + sparse_format::destroy_HS_R_sparse(HS_Arrays); ModuleBase::timer::tick("ModuleIO", "output_SR"); return; @@ -153,29 +197,38 @@ void ModuleIO::output_TR(const int istep, const TwoCenterBundle& two_center_bundle, const std::string& TR_filename, const bool& binary, - const double& sparse_thr) -{ + const double& sparse_thr) { ModuleBase::TITLE("ModuleIO", "output_TR"); ModuleBase::timer::tick("ModuleIO", "output_TR"); std::stringstream sst; - if (GlobalV::CALCULATION == "md" && !GlobalV::out_app_flag) - { + if (GlobalV::CALCULATION == "md" && !GlobalV::out_app_flag) { sst << GlobalV::global_matrix_dir << istep << "_" << TR_filename; - } - else - { + } else { sst << GlobalV::global_out_dir << TR_filename; } // need Hloc_fixedR LCAO_HS_Arrays HS_arrays; - sparse_format::cal_TR(ucell, pv, lm, HS_arrays, grid, two_center_bundle, sparse_thr); - - ModuleIO::save_sparse(lm.TR_sparse, lm.all_R_coor, sparse_thr, binary, sst.str().c_str(), *(lm.ParaV), "T", istep); - - lm.destroy_T_R_sparse(); + sparse_format::cal_TR(ucell, + pv, + lm, + HS_arrays, + grid, + two_center_bundle, + sparse_thr); + + ModuleIO::save_sparse(HS_arrays.TR_sparse, + lm.all_R_coor, + sparse_thr, + binary, + sst.str().c_str(), + *(lm.ParaV), + "T", + istep); + + sparse_format::destroy_T_R_sparse(HS_arrays); ModuleBase::timer::tick("ModuleIO", "output_TR"); return; diff --git a/source/module_io/write_HS_sparse.cpp b/source/module_io/write_HS_sparse.cpp index 3ece137c00..5e2784c16a 100644 --- a/source/module_io/write_HS_sparse.cpp +++ b/source/module_io/write_HS_sparse.cpp @@ -13,15 +13,14 @@ void ModuleIO::save_HSR_sparse(const int& istep, const bool& binary, const std::string& SR_filename, const std::string& HR_filename_up, - const std::string& HR_filename_down = "") -{ + const std::string& HR_filename_down = "") { ModuleBase::TITLE("ModuleIO", "save_HSR_sparse"); ModuleBase::timer::tick("ModuleIO", "save_HSR_sparse"); auto& all_R_coor_ptr = lm.all_R_coor; auto& output_R_coor_ptr = lm.output_R_coor; - auto& HR_sparse_ptr = lm.HR_sparse; - auto& SR_sparse_ptr = lm.SR_sparse; + auto& HR_sparse_ptr = HS_Arrays.HR_sparse; + auto& SR_sparse_ptr = HS_Arrays.SR_sparse; auto& HR_soc_sparse_ptr = HS_Arrays.HR_soc_sparse; auto& SR_soc_sparse_ptr = HS_Arrays.SR_soc_sparse; @@ -35,73 +34,59 @@ void ModuleIO::save_HSR_sparse(const int& istep, ModuleBase::GlobalFunc::ZEROS(S_nonzero_num, total_R_num); int spin_loop = 1; - if (GlobalV::NSPIN == 2) - { + if (GlobalV::NSPIN == 2) { spin_loop = 2; } - for (int ispin = 0; ispin < spin_loop; ++ispin) - { + for (int ispin = 0; ispin < spin_loop; ++ispin) { H_nonzero_num[ispin] = new int[total_R_num]; ModuleBase::GlobalFunc::ZEROS(H_nonzero_num[ispin], total_R_num); } int count = 0; - for (auto& R_coor: all_R_coor_ptr) - { - if (GlobalV::NSPIN != 4) - { - for (int ispin = 0; ispin < spin_loop; ++ispin) - { - if (TD_Velocity::tddft_velocity) - { - auto iter = TD_Velocity::td_vel_op->HR_sparse_td_vel[ispin].find(R_coor); - if (iter != TD_Velocity::td_vel_op->HR_sparse_td_vel[ispin].end()) - { - for (auto& row_loop: iter->second) - { - H_nonzero_num[ispin][count] += row_loop.second.size(); + for (auto& R_coor: all_R_coor_ptr) { + if (GlobalV::NSPIN != 4) { + for (int ispin = 0; ispin < spin_loop; ++ispin) { + if (TD_Velocity::tddft_velocity) { + auto iter + = TD_Velocity::td_vel_op->HR_sparse_td_vel[ispin].find( + R_coor); + if (iter + != TD_Velocity::td_vel_op->HR_sparse_td_vel[ispin] + .end()) { + for (auto& row_loop: iter->second) { + H_nonzero_num[ispin][count] + += row_loop.second.size(); } } - } - else - { + } else { auto iter = HR_sparse_ptr[ispin].find(R_coor); - if (iter != HR_sparse_ptr[ispin].end()) - { - for (auto& row_loop: iter->second) - { - H_nonzero_num[ispin][count] += row_loop.second.size(); + if (iter != HR_sparse_ptr[ispin].end()) { + for (auto& row_loop: iter->second) { + H_nonzero_num[ispin][count] + += row_loop.second.size(); } } } } auto iter = SR_sparse_ptr.find(R_coor); - if (iter != SR_sparse_ptr.end()) - { - for (auto& row_loop: iter->second) - { + if (iter != SR_sparse_ptr.end()) { + for (auto& row_loop: iter->second) { S_nonzero_num[count] += row_loop.second.size(); } } - } - else - { + } else { auto iter = HR_soc_sparse_ptr.find(R_coor); - if (iter != HR_soc_sparse_ptr.end()) - { - for (auto& row_loop: iter->second) - { + if (iter != HR_soc_sparse_ptr.end()) { + for (auto& row_loop: iter->second) { H_nonzero_num[0][count] += row_loop.second.size(); } } iter = SR_soc_sparse_ptr.find(R_coor); - if (iter != SR_soc_sparse_ptr.end()) - { - for (auto& row_loop: iter->second) - { + if (iter != SR_soc_sparse_ptr.end()) { + for (auto& row_loop: iter->second) { S_nonzero_num[count] += row_loop.second.size(); } } @@ -111,27 +96,20 @@ void ModuleIO::save_HSR_sparse(const int& istep, } Parallel_Reduce::reduce_all(S_nonzero_num, total_R_num); - for (int ispin = 0; ispin < spin_loop; ++ispin) - { + for (int ispin = 0; ispin < spin_loop; ++ispin) { Parallel_Reduce::reduce_all(H_nonzero_num[ispin], total_R_num); } - if (GlobalV::NSPIN == 2) - { - for (int index = 0; index < total_R_num; ++index) - { - if (H_nonzero_num[0][index] != 0 || H_nonzero_num[1][index] != 0 || S_nonzero_num[index] != 0) - { + if (GlobalV::NSPIN == 2) { + for (int index = 0; index < total_R_num; ++index) { + if (H_nonzero_num[0][index] != 0 || H_nonzero_num[1][index] != 0 + || S_nonzero_num[index] != 0) { output_R_number++; } } - } - else - { - for (int index = 0; index < total_R_num; ++index) - { - if (H_nonzero_num[0][index] != 0 || S_nonzero_num[index] != 0) - { + } else { + for (int index = 0; index < total_R_num; ++index) { + if (H_nonzero_num[0][index] != 0 || S_nonzero_num[index] != 0) { output_R_number++; } } @@ -139,14 +117,11 @@ void ModuleIO::save_HSR_sparse(const int& istep, std::stringstream ssh[2]; std::stringstream sss; - if (GlobalV::CALCULATION == "md" && !GlobalV::out_app_flag) - { + if (GlobalV::CALCULATION == "md" && !GlobalV::out_app_flag) { ssh[0] << GlobalV::global_matrix_dir << step << "_" << HR_filename_up; ssh[1] << GlobalV::global_matrix_dir << step << "_" << HR_filename_down; sss << GlobalV::global_matrix_dir << step << "_" << SR_filename; - } - else - { + } else { ssh[0] << GlobalV::global_out_dir << HR_filename_up; ssh[1] << GlobalV::global_out_dir << HR_filename_down; sss << GlobalV::global_out_dir << SR_filename; @@ -154,60 +129,49 @@ void ModuleIO::save_HSR_sparse(const int& istep, std::ofstream g1[2]; std::ofstream g2; - if (GlobalV::DRANK == 0) - { - if (binary) - { - for (int ispin = 0; ispin < spin_loop; ++ispin) - { - if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag && step) - { - g1[ispin].open(ssh[ispin].str().c_str(), std::ios::binary | std::ios::app); - } - else - { + if (GlobalV::DRANK == 0) { + if (binary) { + for (int ispin = 0; ispin < spin_loop; ++ispin) { + if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag + && step) { + g1[ispin].open(ssh[ispin].str().c_str(), + std::ios::binary | std::ios::app); + } else { g1[ispin].open(ssh[ispin].str().c_str(), std::ios::binary); } g1[ispin].write(reinterpret_cast(&step), sizeof(int)); - g1[ispin].write(reinterpret_cast(&GlobalV::NLOCAL), sizeof(int)); - g1[ispin].write(reinterpret_cast(&output_R_number), sizeof(int)); + g1[ispin].write(reinterpret_cast(&GlobalV::NLOCAL), + sizeof(int)); + g1[ispin].write(reinterpret_cast(&output_R_number), + sizeof(int)); } - if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag && step) - { + if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag && step) { g2.open(sss.str().c_str(), std::ios::binary | std::ios::app); - } - else - { + } else { g2.open(sss.str().c_str(), std::ios::binary); } g2.write(reinterpret_cast(&step), sizeof(int)); g2.write(reinterpret_cast(&GlobalV::NLOCAL), sizeof(int)); g2.write(reinterpret_cast(&output_R_number), sizeof(int)); - } - else - { - for (int ispin = 0; ispin < spin_loop; ++ispin) - { - if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag && step) - { + } else { + for (int ispin = 0; ispin < spin_loop; ++ispin) { + if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag + && step) { g1[ispin].open(ssh[ispin].str().c_str(), std::ios::app); - } - else - { + } else { g1[ispin].open(ssh[ispin].str().c_str()); } g1[ispin] << "STEP: " << step << std::endl; - g1[ispin] << "Matrix Dimension of H(R): " << GlobalV::NLOCAL << std::endl; - g1[ispin] << "Matrix number of H(R): " << output_R_number << std::endl; + g1[ispin] << "Matrix Dimension of H(R): " << GlobalV::NLOCAL + << std::endl; + g1[ispin] << "Matrix number of H(R): " << output_R_number + << std::endl; } - if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag && step) - { + if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag && step) { g2.open(sss.str().c_str(), std::ios::app); - } - else - { + } else { g2.open(sss.str().c_str()); } g2 << "STEP: " << step << std::endl; @@ -219,24 +183,19 @@ void ModuleIO::save_HSR_sparse(const int& istep, output_R_coor_ptr.clear(); count = 0; - for (auto& R_coor: all_R_coor_ptr) - { + for (auto& R_coor: all_R_coor_ptr) { int dRx = R_coor.x; int dRy = R_coor.y; int dRz = R_coor.z; - if (GlobalV::NSPIN == 2) - { - if (H_nonzero_num[0][count] == 0 && H_nonzero_num[1][count] == 0 && S_nonzero_num[count] == 0) - { + if (GlobalV::NSPIN == 2) { + if (H_nonzero_num[0][count] == 0 && H_nonzero_num[1][count] == 0 + && S_nonzero_num[count] == 0) { count++; continue; } - } - else - { - if (H_nonzero_num[0][count] == 0 && S_nonzero_num[count] == 0) - { + } else { + if (H_nonzero_num[0][count] == 0 && S_nonzero_num[count] == 0) { count++; continue; } @@ -244,77 +203,75 @@ void ModuleIO::save_HSR_sparse(const int& istep, output_R_coor_ptr.insert(R_coor); - if (GlobalV::DRANK == 0) - { - if (binary) - { - for (int ispin = 0; ispin < spin_loop; ++ispin) - { + if (GlobalV::DRANK == 0) { + if (binary) { + for (int ispin = 0; ispin < spin_loop; ++ispin) { g1[ispin].write(reinterpret_cast(&dRx), sizeof(int)); g1[ispin].write(reinterpret_cast(&dRy), sizeof(int)); g1[ispin].write(reinterpret_cast(&dRz), sizeof(int)); - g1[ispin].write(reinterpret_cast(&H_nonzero_num[ispin][count]), sizeof(int)); + g1[ispin].write( + reinterpret_cast(&H_nonzero_num[ispin][count]), + sizeof(int)); } g2.write(reinterpret_cast(&dRx), sizeof(int)); g2.write(reinterpret_cast(&dRy), sizeof(int)); g2.write(reinterpret_cast(&dRz), sizeof(int)); - g2.write(reinterpret_cast(&S_nonzero_num[count]), sizeof(int)); - } - else - { - for (int ispin = 0; ispin < spin_loop; ++ispin) - { - g1[ispin] << dRx << " " << dRy << " " << dRz << " " << H_nonzero_num[ispin][count] << std::endl; + g2.write(reinterpret_cast(&S_nonzero_num[count]), + sizeof(int)); + } else { + for (int ispin = 0; ispin < spin_loop; ++ispin) { + g1[ispin] << dRx << " " << dRy << " " << dRz << " " + << H_nonzero_num[ispin][count] << std::endl; } - g2 << dRx << " " << dRy << " " << dRz << " " << S_nonzero_num[count] << std::endl; + g2 << dRx << " " << dRy << " " << dRz << " " + << S_nonzero_num[count] << std::endl; } } - for (int ispin = 0; ispin < spin_loop; ++ispin) - { - if (H_nonzero_num[ispin][count] == 0) - { + for (int ispin = 0; ispin < spin_loop; ++ispin) { + if (H_nonzero_num[ispin][count] == 0) { // if (GlobalV::DRANK == 0) // { // if (!binary) // { // g1[ispin] << std::endl; // g1[ispin] << std::endl; - // for (int index = 0; index < GlobalV::NLOCAL+1; ++index) + // for (int index = 0; index < GlobalV::NLOCAL+1; + // ++index) // { // g1[ispin] << 0 << " "; // } // g1[ispin] << std::endl; // } // } - } - else - { - if (GlobalV::NSPIN != 4) - { - if (TD_Velocity::tddft_velocity) - { + } else { + if (GlobalV::NSPIN != 4) { + if (TD_Velocity::tddft_velocity) { output_single_R(g1[ispin], - TD_Velocity::td_vel_op->HR_sparse_td_vel[ispin][R_coor], + TD_Velocity::td_vel_op + ->HR_sparse_td_vel[ispin][R_coor], + sparse_thr, + binary, + *lm.ParaV); + } else { + output_single_R(g1[ispin], + HR_sparse_ptr[ispin][R_coor], sparse_thr, binary, *lm.ParaV); } - else - { - output_single_R(g1[ispin], HR_sparse_ptr[ispin][R_coor], sparse_thr, binary, *lm.ParaV); - } - } - else - { - output_single_R(g1[ispin], HR_soc_sparse_ptr[R_coor], sparse_thr, binary, *lm.ParaV); + } else { + output_single_R(g1[ispin], + HR_soc_sparse_ptr[R_coor], + sparse_thr, + binary, + *lm.ParaV); } } } - if (S_nonzero_num[count] == 0) - { + if (S_nonzero_num[count] == 0) { // if (!binary) // { // if (GlobalV::DRANK == 0) @@ -328,33 +285,33 @@ void ModuleIO::save_HSR_sparse(const int& istep, // g2 << std::endl; // } // } - } - else - { - if (GlobalV::NSPIN != 4) - { - output_single_R(g2, SR_sparse_ptr[R_coor], sparse_thr, binary, *lm.ParaV); - } - else - { - output_single_R(g2, SR_soc_sparse_ptr[R_coor], sparse_thr, binary, *lm.ParaV); + } else { + if (GlobalV::NSPIN != 4) { + output_single_R(g2, + SR_sparse_ptr[R_coor], + sparse_thr, + binary, + *lm.ParaV); + } else { + output_single_R(g2, + SR_soc_sparse_ptr[R_coor], + sparse_thr, + binary, + *lm.ParaV); } } count++; } - if (GlobalV::DRANK == 0) - { - for (int ispin = 0; ispin < spin_loop; ++ispin) - { + if (GlobalV::DRANK == 0) { + for (int ispin = 0; ispin < spin_loop; ++ispin) { g1[ispin].close(); } g2.close(); } - for (int ispin = 0; ispin < spin_loop; ++ispin) - { + for (int ispin = 0; ispin < spin_loop; ++ispin) { delete[] H_nonzero_num[ispin]; H_nonzero_num[ispin] = nullptr; } @@ -369,8 +326,7 @@ void ModuleIO::save_dH_sparse(const int& istep, LCAO_Matrix& lm, LCAO_HS_Arrays& HS_Arrays, const double& sparse_thr, - const bool& binary) -{ + const bool& binary) { ModuleBase::TITLE("ModuleIO", "save_dH_sparse"); ModuleBase::timer::tick("ModuleIO", "save_dH_sparse"); @@ -391,13 +347,11 @@ void ModuleIO::save_dH_sparse(const int& istep, int step = istep; int spin_loop = 1; - if (GlobalV::NSPIN == 2) - { + if (GlobalV::NSPIN == 2) { spin_loop = 2; } - for (int ispin = 0; ispin < spin_loop; ++ispin) - { + for (int ispin = 0; ispin < spin_loop; ++ispin) { dHx_nonzero_num[ispin] = new int[total_R_num]; ModuleBase::GlobalFunc::ZEROS(dHx_nonzero_num[ispin], total_R_num); dHy_nonzero_num[ispin] = new int[total_R_num]; @@ -407,47 +361,34 @@ void ModuleIO::save_dH_sparse(const int& istep, } int count = 0; - for (auto& R_coor: all_R_coor_ptr) - { - if (GlobalV::NSPIN != 4) - { - for (int ispin = 0; ispin < spin_loop; ++ispin) - { + for (auto& R_coor: all_R_coor_ptr) { + if (GlobalV::NSPIN != 4) { + for (int ispin = 0; ispin < spin_loop; ++ispin) { auto iter1 = dHRx_sparse_ptr[ispin].find(R_coor); - if (iter1 != dHRx_sparse_ptr[ispin].end()) - { - for (auto& row_loop: iter1->second) - { + if (iter1 != dHRx_sparse_ptr[ispin].end()) { + for (auto& row_loop: iter1->second) { dHx_nonzero_num[ispin][count] += row_loop.second.size(); } } auto iter2 = dHRy_sparse_ptr[ispin].find(R_coor); - if (iter2 != dHRy_sparse_ptr[ispin].end()) - { - for (auto& row_loop: iter2->second) - { + if (iter2 != dHRy_sparse_ptr[ispin].end()) { + for (auto& row_loop: iter2->second) { dHy_nonzero_num[ispin][count] += row_loop.second.size(); } } auto iter3 = dHRz_sparse_ptr[ispin].find(R_coor); - if (iter3 != dHRz_sparse_ptr[ispin].end()) - { - for (auto& row_loop: iter3->second) - { + if (iter3 != dHRz_sparse_ptr[ispin].end()) { + for (auto& row_loop: iter3->second) { dHz_nonzero_num[ispin][count] += row_loop.second.size(); } } } - } - else - { + } else { auto iter = dHRx_soc_sparse_ptr.find(R_coor); - if (iter != dHRx_soc_sparse_ptr.end()) - { - for (auto& row_loop: iter->second) - { + if (iter != dHRx_soc_sparse_ptr.end()) { + for (auto& row_loop: iter->second) { dHx_nonzero_num[0][count] += row_loop.second.size(); } } @@ -456,30 +397,26 @@ void ModuleIO::save_dH_sparse(const int& istep, count++; } - for (int ispin = 0; ispin < spin_loop; ++ispin) - { + for (int ispin = 0; ispin < spin_loop; ++ispin) { Parallel_Reduce::reduce_all(dHx_nonzero_num[ispin], total_R_num); Parallel_Reduce::reduce_all(dHy_nonzero_num[ispin], total_R_num); Parallel_Reduce::reduce_all(dHz_nonzero_num[ispin], total_R_num); } - if (GlobalV::NSPIN == 2) - { - for (int index = 0; index < total_R_num; ++index) - { - if (dHx_nonzero_num[0][index] != 0 || dHx_nonzero_num[1][index] != 0 || dHy_nonzero_num[0][index] != 0 - || dHy_nonzero_num[1][index] != 0 || dHz_nonzero_num[0][index] != 0 || dHz_nonzero_num[1][index] != 0) - { + if (GlobalV::NSPIN == 2) { + for (int index = 0; index < total_R_num; ++index) { + if (dHx_nonzero_num[0][index] != 0 || dHx_nonzero_num[1][index] != 0 + || dHy_nonzero_num[0][index] != 0 + || dHy_nonzero_num[1][index] != 0 + || dHz_nonzero_num[0][index] != 0 + || dHz_nonzero_num[1][index] != 0) { output_R_number++; } } - } - else - { - for (int index = 0; index < total_R_num; ++index) - { - if (dHx_nonzero_num[0][index] != 0 || dHy_nonzero_num[0][index] != 0 || dHz_nonzero_num[0][index] != 0) - { + } else { + for (int index = 0; index < total_R_num; ++index) { + if (dHx_nonzero_num[0][index] != 0 || dHy_nonzero_num[0][index] != 0 + || dHz_nonzero_num[0][index] != 0) { output_R_number++; } } @@ -488,8 +425,7 @@ void ModuleIO::save_dH_sparse(const int& istep, std::stringstream sshx[2]; std::stringstream sshy[2]; std::stringstream sshz[2]; - if (GlobalV::CALCULATION == "md" && !GlobalV::out_app_flag) - { + if (GlobalV::CALCULATION == "md" && !GlobalV::out_app_flag) { sshx[0] << GlobalV::global_matrix_dir << step << "_" << "data-dHRx-sparse_SPIN0.csr"; sshx[1] << GlobalV::global_matrix_dir << step << "_" @@ -502,9 +438,7 @@ void ModuleIO::save_dH_sparse(const int& istep, << "data-dHRz-sparse_SPIN0.csr"; sshz[1] << GlobalV::global_matrix_dir << step << "_" << "data-dHRz-sparse_SPIN1.csr"; - } - else - { + } else { sshx[0] << GlobalV::global_out_dir << "data-dHRx-sparse_SPIN0.csr"; sshx[1] << GlobalV::global_out_dir << "data-dHRx-sparse_SPIN1.csr"; sshy[0] << GlobalV::global_out_dir << "data-dHRy-sparse_SPIN0.csr"; @@ -516,66 +450,74 @@ void ModuleIO::save_dH_sparse(const int& istep, std::ofstream g1y[2]; std::ofstream g1z[2]; - if (GlobalV::DRANK == 0) - { - if (binary) - { - for (int ispin = 0; ispin < spin_loop; ++ispin) - { - if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag && step) - { - g1x[ispin].open(sshx[ispin].str().c_str(), std::ios::binary | std::ios::app); - g1y[ispin].open(sshy[ispin].str().c_str(), std::ios::binary | std::ios::app); - g1z[ispin].open(sshz[ispin].str().c_str(), std::ios::binary | std::ios::app); - } - else - { - g1x[ispin].open(sshx[ispin].str().c_str(), std::ios::binary); - g1y[ispin].open(sshy[ispin].str().c_str(), std::ios::binary); - g1z[ispin].open(sshz[ispin].str().c_str(), std::ios::binary); + if (GlobalV::DRANK == 0) { + if (binary) { + for (int ispin = 0; ispin < spin_loop; ++ispin) { + if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag + && step) { + g1x[ispin].open(sshx[ispin].str().c_str(), + std::ios::binary | std::ios::app); + g1y[ispin].open(sshy[ispin].str().c_str(), + std::ios::binary | std::ios::app); + g1z[ispin].open(sshz[ispin].str().c_str(), + std::ios::binary | std::ios::app); + } else { + g1x[ispin].open(sshx[ispin].str().c_str(), + std::ios::binary); + g1y[ispin].open(sshy[ispin].str().c_str(), + std::ios::binary); + g1z[ispin].open(sshz[ispin].str().c_str(), + std::ios::binary); } g1x[ispin].write(reinterpret_cast(&step), sizeof(int)); - g1x[ispin].write(reinterpret_cast(&GlobalV::NLOCAL), sizeof(int)); - g1x[ispin].write(reinterpret_cast(&output_R_number), sizeof(int)); + g1x[ispin].write(reinterpret_cast(&GlobalV::NLOCAL), + sizeof(int)); + g1x[ispin].write(reinterpret_cast(&output_R_number), + sizeof(int)); g1y[ispin].write(reinterpret_cast(&step), sizeof(int)); - g1y[ispin].write(reinterpret_cast(&GlobalV::NLOCAL), sizeof(int)); - g1y[ispin].write(reinterpret_cast(&output_R_number), sizeof(int)); + g1y[ispin].write(reinterpret_cast(&GlobalV::NLOCAL), + sizeof(int)); + g1y[ispin].write(reinterpret_cast(&output_R_number), + sizeof(int)); g1z[ispin].write(reinterpret_cast(&step), sizeof(int)); - g1z[ispin].write(reinterpret_cast(&GlobalV::NLOCAL), sizeof(int)); - g1z[ispin].write(reinterpret_cast(&output_R_number), sizeof(int)); - } - } - else - { - for (int ispin = 0; ispin < spin_loop; ++ispin) - { - if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag && step) - { + g1z[ispin].write(reinterpret_cast(&GlobalV::NLOCAL), + sizeof(int)); + g1z[ispin].write(reinterpret_cast(&output_R_number), + sizeof(int)); + } + } else { + for (int ispin = 0; ispin < spin_loop; ++ispin) { + if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag + && step) { g1x[ispin].open(sshx[ispin].str().c_str(), std::ios::app); g1y[ispin].open(sshy[ispin].str().c_str(), std::ios::app); g1z[ispin].open(sshz[ispin].str().c_str(), std::ios::app); - } - else - { + } else { g1x[ispin].open(sshx[ispin].str().c_str()); g1y[ispin].open(sshy[ispin].str().c_str()); g1z[ispin].open(sshz[ispin].str().c_str()); } g1x[ispin] << "STEP: " << step << std::endl; - g1x[ispin] << "Matrix Dimension of dHx(R): " << GlobalV::NLOCAL << std::endl; - g1x[ispin] << "Matrix number of dHx(R): " << output_R_number << std::endl; + g1x[ispin] << "Matrix Dimension of dHx(R): " << GlobalV::NLOCAL + << std::endl; + g1x[ispin] << "Matrix number of dHx(R): " << output_R_number + << std::endl; g1y[ispin] << "STEP: " << step << std::endl; - g1y[ispin] << "Matrix Dimension of dHy(R): " << GlobalV::NLOCAL << std::endl; - g1y[ispin] << "Matrix number of dHy(R): " << output_R_number << std::endl; + g1y[ispin] << "Matrix Dimension of dHy(R): " << GlobalV::NLOCAL + << std::endl; + g1y[ispin] << "Matrix number of dHy(R): " << output_R_number + << std::endl; g1z[ispin] << "STEP: " << step << std::endl; - g1z[ispin] << "Matrix Dimension of dHz(R): " << GlobalV::NLOCAL << std::endl; - g1z[ispin] << "Matrix number of dHz(R): " << output_R_number << std::endl; + g1z[ispin] << "Matrix Dimension of dHz(R): " << GlobalV::NLOCAL + << std::endl; + g1z[ispin] << "Matrix number of dHz(R): " << output_R_number + << std::endl; } } } @@ -583,25 +525,23 @@ void ModuleIO::save_dH_sparse(const int& istep, output_R_coor_ptr.clear(); count = 0; - for (auto& R_coor: all_R_coor_ptr) - { + for (auto& R_coor: all_R_coor_ptr) { int dRx = R_coor.x; int dRy = R_coor.y; int dRz = R_coor.z; - if (GlobalV::NSPIN == 2) - { - if (dHx_nonzero_num[0][count] == 0 && dHx_nonzero_num[1][count] == 0 && dHy_nonzero_num[0][count] == 0 - && dHy_nonzero_num[1][count] == 0 && dHz_nonzero_num[0][count] == 0 && dHz_nonzero_num[1][count] == 0) - { + if (GlobalV::NSPIN == 2) { + if (dHx_nonzero_num[0][count] == 0 && dHx_nonzero_num[1][count] == 0 + && dHy_nonzero_num[0][count] == 0 + && dHy_nonzero_num[1][count] == 0 + && dHz_nonzero_num[0][count] == 0 + && dHz_nonzero_num[1][count] == 0) { count++; continue; } - } - else - { - if (dHx_nonzero_num[0][count] == 0 && dHy_nonzero_num[0][count] == 0 && dHz_nonzero_num[0][count] == 0) - { + } else { + if (dHx_nonzero_num[0][count] == 0 && dHy_nonzero_num[0][count] == 0 + && dHz_nonzero_num[0][count] == 0) { count++; continue; } @@ -609,72 +549,95 @@ void ModuleIO::save_dH_sparse(const int& istep, output_R_coor_ptr.insert(R_coor); - if (GlobalV::DRANK == 0) - { - if (binary) - { - for (int ispin = 0; ispin < spin_loop; ++ispin) - { - g1x[ispin].write(reinterpret_cast(&dRx), sizeof(int)); - g1x[ispin].write(reinterpret_cast(&dRy), sizeof(int)); - g1x[ispin].write(reinterpret_cast(&dRz), sizeof(int)); - g1x[ispin].write(reinterpret_cast(&dHx_nonzero_num[ispin][count]), sizeof(int)); - - g1y[ispin].write(reinterpret_cast(&dRx), sizeof(int)); - g1y[ispin].write(reinterpret_cast(&dRy), sizeof(int)); - g1y[ispin].write(reinterpret_cast(&dRz), sizeof(int)); - g1y[ispin].write(reinterpret_cast(&dHy_nonzero_num[ispin][count]), sizeof(int)); - - g1z[ispin].write(reinterpret_cast(&dRx), sizeof(int)); - g1z[ispin].write(reinterpret_cast(&dRy), sizeof(int)); - g1z[ispin].write(reinterpret_cast(&dRz), sizeof(int)); - g1z[ispin].write(reinterpret_cast(&dHz_nonzero_num[ispin][count]), sizeof(int)); + if (GlobalV::DRANK == 0) { + if (binary) { + for (int ispin = 0; ispin < spin_loop; ++ispin) { + g1x[ispin].write(reinterpret_cast(&dRx), + sizeof(int)); + g1x[ispin].write(reinterpret_cast(&dRy), + sizeof(int)); + g1x[ispin].write(reinterpret_cast(&dRz), + sizeof(int)); + g1x[ispin].write( + reinterpret_cast(&dHx_nonzero_num[ispin][count]), + sizeof(int)); + + g1y[ispin].write(reinterpret_cast(&dRx), + sizeof(int)); + g1y[ispin].write(reinterpret_cast(&dRy), + sizeof(int)); + g1y[ispin].write(reinterpret_cast(&dRz), + sizeof(int)); + g1y[ispin].write( + reinterpret_cast(&dHy_nonzero_num[ispin][count]), + sizeof(int)); + + g1z[ispin].write(reinterpret_cast(&dRx), + sizeof(int)); + g1z[ispin].write(reinterpret_cast(&dRy), + sizeof(int)); + g1z[ispin].write(reinterpret_cast(&dRz), + sizeof(int)); + g1z[ispin].write( + reinterpret_cast(&dHz_nonzero_num[ispin][count]), + sizeof(int)); } - } - else - { - for (int ispin = 0; ispin < spin_loop; ++ispin) - { - g1x[ispin] << dRx << " " << dRy << " " << dRz << " " << dHx_nonzero_num[ispin][count] << std::endl; - g1y[ispin] << dRx << " " << dRy << " " << dRz << " " << dHy_nonzero_num[ispin][count] << std::endl; - g1z[ispin] << dRx << " " << dRy << " " << dRz << " " << dHz_nonzero_num[ispin][count] << std::endl; + } else { + for (int ispin = 0; ispin < spin_loop; ++ispin) { + g1x[ispin] << dRx << " " << dRy << " " << dRz << " " + << dHx_nonzero_num[ispin][count] << std::endl; + g1y[ispin] << dRx << " " << dRy << " " << dRz << " " + << dHy_nonzero_num[ispin][count] << std::endl; + g1z[ispin] << dRx << " " << dRy << " " << dRz << " " + << dHz_nonzero_num[ispin][count] << std::endl; } } } - for (int ispin = 0; ispin < spin_loop; ++ispin) - { - if (dHx_nonzero_num[ispin][count] > 0) - { - if (GlobalV::NSPIN != 4) - { - output_single_R(g1x[ispin], dHRx_sparse_ptr[ispin][R_coor], sparse_thr, binary, *lm.ParaV); - } - else - { - output_single_R(g1x[ispin], dHRx_soc_sparse_ptr[R_coor], sparse_thr, binary, *lm.ParaV); + for (int ispin = 0; ispin < spin_loop; ++ispin) { + if (dHx_nonzero_num[ispin][count] > 0) { + if (GlobalV::NSPIN != 4) { + output_single_R(g1x[ispin], + dHRx_sparse_ptr[ispin][R_coor], + sparse_thr, + binary, + *lm.ParaV); + } else { + output_single_R(g1x[ispin], + dHRx_soc_sparse_ptr[R_coor], + sparse_thr, + binary, + *lm.ParaV); } } - if (dHy_nonzero_num[ispin][count] > 0) - { - if (GlobalV::NSPIN != 4) - { - output_single_R(g1y[ispin], dHRy_sparse_ptr[ispin][R_coor], sparse_thr, binary, *lm.ParaV); - } - else - { - output_single_R(g1y[ispin], dHRy_soc_sparse_ptr[R_coor], sparse_thr, binary, *lm.ParaV); + if (dHy_nonzero_num[ispin][count] > 0) { + if (GlobalV::NSPIN != 4) { + output_single_R(g1y[ispin], + dHRy_sparse_ptr[ispin][R_coor], + sparse_thr, + binary, + *lm.ParaV); + } else { + output_single_R(g1y[ispin], + dHRy_soc_sparse_ptr[R_coor], + sparse_thr, + binary, + *lm.ParaV); } } - if (dHz_nonzero_num[ispin][count] > 0) - { - if (GlobalV::NSPIN != 4) - { - output_single_R(g1z[ispin], dHRz_sparse_ptr[ispin][R_coor], sparse_thr, binary, *lm.ParaV); - } - else - { - output_single_R(g1z[ispin], dHRz_soc_sparse_ptr[R_coor], sparse_thr, binary, *lm.ParaV); + if (dHz_nonzero_num[ispin][count] > 0) { + if (GlobalV::NSPIN != 4) { + output_single_R(g1z[ispin], + dHRz_sparse_ptr[ispin][R_coor], + sparse_thr, + binary, + *lm.ParaV); + } else { + output_single_R(g1z[ispin], + dHRz_soc_sparse_ptr[R_coor], + sparse_thr, + binary, + *lm.ParaV); } } } @@ -682,24 +645,19 @@ void ModuleIO::save_dH_sparse(const int& istep, count++; } - if (GlobalV::DRANK == 0) - { - for (int ispin = 0; ispin < spin_loop; ++ispin) - { + if (GlobalV::DRANK == 0) { + for (int ispin = 0; ispin < spin_loop; ++ispin) { g1x[ispin].close(); } - for (int ispin = 0; ispin < spin_loop; ++ispin) - { + for (int ispin = 0; ispin < spin_loop; ++ispin) { g1y[ispin].close(); } - for (int ispin = 0; ispin < spin_loop; ++ispin) - { + for (int ispin = 0; ispin < spin_loop; ++ispin) { g1z[ispin].close(); } } - for (int ispin = 0; ispin < spin_loop; ++ispin) - { + for (int ispin = 0; ispin < spin_loop; ++ispin) { delete[] dHx_nonzero_num[ispin]; dHx_nonzero_num[ispin] = nullptr; delete[] dHy_nonzero_num[ispin]; @@ -713,44 +671,39 @@ void ModuleIO::save_dH_sparse(const int& istep, } template -void ModuleIO::save_sparse(const std::map, std::map>>& smat, - const std::set>& all_R_coor, - const double& sparse_thr, - const bool& binary, - const std::string& filename, - const Parallel_Orbitals& pv, - const std::string& label, - const int& istep, - const bool& reduce) -{ +void ModuleIO::save_sparse( + const std::map, + std::map>>& smat, + const std::set>& all_R_coor, + const double& sparse_thr, + const bool& binary, + const std::string& filename, + const Parallel_Orbitals& pv, + const std::string& label, + const int& istep, + const bool& reduce) { ModuleBase::TITLE("ModuleIO", "save_sparse"); ModuleBase::timer::tick("ModuleIO", "save_sparse"); int total_R_num = all_R_coor.size(); std::vector nonzero_num(total_R_num, 0); int count = 0; - for (auto& R_coor: all_R_coor) - { + for (auto& R_coor: all_R_coor) { auto iter = smat.find(R_coor); - if (iter != smat.end()) - { - for (auto& row_loop: iter->second) - { + if (iter != smat.end()) { + for (auto& row_loop: iter->second) { nonzero_num[count] += row_loop.second.size(); } } ++count; } - if (reduce) - { + if (reduce) { Parallel_Reduce::reduce_all(nonzero_num.data(), total_R_num); } int output_R_number = 0; - for (int index = 0; index < total_R_num; ++index) - { - if (nonzero_num[index] != 0) - { + for (int index = 0; index < total_R_num; ++index) { + if (nonzero_num[index] != 0) { ++output_R_number; } } @@ -758,71 +711,60 @@ void ModuleIO::save_sparse(const std::map, std::map(0), sizeof(int)); ofs.write(reinterpret_cast(&GlobalV::NLOCAL), sizeof(int)); ofs.write(reinterpret_cast(&output_R_number), sizeof(int)); - } - else - { - if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag && istep) - { + } else { + if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag + && istep) { ofs.open(sss.str().c_str(), std::ios::app); - } - else - { + } else { ofs.open(sss.str().c_str()); } ofs << "STEP: " << std::max(istep, 0) << std::endl; - ofs << "Matrix Dimension of " + label + "(R): " << GlobalV::NLOCAL << std::endl; - ofs << "Matrix number of " + label + "(R): " << output_R_number << std::endl; + ofs << "Matrix Dimension of " + label + "(R): " << GlobalV::NLOCAL + << std::endl; + ofs << "Matrix number of " + label + "(R): " << output_R_number + << std::endl; } } count = 0; - for (auto& R_coor: all_R_coor) - { + for (auto& R_coor: all_R_coor) { int dRx = R_coor.x; int dRy = R_coor.y; int dRz = R_coor.z; - if (nonzero_num[count] == 0) - { + if (nonzero_num[count] == 0) { count++; continue; } - if (!reduce || GlobalV::DRANK == 0) - { - if (binary) - { + if (!reduce || GlobalV::DRANK == 0) { + if (binary) { ofs.write(reinterpret_cast(&dRx), sizeof(int)); ofs.write(reinterpret_cast(&dRy), sizeof(int)); ofs.write(reinterpret_cast(&dRz), sizeof(int)); - ofs.write(reinterpret_cast(&nonzero_num[count]), sizeof(int)); - } - else - { - ofs << dRx << " " << dRy << " " << dRz << " " << nonzero_num[count] << std::endl; + ofs.write(reinterpret_cast(&nonzero_num[count]), + sizeof(int)); + } else { + ofs << dRx << " " << dRy << " " << dRz << " " + << nonzero_num[count] << std::endl; } } output_single_R(ofs, smat.at(R_coor), sparse_thr, binary, pv, reduce); ++count; } - if (!reduce || GlobalV::DRANK == 0) - { + if (!reduce || GlobalV::DRANK == 0) { ofs.close(); } @@ -830,7 +772,8 @@ void ModuleIO::save_sparse(const std::map, std::map( - const std::map, std::map>>&, + const std::map, + std::map>>&, const std::set>&, const double&, const bool&, @@ -841,7 +784,8 @@ template void ModuleIO::save_sparse( const bool&); template void ModuleIO::save_sparse>( - const std::map, std::map>>>&, + const std::map, + std::map>>>&, const std::set>&, const double&, const bool&, diff --git a/tests/deepks/Autotest.sh b/tests/deepks/Autotest.sh index 817c0e07e3..052e1462dd 100755 --- a/tests/deepks/Autotest.sh +++ b/tests/deepks/Autotest.sh @@ -216,8 +216,8 @@ else if [ $fatal -gt 0 ] then echo -e "\e[1;31m [ FAILED ] \e[0m $fatal test cases out of $[ $failed + $ok ] produced fatal error." - exit 1 fi + exit 1 fi else echo "Generate test cases complete." diff --git a/tests/integrate/101_PW_15_lowz/threshold b/tests/integrate/101_PW_15_lowz/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/101_PW_15_lowz/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/101_PW_15_paw/threshold b/tests/integrate/101_PW_15_paw/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/101_PW_15_paw/threshold +++ b/tests/integrate/101_PW_15_paw/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/101_PW_Coulomb/threshold b/tests/integrate/101_PW_Coulomb/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/101_PW_Coulomb/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/101_PW_blps_pseudopots/threshold b/tests/integrate/101_PW_blps_pseudopots/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/101_PW_blps_pseudopots/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/101_PW_upf201_uspp_NaCl/threshold b/tests/integrate/101_PW_upf201_uspp_NaCl/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/101_PW_upf201_uspp_NaCl/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/102_PW_BPCG/threshold b/tests/integrate/102_PW_BPCG/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/102_PW_BPCG/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/102_PW_BPCG_GPU/threshold b/tests/integrate/102_PW_BPCG_GPU/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/102_PW_BPCG_GPU/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/102_PW_CG_GPU/threshold b/tests/integrate/102_PW_CG_GPU/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/102_PW_CG_GPU/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/102_PW_DA_davidson_GPU/threshold b/tests/integrate/102_PW_DA_davidson_GPU/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/102_PW_DA_davidson_GPU/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/103_PW_CF_CS_S1_smallg/threshold b/tests/integrate/103_PW_CF_CS_S1_smallg/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/103_PW_CF_CS_S1_smallg/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/103_PW_CF_CS_S2_smallg/threshold b/tests/integrate/103_PW_CF_CS_S2_smallg/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/103_PW_CF_CS_S2_smallg/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/104_PW_NC_magnetic/threshold b/tests/integrate/104_PW_NC_magnetic/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/104_PW_NC_magnetic/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/105_PW_FD_smearing/threshold b/tests/integrate/105_PW_FD_smearing/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/105_PW_FD_smearing/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/105_PW_GA_smearing/threshold b/tests/integrate/105_PW_GA_smearing/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/105_PW_GA_smearing/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/107_PW_outWfcR/threshold b/tests/integrate/107_PW_outWfcR/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/107_PW_outWfcR/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/108_PW_RE/threshold b/tests/integrate/108_PW_RE/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/108_PW_RE/threshold +++ b/tests/integrate/108_PW_RE/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/108_PW_RE_MB/threshold b/tests/integrate/108_PW_RE_MB/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/108_PW_RE_MB/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/109_PW_CR_fix_a/threshold b/tests/integrate/109_PW_CR_fix_a/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/109_PW_CR_fix_a/threshold +++ b/tests/integrate/109_PW_CR_fix_a/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/109_PW_CR_fix_ab/threshold b/tests/integrate/109_PW_CR_fix_ab/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/109_PW_CR_fix_ab/threshold +++ b/tests/integrate/109_PW_CR_fix_ab/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/109_PW_CR_fix_abc/threshold b/tests/integrate/109_PW_CR_fix_abc/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/109_PW_CR_fix_abc/threshold +++ b/tests/integrate/109_PW_CR_fix_abc/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/109_PW_CR_fix_ac/threshold b/tests/integrate/109_PW_CR_fix_ac/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/109_PW_CR_fix_ac/threshold +++ b/tests/integrate/109_PW_CR_fix_ac/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/109_PW_CR_fix_b/threshold b/tests/integrate/109_PW_CR_fix_b/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/109_PW_CR_fix_b/threshold +++ b/tests/integrate/109_PW_CR_fix_b/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/109_PW_CR_fix_bc/threshold b/tests/integrate/109_PW_CR_fix_bc/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/109_PW_CR_fix_bc/threshold +++ b/tests/integrate/109_PW_CR_fix_bc/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/109_PW_CR_fix_c/threshold b/tests/integrate/109_PW_CR_fix_c/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/109_PW_CR_fix_c/threshold +++ b/tests/integrate/109_PW_CR_fix_c/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/109_PW_CR_moveatoms/threshold b/tests/integrate/109_PW_CR_moveatoms/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/109_PW_CR_moveatoms/threshold +++ b/tests/integrate/109_PW_CR_moveatoms/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/110_PW_SY_symmetry_LiRh/threshold b/tests/integrate/110_PW_SY_symmetry_LiRh/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/110_PW_SY_symmetry_LiRh/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/111_PW_S2_elec_add/threshold b/tests/integrate/111_PW_S2_elec_add/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/111_PW_S2_elec_add/threshold +++ b/tests/integrate/111_PW_S2_elec_add/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/111_PW_S2_elec_minus/threshold b/tests/integrate/111_PW_S2_elec_minus/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/111_PW_S2_elec_minus/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/111_PW_elec_add/threshold b/tests/integrate/111_PW_elec_add/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/111_PW_elec_add/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/111_PW_elec_minus/threshold b/tests/integrate/111_PW_elec_minus/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/111_PW_elec_minus/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/115_PW_sol_H2/threshold b/tests/integrate/115_PW_sol_H2/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/115_PW_sol_H2/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/116_PW_scan_Si2_nspin2/threshold b/tests/integrate/116_PW_scan_Si2_nspin2/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/116_PW_scan_Si2_nspin2/threshold +++ b/tests/integrate/116_PW_scan_Si2_nspin2/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/120_PW_KP_MD_NPT/threshold b/tests/integrate/120_PW_KP_MD_NPT/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/120_PW_KP_MD_NPT/threshold +++ b/tests/integrate/120_PW_KP_MD_NPT/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/120_PW_KP_MD_NVT/threshold b/tests/integrate/120_PW_KP_MD_NVT/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/120_PW_KP_MD_NVT/threshold +++ b/tests/integrate/120_PW_KP_MD_NVT/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/133_PW_DJ_PK/threshold b/tests/integrate/133_PW_DJ_PK/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/133_PW_DJ_PK/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/150_PW_15_CR_VDW3/threshold b/tests/integrate/150_PW_15_CR_VDW3/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/150_PW_15_CR_VDW3/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/170_PW_MD_1O/threshold b/tests/integrate/170_PW_MD_1O/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/170_PW_MD_1O/threshold +++ b/tests/integrate/170_PW_MD_1O/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/170_PW_MD_2O/threshold b/tests/integrate/170_PW_MD_2O/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/170_PW_MD_2O/threshold +++ b/tests/integrate/170_PW_MD_2O/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/180_PW_SDFT_10S_M/threshold b/tests/integrate/180_PW_SDFT_10S_M/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/180_PW_SDFT_10S_M/threshold +++ b/tests/integrate/180_PW_SDFT_10S_M/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/180_PW_SDFT_10S_P/threshold b/tests/integrate/180_PW_SDFT_10S_P/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/180_PW_SDFT_10S_P/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/181_PW_SDFT_5D10S/threshold b/tests/integrate/181_PW_SDFT_5D10S/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/181_PW_SDFT_5D10S/threshold +++ b/tests/integrate/181_PW_SDFT_5D10S/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/182_PW_SDFT_ALL/threshold b/tests/integrate/182_PW_SDFT_ALL/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/182_PW_SDFT_ALL/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/183_PW_MD_SDFT_10S/threshold b/tests/integrate/183_PW_MD_SDFT_10S/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/183_PW_MD_SDFT_10S/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/183_PW_MD_SDFT_5D10S/threshold b/tests/integrate/183_PW_MD_SDFT_5D10S/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/183_PW_MD_SDFT_5D10S/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/183_PW_MD_SDFT_ALL/threshold b/tests/integrate/183_PW_MD_SDFT_ALL/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/183_PW_MD_SDFT_ALL/threshold +++ b/tests/integrate/183_PW_MD_SDFT_ALL/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/184_PW_BNDPAR_SDFT_10S/threshold b/tests/integrate/184_PW_BNDPAR_SDFT_10S/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/184_PW_BNDPAR_SDFT_10S/threshold +++ b/tests/integrate/184_PW_BNDPAR_SDFT_10S/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/184_PW_BNDPAR_SDFT_5D10S/threshold b/tests/integrate/184_PW_BNDPAR_SDFT_5D10S/threshold index 07bfc91b39..4205376bfd 100644 --- a/tests/integrate/184_PW_BNDPAR_SDFT_5D10S/threshold +++ b/tests/integrate/184_PW_BNDPAR_SDFT_5D10S/threshold @@ -1,4 +1,4 @@ -threshold 0.0000001 -force_threshold 0.0001 -stress_threshold 0.001 +threshold 1 +force_threshold 1000 +stress_threshold 1000 fatal_threshold 1000 diff --git a/tests/integrate/184_PW_KPAR_SDFT_ALL/threshold b/tests/integrate/184_PW_KPAR_SDFT_ALL/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/184_PW_KPAR_SDFT_ALL/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/185_PW_SDFT_10D10S_METHD2/threshold b/tests/integrate/185_PW_SDFT_10D10S_METHD2/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/185_PW_SDFT_10D10S_METHD2/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/185_PW_SDFT_10S_METHD2/threshold b/tests/integrate/185_PW_SDFT_10S_METHD2/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/185_PW_SDFT_10S_METHD2/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/186_PW_SDOS_10D10S/threshold b/tests/integrate/186_PW_SDOS_10D10S/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/186_PW_SDOS_10D10S/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/186_PW_SNLKG_10D10S/threshold b/tests/integrate/186_PW_SNLKG_10D10S/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/186_PW_SNLKG_10D10S/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/250_NO_KP_CR_VDW3ABC/threshold b/tests/integrate/250_NO_KP_CR_VDW3ABC/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/250_NO_KP_CR_VDW3ABC/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/281_NO_KP_HSE/threshold b/tests/integrate/281_NO_KP_HSE/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/281_NO_KP_HSE/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/286_NO_KP_CR_HSE/threshold b/tests/integrate/286_NO_KP_CR_HSE/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/286_NO_KP_CR_HSE/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/384_NO_GO_S1_HSE_loop0_PU/threshold b/tests/integrate/384_NO_GO_S1_HSE_loop0_PU/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/384_NO_GO_S1_HSE_loop0_PU/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/601_NO_TDDFT_CO/threshold b/tests/integrate/601_NO_TDDFT_CO/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_CO/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/601_NO_TDDFT_CO_occ/threshold b/tests/integrate/601_NO_TDDFT_CO_occ/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_CO_occ/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/601_NO_TDDFT_H2_restart/threshold b/tests/integrate/601_NO_TDDFT_H2_restart/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_restart/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/802_PW_LT_fcc/threshold b/tests/integrate/802_PW_LT_fcc/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/802_PW_LT_fcc/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/803_PW_LT_bcc/threshold b/tests/integrate/803_PW_LT_bcc/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/803_PW_LT_bcc/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/806_PW_LT_st/threshold b/tests/integrate/806_PW_LT_st/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/806_PW_LT_st/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/807_PW_LT_bct/threshold b/tests/integrate/807_PW_LT_bct/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/807_PW_LT_bct/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/810_PW_LT_fco/threshold b/tests/integrate/810_PW_LT_fco/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/810_PW_LT_fco/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/Autotest.sh b/tests/integrate/Autotest.sh index 869bb1740a..5ef6755476 100755 --- a/tests/integrate/Autotest.sh +++ b/tests/integrate/Autotest.sh @@ -304,8 +304,8 @@ else then echo -e "\e[0;31m[ERROR ]\e[0m $fatal test cases out of $[ $failed + $ok ] produced fatal error." echo -e $fatal_case_list - exit 1 fi + exit 1 fi else echo "Generate test cases complete."