diff --git a/source/source_estate/module_dm/cal_dm_psi.cpp b/source/source_estate/module_dm/cal_dm_psi.cpp index 5ad2144c5bc..cd1e0402625 100644 --- a/source/source_estate/module_dm/cal_dm_psi.cpp +++ b/source/source_estate/module_dm/cal_dm_psi.cpp @@ -24,6 +24,9 @@ void cal_dm_psi(const Parallel_Orbitals* ParaV, // dm = wfc.T * wg * wfc.conj() // dm[is](iw1,iw2) = \sum_{ib} wfc[is](ib,iw1).T * wg(is,ib) * wfc[is](ib,iw2).conj() + // i.e. DMK(μ,ν) = \sum_ib wg_ib * C*_{μ,ib} * C_{ν,ib} = (C f C†)^T = D_std^T. + // Consumers that read DMK elements explicitly by (μ,ν) must follow this + // transposed convention, not treat DMK as C f C† itself. for (int ik = 0; ik < wfc.get_nk(); ++ik) { diff --git a/source/source_estate/module_dm/density_matrix.cpp b/source/source_estate/module_dm/density_matrix.cpp index 586010844bb..f2a15bdaab8 100644 --- a/source/source_estate/module_dm/density_matrix.cpp +++ b/source/source_estate/module_dm/density_matrix.cpp @@ -110,13 +110,14 @@ void DensityMatrix_Tools::cal_DMR( for(int ik = 0; ik < dm._nk; ++ik) { if(ik_in >= 0 && ik_in != ik) { continue; } - // cal k_phase - // if TK==std::complex, kphase is e^{ikR} + // Inverse Fourier transform: D(R) = sum_k D(k) * exp(-i*k*R) + // k-point weights are embedded in DMK, so there is no 1/Nk + // prefactor. const ModuleBase::Vector3 dR(R_index[0], R_index[1], R_index[2]); const double arg = (dm._kvec_d[ik] * dR) * ModuleBase::TWO_PI; double sinp, cosp; ModuleBase::libm::sincos(arg, &sinp, &cosp); - kphase_vec[ik][iR] = TK(cosp, sinp); + kphase_vec[ik][iR] = TK(cosp, -sinp); } } @@ -265,13 +266,14 @@ void DensityMatrix_Tools::cal_DMR_td( for(int ik = 0; ik < dm._nk; ++ik) { if(ik_in >= 0 && ik_in != ik) { continue; } - // cal k_phase - // if TK==std::complex, kphase is e^{ikR} + // Inverse Fourier transform: D(R) = sum_k D(k) * exp(-i*k*R) + // k-point weights are embedded in DMK, so there is no 1/Nk + // prefactor. const ModuleBase::Vector3 dR(R_index[0], R_index[1], R_index[2]); const double arg = (dm._kvec_d[ik] * dR) * ModuleBase::TWO_PI; double sinp, cosp; ModuleBase::libm::sincos(arg, &sinp, &cosp); - kphase_vec[ik][iR] = TK(cosp, sinp); + kphase_vec[ik][iR] = TK(cosp, -sinp); if(PARAM.inp.td_stype==2) { //phase for hybrid gauge tddft @@ -422,13 +424,14 @@ void DensityMatrix_Tools::cal_DMR_full( for(int ik = 0; ik < dm._nk; ++ik) { if(ik_in >= 0 && ik_in != ik) { continue; } - // cal k_phase - // if TK==std::complex, kphase is e^{ikR} + // Inverse Fourier transform: D(R) = sum_k D(k) * exp(-i*k*R) + // Phase factor: exp(-i*k*R) = cos(k·R) - i*sin(k·R) + // k-point weights are embedded in DMK, so there is no 1/Nk prefactor. const ModuleBase::Vector3 dR(R_index[0], R_index[1], R_index[2]); const double arg = (dm._kvec_d[ik] * dR) * ModuleBase::TWO_PI; double sinp, cosp; ModuleBase::libm::sincos(arg, &sinp, &cosp); - kphase_vec[ik][iR] = TK(cosp, sinp); + kphase_vec[ik][iR] = TK(cosp, -sinp); } } diff --git a/source/source_estate/module_dm/density_matrix.h b/source/source_estate/module_dm/density_matrix.h index 390f1e7884b..4c1513868bd 100644 --- a/source/source_estate/module_dm/density_matrix.h +++ b/source/source_estate/module_dm/density_matrix.h @@ -9,6 +9,26 @@ namespace elecstate { +// --------------------------------------------------------------------------- +// Density-matrix conventions (DMK/DMR), frozen by the DM/DMK/DMR chain. +// +// 1. DMK(μ,ν;k) = Σ_n f_nk C*_{μn}(k) C_{νn}(k) = (C f C†)^T = D_std^T. +// μ = row (bra) orbital, ν = column (ket) orbital. Stored as the transpose +// of the textbook D_std (cal_dm_psi.cpp); for Hermitian D this equals the +// conjugate. Consumers that read DMK elements explicitly by (μ,ν) must +// follow this convention. +// 2. D(R) = Σ_k e^{-ikR} DMK(k), k·R in direct coordinates with 2π. +// k-point weights w_k are embedded in DMK (wg = w_k*occ), so cal_DMR must +// NOT multiply by 1/Nk again. The forward pair used by operators is +// O(k) = Σ_R e^{+ikR} O(R) (folding_HR). +// 3. Closed-trace protection: a contraction Σ_{μνR} D(μ,ν;R)·O(μ,ν;R) equals +// Σ_k w_k·Re Tr(DMK(k)·O(k)) iff +// (a) the contraction is a closed Frobenius trace over all orbitals, +// (b) the paired operator is Hermitian per k (the HContainer stores +// (iat1,iat2,R) and (iat2,iat1,-R) with O(-R) = O(R)†), +// (c) only the real part is kept. +// --------------------------------------------------------------------------- + /** * @brief DensityMatrix Class * = for Gamma-only calculation diff --git a/source/source_estate/module_dm/test/CMakeLists.txt b/source/source_estate/module_dm/test/CMakeLists.txt index 02812bfa7db..ce81b22a261 100644 --- a/source/source_estate/module_dm/test/CMakeLists.txt +++ b/source/source_estate/module_dm/test/CMakeLists.txt @@ -46,5 +46,6 @@ AddTest( ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/base_matrix.cpp ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/hcontainer.cpp ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/atom_pair.cpp + ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/func_folding.cpp ${ABACUS_SOURCE_DIR}/source_basis/module_ao/parallel_orbitals.cpp ) diff --git a/source/source_estate/module_dm/test/test_cal_dm_r.cpp b/source/source_estate/module_dm/test/test_cal_dm_r.cpp index c150690d26d..ed2133907fb 100644 --- a/source/source_estate/module_dm/test/test_cal_dm_r.cpp +++ b/source/source_estate/module_dm/test/test_cal_dm_r.cpp @@ -1,11 +1,61 @@ #include +#include +#include +#include #include "gmock/gmock.h" #include "gtest/gtest.h" #include "source_estate/module_dm/density_matrix.h" #include "source_hamilt/module_hcontainer/hcontainer.h" +#include "source_hamilt/module_hcontainer/hcontainer_funcs.h" #include "source_cell/klist.h" +namespace +{ +// Fill a Hermitian matrix H (row-major, nw*nw) with deterministic pseudo-random +// values. If real_only is true, H is real symmetric (used for k-points that are +// their own time-reversal partner). +void fill_hermitian(std::vector>& H, const int nw, const unsigned seed, const bool real_only = false) +{ + H.assign(nw * nw, std::complex(0.0, 0.0)); + for (int i = 0; i < nw; ++i) + { + for (int j = i; j < nw; ++j) + { + const double re = 0.5 * std::sin(double(i * 7 + j * 13 + seed * 3)); + const double im = (i == j || real_only) ? 0.0 : 0.5 * std::cos(double(i * 11 + j * 5 + seed * 7)); + H[i * nw + j] = std::complex(re, im); + H[j * nw + i] = std::complex(re, -im); + } + } +} + +// Build an HContainer with the given (iat1, iat2, R) entries. +// If paired is true, both (iat1,iat2,R) and (iat2,iat1,-R) are inserted. +template +hamilt::HContainer build_r_container( + const Parallel_Orbitals* paraV, + const std::vector>& pairs, + const std::vector>& r_list, + const bool paired) +{ + hamilt::HContainer hR(paraV); + for (const auto& pr: pairs) + { + for (const auto& R: r_list) + { + hR.insert_pair(hamilt::AtomPair(pr.first, pr.second, R, paraV)); + if (paired) + { + hR.insert_pair(hamilt::AtomPair(pr.second, pr.first, -R, paraV)); + } + } + } + hR.allocate(nullptr, true); + return hR; +} +} // namespace + /************************************************ * unit test of DensityMatrix constructor ***********************************************/ @@ -338,6 +388,326 @@ TEST_F(DMTest, cal_DMR_blas_complex) delete kv; } +// T1: Fourier round-trip consistency between cal_DMR (e^{-ikR}) and +// folding_HR (e^{+ikR}). This is the test that locks the Fourier sign of the +// inverse transform: with the k-phase flipped (density_matrix.cpp: -sinp -> +// +sinp), D~(k) = DMK(-k) = DMK(k)† and the assertions below turn red. +TEST_F(DMTest, T1_fourier_round_trip) +{ + // k-grid {0, 1/4, 1/2, 3/4} along x: contains non-Gamma k-points and k and + // -k as distinct points (a {0, 1/2} grid is blind to the sign because + // -k == k there). + const int nk = 4; + const int nR = 4; + std::vector> kvec_d(nk); + for (int ik = 0; ik < nk; ++ik) + { + kvec_d[ik] = ModuleBase::Vector3(0.25 * ik, 0.0, 0.0); + } + + // Time-reversal constrained Hermitian DMK per k-point: k=0 and k=1/2 are + // their own TR partner (real symmetric), k=1/4 and k=3/4 are conjugate + // pairs. This makes D(R) real, so the real DMR path (which stores + // Re[D(R)]) is lossless and the full complex round trip is exact. + const int nw = test_nw; + std::vector>> M(nk); + fill_hermitian(M[0], nw, 1, true); + fill_hermitian(M[2], nw, 2, true); + fill_hermitian(M[1], nw, 3, false); + M[3].assign(nw * nw, std::complex(0.0, 0.0)); + for (int i = 0; i < nw * nw; ++i) + { + M[3][i] = std::conj(M[1][i]); + } + + // ---- real DMR path (nspin<4): DensityMatrix ---- + { + elecstate::DensityMatrix, double> DM(paraV, 1, kvec_d, nk); + // set_DMK stores the transposed element (see density_matrix.h): to set + // the element read by cal_DMR as (mu,nu), pass (nu,mu) to the setter. + for (int ik = 0; ik < nk; ++ik) + { + for (int i = 0; i < nw; ++i) + { + for (int j = 0; j < nw; ++j) + { + DM.set_DMK(1, ik, j, i, M[ik][i * nw + j]); + } + } + } + // synthetic HContainer: pair (0,0) with the complete set of R + // representatives {0,1,2,3} of the k-grid period + hamilt::HContainer hR(paraV); + for (int rx = 0; rx < nR; ++rx) + { + hR.insert_pair(hamilt::AtomPair(0, 0, rx, 0, 0, paraV)); + } + hR.allocate(nullptr, true); + DM.init_DMR(hR); + DM.cal_DMR(); + + for (int ik = 0; ik < nk; ++ik) + { + std::vector> hk(paraV->nrow * paraV->ncol, std::complex(0.0, 0.0)); + hamilt::folding_HR(*DM.get_DMR_pointer(1), hk.data(), kvec_d[ik], paraV->ncol, 0); + for (int i = 0; i < nw; ++i) + { + for (int j = 0; j < nw; ++j) + { + // folding_HR has no weight normalization: on a complete + // R-representative set D~(k) = nR * DMK(k) (in production + // the k weights are embedded in DMK). + const std::complex expect = double(nR) * M[ik][i * nw + j]; + EXPECT_NEAR(hk[i * paraV->ncol + j].real(), expect.real(), 1e-10) + << "real DMR path, ik=" << ik << " mu=" << i << " nu=" << j; + EXPECT_NEAR(hk[i * paraV->ncol + j].imag(), expect.imag(), 1e-10) + << "real DMR path, ik=" << ik << " mu=" << i << " nu=" << j; + } + } + } + } + + // ---- complex DMR path (nspin=4 via cal_DMR_full) ---- + { + elecstate::DensityMatrix, double> DM(paraV, 4, kvec_d, nk); + for (int ik = 0; ik < nk; ++ik) + { + for (int i = 0; i < nw; ++i) + { + for (int j = 0; j < nw; ++j) + { + DM.set_DMK(1, ik, j, i, M[ik][i * nw + j]); + } + } + } + hamilt::HContainer> hR(paraV); + for (int rx = 0; rx < nR; ++rx) + { + hR.insert_pair(hamilt::AtomPair>(0, 0, rx, 0, 0, paraV)); + } + hR.allocate(nullptr, true); + DM.cal_DMR_full(&hR); + + for (int ik = 0; ik < nk; ++ik) + { + std::vector> hk(paraV->nrow * paraV->ncol, std::complex(0.0, 0.0)); + hamilt::folding_HR(hR, hk.data(), kvec_d[ik], paraV->ncol, 0); + for (int i = 0; i < nw; ++i) + { + for (int j = 0; j < nw; ++j) + { + const std::complex expect = double(nR) * M[ik][i * nw + j]; + EXPECT_NEAR(hk[i * paraV->ncol + j].real(), expect.real(), 1e-10) + << "complex DMR path, ik=" << ik << " mu=" << i << " nu=" << j; + EXPECT_NEAR(hk[i * paraV->ncol + j].imag(), expect.imag(), 1e-10) + << "complex DMR path, ik=" << ik << " mu=" << i << " nu=" << j; + } + } + } + } +} + +// T2: DMR Hermiticity. For a Hermitian DMK and a paired HContainer +// ((iat1,iat2,R) and (iat2,iat1,-R)), cal_DMR must satisfy +// D(iat2,iat1,-R) = D(iat1,iat2,R)^T (real DMR, nspin<4) +// D(iat2,iat1,-R) = D(iat1,iat2,R)^\dagger (complex DMR, nspin=4) +TEST_F(DMTest, T2_dmr_hermiticity) +{ + // two k-points, both non-Gamma, k and -k distinct + std::vector> kvec_d = { + ModuleBase::Vector3(0.25, 0.0, 0.0), + ModuleBase::Vector3(0.75, 0.0, 0.0)}; + const int nk = 2; + const int nw = test_nw; + + // Hermitian DMK (arbitrary, no TRS constraint needed for Hermiticity) + std::vector>> M(nk); + for (int ik = 0; ik < nk; ++ik) + { + fill_hermitian(M[ik], nw, 10 + ik); + } + + // pairs: cross (0,1)/(1,0) and self (0,0); R: 3 non-zero + origin + const std::vector> pairs = {{0, 1}, {0, 0}}; + const std::vector> r_list = { + ModuleBase::Vector3(0, 0, 0), + ModuleBase::Vector3(1, 0, 0), + ModuleBase::Vector3(0, 1, 0), + ModuleBase::Vector3(0, 0, 1)}; + + // ---- real DMR path ---- + { + elecstate::DensityMatrix, double> DM(paraV, 1, kvec_d, nk); + for (int ik = 0; ik < nk; ++ik) + { + for (int i = 0; i < nw; ++i) + { + for (int j = 0; j < nw; ++j) + { + DM.set_DMK(1, ik, j, i, M[ik][i * nw + j]); + } + } + } + hamilt::HContainer hR = build_r_container(paraV, pairs, r_list, true); + DM.init_DMR(hR); + DM.cal_DMR(); + const hamilt::HContainer* dmr = DM.get_DMR_pointer(1); + for (int iap = 0; iap < dmr->size_atom_pairs(); ++iap) + { + const hamilt::AtomPair& ap = dmr->get_atom_pair(iap); + const int iat1 = ap.get_atom_i(); + const int iat2 = ap.get_atom_j(); + for (int ir = 0; ir < ap.get_R_size(); ++ir) + { + const ModuleBase::Vector3 R = ap.get_R_index(ir); + const hamilt::BaseMatrix* mat = ap.find_matrix(R); + const hamilt::BaseMatrix* mat_rev = dmr->find_matrix(iat2, iat1, -R.x, -R.y, -R.z); + ASSERT_NE(mat_rev, nullptr) << "missing reverse pair (" << iat2 << "," << iat1 << "," << -R.x << "," << -R.y << "," << -R.z << ")"; + const int row_size = ap.get_row_size(); + const int col_size = ap.get_col_size(); + for (int i = 0; i < row_size; ++i) + { + for (int j = 0; j < col_size; ++j) + { + EXPECT_NEAR(mat_rev->get_pointer()[j * col_size + i], + mat->get_pointer()[i * col_size + j], + 1e-10) + << "real DMR hermiticity: iat1=" << iat1 << " iat2=" << iat2 + << " R=(" << R.x << "," << R.y << "," << R.z << ")"; + } + } + } + } + } + + // ---- complex DMR path (cal_DMR_full) ---- + { + elecstate::DensityMatrix, double> DM(paraV, 4, kvec_d, nk); + for (int ik = 0; ik < nk; ++ik) + { + for (int i = 0; i < nw; ++i) + { + for (int j = 0; j < nw; ++j) + { + DM.set_DMK(1, ik, j, i, M[ik][i * nw + j]); + } + } + } + hamilt::HContainer> hR = build_r_container>(paraV, pairs, r_list, true); + DM.cal_DMR_full(&hR); + for (int iap = 0; iap < hR.size_atom_pairs(); ++iap) + { + const hamilt::AtomPair>& ap = hR.get_atom_pair(iap); + const int iat1 = ap.get_atom_i(); + const int iat2 = ap.get_atom_j(); + for (int ir = 0; ir < ap.get_R_size(); ++ir) + { + const ModuleBase::Vector3 R = ap.get_R_index(ir); + const hamilt::BaseMatrix>* mat = ap.find_matrix(R); + const hamilt::BaseMatrix>* mat_rev = hR.find_matrix(iat2, iat1, -R.x, -R.y, -R.z); + ASSERT_NE(mat_rev, nullptr) << "missing reverse pair (" << iat2 << "," << iat1 << "," << -R.x << "," << -R.y << "," << -R.z << ")"; + const int row_size = ap.get_row_size(); + const int col_size = ap.get_col_size(); + for (int i = 0; i < row_size; ++i) + { + for (int j = 0; j < col_size; ++j) + { + EXPECT_NEAR(mat_rev->get_pointer()[j * col_size + i].real(), + mat->get_pointer()[i * col_size + j].real(), + 1e-10) + << "complex DMR hermiticity real: iat1=" << iat1 << " iat2=" << iat2 + << " R=(" << R.x << "," << R.y << "," << R.z << ")"; + EXPECT_NEAR(mat_rev->get_pointer()[j * col_size + i].imag(), + -mat->get_pointer()[i * col_size + j].imag(), + 1e-10) + << "complex DMR hermiticity imag: iat1=" << iat1 << " iat2=" << iat2 + << " R=(" << R.x << "," << R.y << "," << R.z << ")"; + } + } + } + } + } +} + +// T8: full-direction pairing storage guard. Every (iat1,iat2,R) block of the +// DMR HContainer built by init_DMR must have a (iat2,iat1,-R) counterpart. +// gint_rho and the force/stress paths rely on this pairing (closed-trace +// protection). A future half-set storage optimization will turn this test +// red. +TEST_F(DMTest, T8_full_direction_pairing_guard) +{ + // multi-k DM (TK=complex) keeps the R structure of init_DMR(Record_adj); + // the TK=double path collapses to gamma-only by design. + const std::vector> kvec = {ModuleBase::Vector3(0.1, 0.0, 0.0)}; + elecstate::DensityMatrix, double> DM(paraV, 1, kvec, 1); + + // synthetic full symmetric neighbor list (as produced by Record_adj::cal_adj) + const int nat = test_size; + Record_adj ra; + ra.na_each = new int[nat]; + ra.info = new int**[nat]; + const std::vector> r_list = { + ModuleBase::Vector3(0, 0, 0), + ModuleBase::Vector3(1, 0, 0), + ModuleBase::Vector3(0, 1, 0), + ModuleBase::Vector3(0, 0, 1)}; + std::vector>> lists(nat); + for (int iat1 = 0; iat1 < nat; ++iat1) + { + for (int iat2 = 0; iat2 < nat; ++iat2) + { + for (const auto& R: r_list) + { + // both directions: (iat1,iat2,R) and (iat2,iat1,-R) + lists[iat1].push_back({R.x, R.y, R.z, 0, iat2}); + lists[iat2].push_back({-R.x, -R.y, -R.z, 0, iat1}); + } + } + } + for (int iat = 0; iat < nat; ++iat) + { + ra.na_each[iat] = lists[iat].size(); + ra.info[iat] = new int*[ra.na_each[iat]]; + for (int ad = 0; ad < ra.na_each[iat]; ++ad) + { + ra.info[iat][ad] = new int[5]; + for (int k = 0; k < 5; ++k) + { + ra.info[iat][ad][k] = lists[iat][ad][k]; + } + } + } + + DM.init_DMR(ra, &ucell); + + const hamilt::HContainer* dmr = DM.get_DMR_pointer(1); + EXPECT_GT(dmr->size_atom_pairs(), 0); + for (int iap = 0; iap < dmr->size_atom_pairs(); ++iap) + { + const hamilt::AtomPair& ap = dmr->get_atom_pair(iap); + const int iat1 = ap.get_atom_i(); + const int iat2 = ap.get_atom_j(); + for (int ir = 0; ir < ap.get_R_size(); ++ir) + { + const ModuleBase::Vector3 R = ap.get_R_index(ir); + EXPECT_NE(dmr->find_matrix(iat2, iat1, -R.x, -R.y, -R.z), nullptr) + << "pairing guard: missing (iat1,iat2,R)=(" << iat1 << "," << iat2 << "," + << R.x << "," << R.y << "," << R.z << ") reverse"; + } + } + + for (int iat = 0; iat < nat; ++iat) + { + for (int ad = 0; ad < ra.na_each[iat]; ++ad) + { + delete[] ra.info[iat][ad]; + } + delete[] ra.info[iat]; + } + delete[] ra.info; + delete[] ra.na_each; +} + int main(int argc, char** argv) { #ifdef __MPI diff --git a/source/source_io/module_mulliken/output_mulliken.cpp b/source/source_io/module_mulliken/output_mulliken.cpp index 3594bf1d04f..6545b2b45e4 100644 --- a/source/source_io/module_mulliken/output_mulliken.cpp +++ b/source/source_io/module_mulliken/output_mulliken.cpp @@ -444,7 +444,7 @@ void Output_Mulliken::collect_MW(ModuleBase::matrix& MecMulP, const ModuleBa const int ic = this->ParaV_->global2local_col(k2); // note that mud is column major MecMulP(1, j) += mud(ic, ir).real(); - // M_y = i(M_{up,down} - M_{down,up}) = -(M_{up,down} - M_{down,up}).imag() + // M_y = Re[i*(M_updown - M_downup)] = Im(M_downup) - Im(M_updown) MecMulP(2, j) -= mud(ic, ir).imag(); } if (this->ParaV_->in_this_processor(k2, k1)) @@ -452,7 +452,7 @@ void Output_Mulliken::collect_MW(ModuleBase::matrix& MecMulP, const ModuleBa const int ir = this->ParaV_->global2local_row(k2); const int ic = this->ParaV_->global2local_col(k1); MecMulP(1, j) += mud(ic, ir).real(); - // M_y = i(M_{up,down} - M_{down,up}) = -(M_{up,down} - M_{down,up}).imag() + // M_y = Re[i*(M_updown - M_downup)] = Im(M_downup) - Im(M_updown) MecMulP(2, j) += mud(ic, ir).imag(); } if (this->ParaV_->in_this_processor(k2, k2)) diff --git a/source/source_lcao/module_lr/dm_trans/dmr_complex.cpp b/source/source_lcao/module_lr/dm_trans/dmr_complex.cpp index 0b65bc610d8..990cf04b148 100644 --- a/source/source_lcao/module_lr/dm_trans/dmr_complex.cpp +++ b/source/source_lcao/module_lr/dm_trans/dmr_complex.cpp @@ -43,13 +43,14 @@ namespace elecstate for (int ik = 0; ik < this->_nk; ++ik) { if (ik_in >= 0 && ik_in != ik) continue; - // cal k_phase - // if TK==std::complex, kphase is e^{ikR} + // Inverse Fourier transform: D(R) = sum_k D(k) * exp(-i*k*R) + // Phase factor: exp(-i*k*R) = cos(k·R) - i*sin(k·R) + // k-point weights are embedded in DMK, so there is no 1/Nk prefactor. const ModuleBase::Vector3 dR(r_index[0], r_index[1], r_index[2]); const double arg = (this->_kvec_d[ik] * dR) * ModuleBase::TWO_PI; double sinp = 0.0, cosp = 0.0; ModuleBase::libm::sincos(arg, &sinp, &cosp); - const std::complex kphase = std::complex(cosp, sinp); + const std::complex kphase = std::complex(cosp, -sinp); // set DMR element std::complex* tmp_DMR_pointer = tmp_matrix->get_pointer(); const std::complex* tmp_DMK_pointer @@ -79,4 +80,4 @@ namespace elecstate ModuleBase::timer::end("DensityMatrix", "cal_DMR"); } // template class DensityMatrix, std::complex>; -} \ No newline at end of file +} diff --git a/source/source_lcao/module_lr/utils/exciton_plotter.cpp b/source/source_lcao/module_lr/utils/exciton_plotter.cpp index 5c121507d26..f8021eeb531 100644 --- a/source/source_lcao/module_lr/utils/exciton_plotter.cpp +++ b/source/source_lcao/module_lr/utils/exciton_plotter.cpp @@ -624,7 +624,7 @@ std::vector>> ExcitonPlotter::build_conditio { const Complex fixed_wfc = std::conj( this->orb_eval_ - .eval_wfc_bloch(r_fix, ik, io, this->psi_ks_vec[0], this->ucell, this->kv.kvec_d[ik])); + .template eval_wfc_bloch(r_fix, ik, io, this->psi_ks_vec[0], this->ucell, this->kv.kvec_d[ik])); for (int iv = 0; iv < nvirt; ++iv) { mixing[ik][iv] += this->X[offset_b + x_start + iv + io * nvirt] * fixed_wfc; @@ -635,7 +635,7 @@ std::vector>> ExcitonPlotter::build_conditio { for (int iv = 0; iv < nvirt; ++iv) { - const Complex fixed_wfc = this->orb_eval_.eval_wfc_bloch(r_fix, + const Complex fixed_wfc = this->orb_eval_.template eval_wfc_bloch(r_fix, ik, nocc + iv, this->psi_ks_vec[0], diff --git a/source/source_lcao/module_operator_lcao/operator_fs_utils.cpp b/source/source_lcao/module_operator_lcao/operator_fs_utils.cpp index 34ccc7c4ee2..68e7cc7eec5 100644 --- a/source/source_lcao/module_operator_lcao/operator_fs_utils.cpp +++ b/source/source_lcao/module_operator_lcao/operator_fs_utils.cpp @@ -18,7 +18,10 @@ void finalize_force_stress( #ifdef __MPI Parallel_Reduce::reduce_all(force.c, force.nr * force.nc); #endif - // Apply factor of 2 for Hermitian matrix + // force_factor is 1.0 because the force loop iterates the full R set + // ((iat1,iat2,R) and (iat2,iat1,-R) are each visited once, see + // cal_force_stress_2center in operator_fs_utils.hpp). A half-set + // iteration would require force_factor = 2. for (int i = 0; i < force.nr * force.nc; i++) { force.c[i] *= force_factor; diff --git a/source/source_lcao/module_operator_lcao/operator_fs_utils.hpp b/source/source_lcao/module_operator_lcao/operator_fs_utils.hpp index 7720903fa0b..c6a2816a967 100644 --- a/source/source_lcao/module_operator_lcao/operator_fs_utils.hpp +++ b/source/source_lcao/module_operator_lcao/operator_fs_utils.hpp @@ -145,7 +145,11 @@ void cal_force_stress_2center( // Calculate force contribution with compile-time sign if (cal_force) { - // Factor of 2 for Hermitian matrix will be applied later + // Full R-set summation: both (iat1,iat2,R) and + // (iat2,iat1,-R) pairs are visited once, so the + // finalize step uses factor=1.0 (no factor 2). + // If a half-set iteration is introduced in the + // future, restore factor=2 in finalize_force_stress. for (int i = 0; i < 3; i++) { force_tmp1[i] += ForceSign * dm_current * olm[i + 1]; diff --git a/source/source_lcao/module_operator_lcao/test/CMakeLists.txt b/source/source_lcao/module_operator_lcao/test/CMakeLists.txt index 4b988f029fc..fd249b505dd 100644 --- a/source/source_lcao/module_operator_lcao/test/CMakeLists.txt +++ b/source/source_lcao/module_operator_lcao/test/CMakeLists.txt @@ -90,6 +90,18 @@ AddTest( tmp_mocks.cpp ../../../source_hamilt/operator.cpp ) +AddTest( + TARGET MODULE_LCAO_operator_dm_trace_test + LIBS parameter psi base device container + SOURCES test_dm_trace.cpp ../../../source_estate/module_dm/density_matrix.cpp + ../../../source_estate/module_dm/density_matrix_io.cpp + ../../../source_hamilt/module_hcontainer/func_folding.cpp + ../../../source_hamilt/module_hcontainer/base_matrix.cpp ../../../source_hamilt/module_hcontainer/hcontainer.cpp ../../../source_hamilt/module_hcontainer/atom_pair.cpp + ../../../source_basis/module_ao/parallel_orbitals.cpp + ../../../source_basis/module_ao/orb_atomic_lm.cpp + tmp_mocks.cpp ../../../source_hamilt/operator.cpp +) + install(FILES parallel_operator_tests.sh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) find_program(BASH bash) add_test(NAME MODULE_LCAO_operators_para_test diff --git a/source/source_lcao/module_operator_lcao/test/test_dm_trace.cpp b/source/source_lcao/module_operator_lcao/test/test_dm_trace.cpp new file mode 100644 index 00000000000..688e2656f06 --- /dev/null +++ b/source/source_lcao/module_operator_lcao/test/test_dm_trace.cpp @@ -0,0 +1,353 @@ +#include +#include +#include +#include + +#include "gtest/gtest.h" +#include "source_estate/module_dm/density_matrix.h" +#include "source_hamilt/module_hcontainer/hcontainer_funcs.h" +#include "source_basis/module_ao/parallel_orbitals.h" + +/************************************************ + * T3: closed Frobenius trace equivalence + * + * For a Hermitian DMK(k), a real DMR obtained via cal_DMR (e^{-ikR}), and a + * paired Hermitian operator O(R) (O(iat2,iat1,-R) = O(iat1,iat2,R)^\dagger), + * the real-space contraction + * + * W_R = \sum_{mu,nu,R} D(mu,nu;R) * O(mu,nu;R) + * + * must equal the k-space contraction + * + * W_k = \sum_k Re Tr( DMK(k) * O(k) ), O(k) = folding_HR(O(R), k). + * + * This holds even when O(R) is NOT symmetric within the same R (the key case + * for overlap/kinetic derivative operators): the protection comes from the + * closed trace + per-k Hermiticity of O(k), not from same-R symmetry. + ************************************************/ + +namespace +{ +// value builders for T = double or std::complex +template +T mk_value(const double re, const double im); +template <> +double mk_value(const double re, const double) +{ + return re; +} +template <> +std::complex mk_value>(const double re, const double im) +{ + return std::complex(re, im); +} + +// conjugate transpose helper: identity for double, conj for complex +template +T conj_val(const T& v); +template <> +double conj_val(const double& v) +{ + return v; +} +template <> +std::complex conj_val>(const std::complex& v) +{ + return std::conj(v); +} + +// Hermitian matrix H (row-major, n*n), deterministic pseudo-random values +void fill_hermitian(std::vector>& H, const int n, const unsigned seed) +{ + H.assign(n * n, std::complex(0.0, 0.0)); + for (int i = 0; i < n; ++i) + { + for (int j = i; j < n; ++j) + { + const double re = 0.5 * std::sin(double(i * 7 + j * 13 + seed * 3)); + const double im = (i == j) ? 0.0 : 0.5 * std::cos(double(i * 11 + j * 5 + seed * 7)); + H[i * n + j] = std::complex(re, im); + H[j * n + i] = std::complex(re, -im); + } + } +} + +// HContainer with (iat1,iat2,R); if paired, also (iat2,iat1,-R) +template +hamilt::HContainer build_r_container( + const Parallel_Orbitals* paraV, + const std::vector>& pairs, + const std::vector>& r_list, + const bool paired) +{ + hamilt::HContainer hR(paraV); + for (const auto& pr: pairs) + { + for (const auto& R: r_list) + { + hR.insert_pair(hamilt::AtomPair(pr.first, pr.second, R, paraV)); + if (paired) + { + hR.insert_pair(hamilt::AtomPair(pr.second, pr.first, -R, paraV)); + } + } + } + hR.allocate(nullptr, true); + return hR; +} + +// Build a paired Hermitian operator O(R) on the container: +// - forward blocks (iat1,iat2,R) filled with deterministic values, +// symmetric within R if symmetric_in_R is true, +// - reverse blocks (iat2,iat1,-R) set to the transpose (real) or conjugate +// transpose (complex) of the forward blocks. +template +hamilt::HContainer build_operator_container( + const Parallel_Orbitals* paraV, + const std::vector>& pairs, + const std::vector>& r_list, + const unsigned seed, + const bool symmetric_in_R) +{ + auto hR = build_r_container(paraV, pairs, r_list, false); + for (int iap = 0; iap < hR.size_atom_pairs(); ++iap) + { + const hamilt::AtomPair& ap = hR.get_atom_pair(iap); + const int row_size = ap.get_row_size(); + const int col_size = ap.get_col_size(); + const bool self_pair = (ap.get_atom_i() == ap.get_atom_j()); + for (int ir = 0; ir < ap.get_R_size(); ++ir) + { + const ModuleBase::Vector3 R = ap.get_R_index(ir); + T* mat = ap.get_pointer(ir); + for (int i = 0; i < row_size; ++i) + { + for (int j = 0; j < col_size; ++j) + { + double re = 0.3 * std::sin(double(i * 3 + j * 5 + iap * 7 + ir * 11 + seed)); + double im = 0.3 * std::cos(double(i * 5 + j * 3 + iap * 11 + ir * 7 + seed)); + mat[i * col_size + j] = mk_value(re, im); + // the R=(0,0,0) self block must be symmetric (real) / + // Hermitian (complex): its reverse partner is itself + const bool force_hermitian = (self_pair && R.x == 0 && R.y == 0 && R.z == 0); + if (symmetric_in_R || force_hermitian) + { + mat[j * col_size + i] = conj_val(mat[i * col_size + j]); + } + } + } + } + } + // insert reverse blocks O(iat2,iat1,-R) = O(iat1,iat2,R)^T (real) or ^\dagger + for (const auto& pr: pairs) + { + for (const auto& R: r_list) + { + // the R=(0,0,0) self block is its own reverse partner + if (pr.first == pr.second && R.x == 0 && R.y == 0 && R.z == 0) + { + continue; + } + const hamilt::BaseMatrix* fwd = hR.find_matrix(pr.first, pr.second, R); + EXPECT_NE(fwd, nullptr) << "forward block missing for (" << pr.first << "," << pr.second << "," << R.x << "," << R.y << "," << R.z << ")"; + if (fwd == nullptr) + { + continue; + } + const int row_size = fwd->get_row_size(); + const int col_size = fwd->get_col_size(); + hamilt::AtomPair ap_rev(pr.second, pr.first, -R, paraV); + hamilt::BaseMatrix& rev_mat = ap_rev.get_HR_values(-R.x, -R.y, -R.z); + rev_mat.allocate(nullptr, true); + T* rev = rev_mat.get_pointer(); + for (int i = 0; i < row_size; ++i) + { + for (int j = 0; j < col_size; ++j) + { + const T fwd_val = fwd->get_pointer()[i * col_size + j]; + rev[j * col_size + i] = conj_val(fwd_val); + } + } + hR.insert_pair(ap_rev); + } + } + hR.allocate(nullptr, false); + return hR; +} + +// W_R = sum over container (pairs, R) of sum_{mu,nu} D(mu,nu;R) * O(mu,nu;R) +template +std::complex compute_WR(const hamilt::HContainer& dmr, const hamilt::HContainer& hr) +{ + std::complex wr(0.0, 0.0); + for (int iap = 0; iap < dmr.size_atom_pairs(); ++iap) + { + const hamilt::AtomPair& ap = dmr.get_atom_pair(iap); + const int iat1 = ap.get_atom_i(); + const int iat2 = ap.get_atom_j(); + const int row_size = ap.get_row_size(); + const int col_size = ap.get_col_size(); + for (int ir = 0; ir < ap.get_R_size(); ++ir) + { + const ModuleBase::Vector3 R = ap.get_R_index(ir); + const hamilt::BaseMatrix* dmat = ap.find_matrix(R); + const hamilt::BaseMatrix* omat = hr.find_matrix(iat1, iat2, R); + EXPECT_NE(omat, nullptr) << "O missing for (" << iat1 << "," << iat2 << "," << R.x << "," << R.y << "," << R.z << ")"; + if (omat == nullptr) + { + continue; + } + for (int i = 0; i < row_size; ++i) + { + for (int j = 0; j < col_size; ++j) + { + wr += dmat->get_pointer()[i * col_size + j] * omat->get_pointer()[i * col_size + j]; + } + } + } + } + return wr; +} + +// W_k = sum_k Re Tr( DMK(k) * folding_HR(O(R), k) ) +template +double compute_Wk(const std::vector>>& M, + const int nk, + const hamilt::HContainer& hr, + const std::vector>& kvec_d, + const int ntot) +{ + double wk = 0.0; + for (int ik = 0; ik < nk; ++ik) + { + std::vector> ok(ntot * ntot, std::complex(0.0, 0.0)); + hamilt::folding_HR(hr, ok.data(), kvec_d[ik], ntot, 0); + std::complex tr(0.0, 0.0); + for (int mu = 0; mu < ntot; ++mu) + { + for (int nu = 0; nu < ntot; ++nu) + { + tr += M[ik][mu * ntot + nu] * ok[nu * ntot + mu]; + } + } + // Tr(DMK(k) O(k)) is real for Hermitian DMK(k) and Hermitian O(k) + EXPECT_NEAR(tr.imag(), 0.0, 1e-10) << "trace not real at ik=" << ik; + wk += tr.real(); + } + return wk; +} +} // namespace + +class DMTraceTest : public ::testing::Test +{ + protected: + enum : int + { + nat = 2, + nw = 3, // orbitals per atom + ntot = nat * nw + }; + + Parallel_Orbitals* paraV; + + void SetUp() override + { + paraV = new Parallel_Orbitals(); + const int iat2iwt[nat] = {0, nw}; + paraV->init(ntot, ntot, 2, MPI_COMM_WORLD); + paraV->set_atomic_trace(iat2iwt, nat, ntot); + } + + void TearDown() override + { + delete paraV; + } +}; + +TEST_F(DMTraceTest, T3_closed_trace_equivalence) +{ + const int nk = 3; + const std::vector> kvec_d = { + ModuleBase::Vector3(0.1, 0.0, 0.0), + ModuleBase::Vector3(0.2, 0.2, 0.0), + ModuleBase::Vector3(0.3, 0.1, 0.2)}; + + // Hermitian DMK per k + std::vector>> M(nk); + for (int ik = 0; ik < nk; ++ik) + { + fill_hermitian(M[ik], ntot, ik); + } + + elecstate::DensityMatrix, double> DM(paraV, 1, kvec_d, nk); + for (int ik = 0; ik < nk; ++ik) + { + for (int i = 0; i < ntot; ++i) + { + for (int j = 0; j < ntot; ++j) + { + // set_DMK stores the transposed element (see density_matrix.h) + DM.set_DMK(1, ik, j, i, M[ik][i * ntot + j]); + } + } + } + + // paired DMR container: cross pair (0,1) and self pairs + const std::vector> pairs = {{0, 1}, {0, 0}, {1, 1}}; + const std::vector> r_list = { + ModuleBase::Vector3(0, 0, 0), + ModuleBase::Vector3(1, 0, 0), + ModuleBase::Vector3(0, 1, 0), + ModuleBase::Vector3(0, 0, 1)}; + + hamilt::HContainer dmr = build_r_container(paraV, pairs, r_list, true); + DM.init_DMR(dmr); + DM.cal_DMR(); + + // case 1: real symmetric O(R) (overlap-like) + { + auto hr = build_operator_container(paraV, pairs, r_list, 1, true); + const std::complex wr = compute_WR(dmr, hr); + EXPECT_NEAR(wr.imag(), 0.0, 1e-10) << "case 1: W_R should be real"; + const double wk = compute_Wk(M, nk, hr, kvec_d, ntot); + EXPECT_NEAR(wr.real(), wk, 1e-10) << "case 1: symmetric real O(R)"; + } + // case 2: non-symmetric real O(R) within the same R (derivative-like). + // The equivalence must still hold: protection comes from the closed trace + // + per-k Hermiticity of O(k), NOT from same-R symmetry. + { + auto hr = build_operator_container(paraV, pairs, r_list, 2, false); + const std::complex wr = compute_WR(dmr, hr); + EXPECT_NEAR(wr.imag(), 0.0, 1e-10) << "case 2: W_R should be real"; + const double wk = compute_Wk(M, nk, hr, kvec_d, ntot); + EXPECT_NEAR(wr.real(), wk, 1e-10) << "case 2: non-symmetric real O(R)"; + } + // case 3: complex Hermitian O(R) (multi-k) + { + auto hr = build_operator_container>(paraV, pairs, r_list, 3, false); + const std::complex wr = compute_WR(dmr, hr); + EXPECT_NEAR(wr.imag(), 0.0, 1e-10) << "case 3: W_R should be real"; + const double wk = compute_Wk(M, nk, hr, kvec_d, ntot); + EXPECT_NEAR(wr.real(), wk, 1e-10) << "case 3: complex Hermitian O(R)"; + } +} + +// Stub required by the DensityMatrix member functions instantiated from +// density_matrix_io.cpp (Record_adj::cal_adj is not linked into this +// unit-test target; see also tmp_mocks.cpp). +#include "source_lcao/record_adj.h" +Record_adj::Record_adj() {} +Record_adj::~Record_adj() {} + +int main(int argc, char** argv) +{ +#ifdef __MPI + MPI_Init(&argc, &argv); +#endif + testing::InitGoogleTest(&argc, argv); + int result = RUN_ALL_TESTS(); +#ifdef __MPI + MPI_Finalize(); +#endif + return result; +} diff --git a/tests/05_rtTDDFT/18_NO_hyb_TDDFT/check_extra.py b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/check_extra.py new file mode 100644 index 00000000000..72f18e9266e --- /dev/null +++ b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/check_extra.py @@ -0,0 +1,205 @@ +#!/usr/bin/env python3 +"""T5: rt-TDDFT current R-pairing / phase checks (multi-k, hybrid gauge Si). + +Appended to result.out / result.ref by the Autotest.sh check_extra.py hook. +Each key is a pass/fail flag (0.0 = pass, 1.0 = fail); the measured values are +written to stderr for diagnostics. All sub-runs use a 4x1x1 Gamma-centred mesh, +which contains the non-time-reversal-invariant pair k=(1/4,0,0) / -k=(3/4,0,0). + +Checks +------ +1. zero_field_current_max: max |j(t)| with td_vext=0. Si is time-reversal + symmetric, so the total current must vanish to ~1e-8 (the +/-k pairs cancel + at the level of the total). A nonzero value indicates a broken R pairing / + phase convention in the DMR -> current contraction. +2. current_path_max_dev: max |j_full - sum_k j_k| over steps/directions with a + finite field (td_vext=1). The full path contracts the all-k DMR + (cal_DMR()); the per-k path contracts each cal_DMR(ik) and sums the + k-weighted contributions. The two constructions must agree to ~1e-10. +3. current_pmk_cancel_max: max |j(k) + j(-k)| over the +/-k pair and steps. + Each per-k current is not physical alone (single-k DMR is not + Hermitian-paired), but time reversal gives j(-k) = -j(k); the residual + must be ~1e-8. + +NOTE on sensitivity: the current is a closed Frobenius trace with a full R sum, +so it is protected against the e^{-ikR}/e^{+ikR} DMR phase choice (flipping the +sign leaves all keys at roundoff). These checks guard against gross +pairing/phase breakage (e.g. a mismatched +phase_hybrid in the hybrid gauge). The DMR sign sentinel is +DMTest.T1_fourier_round_trip. + +Usage: check_extra.py [abacus_bin] [np] [case_dir] +""" + +import os +import re +import shutil +import subprocess +import sys +import tempfile + +# ABACUS integration-test convention: OMP_NUM_THREADS=1 keeps the SCF and the +# two DMR current paths bit-deterministic (multi-threaded summation order makes +# the full-k and per-k paths disagree at ~1e-8 instead of roundoff). +os.environ.setdefault("OMP_NUM_THREADS", "1") + +KPT_TEXT = "K_POINTS\n0\nGamma\n4 1 1 0 0 0\n" +TOL_PATH = 1e-10 +TOL_ZERO = 1e-8 +TOL_PMK = 1e-8 + + +def fail(msg): + print("current_path_max_dev 1.0", flush=True) + print("zero_field_current_max 1.0", flush=True) + print("current_pmk_cancel_max 1.0", flush=True) + sys.stderr.write("check_extra.py: " + msg + "\n") + sys.exit(1) + + +def get_input_value(text, key, default=None): + for line in text.splitlines(): + line = line.strip() + if line.startswith("#"): + continue + parts = line.split() + if parts and parts[0] == key: + return parts[1] + return default + + +def set_input_value(text, key, value): + if re.search(rf"^\s*{key}\s", text, flags=re.M): + return re.sub(rf"^\s*{key}\s+\S+", f"{key} {value}", text, flags=re.M) + return text + f"\n{key} {value}\n" + + +def read_current_rows(path): + rows = [] + for line in open(path): + parts = line.split() + if len(parts) == 4: + rows.append([float(parts[1]), float(parts[2]), float(parts[3])]) + return rows + + +def read_kpoints(running_log): + """Parse the K-point list (direct coords) from KPT.info in the log dir.""" + ks = [] + in_list = False + for line in open(running_log): + if "K-POINTS DIRECT COORDINATES" in line: + in_list = True + continue + if in_list: + parts = line.split() + if len(parts) == 5 and parts[0].isdigit(): + ks.append([float(parts[1]), float(parts[2]), float(parts[3])]) + elif ks and len(parts) != 5: + break + return ks + + +def run_case(case_dir, workdir, abacus, np_, mods, kpt_text=KPT_TEXT): + """Copy the case into workdir, apply INPUT modifications + KPT, run ABACUS, + and return (current_tot rows, {ik: rows} for per-k currents, k vectors).""" + os.makedirs(workdir, exist_ok=True) + inp = open(os.path.join(case_dir, "INPUT")).read() + inp = inp.replace("../../PP_ORB", os.path.abspath(os.path.join(case_dir, "../../PP_ORB"))) + for key, value in mods.items(): + inp = set_input_value(inp, key, value) + with open(os.path.join(workdir, "INPUT"), "w") as f: + f.write(inp) + for name in ("STRU",): + shutil.copy(os.path.join(case_dir, name), os.path.join(workdir, name)) + with open(os.path.join(workdir, "KPT"), "w") as f: + f.write(kpt_text) + with open(os.path.join(workdir, "run.log"), "w") as f: + ret = subprocess.run(["mpirun", "-np", str(np_), abacus], cwd=workdir, + stdout=f, stderr=subprocess.STDOUT) + if ret.returncode != 0: + fail("ABACUS run failed in " + workdir) + out = os.path.join(workdir, "OUT.autotest") + tot = os.path.join(out, "current_tot.txt") + if not os.path.isfile(tot): + fail("current_tot.txt missing in " + workdir) + tot_rows = read_current_rows(tot) + if not tot_rows: + fail("current_tot.txt empty in " + workdir) + perk = {} + for name in sorted(os.listdir(out)): + m = re.fullmatch(r"current_s(\d+)k(\d+)\.txt", name) + if m: + ik = int(m.group(2)) + perk[ik] = read_current_rows(os.path.join(out, name)) + ks = read_kpoints(os.path.join(out, "KPT.info")) + if not ks: + fail("could not parse K-point list from KPT.info in " + workdir) + return tot_rows, perk, ks + + +def max_diff(a_rows, b_rows): + if len(a_rows) != len(b_rows): + fail("step count mismatch between current paths") + return max(abs(x - y) + for ra, rb in zip(a_rows, b_rows) + for x, y in zip(ra, rb)) + + +def main(): + abacus = sys.argv[1] if len(sys.argv) > 1 else "abacus" + np_ = int(sys.argv[2]) if len(sys.argv) > 2 else 4 + case_dir = os.path.abspath(sys.argv[3]) if len(sys.argv) > 3 else os.getcwd() + + with tempfile.TemporaryDirectory(prefix="t5_current_") as tmp: + # finite field: nonzero current; compare full-k DMR path vs per-k DMR path + j_full, _, _ = run_case(case_dir, os.path.join(tmp, "full"), abacus, np_, + {"out_current": "1", "out_current_k": "0", "td_vext": "1"}) + j_eachk, perk, ks = run_case(case_dir, os.path.join(tmp, "eachk"), abacus, np_, + {"out_current": "1", "out_current_k": "1", "td_vext": "1"}) + path_dev = max_diff(j_full, j_eachk) + + # sum of per-k files must reproduce the per-k-path total + sum_dev = 0.0 + if perk: + nstep = len(next(iter(perk.values()))) + acc = [[0.0, 0.0, 0.0] for _ in range(nstep)] + for rows in perk.values(): + for step, row in enumerate(rows): + if step < nstep: + for d in range(3): + acc[step][d] += row[d] + for step, row in enumerate(acc): + if step < len(j_eachk): + sum_dev = max(sum_dev, max(abs(row[d] - j_eachk[step][d]) for d in range(3))) + if sum_dev > TOL_PATH: + fail("per-k file sum deviates from current_tot (per-k path): %.3e" % sum_dev) + + # +/-k cancellation on the per-k currents + pmk_dev = 0.0 + nk = len(ks) + for i in range(nk): + for j in range(i + 1, nk): + if sum(abs(ks[i][d] + ks[j][d]) for d in range(3)) < 1e-10: + ri = perk.get(i + 1) + rj = perk.get(j + 1) + if ri is None or rj is None: + fail("missing per-k current file for k pair %d/%d" % (i + 1, j + 1)) + pmk_dev = max(pmk_dev, max(abs(a + b) + for ra, rb in zip(ri, rj) + for a, b in zip(ra, rb))) + + # zero external field: time-reversal symmetric ground state + j_zero, _, _ = run_case(case_dir, os.path.join(tmp, "zero"), abacus, np_, + {"out_current": "1", "out_current_k": "0", "td_vext": "0"}) + zero_dev = max(abs(v) for row in j_zero for v in row) + + sys.stderr.write("check_extra.py: path_dev=%.3e sum_dev=%.3e pmk_dev=%.3e zero_dev=%.3e\n" + % (path_dev, sum_dev, pmk_dev, zero_dev)) + print("current_path_max_dev %.1f" % (1.0 if path_dev > TOL_PATH else 0.0), flush=True) + print("zero_field_current_max %.1f" % (1.0 if zero_dev > TOL_ZERO else 0.0), flush=True) + print("current_pmk_cancel_max %.1f" % (1.0 if pmk_dev > TOL_PMK else 0.0), flush=True) + + +if __name__ == "__main__": + main() diff --git a/tests/05_rtTDDFT/18_NO_hyb_TDDFT/result.ref b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/result.ref index fe01bf77eb2..b994db43b73 100644 --- a/tests/05_rtTDDFT/18_NO_hyb_TDDFT/result.ref +++ b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/result.ref @@ -1,4 +1,7 @@ -etotref -202.3030324186811 -etotperatomref -101.1515162093 +etotref -202.3030324194245 +etotperatomref -101.1515162097 CompareCurrent_pass 0 -totaltimeref 5.26 +totaltimeref 5.18 +current_path_max_dev 0.0 +zero_field_current_max 0.0 +current_pmk_cancel_max 0.0 diff --git a/tests/integrate/240_NO_KP_15_SO_FD/INPUT b/tests/integrate/240_NO_KP_15_SO_FD/INPUT new file mode 100644 index 00000000000..411632fa4e3 --- /dev/null +++ b/tests/integrate/240_NO_KP_15_SO_FD/INPUT @@ -0,0 +1,23 @@ +INPUT_PARAMETERS +#Parameters (General) +suffix autotest +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB + +#Parameters (Electronic) +calculation scf +basis_type lcao +gamma_only 0 +nspin 4 +lspinorb 1 +ecutwfc 100 +scf_thr 1e-8 +scf_nmax 100 + +ks_solver scalapack_gvx + +mixing_type pulay +mixing_beta 0.7 + +cal_force 1 +cal_stress 0 diff --git a/tests/integrate/240_NO_KP_15_SO_FD/KPT b/tests/integrate/240_NO_KP_15_SO_FD/KPT new file mode 100644 index 00000000000..28006d5e2df --- /dev/null +++ b/tests/integrate/240_NO_KP_15_SO_FD/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +2 2 2 0 0 0 diff --git a/tests/integrate/240_NO_KP_15_SO_FD/README b/tests/integrate/240_NO_KP_15_SO_FD/README new file mode 100644 index 00000000000..0adec89c411 --- /dev/null +++ b/tests/integrate/240_NO_KP_15_SO_FD/README @@ -0,0 +1,8 @@ +T4: SOC nonlocal force finite-difference validation (GaAs, nspin=4, SOC, +multi-k). check_extra.py compares the analytic force with a central FD of the +total energy (delta=1e-3 Bohr) for every atom/direction, and the acoustic sum. +Keys: fd_force_pass, acoustic_sum_pass (0.0 = pass, 1.0 = fail; physical +threshold 1e-4 eV/Bohr; measured max_dev ~5.8e-5). The DMR Fourier sign +sentinel is the unit test DMTest.T1_fourier_round_trip; this case guards +FD-consistency of the analytic SOC force. ecutwfc=100 is required: with a +coarse charge grid the egg-box noise dominates the FD. diff --git a/tests/integrate/240_NO_KP_15_SO_FD/STRU b/tests/integrate/240_NO_KP_15_SO_FD/STRU new file mode 100644 index 00000000000..519e33b46a1 --- /dev/null +++ b/tests/integrate/240_NO_KP_15_SO_FD/STRU @@ -0,0 +1,28 @@ +ATOMIC_SPECIES +As 1 As_ONCV_PBE_FR-1.1.upf upf201 +Ga 1 Ga_ONCV_PBE_FR-1.0.upf.txt upf201 + +LATTICE_CONSTANT +1 + +NUMERICAL_ORBITAL +As_gga_8au_60Ry_2s2p1d.orb +Ga_gga_9au_60Ry_2s2p2d.orb + +LATTICE_VECTORS +5.34197 5.34197 0.0 +0.0 5.34197 5.34197 +5.34197 0.0 5.34197 + +ATOMIC_POSITIONS +Direct + +As +0 +1 +0.2600000 0.2400000 0.245000000 0 0 0 + +Ga +0 +1 +0.0300000 -0.0200000 0.050000000 0 0 0 diff --git a/tests/integrate/240_NO_KP_15_SO_FD/check_extra.py b/tests/integrate/240_NO_KP_15_SO_FD/check_extra.py new file mode 100644 index 00000000000..623f74af1b7 --- /dev/null +++ b/tests/integrate/240_NO_KP_15_SO_FD/check_extra.py @@ -0,0 +1,362 @@ +#!/usr/bin/env python3 +"""T4: SOC nonlocal force finite-difference validation. + +Compares the analytic LCAO force against a central finite difference of the +total energy for every atom and Cartesian direction: + + F_fd = -(E(+delta) - E(-delta)) / (2*delta), delta = 1e-3 Bohr + +and reports two pass/fail keys for result.out (0.0 = pass, 1.0 = fail): + + fd_force_pass max |F_analytic - F_fd| over all atoms/directions < 1e-4 eV/Bohr + acoustic_sum_pass max |sum_i F_analytic(i,d)| over directions d < 1e-4 eV/Bohr + +The acoustic sum (translational invariance) is the first quantity that breaks +if the R pairing of the nonlocal force is wrong; the per-component comparison +prevents accidental cancellations between force terms. The measured max_dev is +written to stderr for diagnostics (reference ~5.8e-5 eV/Bohr on this case, +deterministic at scf_thr=1e-8); comparing raw FD metrics run-to-run at 1e-6 +would be machine-fragile, so the harness comparison uses the pass/fail flags. + +NOTE on sensitivity: this non-magnetic GaAs SOC case exercises the nonlocal +SOC force path, but its force is (numerically) identical under the +e^{-ikR}/e^{+ikR} DMR phase choice (closed-trace protection). The sentinel +for the DMR Fourier sign itself is the unit test DMTest.T1_fourier_round_trip +(source/source_estate/module_dm/test/test_cal_dm_R.cpp), which turns red when +the sign is flipped. This case guards FD-consistency of the analytic force, +including the SOC nonlocal contribution. + +Usage: check_extra.py [abacus_bin] [np] [case_dir] +Keys are printed to stdout so the integrate harness (Autotest.sh) appends them +to result.out / result.ref. +""" + +import math +import os +import re +import shutil +import subprocess +import sys +import tempfile + +BOHR_TO_ANG = 0.529177210903 + + +def fail(msg): + print("fd_force_pass 1.0", flush=True) + print("acoustic_sum_pass 1.0", flush=True) + sys.stderr.write("fd_force.py: " + msg + "\n") + sys.exit(1) + + +def get_input_value(text, key, default=None): + for line in text.splitlines(): + line = line.strip() + if line.startswith("#"): + continue + parts = line.split() + if parts and parts[0] == key: + return parts[1] + return default + + +def parse_stru(text): + """Return (lattice_matrix, atom_list) from an ABACUS STRU file. + + lattice_matrix: 3x3 rows in Bohr. + atom_list: list of dicts {label, idx, frac (3 floats), extra (rest of line)}. + """ + lines = text.splitlines() + i = 0 + n = len(lines) + + def next_data(): + nonlocal i + while i < n: + ln = lines[i].strip() + i += 1 + if ln and not ln.startswith("#"): + return ln + return None + + # ATOMIC_SPECIES block + while True: + ln = next_data() + if ln is None: + raise RuntimeError("ATOMIC_SPECIES not found") + if ln.startswith("ATOMIC_SPECIES"): + break + lat_const = 1.0 + while True: + ln = next_data() + if ln is None: + raise RuntimeError("LATTICE_CONSTANT not found") + if ln.startswith("LATTICE_CONSTANT"): + lat_const = float(next_data()) + break + + # skip NUMERICAL_ORBITAL if present + while True: + ln = next_data() + if ln is None: + raise RuntimeError("LATTICE_VECTORS not found") + if ln.startswith("LATTICE_VECTORS"): + break + latvec = [] + for _ in range(3): + latvec.append([float(x) for x in next_data().split()[:3]]) + lat = [[lat_const * v for v in row] for row in latvec] + + while True: + ln = next_data() + if ln is None: + raise RuntimeError("ATOMIC_POSITIONS not found") + if ln.startswith("ATOMIC_POSITIONS"): + break + coord_type = next_data() + + atoms = [] + while True: + ln = next_data() + if ln is None: + break + if ln.startswith("#") or ln == "": + continue + label = ln.strip() + if next_data() is None: + break # magnetism line + nat_line = next_data() + if nat_line is None: + break + nat = int(nat_line) + for ia in range(nat): + pos = next_data() + if pos is None: + raise RuntimeError("unexpected end of STRU atom block") + fields = pos.split() + frac = [float(fields[0]), float(fields[1]), float(fields[2])] + atoms.append({"label": label, "idx": ia, "frac": frac, "extra": fields[3:]}) + return lat, atoms + + +def rewrite_stru(text, lat, atoms, delta_vec_bohr, target_idx): + """Return STRU text with atom positions shifted by delta (Cartesian, Bohr).""" + inv = invert3(lat) + # fraction-coordinate shift dfrac = delta . A^{-1} (row-vector convention: + # Cartesian delta = dfrac . A). A is not symmetric for non-orthogonal + # lattices, so the index order below matters. + dfrac = [0.0, 0.0, 0.0] + for a in range(3): + for b in range(3): + dfrac[a] += delta_vec_bohr[b] * inv[b][a] + + lines = text.splitlines() + out = [] + i = 0 + n = len(lines) + atom_iter = iter(atoms) + shifted = {} + for idx, atom in enumerate(atoms): + if idx == target_idx: + shifted[id(atom)] = [atom["frac"][a] + dfrac[a] for a in range(3)] + else: + shifted[id(atom)] = list(atom["frac"]) + + # copy everything up to and including the ATOMIC_POSITIONS header line + while i < n: + out.append(lines[i]) + if lines[i].strip().startswith("ATOMIC_POSITIONS"): + i += 1 + if i < n: + out.append(lines[i]) # Direct/Cartesian + i += 1 + break + i += 1 + + # element blocks: label, magnetism, nat, nat coordinate lines + while i < n: + ln = lines[i] + stripped = ln.strip() + if stripped == "" or stripped.startswith("#"): + out.append(ln) + i += 1 + continue + label = stripped + out.append(ln) + i += 1 + # magnetism line + while i < n and (lines[i].strip() == "" or lines[i].strip().startswith("#")): + out.append(lines[i]) + i += 1 + if i >= n: + break + out.append(lines[i]) + i += 1 + # atom count line + while i < n and (lines[i].strip() == "" or lines[i].strip().startswith("#")): + out.append(lines[i]) + i += 1 + if i >= n: + break + out.append(lines[i]) + nat = int(lines[i].strip()) + i += 1 + for ia in range(nat): + while i < n and (lines[i].strip() == "" or lines[i].strip().startswith("#")): + out.append(lines[i]) + i += 1 + if i >= n: + break + pos = lines[i].split() + target = None + for atom in atoms: + if atom["label"] == label and atom["idx"] == ia: + target = atom + break + if target is None: + raise RuntimeError("atom not found: %s %d" % (label, ia)) + frac = shifted[id(target)] + fields = ["%14.9f" % v for v in frac] + fields += pos[3:] + out.append(" ".join(fields)) + i += 1 + return "\n".join(out) + + +def invert3(m): + a, b, c = m[0][0], m[0][1], m[0][2] + d, e, f = m[1][0], m[1][1], m[1][2] + g, h, k = m[2][0], m[2][1], m[2][2] + det = a * (e * k - f * h) - b * (d * k - f * g) + c * (d * h - e * g) + if abs(det) < 1e-30: + raise RuntimeError("singular lattice matrix") + inv = [ + [(e * k - f * h) / det, (c * h - b * k) / det, (b * f - c * e) / det], + [(f * g - d * k) / det, (a * k - c * g) / det, (c * d - a * f) / det], + [(d * h - e * g) / det, (b * g - a * h) / det, (a * e - b * d) / det], + ] + return inv + + +def run_scf(workdir, abacus, np, tag): + """Run ABACUS SCF in workdir, return final total energy in eV.""" + log = os.path.join(workdir, "run_%s.log" % tag) + with open(log, "w") as f: + subprocess.run(["mpirun", "-np", str(np), abacus], + cwd=workdir, stdout=f, stderr=subprocess.STDOUT) + # find the OUT. directory + inp = open(os.path.join(workdir, "INPUT")).read() + suffix = get_input_value(inp, "suffix", "autotest") + outdir = os.path.join(workdir, "OUT." + suffix) + scf = os.path.join(outdir, "running_scf.log") + if not os.path.isfile(scf): + raise RuntimeError("missing " + scf) + text = open(scf).read() + m = re.search(r"!FINAL_ETOT_IS\s+([-0-9.eE+]+)", text) + if not m: + raise RuntimeError("FINAL_ETOT_IS not found in " + scf) + return float(m.group(1)) + + +def read_analytic_force(case_dir, suffix): + scf = os.path.join(case_dir, "OUT." + suffix, "running_scf.log") + if not os.path.isfile(scf): + return None + text = open(scf).read() + m = re.search(r"#TOTAL-FORCE.*?\n(.*?)(?:\n#|$)", text, re.S) + if not m: + return None + block = m.group(1) + forces = [] + for line in block.splitlines(): + parts = line.split() + if len(parts) == 4 and re.match(r"^[A-Za-z]+\d+$", parts[0]): + forces.append([float(parts[1]), float(parts[2]), float(parts[3])]) + return forces + + +def main(): + abacus = sys.argv[1] if len(sys.argv) > 1 else "abacus" + np_ = int(sys.argv[2]) if len(sys.argv) > 2 else 4 + case_dir = os.path.abspath(sys.argv[3]) if len(sys.argv) > 3 else os.getcwd() + delta = 1.0e-3 # Bohr + + inp_path = os.path.join(case_dir, "INPUT") + stru_path = os.path.join(case_dir, "STRU") + if not (os.path.isfile(inp_path) and os.path.isfile(stru_path)): + fail("INPUT/STRU missing in " + case_dir) + + inp_text = open(inp_path).read() + suffix = get_input_value(inp_text, "suffix", "autotest") + pseudo_dir = get_input_value(inp_text, "pseudo_dir", "../../PP_ORB") + orbital_dir = get_input_value(inp_text, "orbital_dir", "../../PP_ORB") + pseudo_abs = os.path.abspath(os.path.join(case_dir, pseudo_dir)) + orbital_abs = os.path.abspath(os.path.join(case_dir, orbital_dir)) + + stru_text = open(stru_path).read() + lat, atoms = parse_stru(stru_text) + natom = len(atoms) + if natom == 0: + fail("no atoms in STRU") + + # analytic force from the base run (already produced by the harness) or run + # a fresh base SCF if missing. + base_forces = read_analytic_force(case_dir, suffix) + if base_forces is None: + shutil.copy(inp_path, case_dir + ".INPUT.bak") + try: + shutil.copy(stru_path, case_dir + ".STRU.bak") + run_scf(case_dir, abacus, np_, "base") + base_forces = read_analytic_force(case_dir, suffix) + finally: + shutil.move(case_dir + ".INPUT.bak", inp_path) + shutil.move(case_dir + ".STRU.bak", stru_path) + if base_forces is None or len(base_forces) != natom: + fail("could not read analytic forces (%d atoms)" % natom) + + fd_force = [[0.0, 0.0, 0.0] for _ in range(natom)] + with tempfile.TemporaryDirectory(prefix="t4_fd_") as tmp: + for ia in range(natom): + for d in range(3): + shift = [0.0, 0.0, 0.0] + shift[d] = delta + energies = [] + for sgn in (1.0, -1.0): + wd = os.path.join(tmp, "a%d_d%d_%s" % (ia, d, "+" if sgn > 0 else "-")) + os.makedirs(wd) + inp_mod = re.sub(r"^\s*pseudo_dir.*$", "pseudo_dir " + pseudo_abs, + inp_text, flags=re.M) + inp_mod = re.sub(r"^\s*orbital_dir.*$", "orbital_dir " + orbital_abs, + inp_mod, flags=re.M) + stru_mod = rewrite_stru(stru_text, lat, atoms, + [sgn * v for v in shift], ia) + with open(os.path.join(wd, "INPUT"), "w") as f: + f.write(inp_mod) + with open(os.path.join(wd, "STRU"), "w") as f: + f.write(stru_mod) + shutil.copy(os.path.join(case_dir, "KPT"), os.path.join(wd, "KPT")) + energies.append(run_scf(wd, abacus, np_, "fd")) + # F = -dE/dx: FD force = -(E(+d)-E(-d)) / (2*d) + fd_force[ia][d] = (energies[1] - energies[0]) / (2.0 * delta) + + # convert analytic force eV/Angstrom -> eV/Bohr (1 Bohr = 0.529177 Angstrom, + # so dE/dx[Bohr] = dE/dx[Ang] * 0.529177) + max_dev = 0.0 + acoustic = [0.0, 0.0, 0.0] + for ia in range(natom): + for d in range(3): + fa = base_forces[ia][d] * BOHR_TO_ANG + dev = abs(fa - fd_force[ia][d]) + max_dev = max(max_dev, dev) + acoustic[d] += fa + acoustic_max = max(abs(v) for v in acoustic) + + print("fd_force_pass %.1f" % (1.0 if max_dev > 1e-4 else 0.0), flush=True) + print("acoustic_sum_pass %.1f" % (1.0 if acoustic_max > 1e-4 else 0.0), flush=True) + sys.stderr.write( + "fd_force.py: max_dev=%.3e eV/Bohr acoustic=%.3e eV/Bohr (threshold 1e-4)\n" + % (max_dev, acoustic_max)) + + +if __name__ == "__main__": + main() diff --git a/tests/integrate/240_NO_KP_15_SO_FD/result.ref b/tests/integrate/240_NO_KP_15_SO_FD/result.ref new file mode 100644 index 00000000000..0e6ba7a1a88 --- /dev/null +++ b/tests/integrate/240_NO_KP_15_SO_FD/result.ref @@ -0,0 +1,6 @@ +etotref -1955.770906167134 +etotperatomref -977.8854530836 +totalforceref 10.072076 +totaltimeref 107.04 +fd_force_pass 0.0 +acoustic_sum_pass 0.0 diff --git a/tests/integrate/240_NO_KP_15_SO_FD/threshold b/tests/integrate/240_NO_KP_15_SO_FD/threshold new file mode 100644 index 00000000000..8fcfcebb260 --- /dev/null +++ b/tests/integrate/240_NO_KP_15_SO_FD/threshold @@ -0,0 +1 @@ +threshold 0.0001 diff --git a/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/INPUT b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/INPUT new file mode 100644 index 00000000000..66e00a5767b --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/INPUT @@ -0,0 +1,38 @@ +INPUT_PARAMETERS +suffix autotest +nbands 40 + +calculation scf +ecutwfc 10 +scf_thr 1.0e-4 +scf_nmax 500 +out_chg 0 + +smearing_method gaussian +smearing_sigma 0.01 + +cal_force 1 +cal_stress 1 + +mixing_type broyden +mixing_beta 0.2 +mixing_restart 5e-3 +mixing_dmr 1 +mixing_ndim 15 + +ks_solver scalapack_gvx +basis_type lcao +gamma_only 0 +symmetry 0 +nspin 2 + +#Parameter DFT+U +dft_plus_u 1 +orbital_corr 2 -1 +hubbard_u 5.0 0.0 +uramping 2.5 +onsite_radius 5.0 +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB + +out_mul 1 diff --git a/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/KPT b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/KPT new file mode 100644 index 00000000000..e769af76382 --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +2 1 1 0 0 0 diff --git a/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/README b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/README new file mode 100644 index 00000000000..b67e80fe522 --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/README @@ -0,0 +1,9 @@ +T7: nspin=2 collinear AFM spin-channel check (LCAO, DFT+U, FeO, uramping). +out_mul=1 produces OUT.autotest/mulliken.txt; check_extra.py asserts the +spin-channel sum rule (spin1+spin2 == total charge), the per-atom charge sum +rule, opposite-sign and comparable-magnitude Fe site moments, small O moments, +and consistency of the Mulliken cell sum with the rho cell magnetization. +Keys: t7_converged_pass, t7_spin_sum_pass, t7_charge_sum_pass, +t7_fe_opposite_sign_pass, t7_fe_mag_ratio_pass, t7_o_small_pass, +t7_mul_rho_pass (0.0 = pass, 1.0 = fail). Checking the two spin channels +separately guards against accidental cancellation of channel errors. diff --git a/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/STRU b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/STRU new file mode 100644 index 00000000000..65d8f7c98f2 --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/STRU @@ -0,0 +1,29 @@ +ATOMIC_SPECIES +Fe 1.000 Fe.upf +O 1.000 O_ONCV_PBE-1.0.upf + +NUMERICAL_ORBITAL +Fe_gga_6au_100Ry_4s2p2d1f.orb +O_gga_7au_100Ry_2s2p1d.orb + +LATTICE_CONSTANT +8.190 + +LATTICE_VECTORS + 1.00 0.50 0.50 + 0.50 1.00 0.50 + 0.50 0.50 1.00 +ATOMIC_POSITIONS +Direct + +Fe +0.0 +2 +0.00 0.00 0.00 mag 1.0 +0.51 0.51 0.51 mag -1.0 + +O +0.0 +2 +0.25 0.25 0.25 1 1 1 +0.75 0.75 0.75 1 1 1 diff --git a/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/check_extra.py b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/check_extra.py new file mode 100644 index 00000000000..a1fef714fd2 --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/check_extra.py @@ -0,0 +1,151 @@ +#!/usr/bin/env python3 +"""T7: collinear (nspin=2) AFM spin-channel check (LCAO, DFT+U, FeO). + +The Mulliken output for nspin=2 reports per-atom charge and magnetization for +each spin channel separately (mulliken.txt "Total charge of spin1/spin2" plus +the per-atom "sum lmz" / "total magnetism on atom i"). AFM order must show up +in the channels individually: the two Fe sites must carry opposite-sign +magnetizations of comparable magnitude, the O sites must stay (almost) +non-magnetic, and the spin-channel charges must add up to the total charge. +Checking the channels separately (rather than only the total charge/force) +guards against an accidental cancellation of the two channel errors. + +Each key is a pass/fail flag: 0.0 = pass, 1.0 = fail. Measured values are +written to stderr for diagnostics. + +Usage: check_extra.py [abacus_bin] [np] [case_dir] +""" + +import os +import re +import sys + +KEYS = [ + "t7_converged_pass", + "t7_spin_sum_pass", + "t7_charge_sum_pass", + "t7_fe_opposite_sign_pass", + "t7_fe_mag_ratio_pass", + "t7_o_small_pass", + "t7_mul_rho_pass", +] + +TOL_CHARGE = 1e-3 # electrons +TOL_MUL_RHO = 0.05 # uB: sum(Mulliken per-atom m) vs rho cell magnetization +FE_MIN = 1.0 # uB, Fe sites must be non-vanishing +RATIO_LO, RATIO_HI = 0.3, 3.0 +O_MAX = 0.5 # uB, O sites stay (almost) non-magnetic + + +def fail(msg): + for k in KEYS: + print(k + " 1.0", flush=True) + sys.stderr.write("t7_check.py: " + msg + "\n") + sys.exit(1) + + +def parse_mulliken(path): + text = open(path).read() + total = None + spins = {} + atoms = {} # idx -> {"chg": float, "mag": float} + for line in text.splitlines(): + m = re.match(r"^\s*Total charge\s+([-+0-9.eE]+)\s*$", line) + if m: + total = float(m.group(1)) + continue + m = re.match(r"^\s*Total charge of spin(\d+)\s+([-+0-9.eE]+)", line) + if m: + spins[int(m.group(1))] = float(m.group(2)) + continue + m = re.match(r"^\s*total charge\s+on atom\s+(\d+)\s+([-+0-9.eE]+)", line) + if m: + atoms.setdefault(int(m.group(1)), {})["chg"] = float(m.group(2)) + continue + m = re.match( + r"^\s*total magnetism\s+on atom\s+(\d+)\s+([-+0-9.eE]+)" + r"(?:\s+([-+0-9.eE]+)\s+([-+0-9.eE]+))?", line) + if m: + vals = [float(m.group(i)) for i in (2, 3, 4) if m.group(i)] + atoms.setdefault(int(m.group(1)), {})["mag"] = vals + return total, spins, atoms + + +def parse_rho_mag(path): + text = open(path).read() + if not re.search(r"#SCF IS CONVERGED#", text): + return None, False + lines = [ln for ln in text.splitlines() if "Total magnetism (Bohr mag/cell)" in ln] + if not lines: + return None, True + # nspin=2: scalar cell magnetization "Total magnetism (Bohr mag/cell) = v" + m = re.search(r"=\s*([-+0-9.eE]+)\s*$", lines[-1]) + if m: + return float(m.group(1)), True + return None, True + + +def main(): + case_dir = os.path.abspath(sys.argv[3]) if len(sys.argv) > 3 else os.getcwd() + mul_path = os.path.join(case_dir, "OUT.autotest", "mulliken.txt") + scf_path = os.path.join(case_dir, "OUT.autotest", "running_scf.log") + if not (os.path.isfile(mul_path) and os.path.isfile(scf_path)): + fail("OUT.autotest/mulliken.txt or running_scf.log missing") + + total, spins, atoms = parse_mulliken(mul_path) + if total is None or not atoms: + fail("could not parse mulliken.txt") + rho_mag, converged = parse_rho_mag(scf_path) + + results = {} + results["t7_converged_pass"] = 0.0 if converged else 1.0 + + # 1. spin-channel sum rule: spin1 + spin2 == total charge + if len(spins) == 2: + results["t7_spin_sum_pass"] = 0.0 if abs(spins[1] + spins[2] - total) < TOL_CHARGE else 1.0 + else: + results["t7_spin_sum_pass"] = 1.0 + + # 2. per-atom charge sum rule + qsum = sum(a["chg"] for a in atoms.values()) + results["t7_charge_sum_pass"] = 0.0 if abs(qsum - total) < TOL_CHARGE else 1.0 + + # 3-6. channel-wise AFM physics (index order follows STRU: Fe, Fe, O, O) + mags = [] + for idx in sorted(atoms): + m = atoms[idx].get("mag") + if m is None or len(m) != 1: + fail("nspin=2 scalar magnetism expected for every atom") + mags.append(m[0]) + if len(mags) != 4: + fail("expected 4 atoms (2 Fe, 2 O)") + m_fe1, m_fe2, m_o1, m_o2 = mags + + results["t7_fe_opposite_sign_pass"] = ( + 0.0 if m_fe1 > FE_MIN and m_fe2 < -FE_MIN else 1.0) + ratio = abs(m_fe1) / abs(m_fe2) if abs(m_fe2) > 1e-8 else 1e6 + results["t7_fe_mag_ratio_pass"] = ( + 0.0 if RATIO_LO <= ratio <= RATIO_HI else 1.0) + results["t7_o_small_pass"] = ( + 0.0 if abs(m_o1) < O_MAX and abs(m_o2) < O_MAX else 1.0) + + # 7. Mulliken cell sum vs rho cell magnetization + if rho_mag is None: + results["t7_mul_rho_pass"] = 1.0 + else: + msum = sum(mags) + results["t7_mul_rho_pass"] = 0.0 if abs(msum - rho_mag) < TOL_MUL_RHO else 1.0 + + for k in KEYS: + print(k + " %.1f" % results[k], flush=True) + sys.stderr.write( + "t7_check.py: total=%.6f spin1=%.6f spin2=%.6f sum_atoms=%.6f | " + "m_Fe1=%.6f m_Fe2=%.6f m_O1=%.6f m_O2=%.6f | " + "sum_mulliken=%.6f rho_cell=%.6f converged=%s\n" + % (total, spins.get(1, float("nan")), spins.get(2, float("nan")), qsum, + m_fe1, m_fe2, m_o1, m_o2, sum(mags), + rho_mag if rho_mag is not None else float("nan"), bool(converged))) + + +if __name__ == "__main__": + main() diff --git a/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/result.ref b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/result.ref new file mode 100644 index 00000000000..f855772d278 --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/result.ref @@ -0,0 +1,11 @@ +etotref -7661.6770236327884049 +etotperatomref -1915.4192559082 +totalforceref 21.733614 +totalstressref 6973.087067 +t7_converged_pass 0.0 +t7_spin_sum_pass 0.0 +t7_charge_sum_pass 0.0 +t7_fe_opposite_sign_pass 0.0 +t7_fe_mag_ratio_pass 0.0 +t7_o_small_pass 0.0 +t7_mul_rho_pass 0.0 diff --git a/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/threshold b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/threshold new file mode 100644 index 00000000000..2bb8a636fea --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING_MUL/threshold @@ -0,0 +1,2 @@ +threshold 1e-6 +force_threshold 1e-4 diff --git a/tests/integrate/260_NO_DJ_PK_PU_SO_4K/INPUT b/tests/integrate/260_NO_DJ_PK_PU_SO_4K/INPUT new file mode 100644 index 00000000000..dde2fabee98 --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_SO_4K/INPUT @@ -0,0 +1,44 @@ +INPUT_PARAMETERS +suffix autotest +nbands 40 + +calculation scf +ecutwfc 10 +scf_thr 1.0e-4 +scf_nmax 50 +out_chg 0 + +#init_chg file +#out_dos 1 +#dos_sigma 0.05 +#out_band 1 + +smearing_method gaussian +smearing_sigma 0.01 + +#force_thr_ev 0.01 +#relax_method cg +#relax_bfgs_init 0.5 + +mixing_type pulay +mixing_beta 0.3 +mixing_restart 1e-3 +mixing_dmr 1 +mixing_gg0 1.1 + +ks_solver scalapack_gvx +basis_type lcao +gamma_only 0 +noncolin 1 +lspinorb 1 +nspin 4 +cal_force 1 +cal_stress 0 + +#Parameter DFT+U +dft_plus_u 1 +orbital_corr 2 +hubbard_u 5.0 +onsite_radius 5.0 +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB diff --git a/tests/integrate/260_NO_DJ_PK_PU_SO_4K/KPT b/tests/integrate/260_NO_DJ_PK_PU_SO_4K/KPT new file mode 100644 index 00000000000..b1f8afcf383 --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_SO_4K/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +4 1 1 0 0 0 diff --git a/tests/integrate/260_NO_DJ_PK_PU_SO_4K/STRU b/tests/integrate/260_NO_DJ_PK_PU_SO_4K/STRU new file mode 100644 index 00000000000..91021e0a697 --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_SO_4K/STRU @@ -0,0 +1,22 @@ +ATOMIC_SPECIES +Fe 1.000 Fe.upf + +NUMERICAL_ORBITAL +Fe_gga_6au_100Ry_4s2p2d1f.orb + +LATTICE_CONSTANT +8.190 + +LATTICE_VECTORS + 1.00 0.50 0.50 + 0.50 1.00 0.50 + 0.50 0.50 1.00 +ATOMIC_POSITIONS +Direct + +Fe +0.0 +2 +0.00 0.00 0.00 mag 1.0 1.0 1.0 +0.51 0.51 0.51 mag 1.0 1.0 1.0 + diff --git a/tests/integrate/260_NO_DJ_PK_PU_SO_4K/result.ref b/tests/integrate/260_NO_DJ_PK_PU_SO_4K/result.ref new file mode 100644 index 00000000000..3d7718e9ca0 --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_SO_4K/result.ref @@ -0,0 +1,4 @@ +etotref -6791.4650419973040698 +etotperatomref -3395.7325209987 +totalforceref 20.417878 +totaltimeref 18.45 diff --git a/tests/integrate/260_NO_DJ_PK_PU_SO_4K/threshold b/tests/integrate/260_NO_DJ_PK_PU_SO_4K/threshold new file mode 100644 index 00000000000..2bb8a636fea --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_SO_4K/threshold @@ -0,0 +1,2 @@ +threshold 1e-6 +force_threshold 1e-4 diff --git a/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/INPUT b/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/INPUT new file mode 100644 index 00000000000..c53c09997f6 --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/INPUT @@ -0,0 +1,46 @@ +INPUT_PARAMETERS +suffix autotest +nbands 40 + +calculation scf +ecutwfc 10 +scf_thr 1.0e-4 +scf_nmax 50 +out_chg 0 + +#init_chg file +#out_dos 1 +#dos_sigma 0.05 +#out_band 1 + +smearing_method gaussian +smearing_sigma 0.01 + +#force_thr_ev 0.01 +#relax_method cg +#relax_bfgs_init 0.5 + +mixing_type pulay +mixing_beta 0.3 +mixing_restart 1e-3 +mixing_dmr 1 +mixing_gg0 1.1 + +ks_solver scalapack_gvx +basis_type lcao +gamma_only 0 +noncolin 1 +lspinorb 1 +nspin 4 +cal_force 1 +cal_stress 0 + +#Parameter DFT+U +dft_plus_u 1 +orbital_corr 2 +hubbard_u 5.0 +onsite_radius 5.0 +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB + +out_mul 1 diff --git a/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/KPT b/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/KPT new file mode 100644 index 00000000000..b1f8afcf383 --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +4 1 1 0 0 0 diff --git a/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/README b/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/README new file mode 100644 index 00000000000..6c8d9567ed6 --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/README @@ -0,0 +1,12 @@ +T6: nspin=4 multi-k Mulliken magnetization calibration (LCAO, DFT+U, Fe). +out_mul=1 produces OUT.autotest/mulliken.txt; check_extra.py validates the +Mulliken per-atom populations/magnetizations against the internal charge sum +rule and the rho-integrated cell magnetization (running_scf.log), in magnitude +and direction. Keys: t6_converged_pass, t6_charge_sum_pass, t6_mul_rho_pass, +t6_mul_rho_dir_pass, t6_mz_sign_pass (0.0 = pass, 1.0 = fail). Tolerance for +the Mulliken-vs-rho cell magnetization is 0.1 uB per component; the measured +deviation on this case is ~0.04 uB on |m|~5.5 uB/cell. No anomalous m_y +deviation (R6) is present, so the B-class gemm change (S^T -> S^dagger) is +not executed. DeltaSpin is disabled for LCAO in this repo +(tests/17_DS_DFTU/CASES_CPU.txt); the Gamma-only control converges to a +different magnetic minimum and is covered in the PR validation record. diff --git a/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/STRU b/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/STRU new file mode 100644 index 00000000000..91021e0a697 --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/STRU @@ -0,0 +1,22 @@ +ATOMIC_SPECIES +Fe 1.000 Fe.upf + +NUMERICAL_ORBITAL +Fe_gga_6au_100Ry_4s2p2d1f.orb + +LATTICE_CONSTANT +8.190 + +LATTICE_VECTORS + 1.00 0.50 0.50 + 0.50 1.00 0.50 + 0.50 0.50 1.00 +ATOMIC_POSITIONS +Direct + +Fe +0.0 +2 +0.00 0.00 0.00 mag 1.0 1.0 1.0 +0.51 0.51 0.51 mag 1.0 1.0 1.0 + diff --git a/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/check_extra.py b/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/check_extra.py new file mode 100644 index 00000000000..c4a0eba6862 --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/check_extra.py @@ -0,0 +1,157 @@ +#!/usr/bin/env python3 +"""T6: nspin=4 multi-k Mulliken magnetization calibration (LCAO, DFT+U, Fe). + +Cross-checks the Mulliken per-atom populations/magnetizations +(OUT.autotest/mulliken.txt, computed from DMK x S(k)) against +1. the internal charge sum rule (sum of per-atom charges == Total charge), +2. the rho-integrated cell magnetization + (OUT.autotest/running_scf.log, "Total magnetism (Bohr mag/cell)"), + both in magnitude and direction. + +This calibration targets the nspin=4 multi-k Mulliken m_y path, which uses +S^T instead of S^dagger: an anomalous m_y deviation (order ||Im S(k)||) would +flag the B-class gemm change. Measured deviation on this case (Fe, 4x1x1, +U=5 eV, converged): + |sum Mulliken - rho_cell| = (0.042, 0.004, 0.044) uB on |m| ~ 5.5 uB/cell, + i.e. the small components agree as well as the large ones, so no anomalous + m_y deviation is present and the B-class change is NOT executed. +The Gamma-only control (same cell, Gamma) converges to a different magnetic +minimum and is not part of the automated test; it is covered in the PR +validation record. DeltaSpin is disabled for LCAO in this repo +(tests/17_DS_DFTU/CASES_CPU.txt), so the DeltaSpin third path is unavailable. + +Each key is a pass/fail flag: 0.0 = pass, 1.0 = fail. Measured values are +written to stderr for diagnostics. + +Usage: check_extra.py [abacus_bin] [np] [case_dir] +""" + +import os +import re +import sys + +KEYS = [ + "t6_converged_pass", + "t6_charge_sum_pass", + "t6_mul_rho_pass", + "t6_mul_rho_dir_pass", + "t6_mz_sign_pass", +] + +TOL_CHARGE = 1e-3 # sum rule, electrons +TOL_MUL_RHO = 0.1 # uB per component +TOL_DIR_DEG = 2.0 # degrees between Mulliken-sum and rho vectors +MZ_MIN = 1.0 # uB, both Fe sites must stay ferromagnetic along z + + +def fail(msg): + for k in KEYS: + print(k + " 1.0", flush=True) + sys.stderr.write("t6_check.py: " + msg + "\n") + sys.exit(1) + + +def parse_mulliken(path): + text = open(path).read() + total = None + atoms = {} # idx -> {"chg": float, "mag": [mx,my,mz]} + for line in text.splitlines(): + m = re.match(r"^\s*Total charge\s+([-+0-9.eE]+)\s*$", line) + if m: + total = float(m.group(1)) + continue + m = re.match(r"^\s*total charge\s+on atom\s+(\d+)\s+([-+0-9.eE]+)", line) + if m: + atoms.setdefault(int(m.group(1)), {})["chg"] = float(m.group(2)) + continue + m = re.match( + r"^\s*total magnetism\s+on atom\s+(\d+)\s+([-+0-9.eE]+)" + r"(?:\s+([-+0-9.eE]+)\s+([-+0-9.eE]+))?", line) + if m: + vals = [float(m.group(i)) for i in (2, 3, 4) if m.group(i)] + atoms.setdefault(int(m.group(1)), {})["mag"] = vals + return total, atoms + + +def parse_rho_mag(path): + text = open(path).read() + if not re.search(r"#SCF IS CONVERGED#", text): + return None, None + lines = [ln for ln in text.splitlines() if "Total magnetism (Bohr mag/cell)" in ln] + if not lines: + return None, None + last = lines[-1] + m = re.search(r"=\s*\[\s*([-+0-9.eE]+)\s*,\s*([-+0-9.eE]+)\s*,\s*([-+0-9.eE]+)", last) + if not m: + return None, None + return [float(m.group(i)) for i in (1, 2, 3)], True + + +def main(): + case_dir = os.path.abspath(sys.argv[3]) if len(sys.argv) > 3 else os.getcwd() + mul_path = os.path.join(case_dir, "OUT.autotest", "mulliken.txt") + scf_path = os.path.join(case_dir, "OUT.autotest", "running_scf.log") + if not (os.path.isfile(mul_path) and os.path.isfile(scf_path)): + fail("OUT.autotest/mulliken.txt or running_scf.log missing") + + total, atoms = parse_mulliken(mul_path) + if total is None or not atoms: + fail("could not parse mulliken.txt") + rho_mag, converged = parse_rho_mag(scf_path) + + results = {} + results["t6_converged_pass"] = 0.0 if converged else 1.0 + + # 1. charge sum rule: sum of per-atom charges == Total charge + qsum = sum(a["chg"] for a in atoms.values()) + results["t6_charge_sum_pass"] = 0.0 if abs(qsum - total) < TOL_CHARGE else 1.0 + + # 2. Mulliken cell magnetization vs rho cell magnetization + msum = [0.0, 0.0, 0.0] + for a in atoms.values(): + mag = a.get("mag") + if mag is None or len(mag) != 3: + fail("nspin=4 magnetism vector expected for every atom") + for c in range(3): + msum[c] += mag[c] + if rho_mag is None: + results["t6_mul_rho_pass"] = 1.0 + results["t6_mul_rho_dir_pass"] = 1.0 + else: + max_dev = max(abs(msum[c] - rho_mag[c]) for c in range(3)) + results["t6_mul_rho_pass"] = 0.0 if max_dev < TOL_MUL_RHO else 1.0 + na = sum(v * v for v in msum) ** 0.5 + nb = sum(v * v for v in rho_mag) ** 0.5 + if na < 1e-8 or nb < 1e-8: + results["t6_mul_rho_dir_pass"] = 1.0 + else: + cosang = sum(msum[c] * rho_mag[c] for c in range(3)) / (na * nb) + cosang = max(-1.0, min(1.0, cosang)) + deg = os_acos_deg(cosang) + results["t6_mul_rho_dir_pass"] = 0.0 if deg < TOL_DIR_DEG else 1.0 + + # 3. magnetic state guard: both Fe sites ferromagnetic along z + mz = [a["mag"][2] for a in sorted(atoms.values(), key=lambda a: a["chg"])] + results["t6_mz_sign_pass"] = 0.0 if all(m > MZ_MIN for m in mz) else 1.0 + + for k in KEYS: + print(k + " %.1f" % results[k], flush=True) + sys.stderr.write( + "t6_check.py: total_charge=%.6f sum_atoms=%.6f | " + "sum_mulliken=(%.6f,%.6f,%.6f) rho_cell=(%.6f,%.6f,%.6f) | " + "max_dev=%.6f uB converged=%s\n" + % (total, qsum, msum[0], msum[1], msum[2], + rho_mag[0] if rho_mag else float("nan"), + rho_mag[1] if rho_mag else float("nan"), + rho_mag[2] if rho_mag else float("nan"), + max(abs(msum[c] - (rho_mag[c] if rho_mag else 0.0)) for c in range(3)), + bool(converged))) + + +def os_acos_deg(cosang): + import math + return math.degrees(math.acos(cosang)) + + +if __name__ == "__main__": + main() diff --git a/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/result.ref b/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/result.ref new file mode 100644 index 00000000000..7a3445de33e --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/result.ref @@ -0,0 +1,8 @@ +etotref -6791.4650419973040698 +etotperatomref -3395.7325209987 +totalforceref 20.417878 +t6_converged_pass 0.0 +t6_charge_sum_pass 0.0 +t6_mul_rho_pass 0.0 +t6_mul_rho_dir_pass 0.0 +t6_mz_sign_pass 0.0 diff --git a/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/threshold b/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/threshold new file mode 100644 index 00000000000..2bb8a636fea --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_SO_4K_MUL/threshold @@ -0,0 +1,2 @@ +threshold 1e-6 +force_threshold 1e-4 diff --git a/tests/integrate/Autotest.sh b/tests/integrate/Autotest.sh index ad00e05ec12..be5ba2ed3d2 100755 --- a/tests/integrate/Autotest.sh +++ b/tests/integrate/Autotest.sh @@ -115,8 +115,8 @@ check_out(){ for key in $properties; do if [ $key == "totaltimeref" ]; then - # echo "time=$cal ref=$ref" - break + # totaltimeref is not compared (wall-clock dependent) + continue fi #-------------------------------------------------- @@ -294,6 +294,12 @@ for dir in $testdir; do if test -z $g then bash -e ../../integrate/tools/catch_properties.sh result.out + # optional case-local extra checks: if check_extra.py exists, + # it is run as `python3 check_extra.py ` + # and its stdout (key-value lines) is appended to result.out + if test -e check_extra.py; then + python3 check_extra.py "$abacus" "$np" "$(pwd)" >> result.out + fi if [ $? -ne 0 ]; then echo -e "\e[0;31m [ERROR ] Fatal Error in catch_properties.sh \e[0m" let fatal++ @@ -308,6 +314,9 @@ for dir in $testdir; do fi else bash -e ../../integrate/tools/catch_properties.sh result.ref + if test -e check_extra.py; then + python3 check_extra.py "$abacus" "$np" "$(pwd)" >> result.ref + fi fi fi diff --git a/tests/integrate/CASES_CPU.txt b/tests/integrate/CASES_CPU.txt index e69de29bb2d..5d2afe26bc5 100644 --- a/tests/integrate/CASES_CPU.txt +++ b/tests/integrate/CASES_CPU.txt @@ -0,0 +1,4 @@ +240_NO_KP_15_SO_FD +260_NO_DJ_PK_PU_SO_4K +260_NO_DJ_PK_PU_SO_4K_MUL +260_NO_DJ_PK_PU_AFM_URAMPING_MUL