diff --git a/source/source_base/module_fft/fft_bundle.cpp b/source/source_base/module_fft/fft_bundle.cpp index 2bca85aefcc..925ce187dae 100644 --- a/source/source_base/module_fft/fft_bundle.cpp +++ b/source/source_base/module_fft/fft_bundle.cpp @@ -110,6 +110,28 @@ void FFT_Bundle::initfft(int nx_in, }else{ ModuleBase::WARNING_QUIT("FFT_Bundle", "Please set the device to cpu or gpu or dsp"); } + + if (this->device != "cpu") + { + // The gpu/dsp backends implement only the device-memory 3D FFT + // interface; the host-memory interface (used by post-processing that + // runs on the CPU wavefunction mirror) is served by the extra + // FFT_CPU backends host_fft_float/host_fft_double. + if (this->float_flag) + { +#if defined(__FLOAT_FFTW) + this->host_fft_float = make_unique>(this->fft_mode); + this->host_fft_float->initfft( + nx_in, ny_in, nz_in, lixy_in, rixy_in, ns_in, nplane_in, nproc_in, gamma_only_in, xprime_in); +#endif + } + // The double host backend is needed regardless of the compute precision: + // post-processing (e.g. the Wannier90 output path) always reads the + // std::complex CPU wavefunction mirror, even for precision=single. + this->host_fft_double = make_unique>(this->fft_mode); + this->host_fft_double->initfft( + nx_in, ny_in, nz_in, lixy_in, rixy_in, ns_in, nplane_in, nproc_in, gamma_only_in, xprime_in); + } } void FFT_Bundle::setupFFT() @@ -118,9 +140,17 @@ void FFT_Bundle::setupFFT() { fft_double->setupFFT(); } + if (host_fft_double != nullptr) + { + host_fft_double->setupFFT(); + } if (float_flag) { fft_float->setupFFT(); + if (host_fft_float != nullptr) + { + host_fft_float->setupFFT(); + } } } @@ -130,9 +160,17 @@ void FFT_Bundle::clearFFT() { fft_double->cleanFFT(); } + if (host_fft_double != nullptr) + { + host_fft_double->cleanFFT(); + } if (float_flag) { fft_float->cleanFFT(); + if (host_fft_float != nullptr) + { + host_fft_float->cleanFFT(); + } } } void FFT_Bundle::clear() @@ -142,9 +180,17 @@ void FFT_Bundle::clear() { fft_double->clear(); } + if (host_fft_double != nullptr) + { + host_fft_double->clear(); + } if (float_flag) { fft_float->clear(); + if (host_fft_float != nullptr) + { + host_fft_float->clear(); + } } } @@ -165,67 +211,67 @@ void FFT_Bundle::resource_handler(const int flag) const template <> void FFT_Bundle::fftxyfor(std::complex* in, std::complex* out) const { - fft_float->fftxyfor(in, out); + this->get_host_fft_float()->fftxyfor(in, out); } template <> void FFT_Bundle::fftxyfor(std::complex* in, std::complex* out) const { - fft_double->fftxyfor(in, out); + this->get_host_fft_double()->fftxyfor(in, out); } template <> void FFT_Bundle::fftzfor(std::complex* in, std::complex* out) const { - fft_float->fftzfor(in, out); + this->get_host_fft_float()->fftzfor(in, out); } template <> void FFT_Bundle::fftzfor(std::complex* in, std::complex* out) const { - fft_double->fftzfor(in, out); + this->get_host_fft_double()->fftzfor(in, out); } template <> void FFT_Bundle::fftxybac(std::complex* in, std::complex* out) const { - fft_float->fftxybac(in, out); + this->get_host_fft_float()->fftxybac(in, out); } template <> void FFT_Bundle::fftxybac(std::complex* in, std::complex* out) const { - fft_double->fftxybac(in, out); + this->get_host_fft_double()->fftxybac(in, out); } template <> void FFT_Bundle::fftzbac(std::complex* in, std::complex* out) const { - fft_float->fftzbac(in, out); + this->get_host_fft_float()->fftzbac(in, out); } template <> void FFT_Bundle::fftzbac(std::complex* in, std::complex* out) const { - fft_double->fftzbac(in, out); + this->get_host_fft_double()->fftzbac(in, out); } template <> void FFT_Bundle::fftxyr2c(float* in, std::complex* out) const { - fft_float->fftxyr2c(in, out); + this->get_host_fft_float()->fftxyr2c(in, out); } template <> void FFT_Bundle::fftxyr2c(double* in, std::complex* out) const { - fft_double->fftxyr2c(in, out); + this->get_host_fft_double()->fftxyr2c(in, out); } template <> void FFT_Bundle::fftxyc2r(std::complex* in, float* out) const { - fft_float->fftxyc2r(in, out); + this->get_host_fft_float()->fftxyc2r(in, out); } template <> void FFT_Bundle::fftxyc2r(std::complex* in, double* out) const { - fft_double->fftxyc2r(in, out); + this->get_host_fft_double()->fftxyc2r(in, out); } template <> @@ -258,34 +304,34 @@ void FFT_Bundle::fft3D_backward(std::complex* in, template <> float* FFT_Bundle::get_rspace_data() const { - return fft_float->get_rspace_data(); + return this->get_host_fft_float()->get_rspace_data(); } template <> double* FFT_Bundle::get_rspace_data() const { - return fft_double->get_rspace_data(); + return this->get_host_fft_double()->get_rspace_data(); } template <> std::complex* FFT_Bundle::get_auxr_data() const { - return fft_float->get_auxr_data(); + return this->get_host_fft_float()->get_auxr_data(); } template <> std::complex* FFT_Bundle::get_auxr_data() const { - return fft_double->get_auxr_data(); + return this->get_host_fft_double()->get_auxr_data(); } template <> std::complex* FFT_Bundle::get_auxg_data() const { - return fft_float->get_auxg_data(); + return this->get_host_fft_float()->get_auxg_data(); } template <> std::complex* FFT_Bundle::get_auxg_data() const { - return fft_double->get_auxg_data(); + return this->get_host_fft_double()->get_auxg_data(); } template <> @@ -298,4 +344,4 @@ std::complex* FFT_Bundle::get_auxr_3d_data() const { return fft_double->get_auxr_3d_data(); } -} // namespace ModuleBase \ No newline at end of file +} // namespace ModuleBase diff --git a/source/source_base/module_fft/fft_bundle.h b/source/source_base/module_fft/fft_bundle.h index d82ca09998f..79163e89dc0 100644 --- a/source/source_base/module_fft/fft_bundle.h +++ b/source/source_base/module_fft/fft_bundle.h @@ -210,8 +210,31 @@ class FFT_Bundle int fft_mode = 0; bool float_flag = false; bool double_flag = false; + // Backends selected by `device`: an FFT_CPU instance for "cpu", + // FFT_CUDA/FFT_ROCM for "gpu", FFT_DSP for "dsp". They serve the + // device-memory ("gpu-like") 3D FFT interface; on "cpu" they also serve + // the host-memory ("cpu-like") interface. std::shared_ptr> fft_float = nullptr; std::shared_ptr> fft_double = nullptr; + // Extra CPU backends serving the host-memory ("cpu-like") interface when + // the configured device is not "cpu", since the gpu/dsp backends implement + // only the 3D FFT interface. They keep host post-processing such as the + // Wannier90 interface and Berry phase working, since those run on the CPU + // wavefunction mirror even in GPU calculations. Null when + // device == "cpu", where fft_float/fft_double already serve both + // interfaces. + std::shared_ptr> host_fft_float = nullptr; + std::shared_ptr> host_fft_double = nullptr; + + // The backend that serves the host-memory interface, whichever it is. + FFT_BASE* get_host_fft_float() const + { + return (this->device == "cpu") ? this->fft_float.get() : this->host_fft_float.get(); + } + FFT_BASE* get_host_fft_double() const + { + return (this->device == "cpu") ? this->fft_double.get() : this->host_fft_double.get(); + } std::string device = "cpu"; std::string precision = "double"; diff --git a/tests/01_PW/CASES_GPU.txt b/tests/01_PW/CASES_GPU.txt index fe156d9365b..2cb811a036d 100644 --- a/tests/01_PW/CASES_GPU.txt +++ b/tests/01_PW/CASES_GPU.txt @@ -84,7 +84,7 @@ scf_pbe0_spin1 scf_pbe0_spin2 099_PW_15_SO_avg #100_PW_DJ_SO -#101_PW_W90 +101_PW_W90 102_PW_MD_1O 103_PW_MD_2O #104_PW_Gene_Descriptors diff --git a/tests/integrate/validation_tools/compare_amn.py b/tests/integrate/validation_tools/compare_amn.py new file mode 100644 index 00000000000..75017244943 --- /dev/null +++ b/tests/integrate/validation_tools/compare_amn.py @@ -0,0 +1,194 @@ +#!/usr/bin/env python3 +"""Gauge-invariant comparison of Wannier90 .amn files. + +The overlap matrix elements A_mn(k) = are covariant under the +band gauge: an arbitrary unitary mixing U of the selected bands maps +A -> U^dagger A. The Gram matrix A^dagger A is invariant under this +transformation, so it is compared here instead of the raw elements. This +lets runs whose eigensolvers pick different gauges (e.g. LAPACK on CPU and +cuSOLVER on GPU) be validated against the same reference. + +Usage: + compare_amn.py ref.amn cal.amn [accuracy] + +accuracy is the number of decimal digits: the tolerance is 10^-accuracy +(default 6). Exit code is 0 when the files agree up to the band gauge and +1 otherwise. +""" + +import math +import sys + +MAX_REPORTED = 10 + + +def fail(message): + print("CompareAMN Error: " + message) + sys.exit(1) + + +def is_int(token): + try: + int(token) + return True + except ValueError: + return False + + +def is_float(token): + try: + float(token) + return True + except ValueError: + return False + + +def read_amn(path): + """Return ((nbnd, nk, nwann), A[k][m][n]) from a Wannier90 .amn file. + + Lines before the 'nbnd nk nwann' header (e.g. ABACUS's + 'Created on ...' line) are ignored. Every entry is checked for range, + duplicates and completeness so that a truncated or corrupted output is + reported instead of silently producing a small Gram matrix. + """ + try: + with open(path) as handle: + lines = handle.readlines() + except OSError: + fail("can not find file %s" % path) + + header = None + data = None + filled = None + nentries = 0 + + for lineno, line in enumerate(lines, 1): + tokens = line.split() + if not tokens: + continue + + if header is None: + if len(tokens) == 3 and all(is_int(t) for t in tokens): + nbnd, nk, nwann = (int(t) for t in tokens) + if nbnd <= 0 or nk <= 0 or nwann <= 0: + fail("%s:%d: invalid header '%s'" % (path, lineno, line.strip())) + header = (nbnd, nk, nwann) + data = [[[0j] * nwann for _ in range(nbnd)] for _ in range(nk)] + filled = [[[False] * nwann for _ in range(nbnd)] for _ in range(nk)] + elif len(tokens) == 5 and is_int(tokens[0]) \ + and is_int(tokens[1]) and is_int(tokens[2]) \ + and is_float(tokens[3]) and is_float(tokens[4]): + fail("%s:%d: data line before the header" % (path, lineno)) + # Otherwise this is a preamble line such as "Created on ...". + continue + + if len(tokens) != 5: + fail("%s:%d: expected 5 columns 'm n k Re Im', got %d" + % (path, lineno, len(tokens))) + if not (is_int(tokens[0]) and is_int(tokens[1]) and is_int(tokens[2]) + and is_float(tokens[3]) and is_float(tokens[4])): + fail("%s:%d: malformed data line '%s'" % (path, lineno, line.strip())) + + m = int(tokens[0]) + n = int(tokens[1]) + k = int(tokens[2]) + real = float(tokens[3]) + imag = float(tokens[4]) + if not (math.isfinite(real) and math.isfinite(imag)): + fail("%s:%d: non-finite matrix element" % (path, lineno)) + + nbnd, nk, nwann = header + if not (1 <= m <= nbnd and 1 <= n <= nwann and 1 <= k <= nk): + fail("%s:%d: index out of range (m=%d n=%d k=%d, header %d %d %d)" + % (path, lineno, m, n, k, nbnd, nk, nwann)) + if filled[k - 1][m - 1][n - 1]: + fail("%s:%d: duplicate entry m=%d n=%d k=%d" % (path, lineno, m, n, k)) + data[k - 1][m - 1][n - 1] = complex(real, imag) + filled[k - 1][m - 1][n - 1] = True + nentries += 1 + + if header is None: + fail("%s: no valid header 'nbnd nk nwann' found" % path) + + expected = header[0] * header[1] * header[2] + if nentries != expected: + fail("%s: found %d entries, expected %d" % (path, nentries, expected)) + for k in range(header[1]): + for m in range(header[0]): + for n in range(header[2]): + if not filled[k][m][n]: + fail("%s: missing entry m=%d n=%d k=%d" % (path, m + 1, n + 1, k + 1)) + return header, data + + +def gram_matrix(matrix): + """Return A^dagger A for a (nbnd x nwann) matrix.""" + nbnd = len(matrix) + nwann = len(matrix[0]) + gram = [[0j] * nwann for _ in range(nwann)] + for m in range(nbnd): + for i in range(nwann): + ai_conj = matrix[m][i].conjugate() + for j in range(nwann): + gram[i][j] += ai_conj * matrix[m][j] + return gram + + +def main(): + if len(sys.argv) not in (3, 4): + print(__doc__.strip()) + sys.exit(1) + + ref_path = sys.argv[1] + cal_path = sys.argv[2] + accuracy = 6 + if len(sys.argv) == 4: + try: + accuracy = int(sys.argv[3]) + except ValueError: + fail("accuracy must be an integer, got '%s'" % sys.argv[3]) + if accuracy < 0: + fail("accuracy must be non-negative") + tolerance = 10.0 ** (-accuracy) + + ref_header, ref_data = read_amn(ref_path) + cal_header, cal_data = read_amn(cal_path) + if ref_header != cal_header: + fail("header mismatch: %s has %s, %s has %s" + % (ref_path, ref_header, cal_path, cal_header)) + + nbnd, nk, nwann = ref_header + failures = 0 + reported = 0 + max_diff = 0.0 + for k in range(nk): + ref_gram = gram_matrix(ref_data[k]) + cal_gram = gram_matrix(cal_data[k]) + for i in range(nwann): + for j in range(nwann): + diff_real = abs(ref_gram[i][j].real - cal_gram[i][j].real) + diff_imag = abs(ref_gram[i][j].imag - cal_gram[i][j].imag) + max_diff = max(max_diff, diff_real, diff_imag) + if diff_real > tolerance or diff_imag > tolerance: + failures += 1 + if reported < MAX_REPORTED: + print("Mismatch at k=%d, (A^dagger A)[%d,%d]: " + "ref=(%.12g,%.12g) cal=(%.12g,%.12g) diff=(%.3g,%.3g)" + % (k + 1, i + 1, j + 1, + ref_gram[i][j].real, ref_gram[i][j].imag, + cal_gram[i][j].real, cal_gram[i][j].imag, + diff_real, diff_imag)) + reported += 1 + + if failures: + print("AMN check FAILED: %d Gram-matrix components exceed %g " + "(max diff %g)" % (failures, tolerance, max_diff)) + sys.exit(1) + + print("AMN check PASSED: %d k-points, %d bands, %d Wannier functions; " + "max |d(A^dagger A)| = %.3g (tolerance %g)" + % (nk, nbnd, nwann, max_diff, tolerance)) + + +if __name__ == "__main__": + main() diff --git a/tests/integrate/validation_tools/compare_mmn.py b/tests/integrate/validation_tools/compare_mmn.py new file mode 100644 index 00000000000..96171433d4c --- /dev/null +++ b/tests/integrate/validation_tools/compare_mmn.py @@ -0,0 +1,213 @@ +#!/usr/bin/env python3 +"""Gauge-invariant comparison of Wannier90 .mmn files. + +M_mn(k,b) = is covariant under the band gauge: +M(k,b) -> U_k^dagger M(k,b) U_{k+b}. The basis of the selected bands is +solver-dependent: the CPU (LAPACK) and GPU (cuSOLVER) diagonalizations of +the same subspace differ by a per-band phase at each k. For this diamond +case the gauge recovered from the .amn files is diagonal to ~1e-7 (median +~1e-10), and the elementwise magnitudes of A agree to ~1e-7. Under such a +diagonal gauge |M_mn| is invariant, so the magnitudes are compared here. +The raw complex elements cannot be compared because their arbitrary phases +are solver-dependent. The diagonal-phase picture assumes the selected +bands are non-degenerate; a mixed degenerate subspace would require the +full gauge-invariant content (loop products of M around closed paths). + +File layout (as written by toW90_PW::cal_Mmn): an optional preamble such +as " Created on ...", the header 'nbnd nk nntot', then nk*nntot blocks. +Each block starts with 'ik ikb G1 G2 G3' (ikb is the 1-based neighbour +k-point index from the .nnkp file), followed by nbnd*nbnd lines 'Re Im' +with m varying fastest. + +Usage: + compare_mmn.py ref.mmn cal.mmn [accuracy] + +accuracy is the number of decimal digits: the tolerance is 10^-accuracy +(default 6). Exit code is 0 when the files agree up to the band gauge and +1 otherwise. +""" + +import math +import sys + +MAX_REPORTED = 10 + + +def fail(message): + print("CompareMMN Error: " + message) + sys.exit(1) + + +def is_int(token): + try: + int(token) + return True + except ValueError: + return False + + +def is_float(token): + try: + float(token) + return True + except ValueError: + return False + + +def read_mmn(path): + """Return ((nbnd, nk, nntot), blocks) from a Wannier90 .mmn file. + + blocks maps (ik, ikb) -> (G, matrix) with matrix[m][n] = M_mn. + Lines before the 'nbnd nk nntot' header (e.g. ABACUS's + 'Created on ...' preamble) are ignored. + """ + try: + with open(path) as handle: + lines = handle.readlines() + except OSError: + fail("can not find file %s" % path) + + header = None + lineno = 0 + for index, line in enumerate(lines): + tokens = line.split() + if not tokens: + continue + if len(tokens) == 3 and all(is_int(t) for t in tokens): + nbnd, nk, nntot = (int(t) for t in tokens) + if nbnd <= 0 or nk <= 0 or nntot <= 0: + fail("%s:%d: invalid header '%s'" % (path, index + 1, line.strip())) + header = (nbnd, nk, nntot) + lineno = index + 1 + break + if len(tokens) == 5 and all(is_int(t) for t in tokens): + fail("%s:%d: block header before the nbnd nk nntot header" + % (path, index + 1)) + if len(tokens) == 2 and all(is_float(t) for t in tokens): + fail("%s:%d: data line before the header" % (path, index + 1)) + # Otherwise this is a preamble line such as "Created on ...". + if header is None: + fail("%s: no valid header 'nbnd nk nntot' found" % path) + + nbnd, nk, nntot = header + blocks = {} + counts = [0] * nk + total = nk * nntot + + def next_content_line(describe): + nonlocal lineno + while lineno < len(lines) and not lines[lineno].strip(): + lineno += 1 + if lineno >= len(lines): + fail("%s: unexpected end of file while %s" % (path, describe)) + line = lines[lineno] + lineno += 1 + return line + + for _ in range(total): + line = next_content_line("reading block header %d of %d" + % (len(blocks) + 1, total)) + tokens = line.split() + if len(tokens) != 5 or not all(is_int(t) for t in tokens): + fail("%s:%d: expected block header 'ik ikb G1 G2 G3', got '%s'" + % (path, lineno, line.strip())) + ik, ikb, gx, gy, gz = (int(t) for t in tokens) + if not (1 <= ik <= nk and 1 <= ikb <= nk): + fail("%s:%d: ik/ikb out of range (ik=%d ikb=%d, nk=%d)" + % (path, lineno, ik, ikb, nk)) + if (ik, ikb) in blocks: + fail("%s:%d: duplicate block ik=%d ikb=%d" % (path, lineno, ik, ikb)) + + matrix = [[0j] * nbnd for _ in range(nbnd)] + for n in range(nbnd): + for m in range(nbnd): + line = next_content_line( + "reading M(%d,%d) of block ik=%d ikb=%d" % (m + 1, n + 1, ik, ikb)) + tokens = line.split() + if len(tokens) != 2 or not all(is_float(t) for t in tokens): + fail("%s:%d: expected 'Re Im' in block ik=%d ikb=%d, got '%s'" + % (path, lineno, ik, ikb, line.strip())) + real = float(tokens[0]) + imag = float(tokens[1]) + if not (math.isfinite(real) and math.isfinite(imag)): + fail("%s:%d: non-finite matrix element" % (path, lineno)) + matrix[m][n] = complex(real, imag) + + blocks[(ik, ikb)] = ((gx, gy, gz), matrix) + counts[ik - 1] += 1 + + for index in range(lineno, len(lines)): + if lines[index].strip(): + fail("%s:%d: unexpected trailing content '%s'" + % (path, index + 1, lines[index].strip())) + for ik in range(nk): + if counts[ik] != nntot: + fail("%s: k-point %d has %d blocks, expected %d" + % (path, ik + 1, counts[ik], nntot)) + return header, blocks + + +def compare(ref_path, cal_path, tolerance): + ref_header, ref_blocks = read_mmn(ref_path) + cal_header, cal_blocks = read_mmn(cal_path) + if ref_header != cal_header: + fail("header mismatch: %s has %s, %s has %s" + % (ref_path, ref_header, cal_path, cal_header)) + + missing = sorted(set(ref_blocks) - set(cal_blocks)) + extra = sorted(set(cal_blocks) - set(ref_blocks)) + if missing or extra: + fail("block set mismatch: %d missing (e.g. %s), %d extra (e.g. %s)" + % (len(missing), missing[:1], len(extra), extra[:1])) + + nbnd, nk, nntot = ref_header + failures = 0 + reported = 0 + max_diff = 0.0 + for ik, ikb in sorted(ref_blocks): + ref_g, ref_matrix = ref_blocks[(ik, ikb)] + cal_g, cal_matrix = cal_blocks[(ik, ikb)] + if ref_g != cal_g: + fail("block ik=%d ikb=%d: G vector differs, ref=%s cal=%s" + % (ik, ikb, ref_g, cal_g)) + for m in range(nbnd): + for n in range(nbnd): + diff = abs(abs(ref_matrix[m][n]) - abs(cal_matrix[m][n])) + max_diff = max(max_diff, diff) + if diff > tolerance: + failures += 1 + if reported < MAX_REPORTED: + print("Mismatch at ik=%d ikb=%d, |M[%d,%d]|: " + "ref=%.12g cal=%.12g diff=%.3g" + % (ik, ikb, m + 1, n + 1, + abs(ref_matrix[m][n]), abs(cal_matrix[m][n]), diff)) + reported += 1 + + if failures: + print("MMN check FAILED: %d magnitude components exceed %g " + "(max diff %g)" % (failures, tolerance, max_diff)) + sys.exit(1) + + print("MMN check PASSED: %d k-points, %d bands, %d b-vectors per k-point; " + "max |d|M|| = %.3g (tolerance %g)" + % (nk, nbnd, nntot, max_diff, tolerance)) + + +def main(): + if len(sys.argv) not in (3, 4): + print(__doc__.strip()) + sys.exit(1) + + accuracy = 6 + if len(sys.argv) == 4: + try: + accuracy = int(sys.argv[3]) + except ValueError: + fail("accuracy must be an integer, got '%s'" % sys.argv[3]) + if accuracy < 0: + fail("accuracy must be non-negative") + compare(sys.argv[1], sys.argv[2], 10.0 ** (-accuracy)) + + +if __name__ == "__main__": + main() diff --git a/tests/integrate/validation_tools/props_common.sh b/tests/integrate/validation_tools/props_common.sh index 590f157b0a0..bca1717e6e6 100644 --- a/tests/integrate/validation_tools/props_common.sh +++ b/tests/integrate/validation_tools/props_common.sh @@ -23,6 +23,8 @@ PROPS_TOOLS_DIR=$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd) COMPARE_SCRIPT="$PROPS_TOOLS_DIR/CompareFile.py" +COMPARE_AMN_SCRIPT="$PROPS_TOOLS_DIR/compare_amn.py" +COMPARE_MMN_SCRIPT="$PROPS_TOOLS_DIR/compare_mmn.py" CUBE_TOOL="$PROPS_TOOLS_DIR/cube_tool.py" COLLECT_NPY_MEANS="$PROPS_TOOLS_DIR/collect_npy_means.py" diff --git a/tests/integrate/validation_tools/props_mat.sh b/tests/integrate/validation_tools/props_mat.sh index 3123a8d4fc9..5128c3dd379 100644 --- a/tests/integrate/validation_tools/props_mat.sh +++ b/tests/integrate/validation_tools/props_mat.sh @@ -86,11 +86,13 @@ if ! test -z "$toW90" && [ $toW90 == 1 ]; then mmncal=OUT.autotest/diamond.mmn eigref=diamond.eig eigcal=OUT.autotest/diamond.eig - sed -i '1d' $amncal - sed -i '1d' $mmncal - python3 $COMPARE_SCRIPT $amnref $amncal 1 -abs 8 + # AMN is covariant under the band gauge (A -> U^dagger A), so compare + # the gauge-invariant Gram matrix A^dagger A instead of the raw elements. + python3 $COMPARE_AMN_SCRIPT $amnref $amncal echo "CompareAMN_pass $?" >>$props_result_file - python3 $COMPARE_SCRIPT $mmnref $mmncal 1 -abs 8 + # MMN is covariant as M(k,b) -> U_k^dagger M(k,b) U_{k+b}; CPU and GPU + # differ by per-band phases, under which |M_mn| is invariant. + python3 $COMPARE_MMN_SCRIPT $mmnref $mmncal echo "CompareMMN_pass $?" >>$props_result_file python3 $COMPARE_SCRIPT $eigref $eigcal 8 echo "CompareEIG_pass $?" >>$props_result_file