diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index a03b3c4e51a..9b8a38f5e1b 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -2374,28 +2374,30 @@ - **Type**: String - **Availability**: *[`basis_type`](#basis_type)==pw or ([`basis_type`](#basis_type)==lcao and [`calculation`](#calculation)==get_pchg)* -- **Description**: Specifies the electronic states to calculate the charge densities with state index for, using a space-separated string of 0s and 1s. Each digit in the string corresponds to a state, starting from the first state. A 1 indicates that the charge density should be calculated for that state, while a 0 means the state will be ignored. The parameter allows a compact and flexible notation (similar to ocp_set), for example the syntax 1 4*0 5*1 0 is used to denote the selection of states: 1 means calculate for the first state, 4*0 skips the next four states, 5*1 means calculate for the following five states, and the final 0 skips the next state. It's essential that the total count of states does not exceed the total number of states (nbands); otherwise, it results in an error, and the process exits. The input string must contain only numbers and the asterisk (*) for repetition, ensuring correct format and intention of state selection. The outputs comprise multiple .cube files following the naming convention pchgi[state]s[spin]k[kpoint].cube. +- **Description**: Selects electronic states for partial (band-decomposed) charge-density output using a space-separated string of `0`s and `1`s, where `1` selects a state and `0` skips it. Repetition follows the `ocp_set` syntax, for example `1 4*0 5*1 0`; the expanded list must not exceed `nbands`. Each output represents a complete one-particle state rather than its SCF occupation. The spin degeneracy is 2 for `nspin=1` and 1 for `nspin=2` or `nspin=4`. For `nspin=1`, `s1` contains the charge density. For `nspin=2`, `s1` and `s2` contain the spin-up and spin-down charge densities, respectively. For `nspin=4`, `s1`, `s2`, `s3`, and `s4` respectively contain $\rho_0$, $m_x$, $m_y$, and $m_z$. With `if_separate_k=true`, files are named `pchgi[state]s[component]k[kpoint].cube`; otherwise, the weighted k-point sum is named `pchgi[state]s[component].cube`. + + > Note: Enabling symmetry may produce unintended partial charge densities because of reduced k-point weights and real-space symmetry operations. If the desired symmetry treatment is uncertain, set `symmetry = -1`. Use the same symmetry setting as in the SCF calculation. - **Default**: none ### out_wfc_norm - **Type**: String - **Availability**: *[`basis_type`](#basis_type)==pw or ([`basis_type`](#basis_type)==lcao and [`calculation`](#calculation)==get_wf)* -- **Description**: Specifies the electronic states to calculate the real-space wave function modulus (norm, or known as the envelope function) with state index. The syntax and state selection rules are identical to out_pchg, but the output is the norm of the wave function. The outputs comprise multiple .cube files following the naming convention wfi[state]s[spin]k[kpoint].cube. +- **Description**: Selects electronic states for real-space wavefunction-modulus output using the selection syntax and complete-state normalization of `out_pchg`. For `nspin=1`, `s1` contains the wavefunction modulus. For `nspin=2`, `s1` and `s2` contain the spin-up and spin-down wavefunction moduli, respectively. For `nspin=4`, `s1` contains the total spinor modulus. Files are named `wfi[state]s[spin]k[kpoint].cube`. - **Default**: none ### out_wfc_re_im - **Type**: String - **Availability**: *[`basis_type`](#basis_type)==pw or ([`basis_type`](#basis_type)==lcao and [`calculation`](#calculation)==get_wf)* -- **Description**: Specifies the electronic states to calculate the real and imaginary parts of the wave function with state index. The syntax and state selection rules are identical to out_pchg, but the output contains both the real and imaginary components of the wave function. The outputs comprise multiple .cube files following the naming convention wfi[state]s[spin]k[kpoint][re/im].cube. +- **Description**: Selects electronic states for real-space wavefunction real- and imaginary-part output using the selection syntax and complete-state normalization of `out_pchg`. For `nspin=1`, `s1` contains the wavefunction. For `nspin=2`, `s1` and `s2` contain the spin-up and spin-down wavefunctions, respectively. For `nspin=4`, `s1` and `s2` contain the upper and lower spinor components, respectively. Files are named `wfi[state]s[spin]k[kpoint][re/im].cube`. - **Default**: none ### if_separate_k - **Type**: Boolean - **Availability**: *([`basis_type`](#basis_type)==pw and [`out_pchg`](#out_pchg)!=none) or ([`basis_type`](#basis_type)==lcao and [`calculation`](#calculation)==get_pchg and [`gamma_only`](#gamma_only)==0)* -- **Description**: Specifies whether to write the partial charge densities for all k-points to individual files or merge them. Warning: Enabling symmetry may produce unwanted results due to reduced k-point weights and symmetry operations in real space. Therefore when calculating partial charge densities, if you are not sure what you want exactly, it is strongly recommended to set symmetry = -1. It is noteworthy that your symmetry setting should remain the same as that in the SCF procedure. +- **Description**: Specifies whether to write partial charge densities for individual k-points or merge them. - **Default**: false ### out_elf diff --git a/docs/parameters.yaml b/docs/parameters.yaml index 78b5997b840..fd381df344d 100644 --- a/docs/parameters.yaml +++ b/docs/parameters.yaml @@ -3624,7 +3624,9 @@ parameters: category: Output information type: String description: | - Specifies the electronic states to calculate the charge densities with state index for, using a space-separated string of 0s and 1s. Each digit in the string corresponds to a state, starting from the first state. A 1 indicates that the charge density should be calculated for that state, while a 0 means the state will be ignored. The parameter allows a compact and flexible notation (similar to ocp_set), for example the syntax 1 4*0 5*1 0 is used to denote the selection of states: 1 means calculate for the first state, 4*0 skips the next four states, 5*1 means calculate for the following five states, and the final 0 skips the next state. It's essential that the total count of states does not exceed the total number of states (nbands); otherwise, it results in an error, and the process exits. The input string must contain only numbers and the asterisk (*) for repetition, ensuring correct format and intention of state selection. The outputs comprise multiple .cube files following the naming convention pchgi[state]s[spin]k[kpoint].cube. + Selects electronic states for partial (band-decomposed) charge-density output using a space-separated string of `0`s and `1`s, where `1` selects a state and `0` skips it. Repetition follows the `ocp_set` syntax, for example `1 4*0 5*1 0`; the expanded list must not exceed `nbands`. Each output represents a complete one-particle state rather than its SCF occupation. The spin degeneracy is 2 for `nspin=1` and 1 for `nspin=2` or `nspin=4`. For `nspin=1`, `s1` contains the charge density. For `nspin=2`, `s1` and `s2` contain the spin-up and spin-down charge densities, respectively. For `nspin=4`, `s1`, `s2`, `s3`, and `s4` respectively contain $\rho_0$, $m_x$, $m_y$, and $m_z$. With `if_separate_k=true`, files are named `pchgi[state]s[component]k[kpoint].cube`; otherwise, the weighted k-point sum is named `pchgi[state]s[component].cube`. + + [NOTE] Enabling symmetry may produce unintended partial charge densities because of reduced k-point weights and real-space symmetry operations. If the desired symmetry treatment is uncertain, set `symmetry = -1`. Use the same symmetry setting as in the SCF calculation. default_value: none unit: "" availability: basis_type==pw or (basis_type==lcao and calculation==get_pchg) @@ -3632,7 +3634,7 @@ parameters: category: Output information type: String description: | - Specifies the electronic states to calculate the real-space wave function modulus (norm, or known as the envelope function) with state index. The syntax and state selection rules are identical to out_pchg, but the output is the norm of the wave function. The outputs comprise multiple .cube files following the naming convention wfi[state]s[spin]k[kpoint].cube. + Selects electronic states for real-space wavefunction-modulus output using the selection syntax and complete-state normalization of `out_pchg`. For `nspin=1`, `s1` contains the wavefunction modulus. For `nspin=2`, `s1` and `s2` contain the spin-up and spin-down wavefunction moduli, respectively. For `nspin=4`, `s1` contains the total spinor modulus. Files are named `wfi[state]s[spin]k[kpoint].cube`. default_value: none unit: "" availability: basis_type==pw or (basis_type==lcao and calculation==get_wf) @@ -3640,7 +3642,7 @@ parameters: category: Output information type: String description: | - Specifies the electronic states to calculate the real and imaginary parts of the wave function with state index. The syntax and state selection rules are identical to out_pchg, but the output contains both the real and imaginary components of the wave function. The outputs comprise multiple .cube files following the naming convention wfi[state]s[spin]k[kpoint][re/im].cube. + Selects electronic states for real-space wavefunction real- and imaginary-part output using the selection syntax and complete-state normalization of `out_pchg`. For `nspin=1`, `s1` contains the wavefunction. For `nspin=2`, `s1` and `s2` contain the spin-up and spin-down wavefunctions, respectively. For `nspin=4`, `s1` and `s2` contain the upper and lower spinor components, respectively. Files are named `wfi[state]s[spin]k[kpoint][re/im].cube`. default_value: none unit: "" availability: basis_type==pw or (basis_type==lcao and calculation==get_wf) @@ -3648,7 +3650,7 @@ parameters: category: Output information type: Boolean description: | - Specifies whether to write the partial charge densities for all k-points to individual files or merge them. Warning: Enabling symmetry may produce unwanted results due to reduced k-point weights and symmetry operations in real space. Therefore when calculating partial charge densities, if you are not sure what you want exactly, it is strongly recommended to set symmetry = -1. It is noteworthy that your symmetry setting should remain the same as that in the SCF procedure. + Specifies whether to write partial charge densities for individual k-points or merge them. default_value: "false" unit: "" availability: "(basis_type==pw and out_pchg!=none) or (basis_type==lcao and calculation==get_pchg and gamma_only==0)" diff --git a/source/source_base/parallel_grid.cpp b/source/source_base/parallel_grid.cpp index fb43a2cb707..8d32951c304 100644 --- a/source/source_base/parallel_grid.cpp +++ b/source/source_base/parallel_grid.cpp @@ -1,9 +1,11 @@ #include "parallel_grid.h" + #include "source_base/global_function.h" #include "source_base/global_variable.h" #ifdef __MPI #include "source_base/parallel_comm.h" // use POOL_WORLD + #include #endif @@ -167,6 +169,95 @@ void Parallel_Grid::z_distribution() return; } +void Parallel_Grid::reduce_across_pools(double* data) const +{ +#ifdef __MPI + if (GlobalV::KPAR <= 1) + { + return; + } + + assert(data != nullptr); + if (KP_WORLD != MPI_COMM_NULL) + { + // Equal-sized pools give corresponding ranks identical z-slab layouts, + // so their local buffers can be summed directly without redistribution. + MPI_Allreduce(MPI_IN_PLACE, data, this->nrxx, MPI_DOUBLE, MPI_SUM, KP_WORLD); + return; + } + + // Uneven pool sizes have no KP_WORLD and may assign different z-slabs to + // corresponding ranks. Validate the local distribution before rebuilding + // a common global layout for the cross-pool reduction. + assert(!this->numz.empty()); + assert(GlobalV::MY_POOL >= 0 && GlobalV::MY_POOL < static_cast(this->numz.size())); + assert(GlobalV::RANK_IN_POOL >= 0 && GlobalV::RANK_IN_POOL < static_cast(this->numz[GlobalV::MY_POOL].size())); + assert(this->nczp == this->numz[GlobalV::MY_POOL][GlobalV::RANK_IN_POOL]); + assert(this->nrxx == this->ncxy * this->nczp); + + const int pool_size = this->nproc_in_pool[GlobalV::MY_POOL]; + std::vector receive_counts(pool_size); + std::vector displacements(pool_size); + for (int ip = 0; ip < pool_size; ++ip) + { + receive_counts[ip] = this->numz[GlobalV::MY_POOL][ip] * this->ncxy; + displacements[ip] = this->startz[GlobalV::MY_POOL][ip] * this->ncxy; + } + + std::vector local_data(this->nrxx); + // The allgather below replicates one complete pool grid on every rank in + // that pool. INT_BGROUP then sums all of those replicas, so divide each + // local slab by the pool size to make each pool contribute exactly once. + const double pool_normalization = 1.0 / static_cast(pool_size); + for (int ir = 0; ir < this->nrxx; ++ir) + { + local_data[ir] = data[ir] * pool_normalization; + } + + std::vector pool_data(this->ncxyz); + // Collect the rank-local [xy][local_z] slabs into rank-contiguous blocks. + MPI_Allgatherv(local_data.data(), + this->nrxx, + MPI_DOUBLE, + pool_data.data(), + receive_counts.data(), + displacements.data(), + MPI_DOUBLE, + POOL_WORLD); + + std::vector global_layout(this->ncxyz); + // Convert the rank-contiguous allgather result to the canonical + // [xy][global_z] order required for element-wise reduction across pools. + for (int ip = 0; ip < pool_size; ++ip) + { + const int local_nz = this->numz[GlobalV::MY_POOL][ip]; + const int global_z_start = this->startz[GlobalV::MY_POOL][ip]; + const int gathered_start = global_z_start * this->ncxy; + for (int ixy = 0; ixy < this->ncxy; ++ixy) + { + for (int iz = 0; iz < local_nz; ++iz) + { + global_layout[ixy * this->ncz + global_z_start + iz] = pool_data[gathered_start + ixy * local_nz + iz]; + } + } + } + + MPI_Allreduce(MPI_IN_PLACE, global_layout.data(), this->ncxyz, MPI_DOUBLE, MPI_SUM, INT_BGROUP); + + // Return only the z-slab owned by this rank under its pool's distribution. + const int local_z_start = this->startz[GlobalV::MY_POOL][GlobalV::RANK_IN_POOL]; + for (int ixy = 0; ixy < this->ncxy; ++ixy) + { + for (int iz = 0; iz < this->nczp; ++iz) + { + data[ixy * this->nczp + iz] = global_layout[ixy * this->ncz + local_z_start + iz]; + } + } +#else + (void)data; +#endif +} + #ifdef __MPI void Parallel_Grid::bcast(const double* const data_global, double* data_local, const int& rank, const bool is_sdft) const { diff --git a/source/source_base/parallel_grid.h b/source/source_base/parallel_grid.h index 8e724957ce0..2ca486bf896 100644 --- a/source/source_base/parallel_grid.h +++ b/source/source_base/parallel_grid.h @@ -21,6 +21,16 @@ class Parallel_Grid const int &nczp, const int &nrxx, const int &nbz, const int &bz, const int nprocgroup); + /** + * @brief Sum a distributed real-space grid across k-point pools. + * + * Uses a direct local-slab reduction for equal-sized pools. For uneven + * pools, reconstructs a common global layout before the cross-pool sum. + * + * @param data Local real-space grid data ordered as [x][y][z]. + */ + void reduce_across_pools(double* data) const; + #ifdef __MPI /// @brief Broadcast data from root to all processors. The index order is [x][y][z]. void bcast(const double* const data_global, double* data_local, const int& rank, const bool is_sdft) const; diff --git a/source/source_estate/module_charge/charge.cpp b/source/source_estate/module_charge/charge.cpp index 235b0839e00..5f25aefe06f 100644 --- a/source/source_estate/module_charge/charge.cpp +++ b/source/source_estate/module_charge/charge.cpp @@ -42,10 +42,6 @@ Charge::Charge() Charge::~Charge() { this->destroy(); -#ifdef __MPI - delete[] rec; - delete[] dis; -#endif } void Charge::set_rhopw(ModulePW::PW_Basis* rhopw_in) diff --git a/source/source_estate/module_charge/charge.h b/source/source_estate/module_charge/charge.h index 761aecb1943..785985da40e 100644 --- a/source/source_estate/module_charge/charge.h +++ b/source/source_estate/module_charge/charge.h @@ -124,11 +124,6 @@ class Charge void init_final_scf(); //LiuXh add 20180619 public: - /** - * @brief init some arrays for mpi_inter_pools, rho_mpi - */ - void init_chgmpi(); - /** * @brief Sum rho at different pools (k-point parallelism). * Only used when GlobalV::KPAR > 1 @@ -169,12 +164,6 @@ class Charge bool allocate_rho; bool allocate_rho_final_scf; // LiuXh add 20180606 - -#ifdef __MPI - int *rec = nullptr; //The number of elements each process should receive into the receive buffer. - int *dis = nullptr; //The displacement (relative to recvbuf) for each process in the receive buffer. -#endif - }; #endif // charge diff --git a/source/source_estate/module_charge/charge_init.cpp b/source/source_estate/module_charge/charge_init.cpp index b4bd4f83717..30ab90598f5 100644 --- a/source/source_estate/module_charge/charge_init.cpp +++ b/source/source_estate/module_charge/charge_init.cpp @@ -249,9 +249,6 @@ void Charge::init_rho(const UnitCell& ucell, GlobalC::restart.info_load.load_charge_finish = true; } -#ifdef __MPI - this->init_chgmpi(); -#endif if (PARAM.inp.init_chg == "wfc") { if (wfcpw == nullptr) diff --git a/source/source_estate/module_charge/charge_mpi.cpp b/source/source_estate/module_charge/charge_mpi.cpp index c178471320c..8ef4cce4a32 100644 --- a/source/source_estate/module_charge/charge_mpi.cpp +++ b/source/source_estate/module_charge/charge_mpi.cpp @@ -2,115 +2,20 @@ #include "source_base/global_function.h" #include "source_base/global_variable.h" #include "source_base/parallel_comm.h" -#include "source_base/parallel_reduce.h" #include "source_base/timer.h" #include "source_hamilt/module_xc/xc_functional.h" #include "source_io/module_parameter/parameter.h" #ifdef __MPI -void Charge::init_chgmpi() -{ - if (KP_WORLD == MPI_COMM_NULL) - { - delete[] rec; - rec = new int[GlobalV::NPROC_IN_POOL]; - delete[] dis; - dis = new int[GlobalV::NPROC_IN_POOL]; - - const int ncxy = this->rhopw->nx * this->rhopw->ny; - for (int ip = 0; ip < GlobalV::NPROC_IN_POOL; ip++) - { - rec[ip] = this->rhopw->numz[ip] * ncxy; - dis[ip] = this->rhopw->startz[ip] * ncxy; - } - } -} - void Charge::reduce_diff_pools(double* array_rho) const { ModuleBase::TITLE("Charge", "reduce_diff_pools"); ModuleBase::timer::start("Charge", "reduce_diff_pools"); - if (KP_WORLD != MPI_COMM_NULL) + if (GlobalV::KPAR > 1) { - MPI_Allreduce(MPI_IN_PLACE, array_rho, this->nrxx, MPI_DOUBLE, MPI_SUM, KP_WORLD); + assert(this->pgrid != nullptr); + this->pgrid->reduce_across_pools(array_rho); } - else - { - double* array_tmp = new double[this->rhopw->nxyz]; - double* array_tot = new double[this->rhopw->nxyz]; - double* array_tot_aux = new double[this->rhopw->nxyz]; - //================================== - // Collect the rho in each pool - //================================== - for (int ir = 0; ir < this->rhopw->nrxx; ++ir) - { - array_tmp[ir] = array_rho[ir] / GlobalV::NPROC_IN_POOL; - } - MPI_Allgatherv(array_tmp, this->rhopw->nrxx, MPI_DOUBLE, array_tot, rec, dis, MPI_DOUBLE, POOL_WORLD); - - const int ncxy = this->rhopw->nx * this->rhopw->ny; - for (int ip = 0; ip < GlobalV::NPROC_IN_POOL; ++ip) - { - for (int ir = 0; ir < ncxy; ++ir) - { - for (int iz = 0; iz < this->rhopw->numz[ip]; ++iz) - { - // ------------------------------------------------- - // very carefully with the order of charge density. - // the data (ir,iz) is now in processor 'ip'. - // different POOL has different ordering. - // we want to collect them in each processor - // in a unit format, - // and then reduce among all POOLS to yield - // the correct charge density. - // we know the division of 'z' is indipendent - // in each processor, so the 'unit format' - // must have no relationship with 'z' divide method. - // ------------------------------------------------- - // rot_tot_aux : suitable among all pools. - // (1) the data save along z direction. - // (2) and each element number of group 'z data' - // is 'this->rhopw->nz' - // (3) however, the data rearrange is occured - // between [ this->rhopw->startz[ip], this->rhopw->startz[ip]+this->rhopw->numz[ip] ) - // (4) this->rhopw->startz[ip] + iz yields correct z coordiante. - // ------------------------------------------------- - // rot_tot: suitable for local pool. - // (1) the data save along z direction, only - // in a small distance. - // (2) however, the number of z in each processor - // 'ip' is this->rhopw->numz[ip] - // (3) the index of data increases with the ip, - // so the data on large 'ip' processor must - // have large 'start position', which we label - // this->rhopw->startz[ip] * ncxy. - // ------------------------------------------------- - array_tot_aux[this->rhopw->nz * ir + this->rhopw->startz[ip] + iz] - = array_tot[this->rhopw->numz[ip] * ir + this->rhopw->startz[ip] * ncxy + iz]; - } - } - } - - //================================== - // Reduce all the rho in each cpu - //================================== - MPI_Allreduce(array_tot_aux, array_tot, this->rhopw->nxyz, MPI_DOUBLE, MPI_SUM, INT_BGROUP); - - //===================================== - // Change the order of rho in each cpu - //===================================== - for (int ir = 0; ir < ncxy; ir++) - { - for (int iz = 0; iz < this->rhopw->numz[GlobalV::RANK_IN_POOL]; iz++) - { - array_rho[this->rhopw->numz[GlobalV::RANK_IN_POOL] * ir + iz] - = array_tot[this->rhopw->nz * ir + this->rhopw->startz_current + iz]; - } - } - delete[] array_tot_aux; - delete[] array_tot; - delete[] array_tmp; - } - if(PARAM.globalv.all_ks_run && PARAM.inp.bndpar > 1) + if (PARAM.globalv.all_ks_run && PARAM.inp.bndpar > 1) { MPI_Allreduce(MPI_IN_PLACE, array_rho, this->nrxx, MPI_DOUBLE, MPI_SUM, BP_WORLD); } @@ -120,10 +25,10 @@ void Charge::reduce_diff_pools(double* array_rho) const void Charge::rho_mpi() { ModuleBase::TITLE("Charge", "rho_mpi"); - if (GlobalV::KPAR * PARAM.inp.bndpar <= 1) - { - return; - } + if (GlobalV::KPAR * PARAM.inp.bndpar <= 1) + { + return; + } ModuleBase::timer::start("Charge", "rho_mpi"); for (int is = 0; is < PARAM.inp.nspin; ++is) diff --git a/source/source_estate/module_charge/symm_rho.cpp b/source/source_estate/module_charge/symm_rho.cpp index 6d4009437e1..a3417e96d9b 100644 --- a/source/source_estate/module_charge/symm_rho.cpp +++ b/source/source_estate/module_charge/symm_rho.cpp @@ -1,5 +1,6 @@ #include "symm_rho.h" +#include "source_estate/module_charge/charge.h" #include "source_hamilt/module_xc/xc_functional.h" #include "source_io/module_parameter/parameter.h" @@ -106,7 +107,7 @@ void Symmetry_rho::begin(const int& spin_now, psymmg(rhog[spin_now], rho_basis, symm); rho_basis->recip2real(rhog[spin_now], rho[spin_now]); - if (XC_Functional::get_ked_flag()) + if (XC_Functional::get_ked_flag() && kin_r != nullptr) { // Use std::vector to manage kin_g instead of raw pointer std::vector> kin_g(ngmc); @@ -123,6 +124,14 @@ void Symmetry_rho::begin(const int& spin_now, void Symmetry_rho::begin_soc(const Charge& chr, const ModulePW::PW_Basis* rho_basis, ModuleSymmetry::Symmetry& symm) const +{ + this->begin_soc(chr.rho, chr.rhog, rho_basis, symm); +} + +void Symmetry_rho::begin_soc(double** rho, + std::complex** rhog, + const ModulePW::PW_Basis* rho_basis, + ModuleSymmetry::Symmetry& symm) const { if (ModuleSymmetry::Symmetry::symm_flag != 1) { @@ -136,14 +145,14 @@ void Symmetry_rho::begin_soc(const Charge& chr, // reciprocal space and symmetrized together (rho[1]=rho^x, rho[2]=rho^y, rho[3]=rho^z). for (int is = 1; is < 4; ++is) { - rho_basis->real2recip(chr.rho[is], chr.rhog[is]); + rho_basis->real2recip(rho[is], rhog[is]); } - psymmg_soc(chr.rhog[1], chr.rhog[2], chr.rhog[3], rho_basis, symm); + psymmg_soc(rhog[1], rhog[2], rhog[3], rho_basis, symm); for (int is = 1; is < 4; ++is) { - rho_basis->recip2real(chr.rhog[is], chr.rho[is]); + rho_basis->recip2real(rhog[is], rho[is]); } ModuleBase::timer::end("Symmetry_rho", "begin_soc"); diff --git a/source/source_estate/module_charge/symm_rho.h b/source/source_estate/module_charge/symm_rho.h index 09b6ab5b351..f9d1cb48a89 100644 --- a/source/source_estate/module_charge/symm_rho.h +++ b/source/source_estate/module_charge/symm_rho.h @@ -2,9 +2,10 @@ #define SYMM_RHO_H #include "source_basis/module_pw/pw_basis.h" #include "source_cell/module_symmetry/symmetry.h" -#include "source_estate/module_charge/charge.h" #include "source_base/parallel_grid.h" +class Charge; + class Symmetry_rho { public: @@ -47,6 +48,19 @@ class Symmetry_rho const ModulePW::PW_Basis* pw, ModuleSymmetry::Symmetry& symm) const; + /** + * @brief Symmetrize raw nspin=4 spin-density arrays with coupled spin rotations. + * + * @param rho Real-space density components ordered as rho0, mx, my, mz. + * @param rhog Reciprocal-space work arrays with the same component ordering. + * @param pw Plane-wave basis used for the Fourier transforms. + * @param symm Symmetry operations and spin rotations. + */ + void begin_soc(double** rho, + std::complex** rhog, + const ModulePW::PW_Basis* pw, + ModuleSymmetry::Symmetry& symm) const; + private: // in real space: void psymm(double* rho_part, diff --git a/source/source_estate/test_mpi/charge_mpi_test.cpp b/source/source_estate/test_mpi/charge_mpi_test.cpp index e3214fbf655..3f37329e8d6 100644 --- a/source/source_estate/test_mpi/charge_mpi_test.cpp +++ b/source/source_estate/test_mpi/charge_mpi_test.cpp @@ -1,27 +1,31 @@ -#include "gmock/gmock.h" -#include "gtest/gtest.h" -#define private public #include "source_base/matrix3.h" #include "source_base/parallel_global.h" +#include "source_base/parallel_grid.h" #include "source_estate/module_charge/charge.h" #include "source_hamilt/module_xc/xc_functional.h" #include "source_io/module_parameter/parameter.h" +#include "gmock/gmock.h" +#include "gtest/gtest.h" + bool XC_Functional::ked_flag = false; +void XC_Functional::set_xc_type(const std::string xc_func_in) +{ + ked_flag = xc_func_in == "scan"; +} Charge::Charge() { } Charge::~Charge() { - delete[] rec; - delete[] dis; } auto sum_array = [](const double* v, const int& nv) { double sum = 0; - for (int i = 0; i < nv; ++i) { + for (int i = 0; i < nv; ++i) + { sum += v[i]; -} + } return sum; }; /************************************************ @@ -44,9 +48,11 @@ class ChargeMpiTest : public ::testing::Test std::string output; double lat0 = 4; ModuleBase::Matrix3 latvec; + int world_size = 1; void SetUp() override { charge = new Charge; + MPI_Comm_size(MPI_COMM_WORLD, &world_size); } void TearDown() override { @@ -60,20 +66,23 @@ TEST_F(ChargeMpiTest, reduce_diff_pools1) { GlobalV::KPAR = 2; Parallel_Global::init_pools(GlobalV::NPROC, - GlobalV::MY_RANK, - PARAM.input.bndpar, - GlobalV::KPAR, - GlobalV::NPROC_IN_BNDGROUP, - GlobalV::RANK_IN_BPGROUP, - GlobalV::MY_BNDGROUP, - GlobalV::NPROC_IN_POOL, - GlobalV::RANK_IN_POOL, - GlobalV::MY_POOL); + GlobalV::MY_RANK, + PARAM.inp.bndpar, + GlobalV::KPAR, + GlobalV::NPROC_IN_BNDGROUP, + GlobalV::RANK_IN_BPGROUP, + GlobalV::MY_BNDGROUP, + GlobalV::NPROC_IN_POOL, + GlobalV::RANK_IN_POOL, + GlobalV::MY_POOL); ModulePW::PW_Basis* rhopw = new ModulePW::PW_Basis(); rhopw->initmpi(GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL, POOL_WORLD); rhopw->initgrids(lat0, latvec, 40); rhopw->initparameters(false, 10); rhopw->setuptransform(); + Parallel_Grid pgrid; + pgrid.init(rhopw->nx, rhopw->ny, rhopw->nz, rhopw->nplane, rhopw->nrxx, rhopw->nz, 1, world_size); + charge->pgrid = &pgrid; int nz = rhopw->nz; const int nrxx = rhopw->nrxx; @@ -90,7 +99,6 @@ TEST_F(ChargeMpiTest, reduce_diff_pools1) } double refsum = sum_array(array_rho, nrxx); - charge->init_chgmpi(); charge->reduce_diff_pools(array_rho); double sum = sum_array(array_rho, nrxx); EXPECT_EQ(sum, refsum * GlobalV::KPAR); @@ -107,7 +115,7 @@ TEST_F(ChargeMpiTest, reduce_diff_pools2) GlobalV::KPAR = 3; Parallel_Global::divide_pools(GlobalV::NPROC, GlobalV::MY_RANK, - PARAM.input.bndpar, + PARAM.inp.bndpar, GlobalV::KPAR, GlobalV::NPROC_IN_BNDGROUP, GlobalV::RANK_IN_BPGROUP, @@ -121,6 +129,9 @@ TEST_F(ChargeMpiTest, reduce_diff_pools2) rhopw->initparameters(false, 10); rhopw->setuptransform(); charge->rhopw = rhopw; + Parallel_Grid pgrid; + pgrid.init(rhopw->nx, rhopw->ny, rhopw->nz, rhopw->nplane, rhopw->nrxx, rhopw->nz, 1, world_size); + charge->pgrid = &pgrid; int nz = rhopw->nz; const int nrxx = rhopw->nrxx; @@ -143,7 +154,6 @@ TEST_F(ChargeMpiTest, reduce_diff_pools2) } } - charge->init_chgmpi(); charge->reduce_diff_pools(array_rho); double sum = sum_array(array_rho, nrxx); MPI_Allreduce(MPI_IN_PLACE, &sum, 1, MPI_DOUBLE, MPI_SUM, POOL_WORLD); @@ -162,7 +172,7 @@ TEST_F(ChargeMpiTest, rho_mpi) GlobalV::KPAR = 2; Parallel_Global::divide_pools(GlobalV::NPROC, GlobalV::MY_RANK, - PARAM.input.bndpar, + PARAM.inp.bndpar, GlobalV::KPAR, GlobalV::NPROC_IN_BNDGROUP, GlobalV::RANK_IN_BPGROUP, @@ -176,7 +186,10 @@ TEST_F(ChargeMpiTest, rho_mpi) rhopw->initparameters(false, 10); rhopw->setuptransform(); charge->rhopw = rhopw; - PARAM.input.nspin = 1; + Parallel_Grid pgrid; + pgrid.init(rhopw->nx, rhopw->ny, rhopw->nz, rhopw->nplane, rhopw->nrxx, rhopw->nz, 1, world_size); + charge->pgrid = &pgrid; + ASSERT_EQ(PARAM.inp.nspin, 1); charge->rho = new double*[1]; charge->kin_r = new double*[1]; @@ -187,7 +200,6 @@ TEST_F(ChargeMpiTest, rho_mpi) charge->nrxx = nrxx; charge->rho[0] = new double[nrxx]; charge->kin_r[0] = new double[nrxx]; - charge->init_chgmpi(); charge->rho_mpi(); delete[] charge->rho[0]; @@ -205,15 +217,14 @@ TEST_F(ChargeMpiTest, kin_r_mpi) { if (GlobalV::NPROC >= 2 && GlobalV::NPROC % 2 == 0) { - const bool ked_flag_old = XC_Functional::ked_flag; - XC_Functional::ked_flag = true; - PARAM.input.nspin = 1; - PARAM.input.bndpar = 1; + XC_Functional::set_xc_type("scan"); + ASSERT_EQ(PARAM.inp.nspin, 1); + ASSERT_EQ(PARAM.inp.bndpar, 1); GlobalV::KPAR = 2; Parallel_Global::divide_pools(GlobalV::NPROC, GlobalV::MY_RANK, - PARAM.input.bndpar, + PARAM.inp.bndpar, GlobalV::KPAR, GlobalV::NPROC_IN_BNDGROUP, GlobalV::RANK_IN_BPGROUP, @@ -227,6 +238,9 @@ TEST_F(ChargeMpiTest, kin_r_mpi) rhopw->initparameters(false, 10); rhopw->setuptransform(); charge->rhopw = rhopw; + Parallel_Grid pgrid; + pgrid.init(rhopw->nx, rhopw->ny, rhopw->nz, rhopw->nplane, rhopw->nrxx, rhopw->nz, 1, world_size); + charge->pgrid = &pgrid; const int nz = rhopw->nz; const int nrxx = rhopw->nrxx; @@ -240,13 +254,11 @@ TEST_F(ChargeMpiTest, kin_r_mpi) { for (int iz = 0; iz < nplane; ++iz) { - charge->kin_r[0][nplane * ir + iz] - = (rhopw->startz_current + iz + ir * nz) / double(nxy * nz); + charge->kin_r[0][nplane * ir + iz] = (rhopw->startz_current + iz + ir * nz) / double(nxy * nz); } } const double refsum = sum_array(charge->kin_r[0], nrxx); - charge->init_chgmpi(); charge->kin_r_mpi(); const double sum = sum_array(charge->kin_r[0], nrxx); EXPECT_EQ(sum, refsum * GlobalV::KPAR); @@ -254,7 +266,7 @@ TEST_F(ChargeMpiTest, kin_r_mpi) delete[] charge->kin_r[0]; delete[] charge->kin_r; delete rhopw; - XC_Functional::ked_flag = ked_flag_old; + XC_Functional::set_xc_type("pbe"); } } diff --git a/source/source_io/module_chgpot/get_pchg_pw.h b/source/source_io/module_chgpot/get_pchg_pw.h index e572d26ef91..cf38872cb60 100644 --- a/source/source_io/module_chgpot/get_pchg_pw.h +++ b/source/source_io/module_chgpot/get_pchg_pw.h @@ -1,247 +1,298 @@ #ifndef GET_PCHG_PW_H #define GET_PCHG_PW_H -#include "source_io/module_output/cube_io.h" +#include "source_base/module_container/ATen/core/tensor.h" #include "source_estate/module_charge/symm_rho.h" +#include "source_io/module_output/cube_io.h" namespace ModuleIO { +/** + * @brief Write band-resolved partial charges from plane-wave coefficients. + * + * The selected bands are transformed to the dense real-space grid. Depending on + * @p if_separate_k, the function either writes one cube per k-point or sums the + * k-point contributions, restores symmetry, and writes one cube per spin/charge + * component. For spinors, the four components are charge, m_x, m_y, and m_z. + */ template void get_pchg_pw(const std::vector& out_pchg, - const int nbands, const int nspin, - const int nxyz, - const int chr_ngmc, UnitCell* ucell, const psi::Psi, Device>* kspw_psi, + const ModulePW::PW_Basis* pw_rho, const ModulePW::PW_Basis* pw_rhod, const ModulePW::PW_Basis_K* pw_wfc, - const Device* ctx, const Parallel_Grid& pgrid, const std::string& global_out_dir, const bool if_separate_k, - const K_Vectors& kv, - const int kpar, - const int my_pool, - const Charge* chr) // Charge class is needed for the charge density reduce + const bool noncolin, + const K_Vectors& kv) { - // Get necessary parameters from kv const int nks = kv.get_nks(); // current process pool k-point count const int nkstot = kv.get_nkstot(); // total k-point count + const int nbands = kspw_psi->get_nbands(); + + const int nks_without_spin = nspin == 2 ? nkstot / 2 : nkstot; + const int smooth_nrxx = pw_wfc->nrxx; + const int dense_nrxx = pw_rhod->nrxx; + // Avoid an extra forward and backward FFT unless a distinct dense grid is in use. + const bool needs_interpolation = pw_rhod != pw_rho; - // Loop over k-parallelism - for (int ip = 0; ip < kpar; ++ip) + // Expand the INPUT selection into a fixed-size mask indexed directly by band. + std::vector bands_picked(nbands, 0); + if (static_cast(out_pchg.size()) > nbands) { - if (my_pool != ip) + ModuleBase::WARNING_QUIT("ModuleIO::get_pchg_pw", + "The number of bands specified by `out_pchg` in the " + "INPUT file exceeds `nbands`!"); + } + for (int value: out_pchg) + { + if (value != 0 && value != 1) { - continue; + ModuleBase::WARNING_QUIT("ModuleIO::get_pchg_pw", + "The elements of `out_pchg` must be either 0 or 1. " + "Invalid values found!"); } - - // bands_picked is a vector of 0s and 1s, where 1 means the band is picked to output - std::vector bands_picked(nbands, 0); - - // Check if length of out_pchg is valid - if (static_cast(out_pchg.size()) > nbands) + } + const int length = std::min(static_cast(out_pchg.size()), nbands); + for (int i = 0; i < length; ++i) + { + bands_picked[i] = static_cast(out_pchg[i]); + } + + // Map the wavefunction backend type to the tensor library's device type. + using ContainerDevice = typename ct::PsiToContainer::type; + const ct::DeviceType device_type = ct::DeviceTypeToEnum::value; + const bool is_cpu = device_type == ct::DeviceType::CpuDevice; + const bool is_spinor = nspin == 4; + // Spinor coefficients store the up and down blocks consecutively. + const int npwx = kspw_psi->get_nbasis() / (is_spinor ? 2 : 1); + + // Zero-length tensors avoid allocating buffers that a given backend or spin mode never uses. + ct::Tensor wfcr_up_smooth(ct::DataType::DT_COMPLEX_DOUBLE, device_type, ct::TensorShape({smooth_nrxx})); + ct::Tensor wfcr_down_smooth(ct::DataType::DT_COMPLEX_DOUBLE, device_type, ct::TensorShape({is_spinor ? smooth_nrxx : 0})); + ct::Tensor wfcr_up_smooth_host(ct::DataType::DT_COMPLEX_DOUBLE, ct::DeviceType::CpuDevice, ct::TensorShape({is_cpu ? 0 : smooth_nrxx})); + ct::Tensor wfcr_down_smooth_host(ct::DataType::DT_COMPLEX_DOUBLE, + ct::DeviceType::CpuDevice, + ct::TensorShape({!is_cpu && is_spinor ? smooth_nrxx : 0})); + ct::Tensor wfcr_up_dense_host(ct::DataType::DT_COMPLEX_DOUBLE, + ct::DeviceType::CpuDevice, + ct::TensorShape({needs_interpolation ? dense_nrxx : 0})); + ct::Tensor wfcr_down_dense_host(ct::DataType::DT_COMPLEX_DOUBLE, + ct::DeviceType::CpuDevice, + ct::TensorShape({needs_interpolation && is_spinor ? dense_nrxx : 0})); + ct::Tensor reciprocal_buffer_host(ct::DataType::DT_COMPLEX_DOUBLE, + ct::DeviceType::CpuDevice, + ct::TensorShape({needs_interpolation ? pw_rhod->npw : 0})); + + // Capture the shared bases and scratch buffers by reference. The returned pointer is owned by + // one of the tensor arguments and remains valid until that tensor is reused or destroyed. + auto transform_wfc = [&](const std::complex* coefficients, + const int ik, + ct::Tensor& smooth, + ct::Tensor& smooth_host, + ct::Tensor& dense_host) -> const std::complex* { + // Perform the wavefunction FFT on its native device, then expose host data for cube output. + pw_wfc->template recip_to_real, Device>(coefficients, smooth.data>(), ik); + const std::complex* smooth_data = smooth.data>(); + if (!is_cpu) { - ModuleBase::WARNING_QUIT("ModuleIO::get_pchg_pw", - "The number of bands specified by `out_pchg` in the " - "INPUT file exceeds `nbands`!"); + ct::kernels::synchronize_memory, ct::DEVICE_CPU, ContainerDevice>()( + smooth_host.data>(), + smooth_data, + smooth_nrxx); + smooth_data = smooth_host.data>(); } - - // Check if all elements in bands_picked are 0 or 1 - for (int value: out_pchg) + if (!needs_interpolation) { - if (value != 0 && value != 1) - { - ModuleBase::WARNING_QUIT("ModuleIO::get_pchg_pw", - "The elements of `out_pchg` must be either 0 or 1. " - "Invalid values found!"); - } + return smooth_data; } - // Fill bands_picked with values from out_pchg - // Remaining bands are already set to 0 - int length = std::min(static_cast(out_pchg.size()), nbands); - for (int i = 0; i < length; ++i) + // Zero padding in reciprocal space transfers the smooth-grid field to the dense rho grid. + reciprocal_buffer_host.zero(); + pw_rho->real2recip(smooth_data, reciprocal_buffer_host.data>()); + pw_rhod->recip2real(reciprocal_buffer_host.data>(), dense_host.data>()); + return dense_host.data>(); + }; + + std::vector> rho_band(nspin, std::vector(dense_nrxx)); + // Convert a two-component spinor into the Pauli-basis fields (rho, m_x, m_y, m_z). + // Per-k output overwrites the fields, whereas k-summed output accumulates weighted fields. + auto accumulate_spinor_density + = [&](const std::complex* up, const std::complex* down, const double weight, const bool accumulate) { + for (int ir = 0; ir < dense_nrxx; ++ir) + { + const double up_norm = std::norm(up[ir]); + const double down_norm = std::norm(down[ir]); + const double rho0 = (up_norm + down_norm) * weight; + const double mx = 2.0 * (up[ir].real() * down[ir].real() + up[ir].imag() * down[ir].imag()) * weight; + const double my = 2.0 * (up[ir].real() * down[ir].imag() - down[ir].real() * up[ir].imag()) * weight; + const double mz = (up_norm - down_norm) * weight; + if (accumulate) + { + rho_band[0][ir] += rho0; + rho_band[1][ir] += noncolin ? mx : 0.0; + rho_band[2][ir] += noncolin ? my : 0.0; + rho_band[3][ir] += mz; + } + else + { + rho_band[0][ir] = rho0; + rho_band[1][ir] = noncolin ? mx : 0.0; + rho_band[2][ir] = noncolin ? my : 0.0; + rho_band[3][ir] = mz; + } + } + }; + + for (int ib = 0; ib < nbands; ++ib) + { + if (!bands_picked[ib]) { - // out_pchg rely on function parse_expression - bands_picked[i] = static_cast(out_pchg[i]); + continue; } - // Allocate host memory - std::vector> wfcr(nxyz); - std::vector> rho_band(nspin, std::vector(nxyz)); - - // Allocate device memory - std::complex* wfcr_device = nullptr; - if (!std::is_same::value) + for (int is = 0; is < nspin; ++is) { - base_device::memory::resize_memory_op, Device>()(wfcr_device, nxyz); + std::fill(rho_band[is].begin(), rho_band[is].end(), 0.0); } - for (int ib = 0; ib < nbands; ++ib) + if (if_separate_k) { - // Skip the loop iteration if bands_picked[ib] is 0 - if (!bands_picked[ib]) - { - continue; - } - - for (int is = 0; is < nspin; ++is) - { - std::fill(rho_band[is].begin(), rho_band[is].end(), 0.0); - } - - if (if_separate_k) + // Preserve each Bloch state's contribution; no Brillouin-zone weight is applied here. + for (int ik = 0; ik < nks; ++ik) { - for (int ik = 0; ik < nks; ++ik) + const int ikstot = kv.ik2iktot[ik]; + const int spin_index = kv.isk[ik]; + // In collinear calculations the two spin channels share the same k-point numbering. + const int k_number = ikstot % nks_without_spin + 1; + + kspw_psi->fix_k(ik); + const std::complex* wfcr_up + = transform_wfc(&kspw_psi[0](ib, 0), ik, wfcr_up_smooth, wfcr_up_smooth_host, wfcr_up_dense_host); + const std::complex* wfcr_up_host_data = wfcr_up; + const std::complex* wfcr_down_host_data = nullptr; + if (is_spinor) { - const int ikstot = kv.ik2iktot[ik]; // global k-point index - const int spin_index = kv.isk[ik]; // spin index - const int k_number = ikstot % (nkstot / nspin) + 1; // k-point number, starting from 1 - - kspw_psi->fix_k(ik); - - // FFT on device and copy result back to host - if (std::is_same::value) - { - pw_wfc->recip_to_real(ctx, &kspw_psi[0](ib, 0), wfcr.data(), ik); - } - else - { - pw_wfc->recip_to_real(ctx, &kspw_psi[0](ib, 0), wfcr_device, ik); - - base_device::memory::synchronize_memory_op, - base_device::DEVICE_CPU, - Device>()(wfcr.data(), wfcr_device, nxyz); - } - - // To ensure the normalization of charge density in multi-k calculation (if if_separate_k is true) - double wg_sum_k = 0.0; - if (nspin == 1) - { - wg_sum_k = 2.0; - } - else if (nspin == 2) - { - wg_sum_k = 1.0; - } - else - { - ModuleBase::WARNING_QUIT("ModuleIO::get_pchg_pw", - "Real space partial charge output currently do not support " - "noncollinear polarized calculation (nspin = 4)!"); - } - - double w1 = static_cast(wg_sum_k / ucell->omega); + const std::complex* wfcr_down + = transform_wfc(&kspw_psi[0](ib, npwx), ik, wfcr_down_smooth, wfcr_down_smooth_host, wfcr_down_dense_host); + wfcr_down_host_data = wfcr_down; + } - for (int i = 0; i < nxyz; ++i) + const double spin_degeneracy = nspin == 1 ? 2.0 : 1.0; + const double weight = spin_degeneracy / ucell->omega; + if (is_spinor) + { + accumulate_spinor_density(wfcr_up_host_data, wfcr_down_host_data, weight, false); + } + else + { + for (int ir = 0; ir < dense_nrxx; ++ir) { - rho_band[spin_index][i] = std::norm(wfcr[i]) * w1; + rho_band[spin_index][ir] = std::norm(wfcr_up_host_data[ir]) * weight; } + } + const int component_begin = is_spinor ? 0 : spin_index; + const int component_end = is_spinor ? 4 : spin_index + 1; + for (int component = component_begin; component < component_end; ++component) + { std::stringstream ssc; - ssc << global_out_dir << "pchgi" << ib + 1 << "s" << spin_index + 1 << "k" << k_number << ".cube"; - + ssc << global_out_dir << "pchgi" << ib + 1 << "s" << component + 1 << "k" << k_number << ".cube"; ModuleIO::write_vdata_palgrid(pgrid, - rho_band[spin_index].data(), - spin_index, + rho_band[component].data(), + component, nspin, 0, ssc.str(), 0.0, ucell, 11, - 1, - PARAM.globalv.two_fermi, - true); // reduce_all_pool is true + 0, + false, + true); } } - else + } + else + { + // Form a Brillouin-zone weighted partial density for the selected band. + for (int ik = 0; ik < nks; ++ik) { - for (int ik = 0; ik < nks; ++ik) + const int spin_index = kv.isk[ik]; + + kspw_psi->fix_k(ik); + const std::complex* wfcr_up + = transform_wfc(&kspw_psi[0](ib, 0), ik, wfcr_up_smooth, wfcr_up_smooth_host, wfcr_up_dense_host); + const std::complex* wfcr_up_host_data = wfcr_up; + const std::complex* wfcr_down_host_data = nullptr; + if (is_spinor) { - const int ikstot = kv.ik2iktot[ik]; // global k-point index - const int spin_index = kv.isk[ik]; // spin index - const int k_number = ikstot % (nkstot / nspin) + 1; // k-point number, starting from 1 - - kspw_psi->fix_k(ik); - - // FFT on device and copy result back to host - if (std::is_same::value) - { - pw_wfc->recip_to_real(ctx, &kspw_psi[0](ib, 0), wfcr.data(), ik); - } - else - { - pw_wfc->recip_to_real(ctx, &kspw_psi[0](ib, 0), wfcr_device, ik); - - base_device::memory::synchronize_memory_op, - base_device::DEVICE_CPU, - Device>()(wfcr.data(), wfcr_device, nxyz); - } - - double w1 = static_cast(kv.wk[ik] / ucell->omega); - - for (int i = 0; i < nxyz; ++i) - { - rho_band[spin_index][i] += std::norm(wfcr[i]) * w1; - } + const std::complex* wfcr_down + = transform_wfc(&kspw_psi[0](ib, npwx), ik, wfcr_down_smooth, wfcr_down_smooth_host, wfcr_down_dense_host); + wfcr_down_host_data = wfcr_down; } -#ifdef __MPI - // Reduce the charge density across all pools if kpar > 1 - if (kpar > 1 && chr != nullptr) + const double weight = static_cast(kv.wk[ik] / ucell->omega); + if (is_spinor) { - for (int is = 0; is < nspin; ++is) + accumulate_spinor_density(wfcr_up_host_data, wfcr_down_host_data, weight, true); + } + else + { + for (int ir = 0; ir < dense_nrxx; ++ir) { - chr->reduce_diff_pools(rho_band[is].data()); + rho_band[spin_index][ir] += std::norm(wfcr_up_host_data[ir]) * weight; } } -#endif + } - // Symmetrize the charge density, otherwise the results are incorrect if the symmetry is on - // std::cout << " Symmetrizing band-decomposed charge density..." << std::endl; - Symmetry_rho srho; +#ifdef __MPI + if (kv.para_k.kpar > 1) + { + // Each pool owns only part of the k-point sum; assemble it before symmetrization. for (int is = 0; is < nspin; ++is) { - // Use vector instead of raw pointers - std::vector rho_save_pointers(nspin); - for (int s = 0; s < nspin; ++s) - { - rho_save_pointers[s] = rho_band[s].data(); - } - - std::vector>> rhog(nspin, - std::vector>(chr_ngmc)); - - // Convert vector of vectors to vector of pointers - std::vector*> rhog_pointers(nspin); - for (int s = 0; s < nspin; ++s) - { - rhog_pointers[s] = rhog[s].data(); - } - - srho.begin(is, - rho_save_pointers.data(), - rhog_pointers.data(), - chr_ngmc, - nullptr, - pw_rhod, - ucell->symm); + pgrid.reduce_across_pools(rho_band[is].data()); } + } +#endif + // Symmetry_rho operates on arrays of component pointers and uses reciprocal workspaces. + Symmetry_rho srho; + std::vector rho_save_pointers(nspin); + std::vector>> rhog(nspin, std::vector>(pw_rhod->npw)); + std::vector*> rhog_pointers(nspin); + for (int is = 0; is < nspin; ++is) + { + rho_save_pointers[is] = rho_band[is].data(); + rhog_pointers[is] = rhog[is].data(); + } + if (is_spinor) + { + // Charge and magnetization components obey different spinor symmetry transformations. + srho.begin(0, rho_save_pointers.data(), rhog_pointers.data(), pw_rhod->npw, nullptr, pw_rhod, ucell->symm); + srho.begin_soc(rho_save_pointers.data(), rhog_pointers.data(), pw_rhod, ucell->symm); + } + else + { for (int is = 0; is < nspin; ++is) { - std::stringstream ssc; - ssc << global_out_dir << "pchgi" << ib + 1 << "s" << is + 1 << ".cube"; - - ModuleIO::write_vdata_palgrid(pgrid, rho_band[is].data(), is, nspin, 0, ssc.str(), 0.0, ucell, 11, 1, PARAM.globalv.two_fermi, false); + srho.begin(is, rho_save_pointers.data(), rhog_pointers.data(), pw_rhod->npw, nullptr, pw_rhod, ucell->symm); } - } // else if_separate_k is false - } // end of ib loop over nbands - } // end of ip loop over kpar -} // get_pchg_pw + } + + for (int is = 0; is < nspin; ++is) + { + std::stringstream ssc; + ssc << global_out_dir << "pchgi" << ib + 1 << "s" << is + 1 << ".cube"; + ModuleIO::write_vdata_palgrid(pgrid, rho_band[is].data(), is, nspin, 0, ssc.str(), 0.0, ucell, 11, 0, false, false); + } + } + } +} } // namespace ModuleIO #endif // GET_PCHG_PW_H diff --git a/source/source_io/module_ctrl/ctrl_output_pw.cpp b/source/source_io/module_ctrl/ctrl_output_pw.cpp index bfd06fd76fd..9fd7783e902 100644 --- a/source/source_io/module_ctrl/ctrl_output_pw.cpp +++ b/source/source_io/module_ctrl/ctrl_output_pw.cpp @@ -1,30 +1,30 @@ #include "ctrl_output_pw.h" -#include "../module_wf/write_wfc_pw.h" // use write_wfc_pw -#include "../module_dos/write_dos_pw.h" // use write_dos_pw -#include "../module_wannier/to_wannier90_pw.h" // wannier90 interface -#include "source_pw/module_pwdft/onsite_proj.h" // use projector #include "../module_bessel/numerical_basis.h" #include "../module_bessel/numerical_descriptor.h" +#include "../module_chgpot/get_pchg_pw.h" #include "../module_dos/cal_ldos.h" +#include "../module_dos/write_dos_pw.h" // use write_dos_pw #include "../module_unk/berryphase.h" -#include "source_lcao/module_deltaspin/spin_constrain.h" -#include "source_base/formatter.h" -#include "../module_chgpot/get_pchg_pw.h" +#include "../module_wannier/to_wannier90_pw.h" // wannier90 interface #include "../module_wf/get_wf_pw.h" +#include "../module_wf/write_wfc_pw.h" // use write_wfc_pw +#include "source_base/formatter.h" +#include "source_lcao/module_deltaspin/spin_constrain.h" #include "source_pw/module_pwdft/elecond.h" +#include "source_pw/module_pwdft/onsite_proj.h" // use projector #ifdef __MLALGO #include "../module_ml/write_mlkedf_desc.h" #endif -void ModuleIO::ctrl_iter_pw(const int istep, - const int iter, - const double &conv_esolver, - psi::Psi, base_device::DEVICE_CPU>* psi, - const K_Vectors &kv, - const ModulePW::PW_Basis_K *pw_wfc, - const Input_para& inp) +void ModuleIO::ctrl_iter_pw(const int istep, + const int iter, + const double& conv_esolver, + psi::Psi, base_device::DEVICE_CPU>* psi, + const K_Vectors& kv, + const ModulePW::PW_Basis_K* pw_wfc, + const Input_para& inp) { ModuleBase::TITLE("ModuleIO", "ctrl_iter_pw"); ModuleBase::timer::start("ModuleIO", "ctrl_iter_pw"); @@ -38,11 +38,11 @@ void ModuleIO::ctrl_iter_pw(const int istep, int istep_in = -1; int iter_in = -1; bool out_wfc_flag = false; - if (inp.out_freq_ion>0) // default value of out_freq_ion is 0 + if (inp.out_freq_ion > 0) // default value of out_freq_ion is 0 { if (istep % inp.out_freq_ion == 0) { - if(iter % inp.out_freq_elec == 0 || iter == inp.scf_nmax || conv_esolver) + if (iter % inp.out_freq_elec == 0 || iter == inp.scf_nmax || conv_esolver) { istep_in = istep; iter_in = iter; @@ -50,56 +50,53 @@ void ModuleIO::ctrl_iter_pw(const int istep, } } } - else if(iter == inp.scf_nmax || conv_esolver) + else if (iter == inp.scf_nmax || conv_esolver) { out_wfc_flag = true; } if (out_wfc_flag) { - ModuleIO::write_wfc_pw(istep_in, iter_in, - GlobalV::KPAR, - GlobalV::MY_POOL, - GlobalV::MY_RANK, - inp.nbands, - inp.nspin, - PARAM.globalv.npol, - GlobalV::RANK_IN_POOL, - GlobalV::NPROC_IN_POOL, - inp.out_wfc_pw, - inp.ecutwfc, - PARAM.globalv.global_out_dir, - psi[0], - kv, - pw_wfc, - GlobalV::ofs_running); + ModuleIO::write_wfc_pw(istep_in, + iter_in, + GlobalV::KPAR, + GlobalV::MY_POOL, + GlobalV::MY_RANK, + inp.nbands, + inp.nspin, + PARAM.globalv.npol, + GlobalV::RANK_IN_POOL, + GlobalV::NPROC_IN_POOL, + inp.out_wfc_pw, + inp.ecutwfc, + PARAM.globalv.global_out_dir, + psi[0], + kv, + pw_wfc, + GlobalV::ofs_running); } - ModuleBase::timer::end("ModuleIO", "ctrl_iter_pw"); - return; + ModuleBase::timer::end("ModuleIO", "ctrl_iter_pw"); + return; } - template void ModuleIO::ctrl_scf_pw(const int istep, - UnitCell& ucell, - elecstate::ElecState* pelec, - const Charge &chr, - const K_Vectors &kv, - const ModulePW::PW_Basis_K *pw_wfc, - const ModulePW::PW_Basis *pw_rho, - const ModulePW::PW_Basis *pw_rhod, - const ModulePW::PW_Basis_Big *pw_big, - Setup_Psi_pw &stp, - const Parallel_Grid ¶_grid, - const Input_para& inp) + UnitCell& ucell, + elecstate::ElecState* pelec, + const Charge& chr, + const K_Vectors& kv, + const ModulePW::PW_Basis_K* pw_wfc, + const ModulePW::PW_Basis* pw_rho, + const ModulePW::PW_Basis* pw_rhod, + const ModulePW::PW_Basis_Big* pw_big, + Setup_Psi_pw& stp, + const Parallel_Grid& para_grid, + const Input_para& inp) { ModuleBase::TITLE("ModuleIO", "ctrl_scf_pw"); ModuleBase::timer::start("ModuleIO", "ctrl_scf_pw"); - // Create local ctx for device type deduction - Device* ctx = nullptr; - // Transfer data from device (GPU) to host (CPU) in pw basis stp.copy_d2h(); @@ -113,16 +110,16 @@ void ModuleIO::ctrl_scf_pw(const int istep, int istep_in = -1; // default value of out_freq_ion is 0 - if(inp.out_freq_ion==0) + if (inp.out_freq_ion == 0) { out_dos_tmp = true; } - else if (inp.out_freq_ion>0) + else if (inp.out_freq_ion > 0) { if (istep % inp.out_freq_ion == 0) { out_dos_tmp = true; - istep_in=istep; + istep_in = istep; } else { @@ -141,23 +138,22 @@ void ModuleIO::ctrl_scf_pw(const int istep, out_dos_tmp = false; } - if(out_dos_tmp) + if (out_dos_tmp) { ModuleIO::write_dos_pw(ucell, - pelec->ekb, - pelec->wg, - kv, - inp.nbands, - istep_in, - pelec->eferm, - inp.dos_edelta_ev, - inp.dos_scale, - inp.dos_sigma, - GlobalV::ofs_running); + pelec->ekb, + pelec->wg, + kv, + inp.nbands, + istep_in, + pelec->eferm, + inp.dos_edelta_ev, + inp.dos_scale, + inp.dos_sigma, + GlobalV::ofs_running); } } - //------------------------------------------------------------------ // 5) calculate band-decomposed (partial) charge density in pw basis //------------------------------------------------------------------ @@ -166,29 +162,20 @@ void ModuleIO::ctrl_scf_pw(const int istep, // update psi_d stp.update_psi_d(); - const int nbands = stp.get_nbands(); - const int ngmc = chr.ngmc; - ModuleIO::get_pchg_pw(inp.out_pchg, - nbands, inp.nspin, - pw_rhod->nxyz, - ngmc, &ucell, stp.template get_psi_d(), + pw_rho, pw_rhod, pw_wfc, - ctx, para_grid, PARAM.globalv.global_out_dir, inp.if_separate_k, - kv, - GlobalV::KPAR, - GlobalV::MY_POOL, - &chr); + inp.noncolin, + kv); } - //------------------------------------------------------------------ //! 6) calculate Wannier functions in pw basis //------------------------------------------------------------------ @@ -207,7 +194,6 @@ void ModuleIO::ctrl_scf_pw(const int istep, std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Wannier functions calculation"); } - //------------------------------------------------------------------ //! 7) calculate Berry phase polarization in pw basis //------------------------------------------------------------------ @@ -225,8 +211,7 @@ void ModuleIO::ctrl_scf_pw(const int istep, //------------------------------------------------------------------ if (inp.sc_mag_switch) { - spinconstrain::SpinConstrain>& sc - = spinconstrain::SpinConstrain>::getScInstance(); + spinconstrain::SpinConstrain>& sc = spinconstrain::SpinConstrain>::getScInstance(); sc.cal_mi_pw(); sc.print_Mag_Force(GlobalV::ofs_running); } @@ -246,34 +231,30 @@ void ModuleIO::ctrl_scf_pw(const int istep, } template -void ModuleIO::ctrl_runner_pw(UnitCell& ucell, - elecstate::ElecState* pelec, - ModulePW::PW_Basis_K* pw_wfc, - ModulePW::PW_Basis* pw_rho, - ModulePW::PW_Basis* pw_rhod, - Charge &chr, - K_Vectors &kv, - Setup_Psi_pw &stp, - Structure_Factor &sf, - pseudopot_cell_vnl &ppcell, - surchem &solvent, - Parallel_Grid ¶_grid, - const Input_para& inp) +void ModuleIO::ctrl_runner_pw(UnitCell& ucell, + elecstate::ElecState* pelec, + ModulePW::PW_Basis_K* pw_wfc, + ModulePW::PW_Basis* pw_rho, + ModulePW::PW_Basis* pw_rhod, + Charge& chr, + K_Vectors& kv, + Setup_Psi_pw& stp, + Structure_Factor& sf, + pseudopot_cell_vnl& ppcell, + surchem& solvent, + Parallel_Grid& para_grid, + const Input_para& inp) { ModuleBase::TITLE("ModuleIO", "ctrl_runner_pw"); ModuleBase::timer::start("ModuleIO", "ctrl_runner_pw"); - // Create local ctx for device type deduction - Device* ctx = nullptr; - - //---------------------------------------------------------- - //! 1) Compute LDOS - //---------------------------------------------------------- - if (inp.out_ldos[0]) - { - ModuleIO::cal_ldos_pw(reinterpret_cast>*>(pelec), - stp.psi_cpu[0], para_grid, ucell); - } + //---------------------------------------------------------- + //! 1) Compute LDOS + //---------------------------------------------------------- + if (inp.out_ldos[0]) + { + ModuleIO::cal_ldos_pw(reinterpret_cast>*>(pelec), stp.psi_cpu[0], para_grid, ucell); + } //---------------------------------------------------------- //! 2) Calculate the spillage value, @@ -288,8 +269,7 @@ void ModuleIO::ctrl_runner_pw(UnitCell& ucell, { if (GlobalV::MY_RANK == 0) { - std::cout << "update value: bessel_nao_rcut <- " << std::fixed << inp.bessel_nao_rcuts[i] - << " a.u." << std::endl; + std::cout << "update value: bessel_nao_rcut <- " << std::fixed << inp.bessel_nao_rcuts[i] << " a.u." << std::endl; } Numerical_Basis numerical_basis; numerical_basis.output_overlap(stp.psi_cpu[0], sf, kv, pw_wfc, ucell, i); @@ -307,18 +287,15 @@ void ModuleIO::ctrl_runner_pw(UnitCell& ucell, ModuleIO::get_wf_pw(inp.out_wfc_norm, inp.out_wfc_re_im, - stp.get_nbands(), inp.nspin, - pw_rhod->nxyz, &ucell, stp.template get_psi_d(), pw_wfc, - ctx, + pw_rho, + pw_rhod, para_grid, PARAM.globalv.global_out_dir, - kv, - GlobalV::KPAR, - GlobalV::MY_POOL); + kv); } //---------------------------------------------------------- @@ -328,14 +305,8 @@ void ModuleIO::ctrl_runner_pw(UnitCell& ucell, { using Real = typename GetTypeReal::type; EleCond elec_cond(&ucell, &kv, pelec, pw_wfc, stp.template get_psi_t(), &ppcell); - elec_cond.KG(inp.cond_smear, - inp.cond_fwhm, - inp.cond_wcut, - inp.cond_dw, - inp.cond_dt, - inp.cond_nonlocal, - inp.cond_mgga_vel, - pelec->wg); + elec_cond + .KG(inp.cond_smear, inp.cond_fwhm, inp.cond_wcut, inp.cond_dw, inp.cond_dt, inp.cond_nonlocal, inp.cond_mgga_vel, pelec->wg); } #ifdef __MLALGO @@ -380,129 +351,121 @@ void ModuleIO::ctrl_runner_pw(UnitCell& ucell, } // complex + CPU -template void ModuleIO::ctrl_scf_pw, base_device::DEVICE_CPU>( - const int nstep, - UnitCell& ucell, - elecstate::ElecState* pelec, - const Charge &chr, - const K_Vectors &kv, - const ModulePW::PW_Basis_K *pw_wfc, - const ModulePW::PW_Basis *pw_rho, - const ModulePW::PW_Basis *pw_rhod, - const ModulePW::PW_Basis_Big *pw_big, - Setup_Psi_pw &stp, - const Parallel_Grid ¶_grid, - const Input_para& inp); +template void ModuleIO::ctrl_scf_pw, base_device::DEVICE_CPU>(const int nstep, + UnitCell& ucell, + elecstate::ElecState* pelec, + const Charge& chr, + const K_Vectors& kv, + const ModulePW::PW_Basis_K* pw_wfc, + const ModulePW::PW_Basis* pw_rho, + const ModulePW::PW_Basis* pw_rhod, + const ModulePW::PW_Basis_Big* pw_big, + Setup_Psi_pw& stp, + const Parallel_Grid& para_grid, + const Input_para& inp); // complex + CPU -template void ModuleIO::ctrl_scf_pw, base_device::DEVICE_CPU>( - const int nstep, - UnitCell& ucell, - elecstate::ElecState* pelec, - const Charge &chr, - const K_Vectors &kv, - const ModulePW::PW_Basis_K *pw_wfc, - const ModulePW::PW_Basis *pw_rho, - const ModulePW::PW_Basis *pw_rhod, - const ModulePW::PW_Basis_Big *pw_big, - Setup_Psi_pw &stp, - const Parallel_Grid ¶_grid, - const Input_para& inp); +template void ModuleIO::ctrl_scf_pw, base_device::DEVICE_CPU>(const int nstep, + UnitCell& ucell, + elecstate::ElecState* pelec, + const Charge& chr, + const K_Vectors& kv, + const ModulePW::PW_Basis_K* pw_wfc, + const ModulePW::PW_Basis* pw_rho, + const ModulePW::PW_Basis* pw_rhod, + const ModulePW::PW_Basis_Big* pw_big, + Setup_Psi_pw& stp, + const Parallel_Grid& para_grid, + const Input_para& inp); #if ((defined __CUDA) || (defined __ROCM)) // complex + GPU -template void ModuleIO::ctrl_scf_pw, base_device::DEVICE_GPU>( - const int nstep, - UnitCell& ucell, - elecstate::ElecState* pelec, - const Charge &chr, - const K_Vectors &kv, - const ModulePW::PW_Basis_K *pw_wfc, - const ModulePW::PW_Basis *pw_rho, - const ModulePW::PW_Basis *pw_rhod, - const ModulePW::PW_Basis_Big *pw_big, - Setup_Psi_pw &stp, - const Parallel_Grid ¶_grid, - const Input_para& inp); +template void ModuleIO::ctrl_scf_pw, base_device::DEVICE_GPU>(const int nstep, + UnitCell& ucell, + elecstate::ElecState* pelec, + const Charge& chr, + const K_Vectors& kv, + const ModulePW::PW_Basis_K* pw_wfc, + const ModulePW::PW_Basis* pw_rho, + const ModulePW::PW_Basis* pw_rhod, + const ModulePW::PW_Basis_Big* pw_big, + Setup_Psi_pw& stp, + const Parallel_Grid& para_grid, + const Input_para& inp); // complex + GPU -template void ModuleIO::ctrl_scf_pw, base_device::DEVICE_GPU>( - const int nstep, - UnitCell& ucell, - elecstate::ElecState* pelec, - const Charge &chr, - const K_Vectors &kv, - const ModulePW::PW_Basis_K *pw_wfc, - const ModulePW::PW_Basis *pw_rho, - const ModulePW::PW_Basis *pw_rhod, - const ModulePW::PW_Basis_Big *pw_big, - Setup_Psi_pw &stp, - const Parallel_Grid ¶_grid, - const Input_para& inp); +template void ModuleIO::ctrl_scf_pw, base_device::DEVICE_GPU>(const int nstep, + UnitCell& ucell, + elecstate::ElecState* pelec, + const Charge& chr, + const K_Vectors& kv, + const ModulePW::PW_Basis_K* pw_wfc, + const ModulePW::PW_Basis* pw_rho, + const ModulePW::PW_Basis* pw_rhod, + const ModulePW::PW_Basis_Big* pw_big, + Setup_Psi_pw& stp, + const Parallel_Grid& para_grid, + const Input_para& inp); #endif // complex + CPU -template void ModuleIO::ctrl_runner_pw, base_device::DEVICE_CPU>( - UnitCell& ucell, - elecstate::ElecState* pelec, - ModulePW::PW_Basis_K* pw_wfc, - ModulePW::PW_Basis* pw_rho, - ModulePW::PW_Basis* pw_rhod, - Charge &chr, - K_Vectors &kv, - Setup_Psi_pw &stp, - Structure_Factor &sf, - pseudopot_cell_vnl &ppcell, - surchem &solvent, - Parallel_Grid ¶_grid, - const Input_para& inp); +template void ModuleIO::ctrl_runner_pw, base_device::DEVICE_CPU>(UnitCell& ucell, + elecstate::ElecState* pelec, + ModulePW::PW_Basis_K* pw_wfc, + ModulePW::PW_Basis* pw_rho, + ModulePW::PW_Basis* pw_rhod, + Charge& chr, + K_Vectors& kv, + Setup_Psi_pw& stp, + Structure_Factor& sf, + pseudopot_cell_vnl& ppcell, + surchem& solvent, + Parallel_Grid& para_grid, + const Input_para& inp); // complex + CPU -template void ModuleIO::ctrl_runner_pw, base_device::DEVICE_CPU>( - UnitCell& ucell, - elecstate::ElecState* pelec, - ModulePW::PW_Basis_K* pw_wfc, - ModulePW::PW_Basis* pw_rho, - ModulePW::PW_Basis* pw_rhod, - Charge &chr, - K_Vectors &kv, - Setup_Psi_pw &stp, - Structure_Factor &sf, - pseudopot_cell_vnl &ppcell, - surchem &solvent, - Parallel_Grid ¶_grid, - const Input_para& inp); +template void ModuleIO::ctrl_runner_pw, base_device::DEVICE_CPU>(UnitCell& ucell, + elecstate::ElecState* pelec, + ModulePW::PW_Basis_K* pw_wfc, + ModulePW::PW_Basis* pw_rho, + ModulePW::PW_Basis* pw_rhod, + Charge& chr, + K_Vectors& kv, + Setup_Psi_pw& stp, + Structure_Factor& sf, + pseudopot_cell_vnl& ppcell, + surchem& solvent, + Parallel_Grid& para_grid, + const Input_para& inp); #if ((defined __CUDA) || (defined __ROCM)) // complex + GPU -template void ModuleIO::ctrl_runner_pw, base_device::DEVICE_GPU>( - UnitCell& ucell, - elecstate::ElecState* pelec, - ModulePW::PW_Basis_K* pw_wfc, - ModulePW::PW_Basis* pw_rho, - ModulePW::PW_Basis* pw_rhod, - Charge &chr, - K_Vectors &kv, - Setup_Psi_pw &stp, - Structure_Factor &sf, - pseudopot_cell_vnl &ppcell, - surchem &solvent, - Parallel_Grid ¶_grid, - const Input_para& inp); +template void ModuleIO::ctrl_runner_pw, base_device::DEVICE_GPU>(UnitCell& ucell, + elecstate::ElecState* pelec, + ModulePW::PW_Basis_K* pw_wfc, + ModulePW::PW_Basis* pw_rho, + ModulePW::PW_Basis* pw_rhod, + Charge& chr, + K_Vectors& kv, + Setup_Psi_pw& stp, + Structure_Factor& sf, + pseudopot_cell_vnl& ppcell, + surchem& solvent, + Parallel_Grid& para_grid, + const Input_para& inp); // complex + GPU -template void ModuleIO::ctrl_runner_pw, base_device::DEVICE_GPU>( - UnitCell& ucell, - elecstate::ElecState* pelec, - ModulePW::PW_Basis_K* pw_wfc, - ModulePW::PW_Basis* pw_rho, - ModulePW::PW_Basis* pw_rhod, - Charge &chr, - K_Vectors &kv, - Setup_Psi_pw &stp, - Structure_Factor &sf, - pseudopot_cell_vnl &ppcell, - surchem &solvent, - Parallel_Grid ¶_grid, - const Input_para& inp); +template void ModuleIO::ctrl_runner_pw, base_device::DEVICE_GPU>(UnitCell& ucell, + elecstate::ElecState* pelec, + ModulePW::PW_Basis_K* pw_wfc, + ModulePW::PW_Basis* pw_rho, + ModulePW::PW_Basis* pw_rhod, + Charge& chr, + K_Vectors& kv, + Setup_Psi_pw& stp, + Structure_Factor& sf, + pseudopot_cell_vnl& ppcell, + surchem& solvent, + Parallel_Grid& para_grid, + const Input_para& inp); #endif diff --git a/source/source_io/module_parameter/read_inp_out.cpp b/source/source_io/module_parameter/read_inp_out.cpp index dd5abd63372..7f7113bdb45 100644 --- a/source/source_io/module_parameter/read_inp_out.cpp +++ b/source/source_io/module_parameter/read_inp_out.cpp @@ -1495,7 +1495,9 @@ If EXX(exact exchange) is calculated (i.e. dft_fuctional==hse/hf/pbe0/scan0 or r item.annotation = "specify the bands to be calculated for the partial (band-decomposed) charge densities"; item.category = "Output information"; item.type = "String"; - item.description = R"(Specifies the electronic states to calculate the charge densities with state index for, using a space-separated string of 0s and 1s. Each digit in the string corresponds to a state, starting from the first state. A 1 indicates that the charge density should be calculated for that state, while a 0 means the state will be ignored. The parameter allows a compact and flexible notation (similar to ocp_set), for example the syntax 1 4*0 5*1 0 is used to denote the selection of states: 1 means calculate for the first state, 4*0 skips the next four states, 5*1 means calculate for the following five states, and the final 0 skips the next state. It's essential that the total count of states does not exceed the total number of states (nbands); otherwise, it results in an error, and the process exits. The input string must contain only numbers and the asterisk (*) for repetition, ensuring correct format and intention of state selection. The outputs comprise multiple .cube files following the naming convention pchgi[state]s[spin]k[kpoint].cube.)"; + item.description = R"(Selects electronic states for partial (band-decomposed) charge-density output using a space-separated string of `0`s and `1`s, where `1` selects a state and `0` skips it. Repetition follows the `ocp_set` syntax, for example `1 4*0 5*1 0`; the expanded list must not exceed `nbands`. Each output represents a complete one-particle state rather than its SCF occupation. The spin degeneracy is 2 for `nspin=1` and 1 for `nspin=2` or `nspin=4`. For `nspin=1`, `s1` contains the charge density. For `nspin=2`, `s1` and `s2` contain the spin-up and spin-down charge densities, respectively. For `nspin=4`, `s1`, `s2`, `s3`, and `s4` respectively contain $\rho_0$, $m_x$, $m_y$, and $m_z$. With `if_separate_k=true`, files are named `pchgi[state]s[component]k[kpoint].cube`; otherwise, the weighted k-point sum is named `pchgi[state]s[component].cube`. + +[NOTE] Enabling symmetry may produce unintended partial charge densities because of reduced k-point weights and real-space symmetry operations. If the desired symmetry treatment is uncertain, set `symmetry = -1`. Use the same symmetry setting as in the SCF calculation.)"; item.default_value = "none"; item.unit = ""; item.set_availability("basis_type==pw or (basis_type==lcao and calculation==get_pchg)"); @@ -1515,7 +1517,7 @@ If EXX(exact exchange) is calculated (i.e. dft_fuctional==hse/hf/pbe0/scan0 or r item.annotation = "specify the bands to be calculated for the norm of wavefunctions"; item.category = "Output information"; item.type = "String"; - item.description = "Specifies the electronic states to calculate the real-space wave function modulus (norm, or known as the envelope function) with state index. The syntax and state selection rules are identical to out_pchg, but the output is the norm of the wave function. The outputs comprise multiple .cube files following the naming convention wfi[state]s[spin]k[kpoint].cube."; + item.description = R"(Selects electronic states for real-space wavefunction-modulus output using the selection syntax and complete-state normalization of `out_pchg`. For `nspin=1`, `s1` contains the wavefunction modulus. For `nspin=2`, `s1` and `s2` contain the spin-up and spin-down wavefunction moduli, respectively. For `nspin=4`, `s1` contains the total spinor modulus. Files are named `wfi[state]s[spin]k[kpoint].cube`.)"; item.default_value = "none"; item.unit = ""; item.set_availability("basis_type==pw or (basis_type==lcao and calculation==get_wf)"); @@ -1536,7 +1538,7 @@ If EXX(exact exchange) is calculated (i.e. dft_fuctional==hse/hf/pbe0/scan0 or r item.annotation = "specify the bands to be calculated for the real and imaginary parts of wavefunctions"; item.category = "Output information"; item.type = "String"; - item.description = "Specifies the electronic states to calculate the real and imaginary parts of the wave function with state index. The syntax and state selection rules are identical to out_pchg, but the output contains both the real and imaginary components of the wave function. The outputs comprise multiple .cube files following the naming convention wfi[state]s[spin]k[kpoint][re/im].cube."; + item.description = R"(Selects electronic states for real-space wavefunction real- and imaginary-part output using the selection syntax and complete-state normalization of `out_pchg`. For `nspin=1`, `s1` contains the wavefunction. For `nspin=2`, `s1` and `s2` contain the spin-up and spin-down wavefunctions, respectively. For `nspin=4`, `s1` and `s2` contain the upper and lower spinor components, respectively. Files are named `wfi[state]s[spin]k[kpoint][re/im].cube`.)"; item.default_value = "none"; item.unit = ""; item.set_availability("basis_type==pw or (basis_type==lcao and calculation==get_wf)"); @@ -1558,7 +1560,8 @@ If EXX(exact exchange) is calculated (i.e. dft_fuctional==hse/hf/pbe0/scan0 or r "or merge them"; item.category = "Output information"; item.type = "Boolean"; - item.description = "Specifies whether to write the partial charge densities for all k-points to individual files or merge them. Warning: Enabling symmetry may produce unwanted results due to reduced k-point weights and symmetry operations in real space. Therefore when calculating partial charge densities, if you are not sure what you want exactly, it is strongly recommended to set symmetry = -1. It is noteworthy that your symmetry setting should remain the same as that in the SCF procedure."; + item.description + = "Specifies whether to write partial charge densities for individual k-points or merge them."; item.default_value = "false"; item.unit = ""; item.set_availability("(basis_type==pw and out_pchg!=none) or (basis_type==lcao and calculation==get_pchg and gamma_only==0)"); diff --git a/source/source_io/module_wf/get_wf_pw.h b/source/source_io/module_wf/get_wf_pw.h index 0b9461ea68b..551f7fab42b 100644 --- a/source/source_io/module_wf/get_wf_pw.h +++ b/source/source_io/module_wf/get_wf_pw.h @@ -1,278 +1,274 @@ #ifndef GET_WF_PW_H #define GET_WF_PW_H +#include "source_base/module_container/ATen/core/tensor.h" + namespace ModuleIO { +/** + * @brief Write real-space norms and complex components of selected PW states. + * + * Every selected band and k-point is written independently. Scalar or collinear + * states produce one component, while an nspin=4 spinor produces a combined norm + * cube and separate real/imaginary cubes for its up and down components. + */ template void get_wf_pw(const std::vector& out_wfc_norm, const std::vector& out_wfc_re_im, - const int nbands, const int nspin, - const int nxyz, UnitCell* ucell, const psi::Psi, Device>* kspw_psi, const ModulePW::PW_Basis_K* pw_wfc, - const Device* ctx, + const ModulePW::PW_Basis* pw_rho, + const ModulePW::PW_Basis* pw_rhod, const Parallel_Grid& pgrid, const std::string& global_out_dir, - const K_Vectors& kv, - const int kpar, - const int my_pool) + const K_Vectors& kv) { - // Get necessary parameters from kv const int nks = kv.get_nks(); // current process pool k-point count const int nkstot = kv.get_nkstot(); // total k-point count + const int nbands = kspw_psi->get_nbands(); - // Loop over k-parallelism - for (int ip = 0; ip < kpar; ++ip) - { - if (my_pool != ip) - { - continue; - } + const int nks_without_spin = nspin == 2 ? nkstot / 2 : nkstot; + const int smooth_nrxx = pw_wfc->nrxx; + const int dense_nrxx = pw_rhod->nrxx; + // Avoid an extra forward and backward FFT unless a distinct dense grid is in use. + const bool needs_interpolation = pw_rhod != pw_rho; - // bands_picked is a vector of 0s and 1s, where 1 means the band is picked to output - std::vector bands_picked_norm(nbands, 0); - std::vector bands_picked_re_im(nbands, 0); + // The two INPUT vectors are independent band masks for norm and Re/Im output. + std::vector bands_picked_norm(nbands, 0); + std::vector bands_picked_re_im(nbands, 0); - // Check if length of out_wfc_norm and out_wfc_re_im is valid - if (static_cast(out_wfc_norm.size()) > nbands || static_cast(out_wfc_re_im.size()) > nbands) - { - ModuleBase::WARNING_QUIT("ModuleIO::get_wf_pw", - "The number of bands specified by `out_wfc_norm` or `out_wfc_re_im` in the " - "INPUT file exceeds `nbands`!"); - } + if (static_cast(out_wfc_norm.size()) > nbands || static_cast(out_wfc_re_im.size()) > nbands) + { + ModuleBase::WARNING_QUIT("ModuleIO::get_wf_pw", + "The number of bands specified by `out_wfc_norm` or `out_wfc_re_im` in the " + "INPUT file exceeds `nbands`!"); + } - // Check if all elements in bands_picked are 0 or 1 - for (int value: out_wfc_norm) + for (int value: out_wfc_norm) + { + if (value != 0 && value != 1) { - if (value != 0 && value != 1) - { - ModuleBase::WARNING_QUIT("ModuleIO::get_wf_pw", - "The elements of `out_wfc_norm` must be either 0 or 1. " - "Invalid values found!"); - } + ModuleBase::WARNING_QUIT("ModuleIO::get_wf_pw", + "The elements of `out_wfc_norm` must be either 0 or 1. " + "Invalid values found!"); } - for (int value: out_wfc_re_im) + } + for (int value: out_wfc_re_im) + { + if (value != 0 && value != 1) { - if (value != 0 && value != 1) - { - ModuleBase::WARNING_QUIT("ModuleIO::get_wf_pw", - "The elements of `out_wfc_re_im` must be either 0 or 1. " - "Invalid values found!"); - } + ModuleBase::WARNING_QUIT("ModuleIO::get_wf_pw", + "The elements of `out_wfc_re_im` must be either 0 or 1. " + "Invalid values found!"); } + } + + int length = std::min(static_cast(out_wfc_norm.size()), nbands); + for (int i = 0; i < length; ++i) + { + bands_picked_norm[i] = static_cast(out_wfc_norm[i]); + } + length = std::min(static_cast(out_wfc_re_im.size()), nbands); + for (int i = 0; i < length; ++i) + { + bands_picked_re_im[i] = static_cast(out_wfc_re_im[i]); + } - // Fill bands_picked with values from out_wfc_norm - // Remaining bands are already set to 0 - int length = std::min(static_cast(out_wfc_norm.size()), nbands); - for (int i = 0; i < length; ++i) + // Map the wavefunction backend type to the tensor library's device type. + using ContainerDevice = typename ct::PsiToContainer::type; + const ct::DeviceType device_type = ct::DeviceTypeToEnum::value; + const bool is_cpu = device_type == ct::DeviceType::CpuDevice; + const bool is_spinor = nspin == 4; + // Spinor coefficients store the up and down blocks consecutively. + const int npwx = kspw_psi->get_nbasis() / (is_spinor ? 2 : 1); + + // Zero-length tensors avoid allocating buffers that a given backend or spin mode never uses. + ct::Tensor wfcr_up_smooth(ct::DataType::DT_COMPLEX_DOUBLE, device_type, ct::TensorShape({smooth_nrxx})); + ct::Tensor wfcr_down_smooth(ct::DataType::DT_COMPLEX_DOUBLE, device_type, ct::TensorShape({is_spinor ? smooth_nrxx : 0})); + ct::Tensor wfcr_up_smooth_host(ct::DataType::DT_COMPLEX_DOUBLE, ct::DeviceType::CpuDevice, ct::TensorShape({is_cpu ? 0 : smooth_nrxx})); + ct::Tensor wfcr_down_smooth_host(ct::DataType::DT_COMPLEX_DOUBLE, + ct::DeviceType::CpuDevice, + ct::TensorShape({!is_cpu && is_spinor ? smooth_nrxx : 0})); + ct::Tensor wfcr_up_dense_host(ct::DataType::DT_COMPLEX_DOUBLE, + ct::DeviceType::CpuDevice, + ct::TensorShape({needs_interpolation ? dense_nrxx : 0})); + ct::Tensor wfcr_down_dense_host(ct::DataType::DT_COMPLEX_DOUBLE, + ct::DeviceType::CpuDevice, + ct::TensorShape({needs_interpolation && is_spinor ? dense_nrxx : 0})); + ct::Tensor reciprocal_buffer_host(ct::DataType::DT_COMPLEX_DOUBLE, + ct::DeviceType::CpuDevice, + ct::TensorShape({needs_interpolation ? pw_rhod->npw : 0})); + + // Capture the shared bases and scratch buffers by reference. The returned pointer is owned by + // one of the tensor arguments and remains valid until that tensor is reused or destroyed. + auto transform_wfc = [&](const std::complex* coefficients, + const int ik, + ct::Tensor& smooth, + ct::Tensor& smooth_host, + ct::Tensor& dense_host) -> const std::complex* { + // Perform the wavefunction FFT on its native device, then expose host data for cube output. + pw_wfc->template recip_to_real, Device>(coefficients, smooth.data>(), ik); + const std::complex* smooth_data = smooth.data>(); + if (!is_cpu) { - // out_wfc_norm rely on function parse_expression - bands_picked_norm[i] = static_cast(out_wfc_norm[i]); + ct::kernels::synchronize_memory, ct::DEVICE_CPU, ContainerDevice>()( + smooth_host.data>(), + smooth_data, + smooth_nrxx); + smooth_data = smooth_host.data>(); } - length = std::min(static_cast(out_wfc_re_im.size()), nbands); - for (int i = 0; i < length; ++i) + if (!needs_interpolation) { - bands_picked_re_im[i] = static_cast(out_wfc_re_im[i]); + return smooth_data; } - // Allocate host memory - std::vector> wfcr_norm(nxyz); - std::vector> rho_band_norm(nspin, std::vector(nxyz)); + // Zero padding in reciprocal space transfers the smooth-grid field to the dense rho grid. + reciprocal_buffer_host.zero(); + pw_rho->real2recip(smooth_data, reciprocal_buffer_host.data>()); + pw_rhod->recip2real(reciprocal_buffer_host.data>(), dense_host.data>()); + return dense_host.data>(); + }; - // Allocate device memory - std::complex* wfcr_norm_device = nullptr; - if (!std::is_same::value) + // Norm cubes are phase invariant; for spinors they contain sqrt(|up|^2 + |down|^2). + std::vector> rho_band_norm(nspin, std::vector(dense_nrxx)); + for (int ib = 0; ib < nbands; ++ib) + { + if (!bands_picked_norm[ib]) { - base_device::memory::resize_memory_op, Device>()(wfcr_norm_device, nxyz); + continue; } - for (int ib = 0; ib < nbands; ++ib) + for (int is = 0; is < nspin; ++is) { - // Skip the loop iteration if bands_picked[ib] is 0 - if (!bands_picked_norm[ib]) + std::fill(rho_band_norm[is].begin(), rho_band_norm[is].end(), 0.0); + } + for (int ik = 0; ik < nks; ++ik) + { + const int ikstot = kv.ik2iktot[ik]; + const int spin_index = kv.isk[ik]; + // In collinear calculations the two spin channels share the same k-point numbering. + const int k_number = ikstot % nks_without_spin + 1; + + kspw_psi->fix_k(ik); + const std::complex* wfcr_up + = transform_wfc(&kspw_psi[0](ib, 0), ik, wfcr_up_smooth, wfcr_up_smooth_host, wfcr_up_dense_host); + const std::complex* wfcr_up_host_data = wfcr_up; + const std::complex* wfcr_down_host_data = nullptr; + if (is_spinor) { - continue; + const std::complex* wfcr_down + = transform_wfc(&kspw_psi[0](ib, npwx), ik, wfcr_down_smooth, wfcr_down_smooth_host, wfcr_down_dense_host); + wfcr_down_host_data = wfcr_down; } - for (int is = 0; is < nspin; ++is) + const double spin_degeneracy = nspin == 1 ? 2.0 : 1.0; + const double scale = std::sqrt(spin_degeneracy / ucell->omega); + for (int ir = 0; ir < dense_nrxx; ++ir) { - std::fill(rho_band_norm[is].begin(), rho_band_norm[is].end(), 0.0); + const double norm = is_spinor ? std::sqrt(std::norm(wfcr_up_host_data[ir]) + std::norm(wfcr_down_host_data[ir])) + : std::abs(wfcr_up_host_data[ir]); + rho_band_norm[spin_index][ir] = norm * scale; } - for (int ik = 0; ik < nks; ++ik) - { - const int ikstot = kv.ik2iktot[ik]; // global k-point index - const int spin_index = kv.isk[ik]; // spin index - const int k_number = ikstot % (nkstot / nspin) + 1; // k-point number, starting from 1 - - kspw_psi->fix_k(ik); - - // FFT on device and copy result back to host - if (std::is_same::value) - { - pw_wfc->recip_to_real(ctx, &kspw_psi[0](ib, 0), wfcr_norm.data(), ik); - } - else - { - pw_wfc->recip_to_real(ctx, &kspw_psi[0](ib, 0), wfcr_norm_device, ik); - - base_device::memory::synchronize_memory_op, base_device::DEVICE_CPU, Device>()( - wfcr_norm.data(), - wfcr_norm_device, - nxyz); - } - - // To ensure the normalization of charge density in multi-k calculation - double wg_sum_k = 0.0; - if (nspin == 1) - { - wg_sum_k = 2.0; - } - else if (nspin == 2) - { - wg_sum_k = 1.0; - } - else - { - ModuleBase::WARNING_QUIT("ModuleIO::get_wf_pw", - "Real space wavefunction output currently do not support noncollinear " - "polarized calculation (nspin = 4)!"); - } - - double w1 = static_cast(wg_sum_k / ucell->omega); - - for (int i = 0; i < nxyz; ++i) - { - rho_band_norm[spin_index][i] = std::abs(wfcr_norm[i]) * std::sqrt(w1); - } - std::stringstream ss_file; - ss_file << global_out_dir << "wfi" << ib + 1 << "s" << spin_index + 1 << "k" << k_number << ".cube"; - - ModuleIO::write_vdata_palgrid(pgrid, - rho_band_norm[spin_index].data(), - spin_index, - nspin, - 0, - ss_file.str(), - 0.0, - ucell, - 11, - 1, - PARAM.globalv.two_fermi, - true); // reduce_all_pool is true - } + std::stringstream ss_file; + ss_file << global_out_dir << "wfi" << ib + 1 << "s" << spin_index + 1 << "k" << k_number << ".cube"; + ModuleIO::write_vdata_palgrid(pgrid, + rho_band_norm[spin_index].data(), + spin_index, + nspin, + 0, + ss_file.str(), + 0.0, + ucell, + 11, + 0, + false, + true); } + } - // Allocate host memory - std::vector> wfc_re_im(nxyz); - std::vector> rho_band_re(nspin, std::vector(nxyz)); - std::vector> rho_band_im(nspin, std::vector(nxyz)); - - // Allocate device memory - std::complex* wfc_re_im_device = nullptr; - if (!std::is_same::value) + // Re/Im cubes retain phase information and therefore write both spinor components separately. + std::vector> rho_band_re(nspin, std::vector(dense_nrxx)); + std::vector> rho_band_im(nspin, std::vector(dense_nrxx)); + for (int ib = 0; ib < nbands; ++ib) + { + if (!bands_picked_re_im[ib]) { - base_device::memory::resize_memory_op, Device>()(wfc_re_im_device, nxyz); + continue; } - for (int ib = 0; ib < nbands; ++ib) + for (int is = 0; is < nspin; ++is) { - // Skip the loop iteration if bands_picked[ib] is 0 - if (!bands_picked_re_im[ib]) + std::fill(rho_band_re[is].begin(), rho_band_re[is].end(), 0.0); + std::fill(rho_band_im[is].begin(), rho_band_im[is].end(), 0.0); + } + for (int ik = 0; ik < nks; ++ik) + { + const int ikstot = kv.ik2iktot[ik]; + const int spin_index = kv.isk[ik]; + const int k_number = ikstot % nks_without_spin + 1; + + kspw_psi->fix_k(ik); + const std::complex* wfcr_up + = transform_wfc(&kspw_psi[0](ib, 0), ik, wfcr_up_smooth, wfcr_up_smooth_host, wfcr_up_dense_host); + const std::complex* wfcr_up_host_data = wfcr_up; + const std::complex* wfcr_down_host_data = nullptr; + if (is_spinor) { - continue; + const std::complex* wfcr_down + = transform_wfc(&kspw_psi[0](ib, npwx), ik, wfcr_down_smooth, wfcr_down_smooth_host, wfcr_down_dense_host); + wfcr_down_host_data = wfcr_down; } - for (int is = 0; is < nspin; ++is) - { - std::fill(rho_band_re[is].begin(), rho_band_re[is].end(), 0.0); - std::fill(rho_band_im[is].begin(), rho_band_im[is].end(), 0.0); - } - for (int ik = 0; ik < nks; ++ik) + const double spin_degeneracy = nspin == 1 ? 2.0 : 1.0; + const double scale = std::sqrt(spin_degeneracy / ucell->omega); + // For scalar/collinear states, select kv.isk[ik]; for spinors, emit both up and down. + const int component_begin = is_spinor ? 0 : spin_index; + const int component_end = is_spinor ? 2 : spin_index + 1; + for (int component = component_begin; component < component_end; ++component) { - const int ikstot = kv.ik2iktot[ik]; // global k-point index - const int spin_index = kv.isk[ik]; // spin index - const int k_number = ikstot % (nkstot / nspin) + 1; // k-point number, starting from 1 - - kspw_psi->fix_k(ik); - - // FFT on device and copy result back to host - if (std::is_same::value) + const std::complex* component_data = is_spinor && component == 1 ? wfcr_down_host_data : wfcr_up_host_data; + for (int ir = 0; ir < dense_nrxx; ++ir) { - pw_wfc->recip_to_real(ctx, &kspw_psi[0](ib, 0), wfc_re_im.data(), ik); - } - else - { - pw_wfc->recip_to_real(ctx, &kspw_psi[0](ib, 0), wfc_re_im_device, ik); - - base_device::memory::synchronize_memory_op, base_device::DEVICE_CPU, Device>()( - wfc_re_im.data(), - wfc_re_im_device, - nxyz); - } - - // To ensure the normalization of charge density in multi-k calculation - double wg_sum_k = 0.0; - if (nspin == 1) - { - wg_sum_k = 2.0; - } - else if (nspin == 2) - { - wg_sum_k = 1.0; - } - else - { - ModuleBase::WARNING_QUIT("ModuleIO::get_wf_pw", - "Real space wavefunction output currently do not support noncollinear " - "polarized calculation (nspin = 4)!"); - } - - double w1 = static_cast(wg_sum_k / ucell->omega); - - for (int i = 0; i < nxyz; ++i) - { - rho_band_re[spin_index][i] = std::real(wfc_re_im[i]) * std::sqrt(w1); - rho_band_im[spin_index][i] = std::imag(wfc_re_im[i]) * std::sqrt(w1); + rho_band_re[component][ir] = std::real(component_data[ir]) * scale; + rho_band_im[component][ir] = std::imag(component_data[ir]) * scale; } std::stringstream ss_real; - ss_real << global_out_dir << "wfi" << ib + 1 << "s" << spin_index + 1 << "k" << k_number << "re.cube"; - + ss_real << global_out_dir << "wfi" << ib + 1 << "s" << component + 1 << "k" << k_number << "re.cube"; ModuleIO::write_vdata_palgrid(pgrid, - rho_band_re[spin_index].data(), - spin_index, + rho_band_re[component].data(), + component, nspin, 0, ss_real.str(), 0.0, ucell, 11, - 1, - PARAM.globalv.two_fermi, - true); // reduce_all_pool is true + 0, + false, + true); std::stringstream ss_imag; - ss_imag << global_out_dir << "wfi" << ib + 1 << "s" << spin_index + 1 << "k" << k_number << "im.cube"; - + ss_imag << global_out_dir << "wfi" << ib + 1 << "s" << component + 1 << "k" << k_number << "im.cube"; ModuleIO::write_vdata_palgrid(pgrid, - rho_band_im[spin_index].data(), - spin_index, + rho_band_im[component].data(), + component, nspin, 0, ss_imag.str(), 0.0, ucell, 11, - 1, - PARAM.globalv.two_fermi, - true); // reduce_all_pool is true + 0, + false, + true); } } } } } // namespace ModuleIO -#endif // GET_WF_PW_H \ No newline at end of file +#endif // GET_WF_PW_H diff --git a/source/source_io/module_wf/read_wf2rho_pw.cpp b/source/source_io/module_wf/read_wf2rho_pw.cpp index 44b4ead083a..b66e519232d 100644 --- a/source/source_io/module_wf/read_wf2rho_pw.cpp +++ b/source/source_io/module_wf/read_wf2rho_pw.cpp @@ -1,54 +1,57 @@ #include "read_wf2rho_pw.h" #include "read_wfc_pw.h" +#include "source_base/module_out/filename.h" #include "source_base/timer.h" +#include "source_estate/kernels/elecstate_op.h" #include "source_estate/module_charge/symm_rho.h" #include "source_io/module_parameter/parameter.h" -#include "source_estate/kernels/elecstate_op.h" -#include "source_base/module_out/filename.h" -void ModuleIO::read_wf2rho_pw( - const ModulePW::PW_Basis_K* pw_wfc, - ModuleSymmetry::Symmetry& symm, - Charge& chg, - const std::string &readin_dir, - const int kpar, - const int my_pool, - const int my_rank, - const int nproc_in_pool, - const int rank_in_pool, - const int nbands, - const int nspin, - const int npol, - const int nkstot, - const std::vector &ik2iktot, - const std::vector &isk, - std::ofstream &ofs_running) +void ModuleIO::read_wf2rho_pw(const ModulePW::PW_Basis_K* pw_wfc, + ModuleSymmetry::Symmetry& symm, + Charge& chg, + const std::string& readin_dir, + const int kpar, + const int my_pool, + const int my_rank, + const int nproc_in_pool, + const int rank_in_pool, + const int nbands, + const int nspin, + const int npol, + const int nkstot, + const std::vector& ik2iktot, + const std::vector& isk, + std::ofstream& ofs_running) { ModuleBase::TITLE("ModuleIO", "read_wf2rho_pw"); ModuleBase::timer::start("ModuleIO", "read_wf2rho_pw"); - ofs_running << " READING WAVE FUNCTIONS" << std::endl; - ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" - ">>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; - ofs_running << " | " - " |" << std::endl; - ofs_running << " | Reading electronic wave functions in plane wave basis set and |" << std::endl; - ofs_running << " | evaluate charge density based on these wave functions |" << std::endl; - ofs_running << " | " - " |" << std::endl; - ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" - ">>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; - - assert(kpar>=1); - assert(my_pool>=0); - assert(my_rank>=0); - assert(nbands>0); - assert(nspin>0); - assert(npol==1 || npol==2); + ofs_running << " READING WAVE FUNCTIONS" << std::endl; + ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" + ">>>>>>>>>>>>>>>>>>>>>>>>>" + << std::endl; + ofs_running << " | " + " |" + << std::endl; + ofs_running << " | Reading electronic wave functions in plane wave basis set and |" << std::endl; + ofs_running << " | evaluate charge density based on these wave functions |" << std::endl; + ofs_running << " | " + " |" + << std::endl; + ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" + ">>>>>>>>>>>>>>>>>>>>>>>>>" + << std::endl; + + assert(kpar >= 1); + assert(my_pool >= 0); + assert(my_rank >= 0); + assert(nbands > 0); + assert(nspin > 0); + assert(npol == 1 || npol == 2); const int nrxx = pw_wfc->nrxx; - assert(nrxx>=0); + assert(nrxx >= 0); for (int is = 0; is < nspin; ++is) { @@ -67,30 +70,30 @@ void ModuleIO::read_wf2rho_pw( std::string filename = readin_dir + "eig.txt"; std::ifstream ifs(filename); - if(!ifs) - { + if (!ifs) + { std::stringstream sss; sss << "Cannot find file " << filename; - ModuleBase::WARNING_QUIT("ModuleIO::read_wf2rho_pw", sss.str()); - } + ModuleBase::WARNING_QUIT("ModuleIO::read_wf2rho_pw", sss.str()); + } else { ofs_running << " Find file containing weights of wave function: " << filename << std::endl; } - std::string useless; - getline(ifs, useless); - getline(ifs, useless); - for (int ik_tot = 0; ik_tot < nkstot; ++ik_tot) - { - ifs >> useless; - getline(ifs, useless); - for (int ib = 0; ib < nbands; ++ib) - { - ifs >> useless >> useless >> wg_tmp(ik_tot, ib); - } - } - } + std::string useless; + getline(ifs, useless); + getline(ifs, useless); + for (int ik_tot = 0; ik_tot < nkstot; ++ik_tot) + { + ifs >> useless; + getline(ifs, useless); + for (int ib = 0; ib < nbands; ++ib) + { + ifs >> useless >> useless >> wg_tmp(ik_tot, ib); + } + } + } #ifdef __MPI MPI_Bcast(wg_tmp.c, nkstot * nbands, MPI_DOUBLE, 0, MPI_COMM_WORLD); @@ -112,14 +115,11 @@ void ModuleIO::read_wf2rho_pw( const bool gamma_only = false; const int istep = -1; - std::string fn = filename_output(readin_dir,"wf","pw",ik,ik2iktot,nspin,nkstot, - out_type,out_app_flag,gamma_only,istep); + std::string fn = filename_output(readin_dir, "wf", "pw", ik, ik2iktot, nspin, nkstot, out_type, out_app_flag, gamma_only, istep); ofs_running << " Reading wave function from file: " << fn << std::endl; - ModuleIO::read_wfc_pw(fn, pw_wfc, - rank_in_pool, nproc_in_pool, nbands, npol, - ik, ikstot, nkstot, wfc_tmp); + ModuleIO::read_wfc_pw(fn, pw_wfc, rank_in_pool, nproc_in_pool, nbands, npol, ik, ikstot, nkstot, wfc_tmp); if (nspin == 4) { @@ -158,23 +158,21 @@ void ModuleIO::read_wf2rho_pw( if (w1 != 0.0) { - base_device::DEVICE_CPU* ctx = nullptr; - elecstate::elecstate_pw_op()(ctx, is, nrxx, nrxx, - w1, chg.rho, rho_tmp.data()); + base_device::DEVICE_CPU* ctx = nullptr; + elecstate::elecstate_pw_op()(ctx, is, nrxx, nrxx, w1, chg.rho, rho_tmp.data()); } } } } #ifdef __MPI - chg.init_chgmpi(); for (int is = 0; is < nspin; ++is) { chg.reduce_diff_pools(chg.rho[is]); } #endif - // Since rho is calculated by psi^2, it is not symmetric. We need to rearrange it. + // Since rho is calculated by psi^2, it is not symmetric. We need to rearrange it. Symmetry_rho srho; for (int is = 0; is < nspin; is++) { diff --git a/source/source_io/test/read_wf2rho_pw_test.cpp b/source/source_io/test/read_wf2rho_pw_test.cpp index 700fbd83e8d..38f586f3a5e 100644 --- a/source/source_io/test/read_wf2rho_pw_test.cpp +++ b/source/source_io/test/read_wf2rho_pw_test.cpp @@ -4,22 +4,22 @@ #undef __LCAO #define private public +#include "source_base/module_out/filename.h" // mohan add 2025-05-17 +#include "source_base/parallel_grid.h" #include "source_cell/klist.h" #include "source_cell/unitcell.h" #include "source_estate/module_charge/charge.h" #include "source_estate/module_charge/symm_rho.h" #include "source_hamilt/module_xc/xc_functional.h" -#include "source_base/parallel_grid.h" +#include "source_io/module_parameter/parameter.h" #include "source_io/module_wf/read_wf2rho_pw.h" #include "source_io/module_wf/write_wfc_pw.h" -#include "source_base/module_out/filename.h" // mohan add 2025-05-17 -#include "source_io/module_parameter/parameter.h" #include "source_psi/psi.h" #ifdef __MPI +#include "mpi.h" #include "source_base/parallel_global.h" #include "source_basis/module_pw/test/test_tool.h" -#include "mpi.h" #endif Charge::Charge() @@ -40,10 +40,18 @@ Magnetism::Magnetism() Magnetism::~Magnetism() { } -SepPot::SepPot(){} -SepPot::~SepPot(){} -Sep_Cell::Sep_Cell() noexcept {} -Sep_Cell::~Sep_Cell() noexcept {} +SepPot::SepPot() +{ +} +SepPot::~SepPot() +{ +} +Sep_Cell::Sep_Cell() noexcept +{ +} +Sep_Cell::~Sep_Cell() noexcept +{ +} int XC_Functional::func_type = 0; bool XC_Functional::ked_flag = false; @@ -53,47 +61,39 @@ Symmetry_rho::Symmetry_rho() Symmetry_rho::~Symmetry_rho() { } -void Symmetry_rho::begin(const int& spin_now, - const Charge& CHR, - const ModulePW::PW_Basis* rho_basis, - ModuleSymmetry::Symmetry& symm) const +void Symmetry_rho::begin(const int& spin_now, const Charge& CHR, const ModulePW::PW_Basis* rho_basis, ModuleSymmetry::Symmetry& symm) const { return; } void cal_ik2iktot(std::vector& ik2iktot, const int& nks, const int& nkstot) { - if(PARAM.inp.kpar==1) - { - for(int ik = 0; ik < nks; ++ik) - { - ik2iktot[ik] = ik; - } + if (PARAM.inp.kpar == 1) + { + for (int ik = 0; ik < nks; ++ik) + { + ik2iktot[ik] = ik; + } } - else if(PARAM.inp.kpar==2) + else if (PARAM.inp.kpar == 2) { - if(GlobalV::MY_POOL==0) - { - for(int ik = 0; ik < nks; ++ik) - { - ik2iktot[ik] = ik; - } - } - else if(GlobalV::MY_POOL==1) - { - for(int ik = 0; ik < nks; ++ik) - { - ik2iktot[ik] = ik+2; // only works for this test - } - } + if (GlobalV::MY_POOL == 0) + { + for (int ik = 0; ik < nks; ++ik) + { + ik2iktot[ik] = ik; + } + } + else if (GlobalV::MY_POOL == 1) + { + for (int ik = 0; ik < nks; ++ik) + { + ik2iktot[ik] = ik + 2; // only works for this test + } + } } } -namespace GlobalC -{ - Parallel_Grid Pgrid; -} // namespace GlobalC - /** * - Tested Functions: * - write_wfc_pw() @@ -109,6 +109,8 @@ class ReadWfcRhoTest : public ::testing::Test psi::Psi>* psi = nullptr; Charge chg; ModuleSymmetry::Symmetry symm; + Parallel_Grid pgrid; + int world_size = 1; virtual void SetUp() { wfcpw = new ModulePW::PW_Basis_K; @@ -116,6 +118,9 @@ class ReadWfcRhoTest : public ::testing::Test kv = new K_Vectors; // output .dat file PARAM.input.out_wfc_pw = 2; +#ifdef __MPI + MPI_Comm_size(MPI_COMM_WORLD, &world_size); +#endif } virtual void TearDown() { @@ -144,8 +149,7 @@ TEST_F(ReadWfcRhoTest, ReadWfcRho) kv->set_nks(nks); kv->isk = {0, 0}; const double shift = my_pool * 0.1; - kv->kvec_d = {ModuleBase::Vector3(shift, shift, shift), - ModuleBase::Vector3(0.5 + shift, 0.5 + shift, 0.5 + shift)}; + kv->kvec_d = {ModuleBase::Vector3(shift, shift, shift), ModuleBase::Vector3(0.5 + shift, 0.5 + shift, 0.5 + shift)}; kv->ik2iktot.resize(nks); cal_ik2iktot(kv->ik2iktot, nks, nkstot); @@ -160,6 +164,7 @@ TEST_F(ReadWfcRhoTest, ReadWfcRho) rhopw->initparameters(false, 80); rhopw->setuptransform(); rhopw->collect_local_pw(); + pgrid.init(rhopw->nx, rhopw->ny, rhopw->nz, rhopw->nplane, rhopw->nrxx, rhopw->nz, 1, world_size); wfcpw->initgrids(5.3233, ModuleBase::Matrix3(-0.5, 0.0, 0.5, 0.0, 0.5, 0.5, -0.5, 0.5, 0.0), 80); wfcpw->initparameters(false, 20, nks, kv->kvec_d.data()); @@ -193,20 +198,18 @@ TEST_F(ReadWfcRhoTest, ReadWfcRho) const int is = 0; // nspin is 1 for (int ik = 0; ik < nkstot; ++ik) - { - ofs << " spin=" << is+1 << " k-point=" - << ik + 1 << "/" << nkstot - << " Cartesian=" << kv->kvec_c[ik].x << " " << kv->kvec_c[ik].y - << " " << kv->kvec_c[ik].z << " (" << kv->ngk[ik] << " plane wave)" << std::endl; + { + ofs << " spin=" << is + 1 << " k-point=" << ik + 1 << "/" << nkstot << " Cartesian=" << kv->kvec_c[ik].x << " " + << kv->kvec_c[ik].y << " " << kv->kvec_c[ik].z << " (" << kv->ngk[ik] << " plane wave)" << std::endl; ofs << std::setprecision(16); ofs << std::setiosflags(std::ios::showpoint); double ekb = -1.23456; // energy for (int ib = 0; ib < nbands; ib++) - { - ofs << " " << ib + 1 << " " << ekb << " " << wg(ik,ib) << std::endl; - } + { + ofs << " " << ib + 1 << " " << ekb << " " << wg(ik, ib) << std::endl; + } ofs << std::endl; } @@ -232,6 +235,7 @@ TEST_F(ReadWfcRhoTest, ReadWfcRho) ModuleBase::GlobalFunc::ZEROS(chg.rho[0], rhopw->nrxx); chg.rhopw = rhopw; chg.nrxx = rhopw->nrxx; + chg.pgrid = &pgrid; //---------------------------------------- // set charge_ref @@ -243,6 +247,7 @@ TEST_F(ReadWfcRhoTest, ReadWfcRho) ModuleBase::GlobalFunc::ZEROS(chg_ref.rho[0], rhopw->nrxx); std::vector> rho_tmp(rhopw->nrxx); chg_ref.nrxx = rhopw->nrxx; + chg_ref.pgrid = &pgrid; for (int ik = 0; ik < nks; ++ik) { @@ -261,41 +266,61 @@ TEST_F(ReadWfcRhoTest, ReadWfcRho) } #ifdef __MPI - chg_ref.init_chgmpi(); chg_ref.reduce_diff_pools(chg_ref.rho[0]); #endif // for spin=1 or 2, npol=1 - const int npol=1; + const int npol = 1; // Write the wave functions to file - const std::string out_dir = "./"; + const std::string out_dir = "./"; // Read the wave functions to charge density std::stringstream ss; ss << "running_log" << GlobalV::MY_RANK << ".txt"; - std::ofstream running_log(ss.str().c_str()); + std::ofstream running_log(ss.str().c_str()); running_log << " rank=" << GlobalV::MY_RANK << std::endl; - const double ecutwfc = 20; // this is a fake number - const int istep = -1; // -1 means ionic iteration number will not appear in file name - const int iter = -1; // -1 means electronic iteration number will not appear in file name - - ModuleIO::write_wfc_pw(istep, iter, - kpar, my_pool, my_rank, nbands, nspin, npol, - GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, - PARAM.input.out_wfc_pw, ecutwfc, out_dir, *psi, *kv, wfcpw, - running_log); - - ModuleIO::read_wf2rho_pw(wfcpw, symm, chg, - out_dir, kpar, my_pool, my_rank, - GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL, - nbands, nspin, npol, - nkstot, kv->ik2iktot, kv->isk, running_log); + const int iter = -1; // -1 means electronic iteration number will not appear in file name + + ModuleIO::write_wfc_pw(istep, + iter, + kpar, + my_pool, + my_rank, + nbands, + nspin, + npol, + GlobalV::RANK_IN_POOL, + GlobalV::NPROC_IN_POOL, + PARAM.input.out_wfc_pw, + ecutwfc, + out_dir, + *psi, + *kv, + wfcpw, + running_log); + + ModuleIO::read_wf2rho_pw(wfcpw, + symm, + chg, + out_dir, + kpar, + my_pool, + my_rank, + GlobalV::NPROC_IN_POOL, + GlobalV::RANK_IN_POOL, + nbands, + nspin, + npol, + nkstot, + kv->ik2iktot, + kv->isk, + running_log); // compare the charge density for (int ir = 0; ir < rhopw->nrxx; ++ir) @@ -303,24 +328,24 @@ TEST_F(ReadWfcRhoTest, ReadWfcRho) EXPECT_NEAR(chg.rho[0][ir], chg_ref.rho[0][ir], 1e-8); } - if (GlobalV::NPROC == 1) - { - EXPECT_NEAR(chg.rho[0][0], 8617.076357957576, 1e-8); - } - else if (GlobalV::NPROC == 4) - { - const std::vector ref = {8207.849135313403, 35.34776105132742, 8207.849135313403, 35.34776105132742}; - EXPECT_NEAR(chg.rho[0][0], ref[GlobalV::MY_RANK], 1e-8); - // for (int ip = 0; ip < GlobalV::NPROC; ++ip) - // { - // if (GlobalV::MY_RANK == ip) - // { - // std::cout.precision(16); - // std::cout << GlobalV::MY_RANK << " " << chg.rho[0][0] << std::endl; - // } - // MPI_Barrier(MPI_COMM_WORLD); - // } - } + if (GlobalV::NPROC == 1) + { + EXPECT_NEAR(chg.rho[0][0], 8617.076357957576, 1e-8); + } + else if (GlobalV::NPROC == 4) + { + const std::vector ref = {8207.849135313403, 35.34776105132742, 8207.849135313403, 35.34776105132742}; + EXPECT_NEAR(chg.rho[0][0], ref[GlobalV::MY_RANK], 1e-8); + // for (int ip = 0; ip < GlobalV::NPROC; ++ip) + // { + // if (GlobalV::MY_RANK == ip) + // { + // std::cout.precision(16); + // std::cout << GlobalV::MY_RANK << " " << chg.rho[0][0] << std::endl; + // } + // MPI_Barrier(MPI_COMM_WORLD); + // } + } delete[] chg.rho; delete[] chg._space_rho; @@ -337,11 +362,11 @@ TEST_F(ReadWfcRhoTest, ReadWfcRho) if (GlobalV::KPAR == 2) { remove("wfs1k3_pw.dat"); - remove("wfs1k4_pw.dat"); - remove("running_log1.txt"); - remove("running_log2.txt"); - remove("running_log3.txt"); - } + remove("wfs1k4_pw.dat"); + remove("running_log1.txt"); + remove("running_log2.txt"); + remove("running_log3.txt"); + } } } diff --git a/tests/01_PW/085_PW_get_pchg/INPUT b/tests/01_PW/085_PW_get_pchg/INPUT index 9b6512f9e3f..2ccdc452305 100644 --- a/tests/01_PW/085_PW_get_pchg/INPUT +++ b/tests/01_PW/085_PW_get_pchg/INPUT @@ -10,6 +10,7 @@ pseudo_dir ../../PP_ORB #Parameters (2.Iteration) ecutwfc 20 +ecutrho 160 scf_thr 1e-9 scf_nmax 100 diff --git a/tests/01_PW/085_PW_get_pchg/result.ref b/tests/01_PW/085_PW_get_pchg/result.ref index f2c31d7ea11..0bbb9a20a62 100644 --- a/tests/01_PW/085_PW_get_pchg/result.ref +++ b/tests/01_PW/085_PW_get_pchg/result.ref @@ -1,5 +1,5 @@ -etotref -197.1405644417785 -etotperatomref -98.5702822209 +etotref -197.1406124279196 +etotperatomref -98.5703062140 pchgi1s1.cube 2 pchgi2s1.cube 2 pchgi3s1.cube 2 diff --git a/tests/01_PW/086_PW_get_wf/INPUT b/tests/01_PW/086_PW_get_wf/INPUT index 6949f9d6970..d49c72a75a4 100644 --- a/tests/01_PW/086_PW_get_wf/INPUT +++ b/tests/01_PW/086_PW_get_wf/INPUT @@ -10,6 +10,7 @@ pseudo_dir ../../PP_ORB #Parameters (2.Iteration) ecutwfc 20 +ecutrho 160 scf_thr 1e-9 scf_nmax 100 @@ -17,6 +18,6 @@ scf_nmax 100 basis_type pw out_wfc_norm 1 3*0 -out_wfc_re_im 0 1 2*0 +out_wfc_re_im 1 3*0 pw_seed 1 diff --git a/tests/01_PW/086_PW_get_wf/README b/tests/01_PW/086_PW_get_wf/README index f663a595df4..257054cc081 100644 --- a/tests/01_PW/086_PW_get_wf/README +++ b/tests/01_PW/086_PW_get_wf/README @@ -1 +1 @@ -test the output of out_wfc_norm and out_wfc_re_im +Test `out_wfc_norm` and `out_wfc_re_im` with the nondegenerate first band. diff --git a/tests/01_PW/086_PW_get_wf/result.ref b/tests/01_PW/086_PW_get_wf/result.ref index 7a951b5f216..01e414180a7 100644 --- a/tests/01_PW/086_PW_get_wf/result.ref +++ b/tests/01_PW/086_PW_get_wf/result.ref @@ -1,9 +1,37 @@ -etotref -197.1405644417785 -etotperatomref -98.5702822209 -wfi1s1k1.cube 22.08506167 -wfi2s1k1im.cube 7.845633773e-07 -wfi2s1k1re.cube 2.250475858e-07 +etotref -197.1406124279197 +etotperatomref -98.5703062140 +wfi1s1k1.cube 22.08506772 +wfi1s1k1_wfc_fp_components 1 +wfi1s1k1_wfc_fp_nx 30 +wfi1s1k1_wfc_fp_ny 30 +wfi1s1k1_wfc_fp_nz 30 +wfi1s1k1_wfc_fp_voxel 9.826000000000e-03 +wfi1s1k1_wfc_fp_rms 8.682498877844e-02 +wfi1s1k1_wfc_fp_power_0 4.862723210811e-01 +wfi1s1k1_wfc_fp_power_1 3.692234187308e-01 +wfi1s1k1_wfc_fp_power_2 9.149814627242e-01 +wfi1s1k1_wfc_fp_power_3 1.855500021326e-01 +wfi1s1k1_wfc_fp_power_4 9.877998067367e-02 +wfi1s1k1_wfc_fp_power_5 1.684602645445e+00 +wfi1s1k1_wfc_fp_power_6 2.661345668253e+00 +wfi1s1k1_wfc_fp_power_7 2.085195698243e-01 +wfi1s1k1_wfc_fp_cross_0_re -4.225603956740e-01 +wfi1s1k1_wfc_fp_cross_0_im 3.139810235665e-02 +wfi1s1k1_wfc_fp_cross_1_re -4.837558034152e-01 +wfi1s1k1_wfc_fp_cross_1_im -3.222001030485e-01 +wfi1s1k1_wfc_fp_cross_2_re -3.447334705779e-01 +wfi1s1k1_wfc_fp_cross_2_im 2.256848391521e-01 +wfi1s1k1_wfc_fp_cross_3_re 1.638891990120e-02 +wfi1s1k1_wfc_fp_cross_3_im -1.343876070519e-01 +wfi1s1k1_wfc_fp_cross_4_re -3.607389513824e-01 +wfi1s1k1_wfc_fp_cross_4_im -1.904532113547e-01 +wfi1s1k1_wfc_fp_cross_5_re -1.663844914021e+00 +wfi1s1k1_wfc_fp_cross_5_im -1.309553380075e+00 +wfi1s1k1_wfc_fp_cross_6_re -7.182870633772e-01 +wfi1s1k1_wfc_fp_cross_6_im -1.975002493235e-01 +wfi1s1k1_wfc_fp_cross_7_re 3.169996514777e-01 +wfi1s1k1_wfc_fp_cross_7_im 3.014160201981e-02 pointgroupref T_d spacegroupref O_h nksibzref 1 -totaltimeref 0.35 +totaltimeref 0.45 diff --git a/tests/01_PW/090_PW_out_pchg_wfc_spinor/INPUT b/tests/01_PW/090_PW_out_pchg_wfc_spinor/INPUT new file mode 100644 index 00000000000..312c8daf5e3 --- /dev/null +++ b/tests/01_PW/090_PW_out_pchg_wfc_spinor/INPUT @@ -0,0 +1,42 @@ +INPUT_PARAMETERS + +init_wfc random +# pw scf non-collinear calculations +basis_type pw +calculation scf +noncolin 1 +nspin 4 +lspinorb 1 +symmetry 0 + +# plane wave and scf +ecutwfc 20 +ecutrho 160 +pw_seed 1 +scf_thr 1.0e-7 +scf_nmax 1 +out_chg 0 +out_pchg 1 +if_separate_k 1 +out_wfc_norm 1 +out_wfc_re_im 1 + +# smearing method +smearing_method gaussian +smearing_sigma 0.02 + +# charge mixing method +mixing_type broyden +mixing_beta 0.2 +mixing_ndim 10 + +# diagonalizaiton +ks_solver dav_subspace +pw_diag_thr 1.0e-12 +diago_smooth_ethr true +pw_diag_ndim 2 + +# directories +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB +suffix autotest diff --git a/tests/01_PW/090_PW_out_pchg_wfc_spinor/KPT b/tests/01_PW/090_PW_out_pchg_wfc_spinor/KPT new file mode 100644 index 00000000000..b5b3bdb1ae2 --- /dev/null +++ b/tests/01_PW/090_PW_out_pchg_wfc_spinor/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 2 0 0 0 diff --git a/tests/01_PW/090_PW_out_pchg_wfc_spinor/README b/tests/01_PW/090_PW_out_pchg_wfc_spinor/README new file mode 100644 index 00000000000..bdf156b3e60 --- /dev/null +++ b/tests/01_PW/090_PW_out_pchg_wfc_spinor/README @@ -0,0 +1 @@ +Test CPU nspin=4 real-space partial-charge and spinor-wavefunction output from a one-step SOC operator probe. diff --git a/tests/01_PW/090_PW_out_pchg_wfc_spinor/STRU b/tests/01_PW/090_PW_out_pchg_wfc_spinor/STRU new file mode 100644 index 00000000000..da4cd00dab3 --- /dev/null +++ b/tests/01_PW/090_PW_out_pchg_wfc_spinor/STRU @@ -0,0 +1,21 @@ +ATOMIC_SPECIES +Fe 1.000 Fe.upf + +NUMERICAL_ORBITAL +Fe_gga_9au_100Ry_4s2p2d1f.orb + +LATTICE_CONSTANT +6 + +LATTICE_VECTORS +1.00 0.00 0.00 +0.00 1.00 0.00 +0.00 0.00 1.00 +ATOMIC_POSITIONS +Direct + +Fe +1.0 +2 +0.00 0.00 0.00 1 1 1 mag 1.0 1.0 1.0 +0.51 0.52 0.53 1 1 1 mag 1.0 1.0 1.0 diff --git a/tests/01_PW/090_PW_out_pchg_wfc_spinor/result.ref b/tests/01_PW/090_PW_out_pchg_wfc_spinor/result.ref new file mode 100644 index 00000000000..be1bec5394f --- /dev/null +++ b/tests/01_PW/090_PW_out_pchg_wfc_spinor/result.ref @@ -0,0 +1,74 @@ +etotref -6448.107384046055 +etotperatomref -3224.0536920230 +pchgi1s1k1.cube 1 +pchgi1s1k2.cube 1 +pchgi1s2k1.cube 0.5773502772 +pchgi1s2k2.cube 0.577350264 +pchgi1s3k1.cube 0.5773502629 +pchgi1s3k2.cube 0.5773502607 +pchgi1s4k1.cube 0.5773502668 +pchgi1s4k2.cube 0.5773502817 +wfi1s1k1.cube 5.895913845 +wfi1s1k2.cube 5.83429784 +wfi1k1_spinor_wfc_fp_components 2 +wfi1k1_spinor_wfc_fp_nx 25 +wfi1k1_spinor_wfc_fp_ny 25 +wfi1k1_spinor_wfc_fp_nz 25 +wfi1k1_spinor_wfc_fp_voxel 1.382400000000e-02 +wfi1k1_spinor_wfc_fp_rms 4.811252243247e-02 +wfi1k1_spinor_wfc_fp_power_0 4.823186182982e-01 +wfi1k1_spinor_wfc_fp_power_1 4.919438297146e+00 +wfi1k1_spinor_wfc_fp_power_2 2.829710784858e-01 +wfi1k1_spinor_wfc_fp_power_3 2.081068044479e+00 +wfi1k1_spinor_wfc_fp_power_4 5.477005767972e-02 +wfi1k1_spinor_wfc_fp_power_5 3.369275024525e-01 +wfi1k1_spinor_wfc_fp_power_6 2.190307287975e+00 +wfi1k1_spinor_wfc_fp_power_7 2.697692375202e-01 +wfi1k1_spinor_wfc_fp_cross_0_re -6.212963591714e-01 +wfi1k1_spinor_wfc_fp_cross_0_im -1.409513219648e+00 +wfi1k1_spinor_wfc_fp_cross_1_re -1.148921954787e+00 +wfi1k1_spinor_wfc_fp_cross_1_im -2.683972844432e-01 +wfi1k1_spinor_wfc_fp_cross_2_re -2.872580266008e-02 +wfi1k1_spinor_wfc_fp_cross_2_im 7.668486794734e-01 +wfi1k1_spinor_wfc_fp_cross_3_re 1.354870440135e-01 +wfi1k1_spinor_wfc_fp_cross_3_im 3.092304605566e-01 +wfi1k1_spinor_wfc_fp_cross_4_re 8.113881396438e-02 +wfi1k1_spinor_wfc_fp_cross_4_im 1.089496746744e-01 +wfi1k1_spinor_wfc_fp_cross_5_re 3.556879566582e-01 +wfi1k1_spinor_wfc_fp_cross_5_im -7.819596163673e-01 +wfi1k1_spinor_wfc_fp_cross_6_re 3.413052048490e-01 +wfi1k1_spinor_wfc_fp_cross_6_im -6.887585093159e-01 +wfi1k1_spinor_wfc_fp_cross_7_re -3.430393018621e-01 +wfi1k1_spinor_wfc_fp_cross_7_im -1.115292037007e-01 +wfi1k2_spinor_wfc_fp_components 2 +wfi1k2_spinor_wfc_fp_nx 25 +wfi1k2_spinor_wfc_fp_ny 25 +wfi1k2_spinor_wfc_fp_nz 25 +wfi1k2_spinor_wfc_fp_voxel 1.382400000000e-02 +wfi1k2_spinor_wfc_fp_rms 4.811252243247e-02 +wfi1k2_spinor_wfc_fp_power_0 1.148625840151e-01 +wfi1k2_spinor_wfc_fp_power_1 2.461537634946e+00 +wfi1k2_spinor_wfc_fp_power_2 1.590561801657e-02 +wfi1k2_spinor_wfc_fp_power_3 2.853259391276e+00 +wfi1k2_spinor_wfc_fp_power_4 8.924986887110e-01 +wfi1k2_spinor_wfc_fp_power_5 6.668603941489e-01 +wfi1k2_spinor_wfc_fp_power_6 3.954039012731e+00 +wfi1k2_spinor_wfc_fp_power_7 2.795818339018e-01 +wfi1k2_spinor_wfc_fp_cross_0_re 5.016579770209e-01 +wfi1k2_spinor_wfc_fp_cross_0_im -1.762891020220e-01 +wfi1k2_spinor_wfc_fp_cross_1_re 1.800957424658e-01 +wfi1k2_spinor_wfc_fp_cross_1_im -8.196219189708e-02 +wfi1k2_spinor_wfc_fp_cross_2_re 1.216731361730e-01 +wfi1k2_spinor_wfc_fp_cross_2_im -1.748670978591e-01 +wfi1k2_spinor_wfc_fp_cross_3_re 5.848582964961e-01 +wfi1k2_spinor_wfc_fp_cross_3_im 1.484746119135e+00 +wfi1k2_spinor_wfc_fp_cross_4_re 4.603614834245e-01 +wfi1k2_spinor_wfc_fp_cross_4_im 6.190632697151e-01 +wfi1k2_spinor_wfc_fp_cross_5_re 1.124163447259e+00 +wfi1k2_spinor_wfc_fp_cross_5_im -1.171771546999e+00 +wfi1k2_spinor_wfc_fp_cross_6_re -1.028088990439e+00 +wfi1k2_spinor_wfc_fp_cross_6_im -2.202510073458e-01 +wfi1k2_spinor_wfc_fp_cross_7_re -1.757685975816e-01 +wfi1k2_spinor_wfc_fp_cross_7_im -3.491263367086e-02 +pw_spinor_cube_identity 0 +totaltimeref 0.82 diff --git a/tests/01_PW/090_PW_VWR/INPUT b/tests/01_PW/091_PW_VWR/INPUT similarity index 100% rename from tests/01_PW/090_PW_VWR/INPUT rename to tests/01_PW/091_PW_VWR/INPUT diff --git a/tests/01_PW/090_PW_VWR/KPT b/tests/01_PW/091_PW_VWR/KPT similarity index 100% rename from tests/01_PW/090_PW_VWR/KPT rename to tests/01_PW/091_PW_VWR/KPT diff --git a/tests/01_PW/090_PW_VWR/README b/tests/01_PW/091_PW_VWR/README similarity index 100% rename from tests/01_PW/090_PW_VWR/README rename to tests/01_PW/091_PW_VWR/README diff --git a/tests/01_PW/090_PW_VWR/STRU b/tests/01_PW/091_PW_VWR/STRU similarity index 100% rename from tests/01_PW/090_PW_VWR/STRU rename to tests/01_PW/091_PW_VWR/STRU diff --git a/tests/01_PW/090_PW_VWR/result.ref b/tests/01_PW/091_PW_VWR/result.ref similarity index 100% rename from tests/01_PW/090_PW_VWR/result.ref rename to tests/01_PW/091_PW_VWR/result.ref diff --git a/tests/01_PW/090_PW_VWR/run.sh b/tests/01_PW/091_PW_VWR/run.sh similarity index 100% rename from tests/01_PW/090_PW_VWR/run.sh rename to tests/01_PW/091_PW_VWR/run.sh diff --git a/tests/01_PW/091_PW_CR_VDW3/INPUT b/tests/01_PW/092_PW_CR_VDW3/INPUT similarity index 100% rename from tests/01_PW/091_PW_CR_VDW3/INPUT rename to tests/01_PW/092_PW_CR_VDW3/INPUT diff --git a/tests/01_PW/091_PW_CR_VDW3/KPT b/tests/01_PW/092_PW_CR_VDW3/KPT similarity index 100% rename from tests/01_PW/091_PW_CR_VDW3/KPT rename to tests/01_PW/092_PW_CR_VDW3/KPT diff --git a/tests/01_PW/091_PW_CR_VDW3/README b/tests/01_PW/092_PW_CR_VDW3/README similarity index 100% rename from tests/01_PW/091_PW_CR_VDW3/README rename to tests/01_PW/092_PW_CR_VDW3/README diff --git a/tests/01_PW/091_PW_CR_VDW3/STRU b/tests/01_PW/092_PW_CR_VDW3/STRU similarity index 100% rename from tests/01_PW/091_PW_CR_VDW3/STRU rename to tests/01_PW/092_PW_CR_VDW3/STRU diff --git a/tests/01_PW/091_PW_CR_VDW3/result.ref b/tests/01_PW/092_PW_CR_VDW3/result.ref similarity index 100% rename from tests/01_PW/091_PW_CR_VDW3/result.ref rename to tests/01_PW/092_PW_CR_VDW3/result.ref diff --git a/tests/01_PW/092_PW_MSST/INPUT b/tests/01_PW/093_PW_MSST/INPUT similarity index 100% rename from tests/01_PW/092_PW_MSST/INPUT rename to tests/01_PW/093_PW_MSST/INPUT diff --git a/tests/01_PW/092_PW_MSST/KPT b/tests/01_PW/093_PW_MSST/KPT similarity index 100% rename from tests/01_PW/092_PW_MSST/KPT rename to tests/01_PW/093_PW_MSST/KPT diff --git a/tests/01_PW/092_PW_MSST/README b/tests/01_PW/093_PW_MSST/README similarity index 100% rename from tests/01_PW/092_PW_MSST/README rename to tests/01_PW/093_PW_MSST/README diff --git a/tests/01_PW/092_PW_MSST/STRU b/tests/01_PW/093_PW_MSST/STRU similarity index 100% rename from tests/01_PW/092_PW_MSST/STRU rename to tests/01_PW/093_PW_MSST/STRU diff --git a/tests/01_PW/092_PW_MSST/result.ref b/tests/01_PW/093_PW_MSST/result.ref similarity index 100% rename from tests/01_PW/092_PW_MSST/result.ref rename to tests/01_PW/093_PW_MSST/result.ref diff --git a/tests/01_PW/093_PW_MSST2/INPUT b/tests/01_PW/094_PW_MSST2/INPUT similarity index 100% rename from tests/01_PW/093_PW_MSST2/INPUT rename to tests/01_PW/094_PW_MSST2/INPUT diff --git a/tests/01_PW/093_PW_MSST2/KPT b/tests/01_PW/094_PW_MSST2/KPT similarity index 100% rename from tests/01_PW/093_PW_MSST2/KPT rename to tests/01_PW/094_PW_MSST2/KPT diff --git a/tests/01_PW/093_PW_MSST2/README b/tests/01_PW/094_PW_MSST2/README similarity index 100% rename from tests/01_PW/093_PW_MSST2/README rename to tests/01_PW/094_PW_MSST2/README diff --git a/tests/01_PW/093_PW_MSST2/STRU b/tests/01_PW/094_PW_MSST2/STRU similarity index 100% rename from tests/01_PW/093_PW_MSST2/STRU rename to tests/01_PW/094_PW_MSST2/STRU diff --git a/tests/01_PW/093_PW_MSST2/result.ref b/tests/01_PW/094_PW_MSST2/result.ref similarity index 100% rename from tests/01_PW/093_PW_MSST2/result.ref rename to tests/01_PW/094_PW_MSST2/result.ref diff --git a/tests/01_PW/094_PW_NPT/INPUT b/tests/01_PW/095_PW_NPT/INPUT similarity index 100% rename from tests/01_PW/094_PW_NPT/INPUT rename to tests/01_PW/095_PW_NPT/INPUT diff --git a/tests/01_PW/094_PW_NPT/KPT b/tests/01_PW/095_PW_NPT/KPT similarity index 100% rename from tests/01_PW/094_PW_NPT/KPT rename to tests/01_PW/095_PW_NPT/KPT diff --git a/tests/01_PW/094_PW_NPT/README b/tests/01_PW/095_PW_NPT/README similarity index 100% rename from tests/01_PW/094_PW_NPT/README rename to tests/01_PW/095_PW_NPT/README diff --git a/tests/01_PW/094_PW_NPT/STRU b/tests/01_PW/095_PW_NPT/STRU similarity index 100% rename from tests/01_PW/094_PW_NPT/STRU rename to tests/01_PW/095_PW_NPT/STRU diff --git a/tests/01_PW/094_PW_NPT/result.ref b/tests/01_PW/095_PW_NPT/result.ref similarity index 100% rename from tests/01_PW/094_PW_NPT/result.ref rename to tests/01_PW/095_PW_NPT/result.ref diff --git a/tests/01_PW/095_PW_NVT/INPUT b/tests/01_PW/096_PW_NVT/INPUT similarity index 100% rename from tests/01_PW/095_PW_NVT/INPUT rename to tests/01_PW/096_PW_NVT/INPUT diff --git a/tests/01_PW/095_PW_NVT/KPT b/tests/01_PW/096_PW_NVT/KPT similarity index 100% rename from tests/01_PW/095_PW_NVT/KPT rename to tests/01_PW/096_PW_NVT/KPT diff --git a/tests/01_PW/095_PW_NVT/README b/tests/01_PW/096_PW_NVT/README similarity index 100% rename from tests/01_PW/095_PW_NVT/README rename to tests/01_PW/096_PW_NVT/README diff --git a/tests/01_PW/095_PW_NVT/Restart_md.txt b/tests/01_PW/096_PW_NVT/Restart_md.txt similarity index 100% rename from tests/01_PW/095_PW_NVT/Restart_md.txt rename to tests/01_PW/096_PW_NVT/Restart_md.txt diff --git a/tests/01_PW/095_PW_NVT/STRU_MD_2 b/tests/01_PW/096_PW_NVT/STRU_MD_2 similarity index 100% rename from tests/01_PW/095_PW_NVT/STRU_MD_2 rename to tests/01_PW/096_PW_NVT/STRU_MD_2 diff --git a/tests/01_PW/095_PW_NVT/result.ref b/tests/01_PW/096_PW_NVT/result.ref similarity index 100% rename from tests/01_PW/095_PW_NVT/result.ref rename to tests/01_PW/096_PW_NVT/result.ref diff --git a/tests/01_PW/096_PW_PBE0/INPUT b/tests/01_PW/097_PW_PBE0/INPUT similarity index 100% rename from tests/01_PW/096_PW_PBE0/INPUT rename to tests/01_PW/097_PW_PBE0/INPUT diff --git a/tests/01_PW/096_PW_PBE0/KPT b/tests/01_PW/097_PW_PBE0/KPT similarity index 100% rename from tests/01_PW/096_PW_PBE0/KPT rename to tests/01_PW/097_PW_PBE0/KPT diff --git a/tests/01_PW/096_PW_PBE0/README b/tests/01_PW/097_PW_PBE0/README similarity index 100% rename from tests/01_PW/096_PW_PBE0/README rename to tests/01_PW/097_PW_PBE0/README diff --git a/tests/01_PW/096_PW_PBE0/STRU b/tests/01_PW/097_PW_PBE0/STRU similarity index 100% rename from tests/01_PW/096_PW_PBE0/STRU rename to tests/01_PW/097_PW_PBE0/STRU diff --git a/tests/01_PW/096_PW_PBE0/result.ref b/tests/01_PW/097_PW_PBE0/result.ref similarity index 100% rename from tests/01_PW/096_PW_PBE0/result.ref rename to tests/01_PW/097_PW_PBE0/result.ref diff --git a/tests/01_PW/096_PW_PBE0_AFM/INPUT b/tests/01_PW/097_PW_PBE0_AFM/INPUT similarity index 100% rename from tests/01_PW/096_PW_PBE0_AFM/INPUT rename to tests/01_PW/097_PW_PBE0_AFM/INPUT diff --git a/tests/01_PW/096_PW_PBE0_AFM/KPT b/tests/01_PW/097_PW_PBE0_AFM/KPT similarity index 100% rename from tests/01_PW/096_PW_PBE0_AFM/KPT rename to tests/01_PW/097_PW_PBE0_AFM/KPT diff --git a/tests/01_PW/096_PW_PBE0_AFM/README b/tests/01_PW/097_PW_PBE0_AFM/README similarity index 100% rename from tests/01_PW/096_PW_PBE0_AFM/README rename to tests/01_PW/097_PW_PBE0_AFM/README diff --git a/tests/01_PW/096_PW_PBE0_AFM/STRU b/tests/01_PW/097_PW_PBE0_AFM/STRU similarity index 100% rename from tests/01_PW/096_PW_PBE0_AFM/STRU rename to tests/01_PW/097_PW_PBE0_AFM/STRU diff --git a/tests/01_PW/096_PW_PBE0_AFM/result.ref b/tests/01_PW/097_PW_PBE0_AFM/result.ref similarity index 100% rename from tests/01_PW/096_PW_PBE0_AFM/result.ref rename to tests/01_PW/097_PW_PBE0_AFM/result.ref diff --git a/tests/01_PW/096_PW_PBE0_AFM/threshold b/tests/01_PW/097_PW_PBE0_AFM/threshold similarity index 100% rename from tests/01_PW/096_PW_PBE0_AFM/threshold rename to tests/01_PW/097_PW_PBE0_AFM/threshold diff --git a/tests/01_PW/096_PW_PBE0_FM/INPUT b/tests/01_PW/097_PW_PBE0_FM/INPUT similarity index 100% rename from tests/01_PW/096_PW_PBE0_FM/INPUT rename to tests/01_PW/097_PW_PBE0_FM/INPUT diff --git a/tests/01_PW/096_PW_PBE0_FM/KPT b/tests/01_PW/097_PW_PBE0_FM/KPT similarity index 100% rename from tests/01_PW/096_PW_PBE0_FM/KPT rename to tests/01_PW/097_PW_PBE0_FM/KPT diff --git a/tests/01_PW/096_PW_PBE0_FM/README b/tests/01_PW/097_PW_PBE0_FM/README similarity index 100% rename from tests/01_PW/096_PW_PBE0_FM/README rename to tests/01_PW/097_PW_PBE0_FM/README diff --git a/tests/01_PW/096_PW_PBE0_FM/STRU b/tests/01_PW/097_PW_PBE0_FM/STRU similarity index 100% rename from tests/01_PW/096_PW_PBE0_FM/STRU rename to tests/01_PW/097_PW_PBE0_FM/STRU diff --git a/tests/01_PW/096_PW_PBE0_FM/result.ref b/tests/01_PW/097_PW_PBE0_FM/result.ref similarity index 100% rename from tests/01_PW/096_PW_PBE0_FM/result.ref rename to tests/01_PW/097_PW_PBE0_FM/result.ref diff --git a/tests/01_PW/096_PW_PBE0_FM/threshold b/tests/01_PW/097_PW_PBE0_FM/threshold similarity index 100% rename from tests/01_PW/096_PW_PBE0_FM/threshold rename to tests/01_PW/097_PW_PBE0_FM/threshold diff --git a/tests/01_PW/098_PW_15_SO_avg/INPUT b/tests/01_PW/099_PW_15_SO_avg/INPUT similarity index 100% rename from tests/01_PW/098_PW_15_SO_avg/INPUT rename to tests/01_PW/099_PW_15_SO_avg/INPUT diff --git a/tests/01_PW/098_PW_15_SO_avg/KPT b/tests/01_PW/099_PW_15_SO_avg/KPT similarity index 100% rename from tests/01_PW/098_PW_15_SO_avg/KPT rename to tests/01_PW/099_PW_15_SO_avg/KPT diff --git a/tests/01_PW/098_PW_15_SO_avg/README b/tests/01_PW/099_PW_15_SO_avg/README similarity index 100% rename from tests/01_PW/098_PW_15_SO_avg/README rename to tests/01_PW/099_PW_15_SO_avg/README diff --git a/tests/01_PW/098_PW_15_SO_avg/STRU b/tests/01_PW/099_PW_15_SO_avg/STRU similarity index 100% rename from tests/01_PW/098_PW_15_SO_avg/STRU rename to tests/01_PW/099_PW_15_SO_avg/STRU diff --git a/tests/01_PW/098_PW_15_SO_avg/result.ref b/tests/01_PW/099_PW_15_SO_avg/result.ref similarity index 100% rename from tests/01_PW/098_PW_15_SO_avg/result.ref rename to tests/01_PW/099_PW_15_SO_avg/result.ref diff --git a/tests/01_PW/099_PW_DJ_SO/INPUT b/tests/01_PW/100_PW_DJ_SO/INPUT similarity index 100% rename from tests/01_PW/099_PW_DJ_SO/INPUT rename to tests/01_PW/100_PW_DJ_SO/INPUT diff --git a/tests/01_PW/099_PW_DJ_SO/KPT b/tests/01_PW/100_PW_DJ_SO/KPT similarity index 100% rename from tests/01_PW/099_PW_DJ_SO/KPT rename to tests/01_PW/100_PW_DJ_SO/KPT diff --git a/tests/01_PW/099_PW_DJ_SO/README b/tests/01_PW/100_PW_DJ_SO/README similarity index 100% rename from tests/01_PW/099_PW_DJ_SO/README rename to tests/01_PW/100_PW_DJ_SO/README diff --git a/tests/01_PW/099_PW_DJ_SO/STRU b/tests/01_PW/100_PW_DJ_SO/STRU similarity index 100% rename from tests/01_PW/099_PW_DJ_SO/STRU rename to tests/01_PW/100_PW_DJ_SO/STRU diff --git a/tests/01_PW/099_PW_DJ_SO/result.ref b/tests/01_PW/100_PW_DJ_SO/result.ref similarity index 100% rename from tests/01_PW/099_PW_DJ_SO/result.ref rename to tests/01_PW/100_PW_DJ_SO/result.ref diff --git a/tests/01_PW/100_PW_W90/INPUT b/tests/01_PW/101_PW_W90/INPUT similarity index 100% rename from tests/01_PW/100_PW_W90/INPUT rename to tests/01_PW/101_PW_W90/INPUT diff --git a/tests/01_PW/100_PW_W90/KPT b/tests/01_PW/101_PW_W90/KPT similarity index 100% rename from tests/01_PW/100_PW_W90/KPT rename to tests/01_PW/101_PW_W90/KPT diff --git a/tests/01_PW/100_PW_W90/STRU b/tests/01_PW/101_PW_W90/STRU similarity index 100% rename from tests/01_PW/100_PW_W90/STRU rename to tests/01_PW/101_PW_W90/STRU diff --git a/tests/01_PW/100_PW_W90/chg.cube b/tests/01_PW/101_PW_W90/chg.cube similarity index 100% rename from tests/01_PW/100_PW_W90/chg.cube rename to tests/01_PW/101_PW_W90/chg.cube diff --git a/tests/01_PW/100_PW_W90/diamond.amn b/tests/01_PW/101_PW_W90/diamond.amn similarity index 100% rename from tests/01_PW/100_PW_W90/diamond.amn rename to tests/01_PW/101_PW_W90/diamond.amn diff --git a/tests/01_PW/100_PW_W90/diamond.eig b/tests/01_PW/101_PW_W90/diamond.eig similarity index 100% rename from tests/01_PW/100_PW_W90/diamond.eig rename to tests/01_PW/101_PW_W90/diamond.eig diff --git a/tests/01_PW/100_PW_W90/diamond.mmn b/tests/01_PW/101_PW_W90/diamond.mmn similarity index 100% rename from tests/01_PW/100_PW_W90/diamond.mmn rename to tests/01_PW/101_PW_W90/diamond.mmn diff --git a/tests/01_PW/100_PW_W90/diamond.nnkp b/tests/01_PW/101_PW_W90/diamond.nnkp similarity index 100% rename from tests/01_PW/100_PW_W90/diamond.nnkp rename to tests/01_PW/101_PW_W90/diamond.nnkp diff --git a/tests/01_PW/100_PW_W90/result.ref b/tests/01_PW/101_PW_W90/result.ref similarity index 100% rename from tests/01_PW/100_PW_W90/result.ref rename to tests/01_PW/101_PW_W90/result.ref diff --git a/tests/01_PW/101_PW_MD_1O/INPUT b/tests/01_PW/102_PW_MD_1O/INPUT similarity index 100% rename from tests/01_PW/101_PW_MD_1O/INPUT rename to tests/01_PW/102_PW_MD_1O/INPUT diff --git a/tests/01_PW/101_PW_MD_1O/KPT b/tests/01_PW/102_PW_MD_1O/KPT similarity index 100% rename from tests/01_PW/101_PW_MD_1O/KPT rename to tests/01_PW/102_PW_MD_1O/KPT diff --git a/tests/01_PW/101_PW_MD_1O/README b/tests/01_PW/102_PW_MD_1O/README similarity index 100% rename from tests/01_PW/101_PW_MD_1O/README rename to tests/01_PW/102_PW_MD_1O/README diff --git a/tests/01_PW/101_PW_MD_1O/STRU b/tests/01_PW/102_PW_MD_1O/STRU similarity index 100% rename from tests/01_PW/101_PW_MD_1O/STRU rename to tests/01_PW/102_PW_MD_1O/STRU diff --git a/tests/01_PW/101_PW_MD_1O/result.ref b/tests/01_PW/102_PW_MD_1O/result.ref similarity index 100% rename from tests/01_PW/101_PW_MD_1O/result.ref rename to tests/01_PW/102_PW_MD_1O/result.ref diff --git a/tests/01_PW/102_PW_MD_2O/INPUT b/tests/01_PW/103_PW_MD_2O/INPUT similarity index 100% rename from tests/01_PW/102_PW_MD_2O/INPUT rename to tests/01_PW/103_PW_MD_2O/INPUT diff --git a/tests/01_PW/102_PW_MD_2O/KPT b/tests/01_PW/103_PW_MD_2O/KPT similarity index 100% rename from tests/01_PW/102_PW_MD_2O/KPT rename to tests/01_PW/103_PW_MD_2O/KPT diff --git a/tests/01_PW/102_PW_MD_2O/README b/tests/01_PW/103_PW_MD_2O/README similarity index 100% rename from tests/01_PW/102_PW_MD_2O/README rename to tests/01_PW/103_PW_MD_2O/README diff --git a/tests/01_PW/102_PW_MD_2O/STRU b/tests/01_PW/103_PW_MD_2O/STRU similarity index 100% rename from tests/01_PW/102_PW_MD_2O/STRU rename to tests/01_PW/103_PW_MD_2O/STRU diff --git a/tests/01_PW/102_PW_MD_2O/result.ref b/tests/01_PW/103_PW_MD_2O/result.ref similarity index 100% rename from tests/01_PW/102_PW_MD_2O/result.ref rename to tests/01_PW/103_PW_MD_2O/result.ref diff --git a/tests/01_PW/103_PW_Gene_Descriptors/INPUT b/tests/01_PW/104_PW_Gene_Descriptors/INPUT similarity index 100% rename from tests/01_PW/103_PW_Gene_Descriptors/INPUT rename to tests/01_PW/104_PW_Gene_Descriptors/INPUT diff --git a/tests/01_PW/103_PW_Gene_Descriptors/KPT b/tests/01_PW/104_PW_Gene_Descriptors/KPT similarity index 100% rename from tests/01_PW/103_PW_Gene_Descriptors/KPT rename to tests/01_PW/104_PW_Gene_Descriptors/KPT diff --git a/tests/01_PW/103_PW_Gene_Descriptors/README b/tests/01_PW/104_PW_Gene_Descriptors/README similarity index 100% rename from tests/01_PW/103_PW_Gene_Descriptors/README rename to tests/01_PW/104_PW_Gene_Descriptors/README diff --git a/tests/01_PW/103_PW_Gene_Descriptors/STRU b/tests/01_PW/104_PW_Gene_Descriptors/STRU similarity index 100% rename from tests/01_PW/103_PW_Gene_Descriptors/STRU rename to tests/01_PW/104_PW_Gene_Descriptors/STRU diff --git a/tests/01_PW/103_PW_Gene_Descriptors/result.ref b/tests/01_PW/104_PW_Gene_Descriptors/result.ref similarity index 100% rename from tests/01_PW/103_PW_Gene_Descriptors/result.ref rename to tests/01_PW/104_PW_Gene_Descriptors/result.ref diff --git a/tests/01_PW/CASES_CPU.txt b/tests/01_PW/CASES_CPU.txt index c30cac71f3a..ee1a83c0c8d 100644 --- a/tests/01_PW/CASES_CPU.txt +++ b/tests/01_PW/CASES_CPU.txt @@ -89,21 +89,22 @@ scf_out_chg_tau 087_PW_get_pchg_kpar 088_PW_get_pchg_sepk 089_PW_get_wf_kpar -090_PW_VWR -091_PW_CR_VDW3 -092_PW_MSST -093_PW_MSST2 -094_PW_NPT -095_PW_NVT -096_PW_PBE0 -096_PW_PBE0_AFM -096_PW_PBE0_FM -098_PW_15_SO_avg -099_PW_DJ_SO -100_PW_W90 -101_PW_MD_1O -102_PW_MD_2O -103_PW_Gene_Descriptors +090_PW_out_pchg_wfc_spinor +091_PW_VWR +092_PW_CR_VDW3 +093_PW_MSST +094_PW_MSST2 +095_PW_NPT +096_PW_NVT +097_PW_PBE0 +097_PW_PBE0_AFM +097_PW_PBE0_FM +099_PW_15_SO_avg +100_PW_DJ_SO +101_PW_W90 +102_PW_MD_1O +103_PW_MD_2O +104_PW_Gene_Descriptors 201_PW_UPF201_Ce_f 202_PW_ONCV_Libxc 204_PW_SY diff --git a/tests/01_PW/CASES_GPU.txt b/tests/01_PW/CASES_GPU.txt index 20f6e94467e..b859eb8f8d3 100644 --- a/tests/01_PW/CASES_GPU.txt +++ b/tests/01_PW/CASES_GPU.txt @@ -89,21 +89,22 @@ scf_out_elf #087_PW_get_pchg_kpar #088_PW_get_pchg_sepk #089_PW_get_wf_kpar -090_PW_VWR -091_PW_CR_VDW3 -#092_PW_MSST -#093_PW_MSST2 -094_PW_NPT -#095_PW_NVT -#096_PW_PBE0 -#096_PW_PBE0_AFM -#096_PW_PBE0_FM -098_PW_15_SO_avg -#099_PW_DJ_SO -#100_PW_W90 -101_PW_MD_1O -102_PW_MD_2O -#103_PW_Gene_Descriptors +#090_PW_out_pchg_wfc_spinor +091_PW_VWR +092_PW_CR_VDW3 +#093_PW_MSST +#094_PW_MSST2 +095_PW_NPT +#096_PW_NVT +#097_PW_PBE0 +#097_PW_PBE0_AFM +#097_PW_PBE0_FM +099_PW_15_SO_avg +#100_PW_DJ_SO +#101_PW_W90 +102_PW_MD_1O +103_PW_MD_2O +#104_PW_Gene_Descriptors #201_PW_UPF201_Ce_f #202_PW_ONCV_Libxc #204_PW_SY diff --git a/tests/02_NAO_Gamma/get_wf_spin2/result.ref b/tests/02_NAO_Gamma/get_wf_spin2/result.ref index ef88fc35f1b..f3e66387a24 100644 --- a/tests/02_NAO_Gamma/get_wf_spin2/result.ref +++ b/tests/02_NAO_Gamma/get_wf_spin2/result.ref @@ -1,7 +1,63 @@ wfi1s1.cube -11.7665927 wfi1s2.cube -11.76659269 -wfi2s1im.cube 0 -wfi2s1re.cube -11.76659269 -wfi2s2im.cube 0 -wfi2s2re.cube 3.042711753e-14 -totaltimeref 0.36 +wfi2s1_wfc_fp_components 1 +wfi2s1_wfc_fp_nx 32 +wfi2s1_wfc_fp_ny 32 +wfi2s1_wfc_fp_nz 32 +wfi2s1_wfc_fp_voxel 3.051757812500e-02 +wfi2s1_wfc_fp_rms 4.469486123737e-02 +wfi2s1_wfc_fp_power_0 2.957659892862e+00 +wfi2s1_wfc_fp_power_1 7.905249614220e-01 +wfi2s1_wfc_fp_power_2 1.206048372416e-01 +wfi2s1_wfc_fp_power_3 9.539590384810e-02 +wfi2s1_wfc_fp_power_4 3.333875429220e-01 +wfi2s1_wfc_fp_power_5 9.053006739075e-01 +wfi2s1_wfc_fp_power_6 1.840172866335e+00 +wfi2s1_wfc_fp_power_7 1.750156213015e+00 +wfi2s1_wfc_fp_cross_0_re -1.153064324734e+00 +wfi2s1_wfc_fp_cross_0_im -1.004264226054e+00 +wfi2s1_wfc_fp_cross_1_re -2.546670396469e-01 +wfi2s1_wfc_fp_cross_1_im -1.746019278968e-01 +wfi2s1_wfc_fp_cross_2_re -2.838664706167e-02 +wfi2s1_wfc_fp_cross_2_im 1.034379317548e-01 +wfi2s1_wfc_fp_cross_3_re -1.736168623413e-01 +wfi2s1_wfc_fp_cross_3_im -4.075525855058e-02 +wfi2s1_wfc_fp_cross_4_re 3.603487181729e-01 +wfi2s1_wfc_fp_cross_4_im 4.146863496558e-01 +wfi2s1_wfc_fp_cross_5_re 5.250645449420e-01 +wfi2s1_wfc_fp_cross_5_im 1.179074620049e+00 +wfi2s1_wfc_fp_cross_6_re 7.227505542364e-01 +wfi2s1_wfc_fp_cross_6_im 1.642626436956e+00 +wfi2s1_wfc_fp_cross_7_re -2.194499194887e+00 +wfi2s1_wfc_fp_cross_7_im 6.004499322329e-01 +wfi2s2_wfc_fp_components 1 +wfi2s2_wfc_fp_nx 32 +wfi2s2_wfc_fp_ny 32 +wfi2s2_wfc_fp_nz 32 +wfi2s2_wfc_fp_voxel 3.051757812500e-02 +wfi2s2_wfc_fp_rms 6.318141313616e-02 +wfi2s2_wfc_fp_power_0 3.867728978007e+00 +wfi2s2_wfc_fp_power_1 1.155630552870e-01 +wfi2s2_wfc_fp_power_2 1.135675280751e+00 +wfi2s2_wfc_fp_power_3 4.764999384103e-02 +wfi2s2_wfc_fp_power_4 1.141080657160e-02 +wfi2s2_wfc_fp_power_5 1.713477244289e-01 +wfi2s2_wfc_fp_power_6 5.860299849474e-01 +wfi2s2_wfc_fp_power_7 5.336133513696e-01 +wfi2s2_wfc_fp_cross_0_re -5.234515002146e-01 +wfi2s2_wfc_fp_cross_0_im -4.158907364245e-01 +wfi2s2_wfc_fp_cross_1_re 8.251748522092e-02 +wfi2s2_wfc_fp_cross_1_im -3.527505774487e-01 +wfi2s2_wfc_fp_cross_2_re -2.067650317031e-01 +wfi2s2_wfc_fp_cross_2_im 1.065980384342e-01 +wfi2s2_wfc_fp_cross_3_re 2.085528361228e-02 +wfi2s2_wfc_fp_cross_3_im 1.042986137535e-02 +wfi2s2_wfc_fp_cross_4_re 1.038410493530e-02 +wfi2s2_wfc_fp_cross_4_im 4.298122967802e-02 +wfi2s2_wfc_fp_cross_5_re 1.454718504451e-01 +wfi2s2_wfc_fp_cross_5_im 2.815188183691e-01 +wfi2s2_wfc_fp_cross_6_re -9.786135454397e-02 +wfi2s2_wfc_fp_cross_6_im 5.505784045508e-01 +wfi2s2_wfc_fp_cross_7_re -1.349122341727e+00 +wfi2s2_wfc_fp_cross_7_im 4.937010524570e-01 +totaltimeref 0.39 diff --git a/tests/03_NAO_multik/get_wf/result.ref b/tests/03_NAO_multik/get_wf/result.ref index 6c0d25977fd..d07a329fc78 100644 --- a/tests/03_NAO_multik/get_wf/result.ref +++ b/tests/03_NAO_multik/get_wf/result.ref @@ -1,19 +1,187 @@ wfi1s1k1.cube -11.77152988 -wfi1s1k1im.cube 0 -wfi1s1k1re.cube -11.77169544 wfi1s1k2.cube -3.753605249 -wfi1s1k2im.cube -6.240965225e-15 -wfi1s1k2re.cube -3.753696262 wfi1s1k3.cube -3.753605249 -wfi1s1k3im.cube 1.793009931e-15 -wfi1s1k3re.cube -3.753696262 wfi1s2k1.cube -11.77152988 -wfi1s2k1im.cube 0 -wfi1s2k1re.cube -15.5253917 wfi1s2k2.cube -3.753605249 -wfi1s2k2im.cube 3.852845828e-14 -wfi1s2k2re.cube -7.507392525 wfi1s2k3.cube -3.753605249 -wfi1s2k3im.cube -2.168256114e-14 -wfi1s2k3re.cube -7.507392525 -totaltimeref 0.46 +wfi1s1k1_wfc_fp_components 1 +wfi1s1k1_wfc_fp_nx 32 +wfi1s1k1_wfc_fp_ny 32 +wfi1s1k1_wfc_fp_nz 32 +wfi1s1k1_wfc_fp_voxel 3.051757812500e-02 +wfi1s1k1_wfc_fp_rms 3.161737034233e-02 +wfi1s1k1_wfc_fp_power_0 6.037105428880e-01 +wfi1s1k1_wfc_fp_power_1 8.656860728773e-01 +wfi1s1k1_wfc_fp_power_2 1.239410044490e+00 +wfi1s1k1_wfc_fp_power_3 3.599593248387e-01 +wfi1s1k1_wfc_fp_power_4 5.379684311148e-01 +wfi1s1k1_wfc_fp_power_5 8.823644979484e-01 +wfi1s1k1_wfc_fp_power_6 1.400547108329e+00 +wfi1s1k1_wfc_fp_power_7 1.484139852679e+00 +wfi1s1k1_wfc_fp_cross_0_re -1.073532352801e-01 +wfi1s1k1_wfc_fp_cross_0_im -7.149119469573e-01 +wfi1s1k1_wfc_fp_cross_1_re -8.771577497186e-01 +wfi1s1k1_wfc_fp_cross_1_im 5.509394669176e-01 +wfi1s1k1_wfc_fp_cross_2_re -6.304007682579e-01 +wfi1s1k1_wfc_fp_cross_2_im 2.207534239659e-01 +wfi1s1k1_wfc_fp_cross_3_re -3.629257360684e-01 +wfi1s1k1_wfc_fp_cross_3_im -2.488607308272e-01 +wfi1s1k1_wfc_fp_cross_4_re 4.602064046585e-01 +wfi1s1k1_wfc_fp_cross_4_im 5.127322008066e-01 +wfi1s1k1_wfc_fp_cross_5_re 3.796007637900e-01 +wfi1s1k1_wfc_fp_cross_5_im 1.044842718415e+00 +wfi1s1k1_wfc_fp_cross_6_re 1.008054980940e+00 +wfi1s1k1_wfc_fp_cross_6_im 1.030743874311e+00 +wfi1s1k1_wfc_fp_cross_7_re -6.305721235833e-01 +wfi1s1k1_wfc_fp_cross_7_im 7.059530247421e-01 +wfi1s1k2_wfc_fp_components 1 +wfi1s1k2_wfc_fp_nx 32 +wfi1s1k2_wfc_fp_ny 32 +wfi1s1k2_wfc_fp_nz 32 +wfi1s1k2_wfc_fp_voxel 3.051757812500e-02 +wfi1s1k2_wfc_fp_rms 2.496130010879e-02 +wfi1s1k2_wfc_fp_power_0 3.012892846934e+00 +wfi1s1k2_wfc_fp_power_1 1.426989648828e+00 +wfi1s1k2_wfc_fp_power_2 5.013489543082e-01 +wfi1s1k2_wfc_fp_power_3 1.524297186710e-01 +wfi1s1k2_wfc_fp_power_4 2.958521226469e+00 +wfi1s1k2_wfc_fp_power_5 5.369474479834e-01 +wfi1s1k2_wfc_fp_power_6 1.243916024098e-01 +wfi1s1k2_wfc_fp_power_7 1.001604632053e-01 +wfi1s1k2_wfc_fp_cross_0_re -1.844628302122e+00 +wfi1s1k2_wfc_fp_cross_0_im -9.469494878885e-01 +wfi1s1k2_wfc_fp_cross_1_re 1.496929111989e-01 +wfi1s1k2_wfc_fp_cross_1_im 8.324733032268e-01 +wfi1s1k2_wfc_fp_cross_2_re -5.825467221961e-02 +wfi1s1k2_wfc_fp_cross_2_im 2.702348482817e-01 +wfi1s1k2_wfc_fp_cross_3_re 6.139181616298e-01 +wfi1s1k2_wfc_fp_cross_3_im 2.721599696025e-01 +wfi1s1k2_wfc_fp_cross_4_re 1.252845461574e+00 +wfi1s1k2_wfc_fp_cross_4_im 1.376541745441e-01 +wfi1s1k2_wfc_fp_cross_5_re -2.042735083501e-01 +wfi1s1k2_wfc_fp_cross_5_im 1.583164149760e-01 +wfi1s1k2_wfc_fp_cross_6_re 6.698104084379e-02 +wfi1s1k2_wfc_fp_cross_6_im 8.928975688005e-02 +wfi1s1k2_wfc_fp_cross_7_re -1.452159174678e-01 +wfi1s1k2_wfc_fp_cross_7_im 5.297972069112e-01 +wfi1s1k3_wfc_fp_components 1 +wfi1s1k3_wfc_fp_nx 32 +wfi1s1k3_wfc_fp_ny 32 +wfi1s1k3_wfc_fp_nz 32 +wfi1s1k3_wfc_fp_voxel 3.051757812500e-02 +wfi1s1k3_wfc_fp_rms 2.495642827501e-02 +wfi1s1k3_wfc_fp_power_0 2.839961928444e+00 +wfi1s1k3_wfc_fp_power_1 1.426303371101e+00 +wfi1s1k3_wfc_fp_power_2 4.951108023002e-01 +wfi1s1k3_wfc_fp_power_3 1.713208293961e-01 +wfi1s1k3_wfc_fp_power_4 2.922711053847e+00 +wfi1s1k3_wfc_fp_power_5 6.332953032800e-01 +wfi1s1k3_wfc_fp_power_6 1.168508755464e-01 +wfi1s1k3_wfc_fp_power_7 1.224522633109e-01 +wfi1s1k3_wfc_fp_cross_0_re -1.769264122843e+00 +wfi1s1k3_wfc_fp_cross_0_im -9.593496421844e-01 +wfi1s1k3_wfc_fp_cross_1_re 1.431878759246e-01 +wfi1s1k3_wfc_fp_cross_1_im 8.280552146914e-01 +wfi1s1k3_wfc_fp_cross_2_re -2.650778645719e-02 +wfi1s1k3_wfc_fp_cross_2_im 2.900347057684e-01 +wfi1s1k3_wfc_fp_cross_3_re 6.291111804591e-01 +wfi1s1k3_wfc_fp_cross_3_im 3.239450639409e-01 +wfi1s1k3_wfc_fp_cross_4_re 1.348413095347e+00 +wfi1s1k3_wfc_fp_cross_4_im 1.808903190944e-01 +wfi1s1k3_wfc_fp_cross_5_re -2.082481263969e-01 +wfi1s1k3_wfc_fp_cross_5_im 1.750252225248e-01 +wfi1s1k3_wfc_fp_cross_6_re 6.226494485112e-02 +wfi1s1k3_wfc_fp_cross_6_im 1.021358449478e-01 +wfi1s1k3_wfc_fp_cross_7_re -1.291285378447e-01 +wfi1s1k3_wfc_fp_cross_7_im 5.754003706715e-01 +wfi1s2k1_wfc_fp_components 1 +wfi1s2k1_wfc_fp_nx 32 +wfi1s2k1_wfc_fp_ny 32 +wfi1s2k1_wfc_fp_nz 32 +wfi1s2k1_wfc_fp_voxel 3.051757812500e-02 +wfi1s2k1_wfc_fp_rms 4.880517600663e-02 +wfi1s2k1_wfc_fp_power_0 1.910695038838e+00 +wfi1s2k1_wfc_fp_power_1 1.103806840098e+00 +wfi1s2k1_wfc_fp_power_2 3.549508359597e-01 +wfi1s2k1_wfc_fp_power_3 3.236357619838e-02 +wfi1s2k1_wfc_fp_power_4 1.895457785219e+00 +wfi1s2k1_wfc_fp_power_5 7.214970541804e-01 +wfi1s2k1_wfc_fp_power_6 9.144806092504e-01 +wfi1s2k1_wfc_fp_power_7 1.025574948373e+00 +wfi1s2k1_wfc_fp_cross_0_re -8.143236483523e-01 +wfi1s2k1_wfc_fp_cross_0_im -1.202462161129e+00 +wfi1s2k1_wfc_fp_cross_1_re -5.297719562035e-01 +wfi1s2k1_wfc_fp_cross_1_im 3.333749166497e-01 +wfi1s2k1_wfc_fp_cross_2_re -1.065973429297e-01 +wfi1s2k1_wfc_fp_cross_2_im -1.115728042995e-02 +wfi1s2k1_wfc_fp_cross_3_re -2.098559975955e-01 +wfi1s2k1_wfc_fp_cross_3_im -1.315456298626e-01 +wfi1s2k1_wfc_fp_cross_4_re 9.959578342628e-01 +wfi1s2k1_wfc_fp_cross_4_im 6.128908554789e-01 +wfi1s2k1_wfc_fp_cross_5_re -3.924305515258e-02 +wfi1s2k1_wfc_fp_cross_5_im 8.113291861517e-01 +wfi1s2k1_wfc_fp_cross_6_re 6.288135404094e-01 +wfi1s2k1_wfc_fp_cross_6_im 7.365201524862e-01 +wfi1s2k1_wfc_fp_cross_7_re -7.897470030735e-01 +wfi1s2k1_wfc_fp_cross_7_im 1.155794374857e+00 +wfi1s2k2_wfc_fp_components 1 +wfi1s2k2_wfc_fp_nx 32 +wfi1s2k2_wfc_fp_ny 32 +wfi1s2k2_wfc_fp_nz 32 +wfi1s2k2_wfc_fp_voxel 3.051757812500e-02 +wfi1s2k2_wfc_fp_rms 4.992260021758e-02 +wfi1s2k2_wfc_fp_power_0 3.012892846933e+00 +wfi1s2k2_wfc_fp_power_1 1.426989648820e+00 +wfi1s2k2_wfc_fp_power_2 5.013489543058e-01 +wfi1s2k2_wfc_fp_power_3 1.524297186715e-01 +wfi1s2k2_wfc_fp_power_4 2.958521226468e+00 +wfi1s2k2_wfc_fp_power_5 5.369474479841e-01 +wfi1s2k2_wfc_fp_power_6 1.243916024103e-01 +wfi1s2k2_wfc_fp_power_7 1.001604632053e-01 +wfi1s2k2_wfc_fp_cross_0_re -1.844628302119e+00 +wfi1s2k2_wfc_fp_cross_0_im -9.469494878813e-01 +wfi1s2k2_wfc_fp_cross_1_re 1.496929111974e-01 +wfi1s2k2_wfc_fp_cross_1_im 8.324733032226e-01 +wfi1s2k2_wfc_fp_cross_2_re -5.825467221896e-02 +wfi1s2k2_wfc_fp_cross_2_im 2.702348482816e-01 +wfi1s2k2_wfc_fp_cross_3_re 6.139181616301e-01 +wfi1s2k2_wfc_fp_cross_3_im 2.721599696043e-01 +wfi1s2k2_wfc_fp_cross_4_re 1.252845461575e+00 +wfi1s2k2_wfc_fp_cross_4_im 1.376541745464e-01 +wfi1s2k2_wfc_fp_cross_5_re -2.042735083496e-01 +wfi1s2k2_wfc_fp_cross_5_im 1.583164149777e-01 +wfi1s2k2_wfc_fp_cross_6_re 6.698104084396e-02 +wfi1s2k2_wfc_fp_cross_6_im 8.928975688022e-02 +wfi1s2k2_wfc_fp_cross_7_re -1.452159174712e-01 +wfi1s2k2_wfc_fp_cross_7_im 5.297972069102e-01 +wfi1s2k3_wfc_fp_components 1 +wfi1s2k3_wfc_fp_nx 32 +wfi1s2k3_wfc_fp_ny 32 +wfi1s2k3_wfc_fp_nz 32 +wfi1s2k3_wfc_fp_voxel 3.051757812500e-02 +wfi1s2k3_wfc_fp_rms 4.991285655002e-02 +wfi1s2k3_wfc_fp_power_0 2.839961928443e+00 +wfi1s2k3_wfc_fp_power_1 1.426303371104e+00 +wfi1s2k3_wfc_fp_power_2 4.951108022983e-01 +wfi1s2k3_wfc_fp_power_3 1.713208293970e-01 +wfi1s2k3_wfc_fp_power_4 2.922711053839e+00 +wfi1s2k3_wfc_fp_power_5 6.332953032791e-01 +wfi1s2k3_wfc_fp_power_6 1.168508755451e-01 +wfi1s2k3_wfc_fp_power_7 1.224522633118e-01 +wfi1s2k3_wfc_fp_cross_0_re -1.769264122845e+00 +wfi1s2k3_wfc_fp_cross_0_im -9.593496421843e-01 +wfi1s2k3_wfc_fp_cross_1_re 1.431878759235e-01 +wfi1s2k3_wfc_fp_cross_1_im 8.280552146908e-01 +wfi1s2k3_wfc_fp_cross_2_re -2.650778645786e-02 +wfi1s2k3_wfc_fp_cross_2_im 2.900347057686e-01 +wfi1s2k3_wfc_fp_cross_3_re 6.291111804610e-01 +wfi1s2k3_wfc_fp_cross_3_im 3.239450639389e-01 +wfi1s2k3_wfc_fp_cross_4_re 1.348413095344e+00 +wfi1s2k3_wfc_fp_cross_4_im 1.808903190981e-01 +wfi1s2k3_wfc_fp_cross_5_re -2.082481263968e-01 +wfi1s2k3_wfc_fp_cross_5_im 1.750252225223e-01 +wfi1s2k3_wfc_fp_cross_6_re 6.226494485243e-02 +wfi1s2k3_wfc_fp_cross_6_im 1.021358449468e-01 +wfi1s2k3_wfc_fp_cross_7_re -1.291285378481e-01 +wfi1s2k3_wfc_fp_cross_7_im 5.754003706729e-01 +totaltimeref 0.73 diff --git a/tests/03_NAO_multik/get_wf0/result.ref b/tests/03_NAO_multik/get_wf0/result.ref index e4d060972bf..766de44d5de 100644 --- a/tests/03_NAO_multik/get_wf0/result.ref +++ b/tests/03_NAO_multik/get_wf0/result.ref @@ -1,19 +1,187 @@ wfi1s1k1.cube -11.77152988 -wfi1s1k1im.cube 0 -wfi1s1k1re.cube -11.77169544 wfi1s1k2.cube -3.754629032 -wfi1s1k2im.cube -5.256421773e-15 -wfi1s1k2re.cube -3.754662163 wfi1s1k3.cube -3.754629032 -wfi1s1k3im.cube -2.592407863e-14 -wfi1s1k3re.cube -3.754662163 wfi1s2k1.cube -11.77152988 -wfi1s2k1im.cube 0 -wfi1s2k1re.cube -15.5263576 wfi1s2k2.cube -3.754629032 -wfi1s2k2im.cube 8.143533574e-14 -wfi1s2k2re.cube -7.509324327 wfi1s2k3.cube -3.754629032 -wfi1s2k3im.cube -3.563679367e-14 -wfi1s2k3re.cube -7.509324327 -totaltimeref 0.22 +wfi1s1k1_wfc_fp_components 1 +wfi1s1k1_wfc_fp_nx 32 +wfi1s1k1_wfc_fp_ny 32 +wfi1s1k1_wfc_fp_nz 32 +wfi1s1k1_wfc_fp_voxel 3.051757812500e-02 +wfi1s1k1_wfc_fp_rms 3.161737034233e-02 +wfi1s1k1_wfc_fp_power_0 6.037105428880e-01 +wfi1s1k1_wfc_fp_power_1 8.656860728773e-01 +wfi1s1k1_wfc_fp_power_2 1.239410044490e+00 +wfi1s1k1_wfc_fp_power_3 3.599593248387e-01 +wfi1s1k1_wfc_fp_power_4 5.379684311148e-01 +wfi1s1k1_wfc_fp_power_5 8.823644979484e-01 +wfi1s1k1_wfc_fp_power_6 1.400547108329e+00 +wfi1s1k1_wfc_fp_power_7 1.484139852679e+00 +wfi1s1k1_wfc_fp_cross_0_re -1.073532352801e-01 +wfi1s1k1_wfc_fp_cross_0_im -7.149119469573e-01 +wfi1s1k1_wfc_fp_cross_1_re -8.771577497186e-01 +wfi1s1k1_wfc_fp_cross_1_im 5.509394669176e-01 +wfi1s1k1_wfc_fp_cross_2_re -6.304007682579e-01 +wfi1s1k1_wfc_fp_cross_2_im 2.207534239659e-01 +wfi1s1k1_wfc_fp_cross_3_re -3.629257360684e-01 +wfi1s1k1_wfc_fp_cross_3_im -2.488607308272e-01 +wfi1s1k1_wfc_fp_cross_4_re 4.602064046585e-01 +wfi1s1k1_wfc_fp_cross_4_im 5.127322008066e-01 +wfi1s1k1_wfc_fp_cross_5_re 3.796007637900e-01 +wfi1s1k1_wfc_fp_cross_5_im 1.044842718415e+00 +wfi1s1k1_wfc_fp_cross_6_re 1.008054980940e+00 +wfi1s1k1_wfc_fp_cross_6_im 1.030743874311e+00 +wfi1s1k1_wfc_fp_cross_7_re -6.305721235833e-01 +wfi1s1k1_wfc_fp_cross_7_im 7.059530247421e-01 +wfi1s1k2_wfc_fp_components 1 +wfi1s1k2_wfc_fp_nx 32 +wfi1s1k2_wfc_fp_ny 32 +wfi1s1k2_wfc_fp_nz 32 +wfi1s1k2_wfc_fp_voxel 3.051757812500e-02 +wfi1s1k2_wfc_fp_rms 2.496013530393e-02 +wfi1s1k2_wfc_fp_power_0 3.012855832058e+00 +wfi1s1k2_wfc_fp_power_1 1.426974221681e+00 +wfi1s1k2_wfc_fp_power_2 5.011415028638e-01 +wfi1s1k2_wfc_fp_power_3 1.523351042733e-01 +wfi1s1k2_wfc_fp_power_4 2.958463203107e+00 +wfi1s1k2_wfc_fp_power_5 5.369897182208e-01 +wfi1s1k2_wfc_fp_power_6 1.244952727688e-01 +wfi1s1k2_wfc_fp_power_7 1.002522268985e-01 +wfi1s1k2_wfc_fp_cross_0_re -1.844514747575e+00 +wfi1s1k2_wfc_fp_cross_0_im -9.471182354745e-01 +wfi1s1k2_wfc_fp_cross_1_re 1.495776213717e-01 +wfi1s1k2_wfc_fp_cross_1_im 8.323115649719e-01 +wfi1s1k2_wfc_fp_cross_2_re -5.816726606766e-02 +wfi1s1k2_wfc_fp_cross_2_im 2.701074087334e-01 +wfi1s1k2_wfc_fp_cross_3_re 6.137144459097e-01 +wfi1s1k2_wfc_fp_cross_3_im 2.720889182159e-01 +wfi1s1k2_wfc_fp_cross_4_re 1.252862399523e+00 +wfi1s1k2_wfc_fp_cross_4_im 1.378409578595e-01 +wfi1s1k2_wfc_fp_cross_5_re -2.043173440442e-01 +wfi1s1k2_wfc_fp_cross_5_im 1.584522147736e-01 +wfi1s1k2_wfc_fp_cross_6_re 6.704845851142e-02 +wfi1s1k2_wfc_fp_cross_6_im 8.936124744350e-02 +wfi1s1k2_wfc_fp_cross_7_re -1.454152274095e-01 +wfi1s1k2_wfc_fp_cross_7_im 5.299999227599e-01 +wfi1s1k3_wfc_fp_components 1 +wfi1s1k3_wfc_fp_nx 32 +wfi1s1k3_wfc_fp_ny 32 +wfi1s1k3_wfc_fp_nz 32 +wfi1s1k3_wfc_fp_voxel 3.051757812500e-02 +wfi1s1k3_wfc_fp_rms 2.495526418846e-02 +wfi1s1k3_wfc_fp_power_0 2.839936070637e+00 +wfi1s1k3_wfc_fp_power_1 1.426287181674e+00 +wfi1s1k3_wfc_fp_power_2 4.949068803198e-01 +wfi1s1k3_wfc_fp_power_3 1.712198731945e-01 +wfi1s1k3_wfc_fp_power_4 2.922648614231e+00 +wfi1s1k3_wfc_fp_power_5 6.333290094545e-01 +wfi1s1k3_wfc_fp_power_6 1.169543141988e-01 +wfi1s1k3_wfc_fp_power_7 1.225502210538e-01 +wfi1s1k3_wfc_fp_cross_0_re -1.769152629960e+00 +wfi1s1k3_wfc_fp_cross_0_im -9.595120563221e-01 +wfi1s1k3_wfc_fp_cross_1_re 1.430722861016e-01 +wfi1s1k3_wfc_fp_cross_1_im 8.278947158135e-01 +wfi1s1k3_wfc_fp_cross_2_re -2.642449496997e-02 +wfi1s1k3_wfc_fp_cross_2_im 2.898959112458e-01 +wfi1s1k3_wfc_fp_cross_3_re 6.289011459759e-01 +wfi1s1k3_wfc_fp_cross_3_im 3.238809560795e-01 +wfi1s1k3_wfc_fp_cross_4_re 1.348410648719e+00 +wfi1s1k3_wfc_fp_cross_4_im 1.810714617389e-01 +wfi1s1k3_wfc_fp_cross_5_re -2.082983475411e-01 +wfi1s1k3_wfc_fp_cross_5_im 1.751638044078e-01 +wfi1s1k3_wfc_fp_cross_6_re 6.232951210537e-02 +wfi1s1k3_wfc_fp_cross_6_im 1.022145243053e-01 +wfi1s1k3_wfc_fp_cross_7_re -1.293195261815e-01 +wfi1s1k3_wfc_fp_cross_7_im 5.755964327403e-01 +wfi1s2k1_wfc_fp_components 1 +wfi1s2k1_wfc_fp_nx 32 +wfi1s2k1_wfc_fp_ny 32 +wfi1s2k1_wfc_fp_nz 32 +wfi1s2k1_wfc_fp_voxel 3.051757812500e-02 +wfi1s2k1_wfc_fp_rms 4.880612460492e-02 +wfi1s2k1_wfc_fp_power_0 1.910494777217e+00 +wfi1s2k1_wfc_fp_power_1 1.103766518881e+00 +wfi1s2k1_wfc_fp_power_2 3.550339036901e-01 +wfi1s2k1_wfc_fp_power_3 3.238831076718e-02 +wfi1s2k1_wfc_fp_power_4 1.895240724960e+00 +wfi1s2k1_wfc_fp_power_5 7.215264876129e-01 +wfi1s2k1_wfc_fp_power_6 9.145785132890e-01 +wfi1s2k1_wfc_fp_power_7 1.025672423339e+00 +wfi1s2k1_wfc_fp_cross_0_re -8.141918868144e-01 +wfi1s2k1_wfc_fp_cross_0_im -1.202427436910e+00 +wfi1s2k1_wfc_fp_cross_1_re -5.298486011761e-01 +wfi1s2k1_wfc_fp_cross_1_im 3.333691584304e-01 +wfi1s2k1_wfc_fp_cross_2_re -1.066517727693e-01 +wfi1s2k1_wfc_fp_cross_2_im -1.115113316105e-02 +wfi1s2k1_wfc_fp_cross_3_re -2.099269568976e-01 +wfi1s2k1_wfc_fp_cross_3_im -1.315838833077e-01 +wfi1s2k1_wfc_fp_cross_4_re 9.958861686017e-01 +wfi1s2k1_wfc_fp_cross_4_im 6.129250546753e-01 +wfi1s2k1_wfc_fp_cross_5_re -3.919134059554e-02 +wfi1s2k1_wfc_fp_cross_5_im 8.113918049636e-01 +wfi1s2k1_wfc_fp_cross_6_re 6.288866395687e-01 +wfi1s2k1_wfc_fp_cross_6_im 7.365864203412e-01 +wfi1s2k1_wfc_fp_cross_7_re -7.897659009530e-01 +wfi1s2k1_wfc_fp_cross_7_im 1.155773173948e+00 +wfi1s2k2_wfc_fp_components 1 +wfi1s2k2_wfc_fp_nx 32 +wfi1s2k2_wfc_fp_ny 32 +wfi1s2k2_wfc_fp_nz 32 +wfi1s2k2_wfc_fp_voxel 3.051757812500e-02 +wfi1s2k2_wfc_fp_rms 4.992027060786e-02 +wfi1s2k2_wfc_fp_power_0 3.012855832062e+00 +wfi1s2k2_wfc_fp_power_1 1.426974221679e+00 +wfi1s2k2_wfc_fp_power_2 5.011415028664e-01 +wfi1s2k2_wfc_fp_power_3 1.523351042729e-01 +wfi1s2k2_wfc_fp_power_4 2.958463203106e+00 +wfi1s2k2_wfc_fp_power_5 5.369897182222e-01 +wfi1s2k2_wfc_fp_power_6 1.244952727698e-01 +wfi1s2k2_wfc_fp_power_7 1.002522268974e-01 +wfi1s2k2_wfc_fp_cross_0_re -1.844514747573e+00 +wfi1s2k2_wfc_fp_cross_0_im -9.471182354768e-01 +wfi1s2k2_wfc_fp_cross_1_re 1.495776213722e-01 +wfi1s2k2_wfc_fp_cross_1_im 8.323115649733e-01 +wfi1s2k2_wfc_fp_cross_2_re -5.816726606655e-02 +wfi1s2k2_wfc_fp_cross_2_im 2.701074087340e-01 +wfi1s2k2_wfc_fp_cross_3_re 6.137144459080e-01 +wfi1s2k2_wfc_fp_cross_3_im 2.720889182175e-01 +wfi1s2k2_wfc_fp_cross_4_re 1.252862399524e+00 +wfi1s2k2_wfc_fp_cross_4_im 1.378409578638e-01 +wfi1s2k2_wfc_fp_cross_5_re -2.043173440455e-01 +wfi1s2k2_wfc_fp_cross_5_im 1.584522147741e-01 +wfi1s2k2_wfc_fp_cross_6_re 6.704845851164e-02 +wfi1s2k2_wfc_fp_cross_6_im 8.936124744307e-02 +wfi1s2k2_wfc_fp_cross_7_re -1.454152274083e-01 +wfi1s2k2_wfc_fp_cross_7_im 5.299999227575e-01 +wfi1s2k3_wfc_fp_components 1 +wfi1s2k3_wfc_fp_nx 32 +wfi1s2k3_wfc_fp_ny 32 +wfi1s2k3_wfc_fp_nz 32 +wfi1s2k3_wfc_fp_voxel 3.051757812500e-02 +wfi1s2k3_wfc_fp_rms 4.991052837692e-02 +wfi1s2k3_wfc_fp_power_0 2.839936070630e+00 +wfi1s2k3_wfc_fp_power_1 1.426287181674e+00 +wfi1s2k3_wfc_fp_power_2 4.949068803192e-01 +wfi1s2k3_wfc_fp_power_3 1.712198731964e-01 +wfi1s2k3_wfc_fp_power_4 2.922648614235e+00 +wfi1s2k3_wfc_fp_power_5 6.333290094572e-01 +wfi1s2k3_wfc_fp_power_6 1.169543141998e-01 +wfi1s2k3_wfc_fp_power_7 1.225502210548e-01 +wfi1s2k3_wfc_fp_cross_0_re -1.769152629960e+00 +wfi1s2k3_wfc_fp_cross_0_im -9.595120563175e-01 +wfi1s2k3_wfc_fp_cross_1_re 1.430722860990e-01 +wfi1s2k3_wfc_fp_cross_1_im 8.278947158134e-01 +wfi1s2k3_wfc_fp_cross_2_re -2.642449496791e-02 +wfi1s2k3_wfc_fp_cross_2_im 2.898959112474e-01 +wfi1s2k3_wfc_fp_cross_3_re 6.289011459773e-01 +wfi1s2k3_wfc_fp_cross_3_im 3.238809560861e-01 +wfi1s2k3_wfc_fp_cross_4_re 1.348410648723e+00 +wfi1s2k3_wfc_fp_cross_4_im 1.810714617368e-01 +wfi1s2k3_wfc_fp_cross_5_re -2.082983475425e-01 +wfi1s2k3_wfc_fp_cross_5_im 1.751638044088e-01 +wfi1s2k3_wfc_fp_cross_6_re 6.232951210595e-02 +wfi1s2k3_wfc_fp_cross_6_im 1.022145243062e-01 +wfi1s2k3_wfc_fp_cross_7_re -1.293195261824e-01 +wfi1s2k3_wfc_fp_cross_7_im 5.755964327419e-01 +totaltimeref 0.76 diff --git a/tests/11_PW_GPU/CASES_GPU.txt b/tests/11_PW_GPU/CASES_GPU.txt index a10f2a682e2..a7a12920b4d 100644 --- a/tests/11_PW_GPU/CASES_GPU.txt +++ b/tests/11_PW_GPU/CASES_GPU.txt @@ -5,3 +5,4 @@ scf_dav scf_dav_sub scf_out_wf scf_out_wf_norm +scf_out_wf_spinor diff --git a/tests/11_PW_GPU/scf_out_wf_spinor/INPUT b/tests/11_PW_GPU/scf_out_wf_spinor/INPUT new file mode 100644 index 00000000000..3e43be181e7 --- /dev/null +++ b/tests/11_PW_GPU/scf_out_wf_spinor/INPUT @@ -0,0 +1,40 @@ +INPUT_PARAMETERS + +suffix autotest +calculation scf +init_wfc random +basis_type pw +device gpu +kpar 2 +noncolin 1 +nspin 4 +lspinorb 1 +symmetry 0 +cal_force 1 +cal_stress 1 + +ecutwfc 20 +ecutrho 160 +pw_seed 1 +scf_thr 1.0e-7 +scf_nmax 1 +out_chg 0 +out_pchg 1 +if_separate_k 1 +out_wfc_norm 1 +out_wfc_re_im 1 + +smearing_method gaussian +smearing_sigma 0.02 + +mixing_type broyden +mixing_beta 0.2 +mixing_ndim 10 + +ks_solver dav_subspace +pw_diag_thr 1.0e-12 +diago_smooth_ethr true +pw_diag_ndim 2 + +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB diff --git a/tests/11_PW_GPU/scf_out_wf_spinor/KPT b/tests/11_PW_GPU/scf_out_wf_spinor/KPT new file mode 100644 index 00000000000..b5b3bdb1ae2 --- /dev/null +++ b/tests/11_PW_GPU/scf_out_wf_spinor/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 2 0 0 0 diff --git a/tests/11_PW_GPU/scf_out_wf_spinor/README b/tests/11_PW_GPU/scf_out_wf_spinor/README new file mode 100644 index 00000000000..ad31ec6e6d1 --- /dev/null +++ b/tests/11_PW_GPU/scf_out_wf_spinor/README @@ -0,0 +1 @@ +GPU PW nspin=4 real-space partial-charge and spinor-wavefunction output from a one-step SOC operator probe with KPAR. diff --git a/tests/11_PW_GPU/scf_out_wf_spinor/STRU b/tests/11_PW_GPU/scf_out_wf_spinor/STRU new file mode 100644 index 00000000000..1565e73eb0d --- /dev/null +++ b/tests/11_PW_GPU/scf_out_wf_spinor/STRU @@ -0,0 +1,22 @@ +ATOMIC_SPECIES +Fe 1.000 Fe.upf + +NUMERICAL_ORBITAL +Fe_gga_9au_100Ry_4s2p2d1f.orb + +LATTICE_CONSTANT +6 + +LATTICE_VECTORS +1.00 0.00 0.00 +0.00 1.00 0.00 +0.00 0.00 1.00 + +ATOMIC_POSITIONS +Direct + +Fe +1.0 +2 +0.00 0.00 0.00 1 1 1 mag 1.0 1.0 1.0 +0.51 0.52 0.53 1 1 1 mag 1.0 1.0 1.0 diff --git a/tests/11_PW_GPU/scf_out_wf_spinor/result.ref b/tests/11_PW_GPU/scf_out_wf_spinor/result.ref new file mode 100644 index 00000000000..2767b554bcb --- /dev/null +++ b/tests/11_PW_GPU/scf_out_wf_spinor/result.ref @@ -0,0 +1,76 @@ +etotref -6448.1073841145562255 +etotperatomref -3224.0536920573 +totalforceref 15.338334 +totalstressref 110032.560580 +pchgi1s1k1.cube 1 +pchgi1s1k2.cube 1 +pchgi1s2k1.cube 0.5773502772 +pchgi1s2k2.cube 0.577350264 +pchgi1s3k1.cube 0.5773502629 +pchgi1s3k2.cube 0.5773502607 +pchgi1s4k1.cube 0.5773502668 +pchgi1s4k2.cube 0.5773502817 +wfi1s1k1.cube 5.895913845 +wfi1s1k2.cube 5.83429784 +wfi1k1_spinor_wfc_fp_components 2 +wfi1k1_spinor_wfc_fp_nx 25 +wfi1k1_spinor_wfc_fp_ny 25 +wfi1k1_spinor_wfc_fp_nz 25 +wfi1k1_spinor_wfc_fp_voxel 1.382400000000e-02 +wfi1k1_spinor_wfc_fp_rms 4.811252243247e-02 +wfi1k1_spinor_wfc_fp_power_0 4.823186139229e-01 +wfi1k1_spinor_wfc_fp_power_1 4.919438312365e+00 +wfi1k1_spinor_wfc_fp_power_2 2.829710799638e-01 +wfi1k1_spinor_wfc_fp_power_3 2.081068044140e+00 +wfi1k1_spinor_wfc_fp_power_4 5.477005787430e-02 +wfi1k1_spinor_wfc_fp_power_5 3.369275028192e-01 +wfi1k1_spinor_wfc_fp_power_6 2.190307306568e+00 +wfi1k1_spinor_wfc_fp_power_7 2.697692347463e-01 +wfi1k1_spinor_wfc_fp_cross_0_re -6.212963474650e-01 +wfi1k1_spinor_wfc_fp_cross_0_im -1.409513219777e+00 +wfi1k1_spinor_wfc_fp_cross_1_re -1.148921960979e+00 +wfi1k1_spinor_wfc_fp_cross_1_im -2.683972795048e-01 +wfi1k1_spinor_wfc_fp_cross_2_re -2.872580418202e-02 +wfi1k1_spinor_wfc_fp_cross_2_im 7.668486813593e-01 +wfi1k1_spinor_wfc_fp_cross_3_re 1.354870378005e-01 +wfi1k1_spinor_wfc_fp_cross_3_im 3.092304639035e-01 +wfi1k1_spinor_wfc_fp_cross_4_re 8.113881655331e-02 +wfi1k1_spinor_wfc_fp_cross_4_im 1.089496731394e-01 +wfi1k1_spinor_wfc_fp_cross_5_re 3.556879613430e-01 +wfi1k1_spinor_wfc_fp_cross_5_im -7.819596187555e-01 +wfi1k1_spinor_wfc_fp_cross_6_re 3.413052046510e-01 +wfi1k1_spinor_wfc_fp_cross_6_im -6.887585086445e-01 +wfi1k1_spinor_wfc_fp_cross_7_re -3.430392992538e-01 +wfi1k1_spinor_wfc_fp_cross_7_im -1.115292004335e-01 +wfi1k2_spinor_wfc_fp_components 2 +wfi1k2_spinor_wfc_fp_nx 25 +wfi1k2_spinor_wfc_fp_ny 25 +wfi1k2_spinor_wfc_fp_nz 25 +wfi1k2_spinor_wfc_fp_voxel 1.382400000000e-02 +wfi1k2_spinor_wfc_fp_rms 4.811252243247e-02 +wfi1k2_spinor_wfc_fp_power_0 1.148625734803e-01 +wfi1k2_spinor_wfc_fp_power_1 2.461537655935e+00 +wfi1k2_spinor_wfc_fp_power_2 1.590562029325e-02 +wfi1k2_spinor_wfc_fp_power_3 2.853259400966e+00 +wfi1k2_spinor_wfc_fp_power_4 8.924986828253e-01 +wfi1k2_spinor_wfc_fp_power_5 6.668603981526e-01 +wfi1k2_spinor_wfc_fp_power_6 3.954039023282e+00 +wfi1k2_spinor_wfc_fp_power_7 2.795818262702e-01 +wfi1k2_spinor_wfc_fp_cross_0_re 5.016579560579e-01 +wfi1k2_spinor_wfc_fp_cross_0_im -1.762890949638e-01 +wfi1k2_spinor_wfc_fp_cross_1_re 1.800957478739e-01 +wfi1k2_spinor_wfc_fp_cross_1_im -8.196221623755e-02 +wfi1k2_spinor_wfc_fp_cross_2_re 1.216731602046e-01 +wfi1k2_spinor_wfc_fp_cross_2_im -1.748671001524e-01 +wfi1k2_spinor_wfc_fp_cross_3_re 5.848583247221e-01 +wfi1k2_spinor_wfc_fp_cross_3_im 1.484746105273e+00 +wfi1k2_spinor_wfc_fp_cross_4_re 4.603614622071e-01 +wfi1k2_spinor_wfc_fp_cross_4_im 6.190632852093e-01 +wfi1k2_spinor_wfc_fp_cross_5_re 1.124163429914e+00 +wfi1k2_spinor_wfc_fp_cross_5_im -1.171771573397e+00 +wfi1k2_spinor_wfc_fp_cross_6_re -1.028088976048e+00 +wfi1k2_spinor_wfc_fp_cross_6_im -2.202510127127e-01 +wfi1k2_spinor_wfc_fp_cross_7_re -1.757685867231e-01 +wfi1k2_spinor_wfc_fp_cross_7_im -3.491263360254e-02 +pw_spinor_cube_identity 0 +totaltimeref 1.35 diff --git a/tests/integrate/tools/catch_properties.sh b/tests/integrate/tools/catch_properties.sh index 15ca5ffb879..7fae3369292 100755 --- a/tests/integrate/tools/catch_properties.sh +++ b/tests/integrate/tools/catch_properties.sh @@ -4,7 +4,7 @@ # this compare script is used in different integrate tests COMPARE_SCRIPT="../../integrate/tools/CompareFile.py" #COMPARE_SCRIPT="../../integrate/tools/compare_file.py" -SUM_CUBE_EXE="python3 ../../integrate/tools/sum_cube.py" +CUBE_TOOL="../../integrate/tools/cube_tool.py" COLLECT_NPY_MEANS="../../integrate/tools/collect_npy_means.py" @@ -711,16 +711,71 @@ fi # Process .cube files if needed if [ "$need_process_cube" = true ]; then cubefiles=$(ls OUT.autotest/ | grep -E '.cube$') - + wavefunction_re_files=() + if [ -z "$cubefiles" ]; then echo "Error: No .cube files found in OUT.autotest/" exit 1 else for cube in $cubefiles; do - total_chg=$($SUM_CUBE_EXE OUT.autotest/$cube) + if [[ "$cube" =~ ^wfi[0-9]+s[0-9]+(k[0-9]+)?re[.]cube$ ]]; then + wavefunction_re_files+=("$cube") + continue + fi + if [[ "$cube" =~ ^wfi[0-9]+s[0-9]+(k[0-9]+)?im[.]cube$ ]]; then + continue + fi + total_chg=$(python3 "$CUBE_TOOL" integrate "OUT.autotest/$cube") echo "$cube $total_chg" >> $1 done fi + + for cube in "${wavefunction_re_files[@]}"; do + if [[ "$cube" =~ ^wfi([0-9]+)s([0-9]+)(k[0-9]+)?re[.]cube$ ]]; then + band=${BASH_REMATCH[1]} + spin=${BASH_REMATCH[2]} + kpoint=${BASH_REMATCH[3]} + state_prefix=${cube%re.cube} + fingerprint_args=( + "OUT.autotest/$cube" + "OUT.autotest/${state_prefix}im.cube" + ) + if [ "$nspin" = "4" ]; then + if [ "$spin" != "1" ]; then + continue + fi + lower_prefix="wfi${band}s2${kpoint}" + fingerprint_args+=( + "OUT.autotest/${lower_prefix}re.cube" + "OUT.autotest/${lower_prefix}im.cube" + ) + result_prefix="wfi${band}${kpoint}_spinor_wfc_fp" + else + result_prefix="${state_prefix}_wfc_fp" + fi + + if fingerprint=$(python3 "$CUBE_TOOL" fingerprint-wfc "${fingerprint_args[@]}"); then + while read -r metric value; do + echo "${result_prefix}_${metric} $value" >> "$1" + done <<< "$fingerprint" + else + echo "Error: Failed to generate wavefunction fingerprint for $state_prefix" + exit 1 + fi + fi + done +fi + +# Check the pointwise Pauli identities when all PW nspin=4 spinor outputs are available. +nspin=$(get_input_key_value "nspin" "INPUT") +if_separate_k=$(get_input_key_value "if_separate_k" "INPUT") +if [ "$nspin" = "4" ] && { [ "$if_separate_k" = "1" ] || [ "$if_separate_k" = "true" ]; } \ + && [ -n "$out_wfc_norm" ] && [ -n "$out_wfc_re_im" ] && [ -n "$out_pchg" ]; then + if python3 "$CUBE_TOOL" check-spinor OUT.autotest; then + echo "pw_spinor_cube_identity 0" >> "$1" + else + echo "pw_spinor_cube_identity 1" >> "$1" + fi fi #-------------------------------------------- diff --git a/tests/integrate/tools/cube_tool.py b/tests/integrate/tools/cube_tool.py new file mode 100755 index 00000000000..833e8443423 --- /dev/null +++ b/tests/integrate/tools/cube_tool.py @@ -0,0 +1,462 @@ +#!/usr/bin/env python3 +"""Utilities for cube files used by ABACUS integration tests. + +This file intentionally keeps the cube parser and all cube-specific regression +checks together. Use one of the following subcommands: + +``integrate CUBE`` + Integrate one scalar cube file over real space. This is suitable for charge + density, partial charge, wavefunction modulus, potential, and other scalar + cube outputs. The command prints one number: ``sum(values) * voxel_volume``. + +``compare-wfc CAL_RE CAL_IM REF_RE REF_IM [CAL_RE CAL_IM REF_RE REF_IM ...]`` + Compare one complex wavefunction state after removing its arbitrary global + U(1) phase. Pass one group of four files for a scalar wavefunction or for an + independently aligned spin channel. Pass two groups for an nspin=4 spinor; + all groups then share the same fitted phase. The command prints the maximum + pointwise complex error after alignment. It does not perform unitary + alignment within a degenerate subspace. + +``fingerprint-wfc RE IM [RE IM ...]`` + Generate a compact, deterministic fingerprint of one complex wavefunction + state. One real/imaginary pair represents a scalar component; two pairs in + upper/lower order represent an nspin=4 spinor. The fingerprint contains + grid metadata, the RMS amplitude, and phase-invariant products of eight + fixed complex projections. It is insensitive to one common U(1) phase but + remains sensitive to spatial changes and relative spinor phases. This is a + regression summary, not a proof of pointwise equality and not an invariant + of unitary rotations within a degenerate band subspace. + + For flattened values ``psi[i]``, let ``E = sum(abs(psi[i])**2)`` and + ``q[j] = sum(w[j,i] * psi[i]) / sqrt(E)``. Eight deterministic SplitMix64 + sketches map each ``(j, i)`` to one of ``{1, -1, 1j, -1j}``. The reported + powers ``abs(q[j])**2`` and cyclic products + ``q[j] * conj(q[(j + 1) % 8])`` are unchanged by one common phase applied + to every component. The constants and flattening order are part of the + reference format and must not be changed without regenerating references. + +``check-spinor DIRECTORY [--tolerance VALUE]`` + Check pointwise Pauli identities among separate-k PW nspin=4 outputs in one + ABACUS output directory. This requires wavefunction modulus, spinor real and + imaginary parts, and all four partial-charge components. It also checks that + their grids agree with each other and, when available, with the dense grid + reported in the running log. Successful checks produce no output. + +Examples: + + python3 cube_tool.py integrate OUT.autotest/charge.cube + python3 cube_tool.py compare-wfc cal_re.cube cal_im.cube ref_re.cube ref_im.cube + python3 cube_tool.py fingerprint-wfc state_re.cube state_im.cube + python3 cube_tool.py check-spinor OUT.autotest --tolerance 1e-8 + +All commands use only the Python standard library. Invalid or missing inputs are +reported on standard error and return a nonzero exit status. +""" + +import argparse +import math +import re +import sys +from pathlib import Path + + +FINGERPRINT_SKETCHES = 8 +UINT64_MASK = (1 << 64) - 1 +SPLITMIX_GOLDEN = 0x9E3779B97F4A7C15 +SPLITMIX_MIX1 = 0xBF58476D1CE4E5B9 +SPLITMIX_MIX2 = 0x94D049BB133111EB +FINGERPRINT_WEIGHTS = (1.0 + 0.0j, -1.0 + 0.0j, 0.0 + 1.0j, 0.0 - 1.0j) + + +def read_cube(path): + """Return the grid shape, voxel volume, and scalar data from a cube file.""" + input_path = Path(path) + with input_path.open("r", encoding="utf-8") as stream: + stream.readline() + stream.readline() + + atom_line = stream.readline().split() + if len(atom_line) < 4: + raise ValueError(f"invalid atom header in {input_path}") + natom = int(atom_line[0]) + + shape = [] + vectors = [] + for _ in range(3): + grid_line = stream.readline().split() + if len(grid_line) < 4: + raise ValueError(f"invalid grid header in {input_path}") + shape.append(int(grid_line[0])) + vectors.append([float(value) for value in grid_line[1:4]]) + + v1, v2, v3 = vectors + val0 = v2[1] * v3[2] - v2[2] * v3[1] + val1 = v2[0] * v3[2] - v2[2] * v3[0] + val2 = v2[0] * v3[1] - v2[1] * v3[0] + volume = abs(v1[0] * val0 - v1[1] * val1 + v1[2] * val2) + + # Preserve the existing ABACUS cube reader behavior for negative atom counts. + for _ in range(max(natom, 0)): + stream.readline() + + count = shape[0] * shape[1] * shape[2] + values = [] + for line in stream: + values.extend(float(value) for value in line.split()) + if len(values) >= count: + break + + if len(values) < count: + raise ValueError(f"expected {count} grid values in {input_path}, found {len(values)}") + return tuple(shape), volume, values[:count] + + +def integrate_cube(path): + _, volume, values = read_cube(path) + return sum(values) * volume + + +def splitmix64(value): + value = (value + SPLITMIX_GOLDEN) & UINT64_MASK + value = ((value ^ (value >> 30)) * SPLITMIX_MIX1) & UINT64_MASK + value = ((value ^ (value >> 27)) * SPLITMIX_MIX2) & UINT64_MASK + return value ^ (value >> 31) + + +def read_wfc_components(paths): + if len(paths) == 0 or len(paths) % 2 != 0: + raise ValueError("expected one or more RE IM cube path pairs") + + common_shape = None + common_volume = None + components = [] + for index in range(0, len(paths), 2): + real_path = paths[index] + imag_path = paths[index + 1] + real_shape, real_volume, real_values = read_cube(real_path) + imag_shape, imag_volume, imag_values = read_cube(imag_path) + if imag_shape != real_shape: + raise ValueError( + f"grid shape mismatch in {imag_path}: expected {real_shape}, found {imag_shape}" + ) + if imag_volume != real_volume: + raise ValueError( + f"voxel volume mismatch in {imag_path}: expected {real_volume}, found {imag_volume}" + ) + if common_shape is None: + common_shape = real_shape + common_volume = real_volume + elif real_shape != common_shape: + raise ValueError( + f"component grid shape mismatch in {real_path}: " + f"expected {common_shape}, found {real_shape}" + ) + elif real_volume != common_volume: + raise ValueError( + f"component voxel volume mismatch in {real_path}: " + f"expected {common_volume}, found {real_volume}" + ) + components.append( + [complex(real, imag) for real, imag in zip(real_values, imag_values)] + ) + return common_shape, common_volume, components + + +def wfc_fingerprint(paths): + shape, volume, components = read_wfc_components(paths) + point_count = sum(len(component) for component in components) + energy = sum( + value.real * value.real + value.imag * value.imag + for component in components + for value in component + ) + if not math.isfinite(energy) or energy == 0.0: + raise ValueError("wavefunction has zero or non-finite norm") + + projections = [0.0j] * FINGERPRINT_SKETCHES + flat_index = 0 + for component in components: + for value in component: + for sketch in range(FINGERPRINT_SKETCHES): + weight_index = splitmix64(flat_index + sketch * SPLITMIX_GOLDEN) & 3 + projections[sketch] += FINGERPRINT_WEIGHTS[weight_index] * value + flat_index += 1 + + norm = math.sqrt(energy) + normalized = [projection / norm for projection in projections] + metrics = [ + ("components", len(components)), + ("nx", shape[0]), + ("ny", shape[1]), + ("nz", shape[2]), + ("voxel", volume), + ("rms", math.sqrt(energy / point_count)), + ] + metrics.extend( + (f"power_{sketch}", abs(normalized[sketch]) ** 2) + for sketch in range(FINGERPRINT_SKETCHES) + ) + for sketch in range(FINGERPRINT_SKETCHES): + next_sketch = (sketch + 1) % FINGERPRINT_SKETCHES + cross = normalized[sketch] * normalized[next_sketch].conjugate() + metrics.append((f"cross_{sketch}_re", cross.real)) + metrics.append((f"cross_{sketch}_im", cross.imag)) + + if any(not math.isfinite(float(value)) for _, value in metrics): + raise ValueError("wavefunction fingerprint contains a non-finite value") + return metrics + + +def format_metric(value): + if isinstance(value, int): + return str(value) + return f"{value:.12e}" + + +def read_wfc_component(calculated_re, calculated_im, reference_re, reference_im): + paths = (calculated_re, calculated_im, reference_re, reference_im) + cubes = [read_cube(path) for path in paths] + shape = cubes[0][0] + for path, cube in zip(paths[1:], cubes[1:]): + if cube[0] != shape: + raise ValueError(f"grid shape mismatch in {path}: expected {shape}, found {cube[0]}") + + calculated = [complex(real, imag) for real, imag in zip(cubes[0][2], cubes[1][2])] + reference = [complex(real, imag) for real, imag in zip(cubes[2][2], cubes[3][2])] + return shape, calculated, reference + + +def phase_aligned_error(file_groups): + calculated_components = [] + reference_components = [] + common_shape = None + for group in file_groups: + shape, calculated, reference = read_wfc_component(*group) + if common_shape is None: + common_shape = shape + elif shape != common_shape: + raise ValueError( + f"spinor component grid shape mismatch in {group[0]}: " + f"expected {common_shape}, found {shape}" + ) + calculated_components.append(calculated) + reference_components.append(reference) + + calculated_norm = sum( + abs(value) ** 2 for component in calculated_components for value in component + ) + reference_norm = sum( + abs(value) ** 2 for component in reference_components for value in component + ) + if not math.isfinite(calculated_norm) or calculated_norm == 0.0: + raise ValueError("calculated wavefunction has zero or non-finite norm") + if not math.isfinite(reference_norm) or reference_norm == 0.0: + raise ValueError("reference wavefunction has zero or non-finite norm") + + overlap = sum( + reference.conjugate() * calculated + for calculated_component, reference_component in zip( + calculated_components, reference_components + ) + for calculated, reference in zip(calculated_component, reference_component) + ) + if ( + not math.isfinite(overlap.real) + or not math.isfinite(overlap.imag) + or abs(overlap) == 0.0 + ): + raise ValueError("wavefunction overlap is zero or non-finite; the global phase is undefined") + + phase = overlap.conjugate() / abs(overlap) + return max( + abs(phase * calculated - reference) + for calculated_component, reference_component in zip( + calculated_components, reference_components + ) + for calculated, reference in zip(calculated_component, reference_component) + ) + + +def parse_wfc_groups(paths): + if len(paths) % 4 != 0: + raise ValueError("expected four cube paths per component: CAL_RE CAL_IM REF_RE REF_IM") + return [paths[index : index + 4] for index in range(0, len(paths), 4)] + + +def read_dense_shape(directory): + pattern = re.compile( + r"fft grid for dense charge/potential\s*=\s*\[\s*([0-9]+),\s*([0-9]+),\s*([0-9]+)\s*\]" + ) + for log_path in sorted(directory.glob("running_*.log")): + with log_path.open("r", encoding="utf-8") as stream: + for line in stream: + match = pattern.search(line) + if match: + return tuple(int(value) for value in match.groups()) + return None + + +def max_abs_difference(actual, expected): + return max(abs(lhs - rhs) for lhs, rhs in zip(actual, expected)) + + +def check_cube_metadata(path): + with path.open("r", encoding="utf-8") as stream: + stream.readline() + metadata = stream.readline() + if "Fermi energy" in metadata: + raise ValueError(f"unexpected Fermi-energy metadata in {path}") + + +def check_spinor_state(directory, band, kpoint, tolerance, expected_shape): + prefix = f"wfi{band}" + pchg_prefix = f"pchgi{band}" + paths = { + "norm": directory / f"{prefix}s1k{kpoint}.cube", + "up_re": directory / f"{prefix}s1k{kpoint}re.cube", + "up_im": directory / f"{prefix}s1k{kpoint}im.cube", + "down_re": directory / f"{prefix}s2k{kpoint}re.cube", + "down_im": directory / f"{prefix}s2k{kpoint}im.cube", + "rho0": directory / f"{pchg_prefix}s1k{kpoint}.cube", + "mx": directory / f"{pchg_prefix}s2k{kpoint}.cube", + "my": directory / f"{pchg_prefix}s3k{kpoint}.cube", + "mz": directory / f"{pchg_prefix}s4k{kpoint}.cube", + } + missing = [str(path) for path in paths.values() if not path.is_file()] + if missing: + raise ValueError("missing spinor output files: " + ", ".join(missing)) + + shape = None + data = {} + for name, path in paths.items(): + check_cube_metadata(path) + current_shape, _, data[name] = read_cube(path) + if shape is None: + shape = current_shape + elif current_shape != shape: + raise ValueError(f"grid shape mismatch in {path}") + if expected_shape is not None and shape != expected_shape: + raise ValueError(f"output grid {shape} does not match dense charge grid {expected_shape}") + + rho0_expected = [] + mx_expected = [] + my_expected = [] + mz_expected = [] + norm_squared = [] + for up_re, up_im, down_re, down_im, norm in zip( + data["up_re"], data["up_im"], data["down_re"], data["down_im"], data["norm"] + ): + up_norm = up_re * up_re + up_im * up_im + down_norm = down_re * down_re + down_im * down_im + rho0_expected.append(up_norm + down_norm) + mx_expected.append(2.0 * (up_re * down_re + up_im * down_im)) + my_expected.append(2.0 * (up_re * down_im - down_re * up_im)) + mz_expected.append(up_norm - down_norm) + norm_squared.append(norm * norm) + + errors = { + "norm": max_abs_difference(norm_squared, rho0_expected), + "rho0": max_abs_difference(data["rho0"], rho0_expected), + "mx": max_abs_difference(data["mx"], mx_expected), + "my": max_abs_difference(data["my"], my_expected), + "mz": max_abs_difference(data["mz"], mz_expected), + } + failed = { + name: error + for name, error in errors.items() + if not math.isfinite(error) or error > tolerance + } + if failed: + details = ", ".join(f"{name}={error:.6e}" for name, error in failed.items()) + raise ValueError(f"band {band}, k-point {kpoint}: pointwise identity failure ({details})") + + +def check_spinor_directory(directory, tolerance): + pattern = re.compile(r"^pchgi([0-9]+)s1k([0-9]+)[.]cube$") + states = [] + for path in directory.iterdir(): + match = pattern.match(path.name) + if match: + states.append((int(match.group(1)), int(match.group(2)))) + + if not states: + raise ValueError(f"no separate-k PW spinor partial-density files found in {directory}") + + expected_shape = read_dense_shape(directory) + for band, kpoint in sorted(states): + check_spinor_state(directory, band, kpoint, tolerance, expected_shape) + + +def run_integrate(args): + print(f"{integrate_cube(args.cube):.10g}") + + +def run_compare_wfc(args): + print(f"{phase_aligned_error(parse_wfc_groups(args.cube_paths)):.10g}") + + +def run_fingerprint_wfc(args): + for name, value in wfc_fingerprint(args.cube_paths): + print(f"{name} {format_metric(value)}") + + +def run_check_spinor(args): + check_spinor_directory(args.directory, args.tolerance) + + +def build_parser(): + parser = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter + ) + subparsers = parser.add_subparsers(dest="command") + + integrate_parser = subparsers.add_parser( + "integrate", help="integrate one scalar cube file" + ) + integrate_parser.add_argument("cube", type=Path) + integrate_parser.set_defaults(handler=run_integrate) + + compare_parser = subparsers.add_parser( + "compare-wfc", help="compare one complex state after global phase alignment" + ) + compare_parser.add_argument( + "cube_paths", + nargs="+", + help="CAL_RE CAL_IM REF_RE REF_IM, repeated for spinor components", + ) + compare_parser.set_defaults(handler=run_compare_wfc) + + fingerprint_parser = subparsers.add_parser( + "fingerprint-wfc", + help="generate a deterministic global-phase-invariant wavefunction fingerprint", + ) + fingerprint_parser.add_argument( + "cube_paths", + nargs="+", + help="RE IM, repeated in component order for spinors", + ) + fingerprint_parser.set_defaults(handler=run_fingerprint_wfc) + + spinor_parser = subparsers.add_parser( + "check-spinor", help="check pointwise PW nspin=4 Pauli identities" + ) + spinor_parser.add_argument("directory", type=Path) + spinor_parser.add_argument("--tolerance", type=float, default=1.0e-8) + spinor_parser.set_defaults(handler=run_check_spinor) + return parser + + +def main(): + parser = build_parser() + args = parser.parse_args() + if not hasattr(args, "handler"): + parser.error("a subcommand is required") + try: + args.handler(args) + except (OSError, ValueError) as error: + print(error, file=sys.stderr) + sys.exit(1) + + +if __name__ == "__main__": + main() diff --git a/tests/integrate/tools/sum_cube.py b/tests/integrate/tools/sum_cube.py deleted file mode 100755 index 3e5657fdf3a..00000000000 --- a/tests/integrate/tools/sum_cube.py +++ /dev/null @@ -1,96 +0,0 @@ -#!/usr/bin/env python3 -import sys - -def main(): - if len(sys.argv) < 2: - print(f"Can't find {sys.argv[1] if len(sys.argv) > 1 else 'file'} !") - sys.exit(1) - - input_file = sys.argv[1] - - try: - with open(input_file, 'r') as inp: - # skip the first two lines - inp.readline() - inp.readline() - - # read the 3rd line: number of atoms + origin coordinates - line = inp.readline().split() - if not line: - return - natom = int(line[0]) - # origin_x = float(line[1]) - # origin_y = float(line[2]) - # origin_z = float(line[3]) - - # read the grid vectors (support non-orthogonal) - line = inp.readline().split() - nx = int(line[0]) - v1 = [float(line[1]), float(line[2]), float(line[3])] - - line = inp.readline().split() - ny = int(line[0]) - v2 = [float(line[1]), float(line[2]), float(line[3])] - - line = inp.readline().split() - nz = int(line[0]) - v3 = [float(line[1]), float(line[2]), float(line[3])] - - # calculate the volume element |v1 · (v2 × v3)| - val0 = v2[1] * v3[2] - v2[2] * v3[1] - val1 = v2[0] * v3[2] - v2[2] * v3[0] - val2 = v2[0] * v3[1] - v2[1] * v3[0] - - volume = abs(v1[0] * val0 - v1[1] * val1 + v1[2] * val2) - - # skip the atom coordinates - # natom can be negative in cube files sometimes? - # C++ code: for (int i = 0; i < natom; ++i) - # If natom is negative, loop doesn't run. - - atoms_to_skip = natom - if atoms_to_skip < 0: - # In some cube formats, negative natom means second line of header contains units or something? - # Standard: "If the number of atoms is negative, that indicates that the file contains input lines... One line for each non-zero E value." - # But the C++ code simple loops < natom. If natom < 0, it skips 0 lines. - # We will mimic C++ behavior assuming it works for the files they have. - atoms_to_skip = 0 # Loop won't run if natom < 0 - - for _ in range(atoms_to_skip): - inp.readline() - - nr = nx * ny * nz - - total_sum = 0.0 - count = 0 - - # Read grid values - # iterate over remaining lines to handle values spread across lines - for line in inp: - parts = line.split() - for part in parts: - total_sum += float(part) - count += 1 - if count >= nr: - break - if count >= nr: - break - - ne = total_sum * volume - # cout default precision is 6, but setprecision(10) changes it. - # Python's default float printing is usually sufficient, but let's use formatting to be sure. - # {:.10g} prints up to 10 significant digits. - print(f"{ne:.10g}") - - except FileNotFoundError: - print(f"Can't find {input_file} !") - sys.exit(1) - except Exception as e: - # C++ doesn't print other errors generally, but let's be safe. - # But to be functionally equivalent, maybe we shouldn't. - # The C++ code crashes or behaves weirdly on bad input. - # Let's just exit 1. - sys.exit(1) - -if __name__ == "__main__": - main()