bond_list parsing and test

Author: Bernadette-Mohr

src/nomad_simulation_parsers/parsers/gromacs/parser.py

@@ -539,17 +539,44 @@ def get_outputs(self) -> list[dict[str, Any]]:
             outputs.append(data)
         return outputs
 
+    def get_bond_list(self) -> np.ndarray | None:
+        """
+        Extract bond list from MDAnalysis interactions.
+        Returns numpy array of shape (n_bonds, 2) with atom indices.
+        """
+        result = None
+        interactions = self.data_object.get_interactions()
+        if not interactions:
+            return result
+
+        # Filter for bond interactions only
+        bonds = []
+        for interaction in interactions:
+            if interaction.get('type') == 'bond':
+                atom_indices = interaction.get('atom_indices')
+                if atom_indices is not None and len(atom_indices) == 2:
+                    bonds.append(list(atom_indices))
+
+        if bonds:
+            result = np.array(bonds, dtype=int)
+
+        return result
+
     def get_configurations(self) -> list[dict[str, Any]]:
         configurations = []
         for n, _ in enumerate(self._trajectory_steps_sampled):
-            configurations.append(
-                dict(
-                    labels=self.get_atom_labels(n),
-                    positions=self.data_object.get_positions(n),
-                    velocities=self.data_object.get_velocities(n),
-                    lattice_vectors=self.data_object.get_lattice_vectors(n),
-                )
+            config = dict(
+                labels=self.get_atom_labels(n),
+                positions=self.data_object.get_positions(n),
+                velocities=self.data_object.get_velocities(n),
+                lattice_vectors=self.data_object.get_lattice_vectors(n),
             )
+            # Add bond_list only to first configuration (topology is time-independent)
+            if n == 0:
+                bond_list = self.get_bond_list()
+                if bond_list is not None:
+                    config['bond_list'] = bond_list
+            configurations.append(config)
         return configurations
 
     def get_force_field_contributions(self) -> list[dict[str, Any]]:

src/nomad_simulation_parsers/schema_packages/gromacs.py

@@ -47,6 +47,7 @@ class AtomicCell(model_system.Representation):
 class ModelSystem(model_system.ModelSystem):
     add_mapping_annotation(model_system.ModelSystem.velocities, TPR_KEY, '.velocities')
     add_mapping_annotation(model_system.ModelSystem.positions, TPR_KEY, '.positions')
+    add_mapping_annotation(model_system.ModelSystem.bond_list, TPR_KEY, '.bond_list')
     add_mapping_annotation(model_system.AtomsState.m_def, TPR_KEY, '.labels')
     add_mapping_annotation(model_system.Representation.m_def, LOG_KEY, '.@')
     add_mapping_annotation(model_system.Representation.m_def, TPR_KEY, '.@')

tests/parsers/test_gromacs_parser.py

@@ -31,6 +31,10 @@ def get_atom_labels(self, idx):
         n = int(np.asarray(self._positions[idx]).shape[0])
         return ['H'] * n
 
+    def get_interactions(self):
+        # Return empty list for stub (no bonds by default)
+        return []
+
 
 @pytest.fixture
 def simple_mdanalysis_parser():
@@ -358,6 +362,148 @@ def test_get_coordinate_save_frequency():
     assert lp.get_coordinate_save_frequency(None) is None
 
 
+def test_get_bond_list():
+    """Test bond list extraction from MDAnalysis interactions."""
+    mdap = gromacs_parser.GromacsMDAnalysisParser()
+
+    # Test with bond interactions
+    mock_interactions = [
+        {'type': 'bond', 'atom_indices': [0, 1], 'atom_labels': ['O', 'H']},
+        {'type': 'bond', 'atom_indices': [0, 2], 'atom_labels': ['O', 'H']},
+        {'type': 'bond', 'atom_indices': [3, 4], 'atom_labels': ['O', 'H']},
+        {'type': 'angle', 'atom_indices': [1, 0, 2], 'atom_labels': ['H', 'O', 'H']},
+    ]
+
+    mdap.data_object = type(
+        'obj', (object,), {'get_interactions': lambda self: mock_interactions}
+    )()
+
+    bond_list = mdap.get_bond_list()
+
+    assert bond_list is not None
+    assert isinstance(bond_list, np.ndarray)
+    assert bond_list.shape == (3, 2)
+    assert np.array_equal(bond_list[0], [0, 1])
+    assert np.array_equal(bond_list[1], [0, 2])
+    assert np.array_equal(bond_list[2], [3, 4])
+
+
+def test_get_bond_list_no_bonds():
+    """Test bond list extraction when no bonds are present."""
+    mdap = gromacs_parser.GromacsMDAnalysisParser()
+
+    # Only angles, no bonds
+    mock_interactions = [
+        {'type': 'angle', 'atom_indices': [0, 1, 2], 'atom_labels': ['H', 'O', 'H']},
+    ]
+
+    mdap.data_object = type(
+        'obj', (object,), {'get_interactions': lambda self: mock_interactions}
+    )()
+
+    bond_list = mdap.get_bond_list()
+    assert bond_list is None
+
+
+def test_get_bond_list_empty_interactions():
+    """Test bond list extraction with empty interactions."""
+    mdap = gromacs_parser.GromacsMDAnalysisParser()
+
+    mdap.data_object = type('obj', (object,), {'get_interactions': lambda self: []})()
+
+    bond_list = mdap.get_bond_list()
+    assert bond_list is None
+
+
+def test_get_bond_list_invalid_bond_indices():
+    """Test bond list extraction filters out invalid bond entries."""
+    mdap = gromacs_parser.GromacsMDAnalysisParser()
+
+    # Mix of valid and invalid bond interactions
+    mock_interactions = [
+        {'type': 'bond', 'atom_indices': [0, 1], 'atom_labels': ['O', 'H']},
+        {'type': 'bond', 'atom_indices': None, 'atom_labels': ['O', 'H']},
+        {'type': 'bond', 'atom_indices': [2, 3, 4], 'atom_labels': ['O', 'H', 'C']},
+        {'type': 'bond', 'atom_indices': [4, 5], 'atom_labels': ['O', 'H']},
+    ]
+
+    mdap.data_object = type(
+        'obj', (object,), {'get_interactions': lambda self: mock_interactions}
+    )()
+
+    bond_list = mdap.get_bond_list()
+
+    assert bond_list is not None
+    assert bond_list.shape == (2, 2)
+    assert np.array_equal(bond_list[0], [0, 1])
+    assert np.array_equal(bond_list[1], [4, 5])
+
+
+def test_integration_bond_list_in_parsed_system():
+    """Test that bond_list is populated in parsed model_system from TPR file."""
+    base = Path(__file__).parent.parent / 'data' / 'gromacs' / 'water'
+    tpr_file = os.path.join(base, 'topol.tpr')
+
+    if not os.path.exists(tpr_file):
+        pytest.skip(f'TPR file not found: {tpr_file}')
+
+    archive = EntryArchive()
+    parser = gromacs_parser.GromacsParser()
+    parser.parse(tpr_file, archive)
+
+    assert archive.data is not None
+    assert len(archive.data.model_system) > 0
+
+    system = archive.data.model_system[0]
+    assert system.bond_list is not None, 'Bond list should be populated from TPR'
+    assert isinstance(system.bond_list, np.ndarray)
+    assert system.bond_list.shape[1] == 2, 'Bond list should have shape (n_bonds, 2)'
+    assert system.bond_list.shape[0] > 0, 'Should have at least one bond'
+    # Water molecules have 2 O-H bonds each (432 bonds for 144 water molecules)
+    assert system.bond_list.shape[0] == 432
+
+
+def test_get_configurations_includes_bond_list():
+    """Test that get_configurations includes bond_list in first frame only."""
+    mdap = gromacs_parser.GromacsMDAnalysisParser()
+    mdap._trajectory_steps_sampled = [0, 1]
+
+    # Mock interactions with bonds
+    mock_interactions = [
+        {'type': 'bond', 'atom_indices': [0, 1], 'atom_labels': ['O', 'H']},
+        {'type': 'bond', 'atom_indices': [0, 2], 'atom_labels': ['O', 'H']},
+    ]
+
+    # Mock data object
+    class MockDataObject:
+        def get_positions(self, idx):
+            return np.zeros((3, 3))
+
+        def get_velocities(self, idx):
+            return np.zeros((3, 3))
+
+        def get_lattice_vectors(self, idx):
+            return np.eye(3)
+
+        def get_atom_labels(self, idx):
+            return ['O', 'H', 'H']
+
+        def get_interactions(self):
+            return mock_interactions
+
+    mdap.data_object = MockDataObject()
+
+    configs = mdap.get_configurations()
+
+    assert len(configs) == 2
+    # First configuration should have bond_list
+    assert 'bond_list' in configs[0]
+    assert isinstance(configs[0]['bond_list'], np.ndarray)
+    assert configs[0]['bond_list'].shape == (2, 2)
+    # Second configuration should NOT have bond_list (topology is time-independent)
+    assert 'bond_list' not in configs[1]
+
+
 def test_get_thermodynamic_ensemble():
     """Test ensemble determination from thermostat and barostat settings."""
     lp = gromacs_parser.GromacsLogParser()