diff --git a/dpdata/formats/qe/traj.py b/dpdata/formats/qe/traj.py index 382d3acf..4d171529 100644 --- a/dpdata/formats/qe/traj.py +++ b/dpdata/formats/qe/traj.py @@ -66,11 +66,30 @@ def convert_celldm(ibrav, celldm): def load_cell_parameters(lines): - blk = load_block(lines, "CELL_PARAMETERS", 3) - ret = [] - for ii in blk: - ret.append([float(jj) for jj in ii.split()[0:3]]) - return np.array(ret) + """Load ``CELL_PARAMETERS`` and convert its vectors to angstrom. + + CP trajectory ``.cel`` files always use atomic units, but the fallback + cell in the QE input file follows the unit declared on the + ``CELL_PARAMETERS`` card. Keeping that distinction here prevents an + angstrom input cell from being converted a second time when no ``.cel`` + file is available. + """ + for idx, line in enumerate(lines): + # Ignore commented examples such as ``!CELL_PARAMETERS {bohr}``, which + # are common in QE inputs and must not shadow the active card below. + card = line.split("!", 1)[0].strip() + if card.upper().startswith("CELL_PARAMETERS"): + blk = lines[idx + 1 : idx + 4] + break + else: + raise ValueError("CELL_PARAMETERS is required when ibrav is 0") + + cell = np.array([[float(value) for value in row.split()[:3]] for row in blk]) + # Preserve the historical Bohr conversion for the default/Bohr spelling. + # An explicit angstrom card, however, is already in dpdata's length unit. + if "angstrom" not in card.lower(): + cell *= length_convert + return cell def load_atom_names(lines, ntypes): @@ -111,8 +130,8 @@ def load_param_file(fname: FileType): if ibrav == 0: cell = load_cell_parameters(lines) else: - cell = convert_celldm(ibrav, celldm) - cell = cell * length_convert + # celldm and cells reconstructed from it are expressed in Bohr. + cell = convert_celldm(ibrav, celldm) * length_convert # print(atom_names) # print(atom_numbs) # print(atom_types) diff --git a/tests/qe.traj/angstrom_no_cel/cp.in b/tests/qe.traj/angstrom_no_cel/cp.in new file mode 100644 index 00000000..72fb8a89 --- /dev/null +++ b/tests/qe.traj/angstrom_no_cel/cp.in @@ -0,0 +1,20 @@ +&CONTROL + calculation = 'cp', + prefix = 'cp', +/ +&SYSTEM + ibrav = 0, + nat = 1, + ntyp = 1, +/ + +CELL_PARAMETERS { angstrom } + 19.7299995422 0.0000000000 0.0000000000 + 0.0000000000 19.7299995422 0.0000000000 + 0.0000000000 0.0000000000 19.7299995422 + +ATOMIC_SPECIES +H 1.00794 H.UPF + +ATOMIC_POSITIONS { bohr } +H 0.0 0.0 0.0 diff --git a/tests/qe.traj/angstrom_no_cel/cp.pos b/tests/qe.traj/angstrom_no_cel/cp.pos new file mode 100644 index 00000000..c1b958f7 --- /dev/null +++ b/tests/qe.traj/angstrom_no_cel/cp.pos @@ -0,0 +1,2 @@ +0 0.0 +0.0 0.0 0.0 diff --git a/tests/qe.traj/bohr_no_cel/cp.in b/tests/qe.traj/bohr_no_cel/cp.in new file mode 100644 index 00000000..88245ccb --- /dev/null +++ b/tests/qe.traj/bohr_no_cel/cp.in @@ -0,0 +1,20 @@ +&CONTROL + calculation = 'cp', + prefix = 'cp', +/ +&SYSTEM + ibrav = 0, + nat = 1, + ntyp = 1, +/ + +CELL_PARAMETERS { bohr } + 2.0 0.0 0.0 + 0.0 2.0 0.0 + 0.0 0.0 2.0 + +ATOMIC_SPECIES +H 1.00794 H.UPF + +ATOMIC_POSITIONS { bohr } +H 0.0 0.0 0.0 diff --git a/tests/qe.traj/bohr_no_cel/cp.pos b/tests/qe.traj/bohr_no_cel/cp.pos new file mode 100644 index 00000000..c1b958f7 --- /dev/null +++ b/tests/qe.traj/bohr_no_cel/cp.pos @@ -0,0 +1,2 @@ +0 0.0 +0.0 0.0 0.0 diff --git a/tests/qe.traj/missing_cell_no_cel/cp.in b/tests/qe.traj/missing_cell_no_cel/cp.in new file mode 100644 index 00000000..6c873b52 --- /dev/null +++ b/tests/qe.traj/missing_cell_no_cel/cp.in @@ -0,0 +1,15 @@ +&CONTROL + calculation = 'cp', + prefix = 'cp', +/ +&SYSTEM + ibrav = 0, + nat = 1, + ntyp = 1, +/ + +ATOMIC_SPECIES +H 1.00794 H.UPF + +ATOMIC_POSITIONS { bohr } +H 0.0 0.0 0.0 diff --git a/tests/test_qe_cp_traj.py b/tests/test_qe_cp_traj.py index a670bd4d..34a0e972 100644 --- a/tests/test_qe_cp_traj.py +++ b/tests/test_qe_cp_traj.py @@ -61,6 +61,43 @@ def setUp(self): self.system = dpdata.LabeledSystem("qe.traj/oh-md", fmt="qe/cp/traj") +class TestCPTRAJInputCellUnits(unittest.TestCase): + def test_angstrom_cell_without_cel_trajectory(self): + # Earlier tests did not expose this regression: the only missing-.cel + # fixture used ibrav/celldm, whose cell is correctly expressed in Bohr. + # Exercise the distinct fallback path where CELL_PARAMETERS already + # stores angstrom values and therefore must not be converted again. + system = dpdata.System( + "qe.traj/angstrom_no_cel/cp", + fmt="qe/cp/traj", + ) + + np.testing.assert_allclose( + system["cells"][0], + np.eye(3) * 19.7299995422, + ) + + def test_bohr_cell_without_cel_trajectory(self): + system = dpdata.System( + "qe.traj/bohr_no_cel/cp", + fmt="qe/cp/traj", + ) + + np.testing.assert_allclose( + system["cells"][0], + np.eye(3) * 2.0 * bohr2ang, + ) + + def test_missing_cell_parameters_for_ibrav_zero_raises(self): + with self.assertRaisesRegex( + ValueError, "CELL_PARAMETERS is required when ibrav is 0" + ): + dpdata.System( + "qe.traj/missing_cell_no_cel/cp", + fmt="qe/cp/traj", + ) + + class TestConverCellDim(unittest.TestCase): def test_case_null(self): cell = convert_celldm(8, [1, 1, 1])