fix(abacus): transform forces for left-handed SCF cells

dpdata/formats/abacus/scf.py

@@ -221,7 +221,7 @@ def get_frame(fname):
         return data
 
     # read STRU file
-    data = get_frame_from_stru(geometry_path_in)
+    data, cell_flipped = get_frame_from_stru(geometry_path_in, return_flipped=True)
     natoms = sum(data["atom_numbs"])
     # should remove spins from STRU file
     if "spins" in data:
@@ -236,6 +236,8 @@ def get_frame(fname):
         magforce = magforce[-1:]
 
     force = get_force(outlines, natoms)
+    if cell_flipped and force is not None:
+        force = -force
     stress = get_stress(outlines)
 
     data["energies"] = np.array(energy)[np.newaxis]

dpdata/formats/abacus/stru.py

@@ -482,13 +482,15 @@ def validate_duplicate_species(atom_names, masses, pp_files, orb_files):
                 )
 
 
-def get_frame_from_stru(stru):
+def get_frame_from_stru(stru, return_flipped=False):
     """Read the ABACUS STRU file and return the dpdata frame.
 
     The description of ABACUS STRU can be found in https://abacus.deepmodeling.com/en/latest/advanced/input_files/stru.html
 
     Args:
         stru (str): path to the ABACUS STRU file.
+        return_flipped (bool): whether to return a flag indicating that the
+            input cell was left-handed and had to be reflected.
 
     Returns
     -------
@@ -529,6 +531,7 @@ def get_frame_from_stru(stru):
     celldm = parse_lattice_constant_block(blocks["LATTICE_CONSTANT"])
     cell = parse_lattice_vectors_block(blocks["LATTICE_VECTORS"])
     cell = np.array(cell) * celldm * bohr2ang
+    cell_flipped = np.linalg.det(cell) < 0
     atom_numbs, coords, move, mags, velocity, sc, lambda_ = parse_pos(
         blocks["ATOMIC_POSITIONS"], atom_names, celldm, cell
     )
@@ -573,6 +576,8 @@ def get_frame_from_stru(stru):
         data["dpks_descriptor"] = dpks_descriptor[0].strip()
     if len(move) > 0:
         data["move"] = np.array([move])
+    if return_flipped:
+        return data, cell_flipped
     return data
 
 

tests/test_abacus_pw_scf.py

@@ -129,6 +129,27 @@ def test_energy(self):
         # ref_energy = -30007.651851226798
         # self.assertAlmostEqual(self.system_h2o.data['energies'][0], ref_energy)
 
+    def test_left_handed_cell_transforms_forces(self):
+        # ABACUS accepts left-handed cells. The reader reflects both the cell
+        # and coordinates to satisfy dpdata's right-handed-cell convention, so
+        # force vectors must receive the same transformation.
+        with tempfile.TemporaryDirectory() as tmp_dir:
+            shutil.copy("abacus.scf/INPUT.ok", os.path.join(tmp_dir, "INPUT"))
+            shutil.copytree("abacus.scf/OUT.ch4", os.path.join(tmp_dir, "OUT.ch4"))
+            with open("abacus.scf/STRU.ch4") as fp:
+                stru = fp.read().replace("0.0 0.0 1", "0.0 0.0 -1")
+            with open(os.path.join(tmp_dir, "STRU.ch4"), "w") as fp:
+                fp.write(stru)
+
+            left_handed = dpdata.LabeledSystem(tmp_dir, fmt="abacus/scf")
+
+            np.testing.assert_allclose(
+                left_handed.data["coords"], -self.system_ch4.data["coords"]
+            )
+            np.testing.assert_allclose(
+                left_handed.data["forces"], -self.system_ch4.data["forces"]
+            )
+
 
 class TestABACUSLabeledOutputFail(unittest.TestCase):
     def setUp(self):