From cd5073f98b2ee30b13715cf9f0710142e7c332c5 Mon Sep 17 00:00:00 2001 From: guhou-hvi <1322077960@qq.com> Date: Sat, 19 Sep 2026 12:58:18 +0800 Subject: [PATCH 1/2] fix(qe): support ibrav 4-14 in convert_celldm and PW/CP readers --- dpdata/formats/qe/scf.py | 19 +-- dpdata/formats/qe/traj.py | 253 +++++++++++++++++++++++++++++++++++--- tests/test_qe_cp_traj.py | 79 +++++++++++- 3 files changed, 317 insertions(+), 34 deletions(-) diff --git a/dpdata/formats/qe/scf.py b/dpdata/formats/qe/scf.py index 341261d22..728197b5f 100755 --- a/dpdata/formats/qe/scf.py +++ b/dpdata/formats/qe/scf.py @@ -8,7 +8,10 @@ from dpdata.utils import open_file from .traj import ( + convert_celldm, kbar2evperang3, + load_celldm, + resolve_celldm, ry2ev, ) from .traj import ( @@ -64,21 +67,9 @@ def get_cell(lines): for ii in blk: ret.append([float(jj) for jj in ii.split()[0:3]]) ret = np.array(ret) - elif ibrav == 1: - a = None - for iline in lines: - line = iline.replace("=", " ").replace(",", "").split() - if len(line) >= 2 and "a" == line[0]: - # print("line = ", line) - a = float(line[1]) - if len(line) >= 2 and "celldm(1)" == line[0]: - a = float(line[1]) * bohr2ang - # print("a = ", a) - if not a: - raise RuntimeError("parameter 'a' or 'celldm(1)' cannot be found.") - ret = np.array([[a, 0.0, 0.0], [0.0, a, 0.0], [0.0, 0.0, a]]) else: - raise RuntimeError("ibrav > 1 not supported yet.") + celldm = resolve_celldm(lines, ibrav, load_celldm(lines)) + ret = convert_celldm(ibrav, celldm) * bohr2ang return ret diff --git a/dpdata/formats/qe/traj.py b/dpdata/formats/qe/traj.py index 954cd386f..9355cd528 100644 --- a/dpdata/formats/qe/traj.py +++ b/dpdata/formats/qe/traj.py @@ -1,7 +1,6 @@ #!/usr/bin/python3 from __future__ import annotations -import warnings from typing import TYPE_CHECKING import numpy as np @@ -55,22 +54,175 @@ def load_block(lines, key, nlines): def convert_celldm(ibrav, celldm): + """Build a 3x3 cell matrix from QE ``ibrav`` and ``celldm(1..6)``. + + Rows are the lattice vectors a1, a2, a3 in the same length unit as + ``celldm[0]`` (Bohr for QE inputs). Formulas follow Quantum ESPRESSO's + ``latgen`` (``Modules/latgen.f90``) so the cell matches QE's own output. + """ + sr2 = np.sqrt(2.0) + sr3 = np.sqrt(3.0) + a = celldm[0] + cell = np.zeros((3, 3)) + if ibrav == 1: - return celldm[0] * np.eye(3) + cell[0, 0] = a + cell[1, 1] = a + cell[2, 2] = a elif ibrav == 2: - return celldm[0] * 0.5 * np.array([[-1, 0, 1], [0, 1, 1], [-1, 1, 0]]) - elif ibrav == 3: - return celldm[0] * 0.5 * np.array([[1, 1, 1], [-1, 1, 1], [-1, -1, 1]]) - elif ibrav == -3: - return celldm[0] * 0.5 * np.array([[-1, 1, 1], [1, -1, 1], [1, 1, -1]]) - else: - warnings.warn( - "unsupported ibrav " - + str(ibrav) - + " if no .cel file, the cell convertion may be wrong. " + t = a / 2.0 + cell[0, 0] = -t + cell[0, 2] = t + cell[1, 1] = t + cell[1, 2] = t + cell[2, 0] = -t + cell[2, 1] = t + elif abs(ibrav) == 3: + t = a / 2.0 + cell[:, :] = t + if ibrav < 0: + cell[0, 0] = -t + cell[1, 1] = -t + cell[2, 2] = -t + else: + cell[1, 0] = -t + cell[2, 0] = -t + cell[2, 1] = -t + elif ibrav == 4: + cell[0, 0] = a + cell[1, 0] = -a / 2.0 + cell[1, 1] = a * sr3 / 2.0 + cell[2, 2] = a * celldm[2] + elif ibrav == 5: + cosab = celldm[3] + term1 = np.sqrt(1.0 + 2.0 * cosab) + term2 = np.sqrt(1.0 - cosab) + cell[1, 1] = sr2 * a * term2 / sr3 + cell[1, 2] = a * term1 / sr3 + cell[0, 0] = a * term2 / sr2 + cell[0, 1] = -cell[0, 0] / sr3 + cell[0, 2] = cell[1, 2] + cell[2, 0] = -cell[0, 0] + cell[2, 1] = cell[0, 1] + cell[2, 2] = cell[1, 2] + elif ibrav == -5: + cosab = celldm[3] + term1 = np.sqrt(1.0 + 2.0 * cosab) + term2 = np.sqrt(1.0 - cosab) + cell[0, 0] = a * (term1 - 2.0 * term2) / 3.0 + cell[0, 1] = a * (term1 + term2) / 3.0 + cell[0, 2] = cell[0, 1] + cell[1, 0] = cell[0, 2] + cell[1, 1] = cell[0, 0] + cell[1, 2] = cell[0, 1] + cell[2, 0] = cell[0, 1] + cell[2, 1] = cell[0, 2] + cell[2, 2] = cell[0, 0] + elif ibrav == 6: + cell[0, 0] = a + cell[1, 1] = a + cell[2, 2] = a * celldm[2] + elif ibrav == 7: + cell[1, 0] = a / 2.0 + cell[1, 1] = cell[1, 0] + cell[1, 2] = celldm[2] * a / 2.0 + cell[0, 0] = cell[1, 0] + cell[0, 1] = -cell[1, 0] + cell[0, 2] = cell[1, 2] + cell[2, 0] = -cell[1, 0] + cell[2, 1] = -cell[1, 0] + cell[2, 2] = cell[1, 2] + elif ibrav == 8: + cell[0, 0] = a + cell[1, 1] = a * celldm[1] + cell[2, 2] = a * celldm[2] + elif ibrav == 9: + cell[0, 0] = a / 2.0 + cell[0, 1] = cell[0, 0] * celldm[1] + cell[1, 0] = -cell[0, 0] + cell[1, 1] = cell[0, 1] + cell[2, 2] = a * celldm[2] + elif ibrav == -9: + cell[0, 0] = a / 2.0 + cell[0, 1] = -cell[0, 0] * celldm[1] + cell[1, 0] = cell[0, 0] + cell[1, 1] = -cell[0, 1] + cell[2, 2] = a * celldm[2] + elif ibrav == 91: + cell[0, 0] = a + cell[1, 1] = a * celldm[1] / 2.0 + cell[1, 2] = -a * celldm[2] / 2.0 + cell[2, 1] = cell[1, 1] + cell[2, 2] = -cell[1, 2] + elif ibrav == 10: + cell[1, 0] = a / 2.0 + cell[1, 1] = cell[1, 0] * celldm[1] + cell[0, 0] = cell[1, 0] + cell[0, 2] = cell[1, 0] * celldm[2] + cell[2, 1] = cell[1, 0] * celldm[1] + cell[2, 2] = cell[0, 2] + elif ibrav == 11: + cell[0, 0] = a / 2.0 + cell[0, 1] = cell[0, 0] * celldm[1] + cell[0, 2] = cell[0, 0] * celldm[2] + cell[1, 0] = -cell[0, 0] + cell[1, 1] = cell[0, 1] + cell[1, 2] = cell[0, 2] + cell[2, 0] = -cell[0, 0] + cell[2, 1] = -cell[0, 1] + cell[2, 2] = cell[0, 2] + elif ibrav == 12: + sen = np.sqrt(1.0 - celldm[3] ** 2) + cell[0, 0] = a + cell[1, 0] = a * celldm[1] * celldm[3] + cell[1, 1] = a * celldm[1] * sen + cell[2, 2] = a * celldm[2] + elif ibrav == -12: + sen = np.sqrt(1.0 - celldm[4] ** 2) + cell[0, 0] = a + cell[1, 1] = a * celldm[1] + cell[2, 0] = a * celldm[2] * celldm[4] + cell[2, 2] = a * celldm[2] * sen + elif ibrav == 13: + sen = np.sqrt(1.0 - celldm[3] ** 2) + cell[0, 0] = a / 2.0 + cell[0, 2] = -cell[0, 0] * celldm[2] + cell[1, 0] = a * celldm[1] * celldm[3] + cell[1, 1] = a * celldm[1] * sen + cell[2, 0] = cell[0, 0] + cell[2, 2] = -cell[0, 2] + elif ibrav == -13: + sen = np.sqrt(1.0 - celldm[4] ** 2) + cell[0, 0] = a / 2.0 + cell[0, 1] = cell[0, 0] * celldm[1] + cell[1, 0] = -cell[0, 0] + cell[1, 1] = cell[0, 1] + cell[2, 0] = a * celldm[2] * celldm[4] + cell[2, 2] = a * celldm[2] * sen + elif ibrav == 14: + singam = np.sqrt(1.0 - celldm[5] ** 2) + term = ( + 1.0 + + 2.0 * celldm[3] * celldm[4] * celldm[5] + - celldm[3] ** 2 + - celldm[4] ** 2 + - celldm[5] ** 2 ) - return np.eye(3) - # raise RuntimeError('unsupported ibrav ' + str(ibrav)) + if term < 0.0: + raise RuntimeError( + "celldm do not make sense, check your data " + f"(ibrav=14 with cosAB={celldm[5]}, cosAC={celldm[4]}, " + f"cosBC={celldm[3]})" + ) + cell[0, 0] = a + cell[1, 0] = a * celldm[1] * celldm[5] + cell[1, 1] = a * celldm[1] * singam + cell[2, 0] = a * celldm[2] * celldm[4] + cell[2, 1] = a * celldm[2] * (celldm[3] - celldm[4] * celldm[5]) / singam + cell[2, 2] = a * celldm[2] * np.sqrt(term / (1.0 - celldm[5] ** 2)) + else: + raise RuntimeError(f"nonexistent bravais lattice {ibrav}") + return cell def load_cell_parameters(lines, lattice_parameter=None): @@ -137,6 +289,78 @@ def load_celldm(lines): return celldm +def _float_after(lines, name): + """Return the QE namelist value assigned to ``name`` in ``&SYSTEM``. + + The regex keeps ``name`` a whole word so ``B`` does not match inside + ``celldm`` or ``cosBC``, and skips comment text and content after the + closing ``/``. + """ + pattern = re.compile( + rf"\b{re.escape(name)}\s*=\s*({_QE_FLOAT_PATTERN})", re.IGNORECASE + ) + in_system = False + for raw_line in lines: + line = raw_line.split("!", 1)[0] + if re.search(r"&SYSTEM\b", line, re.IGNORECASE): + in_system = True + if not in_system: + continue + matched = pattern.search(line) + if matched is not None: + return float(matched.group(1).replace("d", "e").replace("D", "E")) + if "/" in line: + break + return None + + +def abc_to_celldm(ibrav, a=None, b=None, c=None, cosab=None, cosac=None, cosbc=None): + """Convert QE A/B/C (Angstrom) + cosines to ``celldm(1..6)`` in Bohr. + + Mirrors ``abc2celldm`` in QE's ``Modules/latgen.f90``: the angle-cosine + slots (4..6) are ibrav-dependent because each lattice fixes a different + axis pair. + """ + if a is None: + raise RuntimeError( + f"A (or celldm(1)) must be defined when ibrav={ibrav} is nonzero" + ) + if a <= 0: + raise RuntimeError(f"A must be positive (got {a})") + b_ratio = b / a if b is not None else 0.0 + c_ratio = c / a if c is not None else 0.0 + if ibrav in (14, 0): + c4, c5, c6 = cosbc or 0.0, cosac or 0.0, cosab or 0.0 + elif ibrav in (-12, -13): + c4, c5, c6 = 0.0, cosac or 0.0, 0.0 + elif ibrav in (5, -5, 12, 13): + c4, c5, c6 = cosab or 0.0, 0.0, 0.0 + else: + c4 = c5 = c6 = 0.0 + return np.array([a / length_convert, b_ratio, c_ratio, c4, c5, c6]) + + +def resolve_celldm(lines, ibrav, celldm): + """Return a valid ``celldm`` array, deriving it from A/B/C when unset. + + QE accepts either ``celldm(1..6)`` or the crystallographic ``A, B, C, + cosAB, cosAC, cosBC`` form (A/B/C in Angstrom). Both must not be mixed. + """ + if celldm[0] != 0.0: + if _float_after(lines, "A") is not None: + raise ValueError("both A and celldm(1) define the QE lattice parameter") + return celldm + a = _float_after(lines, "A") + if a is None: + return celldm + b = _float_after(lines, "B") + c = _float_after(lines, "C") + cosab = _float_after(lines, "cosAB") + cosac = _float_after(lines, "cosAC") + cosbc = _float_after(lines, "cosBC") + return abc_to_celldm(ibrav, a, b, c, cosab, cosac, cosbc) + + def load_lattice_parameter(lines, celldm): """Return QE's ``alat`` in angstrom, rejecting conflicting definitions.""" a_value = None @@ -188,6 +412,7 @@ def load_param_file(fname: FileType): cell = load_cell_parameters(lines, lattice_parameter) else: # celldm and cells reconstructed from it are expressed in Bohr. + celldm = resolve_celldm(lines, ibrav, celldm) cell = convert_celldm(ibrav, celldm) * length_convert # print(atom_names) # print(atom_numbs) diff --git a/tests/test_qe_cp_traj.py b/tests/test_qe_cp_traj.py index ce932ea9e..c2ba9289f 100644 --- a/tests/test_qe_cp_traj.py +++ b/tests/test_qe_cp_traj.py @@ -138,12 +138,79 @@ def test_missing_cell_parameters_for_ibrav_zero_raises(self): class TestConverCellDim(unittest.TestCase): - def test_case_null(self): - cell = convert_celldm(8, [1, 1, 1]) - ref = np.eye(3) - for ii in range(3): - for jj in range(3): - self.assertAlmostEqual(cell[ii][jj], ref[ii][jj]) + # Reference cells follow latgen in QE's Modules/latgen.f90. They were + # checked against each lattice family's lengths, angles, and volume. + CELLDIM = np.array([4.0, 1.5, 2.0, 0.5, 0.4, 0.6]) + GOLDEN = { + 1: [[4, 0, 0], [0, 4, 0], [0, 0, 4]], + 2: [[-2, 0, 2], [0, 2, 2], [-2, 2, 0]], + 3: [[2, 2, 2], [-2, 2, 2], [-2, -2, 2]], + -3: [[-2, 2, 2], [2, -2, 2], [2, 2, -2]], + 4: [[4, 0, 0], [-2, 2 * np.sqrt(3), 0], [0, 0, 8]], + 5: [ + [2, -2 / np.sqrt(3), 4 * np.sqrt(2 / 3)], + [0, 4 / np.sqrt(3), 4 * np.sqrt(2 / 3)], + [-2, -2 / np.sqrt(3), 4 * np.sqrt(2 / 3)], + ], + -5: [ + [0, 2 * np.sqrt(2), 2 * np.sqrt(2)], + [2 * np.sqrt(2), 0, 2 * np.sqrt(2)], + [2 * np.sqrt(2), 2 * np.sqrt(2), 0], + ], + 6: [[4, 0, 0], [0, 4, 0], [0, 0, 8]], + 7: [[2, -2, 4], [2, 2, 4], [-2, -2, 4]], + 8: [[4, 0, 0], [0, 6, 0], [0, 0, 8]], + 9: [[2, 3, 0], [-2, 3, 0], [0, 0, 8]], + -9: [[2, -3, 0], [2, 3, 0], [0, 0, 8]], + 91: [[4, 0, 0], [0, 3, -4], [0, 3, 4]], + 10: [[2, 0, 4], [2, 3, 0], [0, 3, 4]], + 11: [[2, 3, 4], [-2, 3, 4], [-2, -3, 4]], + 12: [[4, 0, 0], [3, 3 * np.sqrt(3), 0], [0, 0, 8]], + -12: [[4, 0, 0], [0, 6, 0], [3.2, 0, 8 * np.sqrt(1 - 0.4**2)]], + 13: [[2, 0, -4], [3, 3 * np.sqrt(3), 0], [2, 0, 4]], + -13: [[2, 3, 0], [-2, 3, 0], [3.2, 0, 8 * np.sqrt(1 - 0.4**2)]], + 14: [[4, 0, 0], [3.6, 4.8, 0], [3.2, 2.6, 6.8556546004]], + } + + def test_all_ibrav(self): + for ibrav, expected in self.GOLDEN.items(): + with self.subTest(ibrav=ibrav): + cell = convert_celldm(ibrav, self.CELLDIM) + np.testing.assert_allclose(cell, expected, atol=1e-7) + + def test_lengths_angles_and_volume(self): + # (lens, sorted cos of the three pairs, |det|), independent of axis choice. + expected = { + 1: ([4, 4, 4], [0, 0, 0], 64), + 2: ([4 / np.sqrt(2)] * 3, [0.5, 0.5, 0.5], 16), + 4: ([4, 4, 8], [-0.5, 0, 0], 0.5 * np.sqrt(3) * 4 * 4 * 8), + 5: ([4, 4, 4], [0.5, 0.5, 0.5], 45.25483), + 8: ([4, 6, 8], [0, 0, 0], 192), + 12: ([4, 6, 8], [0, 0, 0.5], 192 * np.sqrt(0.75)), + -12: ([4, 6, 8], [0, 0.4, 0], 192 * np.sqrt(1 - 0.4**2)), + 14: ([4, 6, 8], [0.4, 0.5, 0.6], 131.6286), + } + for ibrav, (lens, cos, vol) in expected.items(): + with self.subTest(ibrav=ibrav): + cell = convert_celldm(ibrav, self.CELLDIM) + np.testing.assert_allclose( + np.linalg.norm(cell, axis=1), lens, atol=1e-6 + ) + pairs = [ + np.dot(cell[i], cell[j]) + / np.linalg.norm(cell[i]) + / np.linalg.norm(cell[j]) + for i in range(3) + for j in range(i + 1, 3) + ] + np.testing.assert_allclose(sorted(pairs), sorted(cos), atol=1e-6) + self.assertAlmostEqual(abs(np.linalg.det(cell)), vol, places=4) + + def test_unsupported_ibrav_raises(self): + for ibrav in (0, 15, -1, 99): + with self.subTest(ibrav=ibrav): + with self.assertRaises(RuntimeError): + convert_celldm(ibrav, self.CELLDIM) class TestVirial(unittest.TestCase): From a7a0099e26a57074a40c65ea6f80ca490ad7c105 Mon Sep 17 00:00:00 2001 From: guhou-hvi <1322077960@qq.com> Date: Sat, 19 Sep 2026 13:49:37 +0800 Subject: [PATCH 2/2] fix(qe): validate celldm parameters per ibrav in convert_celldm --- dpdata/formats/qe/traj.py | 107 ++++++++++++++++++++++++++++++++++++-- tests/test_qe_cp_traj.py | 64 ++++++++++++++++++++++- 2 files changed, 166 insertions(+), 5 deletions(-) diff --git a/dpdata/formats/qe/traj.py b/dpdata/formats/qe/traj.py index 9355cd528..e0e68346d 100644 --- a/dpdata/formats/qe/traj.py +++ b/dpdata/formats/qe/traj.py @@ -63,6 +63,8 @@ def convert_celldm(ibrav, celldm): sr2 = np.sqrt(2.0) sr3 = np.sqrt(3.0) a = celldm[0] + if not a > 0.0: + raise RuntimeError(f"wrong celldm(1)={a}: must be positive for ibrav={ibrav}") cell = np.zeros((3, 3)) if ibrav == 1: @@ -89,12 +91,20 @@ def convert_celldm(ibrav, celldm): cell[2, 0] = -t cell[2, 1] = -t elif ibrav == 4: + if not celldm[2] > 0.0: + raise RuntimeError( + f"wrong celldm(3)={celldm[2]} for ibrav=4: c/a must be positive" + ) cell[0, 0] = a cell[1, 0] = -a / 2.0 cell[1, 1] = a * sr3 / 2.0 cell[2, 2] = a * celldm[2] elif ibrav == 5: cosab = celldm[3] + if not -0.5 < cosab < 1.0: + raise RuntimeError( + f"wrong celldm(4)={cosab} for ibrav=5: need -0.5 < cosAB < 1" + ) term1 = np.sqrt(1.0 + 2.0 * cosab) term2 = np.sqrt(1.0 - cosab) cell[1, 1] = sr2 * a * term2 / sr3 @@ -107,6 +117,10 @@ def convert_celldm(ibrav, celldm): cell[2, 2] = cell[1, 2] elif ibrav == -5: cosab = celldm[3] + if not -0.5 < cosab < 1.0: + raise RuntimeError( + f"wrong celldm(4)={cosab} for ibrav=-5: need -0.5 < cosAB < 1" + ) term1 = np.sqrt(1.0 + 2.0 * cosab) term2 = np.sqrt(1.0 - cosab) cell[0, 0] = a * (term1 - 2.0 * term2) / 3.0 @@ -119,10 +133,18 @@ def convert_celldm(ibrav, celldm): cell[2, 1] = cell[0, 2] cell[2, 2] = cell[0, 0] elif ibrav == 6: + if not celldm[2] > 0.0: + raise RuntimeError( + f"wrong celldm(3)={celldm[2]} for ibrav=6: c/a must be positive" + ) cell[0, 0] = a cell[1, 1] = a cell[2, 2] = a * celldm[2] elif ibrav == 7: + if not celldm[2] > 0.0: + raise RuntimeError( + f"wrong celldm(3)={celldm[2]} for ibrav=7: c/a must be positive" + ) cell[1, 0] = a / 2.0 cell[1, 1] = cell[1, 0] cell[1, 2] = celldm[2] * a / 2.0 @@ -133,28 +155,53 @@ def convert_celldm(ibrav, celldm): cell[2, 1] = -cell[1, 0] cell[2, 2] = cell[1, 2] elif ibrav == 8: + if not (celldm[1] > 0.0 and celldm[2] > 0.0): + raise RuntimeError( + f"wrong celldm(2)={celldm[1]} or celldm(3)={celldm[2]} " + f"for ibrav=8: b/a and c/a must be positive" + ) cell[0, 0] = a cell[1, 1] = a * celldm[1] cell[2, 2] = a * celldm[2] elif ibrav == 9: + if not (celldm[1] > 0.0 and celldm[2] > 0.0): + raise RuntimeError( + f"wrong celldm(2)={celldm[1]} or celldm(3)={celldm[2]} " + f"for ibrav=9: b/a and c/a must be positive" + ) cell[0, 0] = a / 2.0 cell[0, 1] = cell[0, 0] * celldm[1] cell[1, 0] = -cell[0, 0] cell[1, 1] = cell[0, 1] cell[2, 2] = a * celldm[2] elif ibrav == -9: + if not (celldm[1] > 0.0 and celldm[2] > 0.0): + raise RuntimeError( + f"wrong celldm(2)={celldm[1]} or celldm(3)={celldm[2]} " + f"for ibrav=-9: b/a and c/a must be positive" + ) cell[0, 0] = a / 2.0 cell[0, 1] = -cell[0, 0] * celldm[1] cell[1, 0] = cell[0, 0] cell[1, 1] = -cell[0, 1] cell[2, 2] = a * celldm[2] elif ibrav == 91: + if not (celldm[1] > 0.0 and celldm[2] > 0.0): + raise RuntimeError( + f"wrong celldm(2)={celldm[1]} or celldm(3)={celldm[2]} " + f"for ibrav=91: b/a and c/a must be positive" + ) cell[0, 0] = a cell[1, 1] = a * celldm[1] / 2.0 cell[1, 2] = -a * celldm[2] / 2.0 cell[2, 1] = cell[1, 1] cell[2, 2] = -cell[1, 2] elif ibrav == 10: + if not (celldm[1] > 0.0 and celldm[2] > 0.0): + raise RuntimeError( + f"wrong celldm(2)={celldm[1]} or celldm(3)={celldm[2]} " + f"for ibrav=10: b/a and c/a must be positive" + ) cell[1, 0] = a / 2.0 cell[1, 1] = cell[1, 0] * celldm[1] cell[0, 0] = cell[1, 0] @@ -162,6 +209,11 @@ def convert_celldm(ibrav, celldm): cell[2, 1] = cell[1, 0] * celldm[1] cell[2, 2] = cell[0, 2] elif ibrav == 11: + if not (celldm[1] > 0.0 and celldm[2] > 0.0): + raise RuntimeError( + f"wrong celldm(2)={celldm[1]} or celldm(3)={celldm[2]} " + f"for ibrav=11: b/a and c/a must be positive" + ) cell[0, 0] = a / 2.0 cell[0, 1] = cell[0, 0] * celldm[1] cell[0, 2] = cell[0, 0] * celldm[2] @@ -172,18 +224,45 @@ def convert_celldm(ibrav, celldm): cell[2, 1] = -cell[0, 1] cell[2, 2] = cell[0, 2] elif ibrav == 12: + if not (celldm[1] > 0.0 and celldm[2] > 0.0): + raise RuntimeError( + f"wrong celldm(2)={celldm[1]} or celldm(3)={celldm[2]} " + f"for ibrav=12: b/a and c/a must be positive" + ) + if not -1.0 < celldm[3] < 1.0: + raise RuntimeError( + f"wrong celldm(4)={celldm[3]} for ibrav=12: |cosAB| < 1 required" + ) sen = np.sqrt(1.0 - celldm[3] ** 2) cell[0, 0] = a cell[1, 0] = a * celldm[1] * celldm[3] cell[1, 1] = a * celldm[1] * sen cell[2, 2] = a * celldm[2] elif ibrav == -12: + if not (celldm[1] > 0.0 and celldm[2] > 0.0): + raise RuntimeError( + f"wrong celldm(2)={celldm[1]} or celldm(3)={celldm[2]} " + f"for ibrav=-12: b/a and c/a must be positive" + ) + if not -1.0 < celldm[4] < 1.0: + raise RuntimeError( + f"wrong celldm(5)={celldm[4]} for ibrav=-12: |cosAC| < 1 required" + ) sen = np.sqrt(1.0 - celldm[4] ** 2) cell[0, 0] = a cell[1, 1] = a * celldm[1] cell[2, 0] = a * celldm[2] * celldm[4] cell[2, 2] = a * celldm[2] * sen elif ibrav == 13: + if not (celldm[1] > 0.0 and celldm[2] > 0.0): + raise RuntimeError( + f"wrong celldm(2)={celldm[1]} or celldm(3)={celldm[2]} " + f"for ibrav=13: b/a and c/a must be positive" + ) + if not -1.0 < celldm[3] < 1.0: + raise RuntimeError( + f"wrong celldm(4)={celldm[3]} for ibrav=13: |cosAB| < 1 required" + ) sen = np.sqrt(1.0 - celldm[3] ** 2) cell[0, 0] = a / 2.0 cell[0, 2] = -cell[0, 0] * celldm[2] @@ -192,6 +271,15 @@ def convert_celldm(ibrav, celldm): cell[2, 0] = cell[0, 0] cell[2, 2] = -cell[0, 2] elif ibrav == -13: + if not (celldm[1] > 0.0 and celldm[2] > 0.0): + raise RuntimeError( + f"wrong celldm(2)={celldm[1]} or celldm(3)={celldm[2]} " + f"for ibrav=-13: b/a and c/a must be positive" + ) + if not -1.0 < celldm[4] < 1.0: + raise RuntimeError( + f"wrong celldm(5)={celldm[4]} for ibrav=-13: |cosAC| < 1 required" + ) sen = np.sqrt(1.0 - celldm[4] ** 2) cell[0, 0] = a / 2.0 cell[0, 1] = cell[0, 0] * celldm[1] @@ -200,6 +288,17 @@ def convert_celldm(ibrav, celldm): cell[2, 0] = a * celldm[2] * celldm[4] cell[2, 2] = a * celldm[2] * sen elif ibrav == 14: + if not (celldm[1] > 0.0 and celldm[2] > 0.0): + raise RuntimeError( + f"wrong celldm(2)={celldm[1]} or celldm(3)={celldm[2]} " + f"for ibrav=14: b/a and c/a must be positive" + ) + for _idx, _key in ((3, "cosBC"), (4, "cosAC"), (5, "cosAB")): + if not -1.0 < celldm[_idx] < 1.0: + raise RuntimeError( + f"wrong celldm({_idx + 1})={celldm[_idx]} for ibrav=14: " + f"|{_key}| < 1 required" + ) singam = np.sqrt(1.0 - celldm[5] ** 2) term = ( 1.0 @@ -208,11 +307,11 @@ def convert_celldm(ibrav, celldm): - celldm[4] ** 2 - celldm[5] ** 2 ) - if term < 0.0: + if term <= 0.0: raise RuntimeError( - "celldm do not make sense, check your data " - f"(ibrav=14 with cosAB={celldm[5]}, cosAC={celldm[4]}, " - f"cosBC={celldm[3]})" + "celldm do not make sense for ibrav=14, check your data " + f"(cosBC={celldm[3]}, cosAC={celldm[4]}, cosAB={celldm[5]} " + "yield a degenerate or invalid cell)" ) cell[0, 0] = a cell[1, 0] = a * celldm[1] * celldm[5] diff --git a/tests/test_qe_cp_traj.py b/tests/test_qe_cp_traj.py index c2ba9289f..33905074e 100644 --- a/tests/test_qe_cp_traj.py +++ b/tests/test_qe_cp_traj.py @@ -169,7 +169,19 @@ class TestConverCellDim(unittest.TestCase): -12: [[4, 0, 0], [0, 6, 0], [3.2, 0, 8 * np.sqrt(1 - 0.4**2)]], 13: [[2, 0, -4], [3, 3 * np.sqrt(3), 0], [2, 0, 4]], -13: [[2, 3, 0], [-2, 3, 0], [3.2, 0, 8 * np.sqrt(1 - 0.4**2)]], - 14: [[4, 0, 0], [3.6, 4.8, 0], [3.2, 2.6, 6.8556546004]], + 14: [ + [4, 0, 0], + [3.6, 4.8, 0], + [ + 3.2, + 2.6, + 8 + * np.sqrt( + (1.0 + 2.0 * 0.5 * 0.4 * 0.6 - 0.5**2 - 0.4**2 - 0.6**2) + / (1.0 - 0.6**2) + ), + ], + ], } def test_all_ibrav(self): @@ -212,6 +224,56 @@ def test_unsupported_ibrav_raises(self): with self.assertRaises(RuntimeError): convert_celldm(ibrav, self.CELLDIM) + def test_invalid_celldm_raises(self): + # (ibrav, celldm index to invalidate, replacement, QE-slot label) + bad_inputs = [ + (1, 0, 0.0, "celldm(1)=0"), + (1, 0, -1.0, "celldm(1)<0"), + (4, 2, 0.0, "ibrav=4 celldm(3)<=0"), + (4, 2, -1.0, "ibrav=4 celldm(3)<0"), + (5, 3, -0.5, "ibrav=5 cosAB=-0.5"), + (5, 3, 1.0, "ibrav=5 cosAB=1"), + (-5, 3, 1.5, "ibrav=-5 cosAB out of range"), + (6, 2, 0.0, "ibrav=6 celldm(3)<=0"), + (7, 2, -0.1, "ibrav=7 celldm(3)<0"), + (8, 1, 0.0, "ibrav=8 celldm(2)=0"), + (8, 2, -1.0, "ibrav=8 celldm(3)<0"), + (9, 1, 0.0, "ibrav=9 celldm(2)=0"), + (-9, 2, 0.0, "ibrav=-9 celldm(3)=0"), + (91, 1, 0.0, "ibrav=91 celldm(2)=0"), + (10, 2, 0.0, "ibrav=10 celldm(3)=0"), + (11, 1, 0.0, "ibrav=11 celldm(2)=0"), + (12, 3, 1.0, "ibrav=12 cosAB=+1 (degenerate sen)"), + (12, 3, -1.0, "ibrav=12 cosAB=-1"), + (12, 1, 0.0, "ibrav=12 celldm(2)=0"), + (-12, 4, 1.0, "ibrav=-12 cosAC=+1"), + (-12, 2, -1.0, "ibrav=-12 celldm(3)<0"), + (13, 3, 1.0, "ibrav=13 cosAB=+1"), + (-13, 4, -1.0, "ibrav=-13 cosAC=-1"), + (14, 1, 0.0, "ibrav=14 celldm(2)=0"), + (14, 2, 0.0, "ibrav=14 celldm(3)=0"), + (14, 3, 1.0, "ibrav=14 cosBC=+1"), + (14, 4, -1.0, "ibrav=14 cosAC=-1"), + (14, 5, 1.0, "ibrav=14 cosAB=+1 (singam=0)"), + ] + for ibrav, idx, value, label in bad_inputs: + with self.subTest(case=label): + celldm = self.CELLDIM.copy() + celldm[idx] = value + with self.assertRaises(RuntimeError): + convert_celldm(ibrav, celldm) + + def test_ibrav_14_degenerate_term_rejects_zero(self): + # cosBC=cosAC=cosAB=0.5 gives 1 + 2*(0.5)**3 - 3*(0.5)**2 = 0.875 > 0 (valid) + # But cosAB=0, cosAC=1-eps, cosBC=1-eps pushes term → 0 exactly at boundary + celldm = self.CELLDIM.copy() + celldm[3] = 0.0 # cosBC + celldm[4] = 1.0 - 1e-9 # cosAC (just inside range) + celldm[5] = 1.0 - 1e-9 # cosAB (just inside range) + # term = 1 + 0 - 0 - (1-1e-9)^2 - (1-1e-9)^2 ≈ -1 + 4e-9 → still rejected + with self.assertRaises(RuntimeError): + convert_celldm(14, celldm) + class TestVirial(unittest.TestCase): def test(self):