From 4a49df45184297d42ade04a519b09f0c9a2e6d4a Mon Sep 17 00:00:00 2001
From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com>
Date: Mon, 24 Aug 2026 09:13:29 +0200
Subject: [PATCH] fix: write QMCFC dummy atom for vel, force and charge output

The QMCFC header always counts one extra atom, but the dummy X row
was only written for position trajectories. Velocity, force and
charge files therefore declared n_atoms+1 atoms while containing
only n_atoms rows per frame, so read_trajectory rejected them as
corrupt. The dummy row is now written for every trajectory type,
matching what the reader strips, so header and body agree.
---
 PQAnalysis/io/traj_file/trajectory_writer.py |  7 +++-
 tests/io/test_trajectoryWriter.py            | 44 ++++++++++++++++++++
 2 files changed, 50 insertions(+), 1 deletion(-)

diff --git a/PQAnalysis/io/traj_file/trajectory_writer.py b/PQAnalysis/io/traj_file/trajectory_writer.py
index e812c9b5..9840cbf7 100644
--- a/PQAnalysis/io/traj_file/trajectory_writer.py
+++ b/PQAnalysis/io/traj_file/trajectory_writer.py
@@ -236,7 +236,7 @@ def _write_xyz(self, xyz: Np2DNumberArray, atoms: List[Atom]) -> None:
             The elements of the frame.
         """
 
-        if self.format == MDEngineFormat.QMCFC and self._type == TrajectoryFormat.XYZ:
+        if self.format == MDEngineFormat.QMCFC:
             print("X   0.0 0.0 0.0", file=self.file)
 
         for i, atom in enumerate(atoms):
@@ -419,6 +419,8 @@ def _write_scalar(
         """
         Writes the charges of the frame to the file.
 
+        If format is 'qmcfc', an additional X 0.0 line is written.
+
         Parameters
         ----------
         scalar : np.array
@@ -427,6 +429,9 @@ def _write_scalar(
             The elements of the frame.
         """
 
+        if self.format == MDEngineFormat.QMCFC:
+            print("X 0.0", file=self.file)
+
         for i, atom in enumerate(atoms):
             print(f"{atom.name} {scalar[i]}", file=self.file)
 
diff --git a/tests/io/test_trajectoryWriter.py b/tests/io/test_trajectoryWriter.py
index e33fec43..e0f5c867 100644
--- a/tests/io/test_trajectoryWriter.py
+++ b/tests/io/test_trajectoryWriter.py
@@ -94,6 +94,50 @@ def test_write_extxyz_roundtrip():
     assert output_frame.cell == frame.cell
 
 
+@pytest.mark.usefixtures("tmpdir")
+@pytest.mark.parametrize("traj_type", ["xyz", "vel", "force", "charge"])
+def test_write_qmcfc_roundtrip(traj_type):
+    atoms = [Atom("h"), Atom("o")]
+    data = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]])
+    charges = np.array([1.0, 2.0])
+    frame = AtomicSystem(
+        atoms=atoms,
+        pos=data,
+        vel=data,
+        forces=data,
+        charges=charges,
+        cell=Cell(10, 10, 10),
+    )
+    traj = Trajectory([frame, frame])
+
+    filename = f"qmcfc_output.{traj_type}"
+    write_trajectory(
+        traj,
+        filename=filename,
+        engine_format="qmcfc",
+        traj_type=traj_type,
+    )
+
+    with open(filename, "r", encoding="utf-8") as file:
+        lines = file.read().splitlines()
+
+    assert lines[0].split()[0] == "3"
+    assert lines[2].split()[0] == "X"
+
+    output = read_trajectory(filename, md_format="qmcfc", traj_format=traj_type)
+
+    assert len(output) == 2
+    assert output[0].atoms == atoms
+    if traj_type == "xyz":
+        assert np.allclose(output[0].pos, data)
+    elif traj_type == "vel":
+        assert np.allclose(output[0].vel, data)
+    elif traj_type == "force":
+        assert np.allclose(output[0].forces, data)
+    else:
+        assert np.allclose(output[0].charges, charges)
+
+
 @pytest.mark.usefixtures("tmpdir")
 def test_write_extxyz_ase_profile_converts_energy_like_units():
     from ase.io import read as ase_read