molecule can live with old-style mapping or with mapping divided into mapping and naming. It will allow us to migrate gradually to new mapping structure.

(cherry picked from commit 453170de7b67f3444052b802602f0fa82a30aac4)

Author: comcon1

src/fairmd/lipids/data/ToyData/Molecules/membrane/BOG/metadata.yaml

@@ -0,0 +1,27 @@
+NMRlipids:
+  id: BOG
+  name: beta-octyl D-glucopyranoside
+  charge: 0
+sameAs:
+  ChEBI: CHEBI:41128
+  pubchem.compound: 62852
+  metabolights: MTBLC41128
+  pdb.ligand: BOG
+  ChEMBL: CHEMBL446037
+bioschema_properties:
+  iupacName: (2R,3S,4S,5R,6R)-2-(hydroxymethyl)-6-octoxyoxane-3,4,5-triol
+  smiles: CCCCCCCCO[C@H]1[C@@H]([C@H]([C@@H]([C@H](O1)CO)O)O)O
+  inChI: InChI=1S/C14H28O6/c1-2-3-4-5-6-7-8-19-14-13(18)12(17)11(16)10(9-15)20-14/h10-18H,2-9H2,1H3/t10-,11-,12+,13-,14-/m1/s1
+  inChIKey: HEGSGKPQLMEBJL-RKQHYHRCSA-N
+  molecularFormula: C14H28O6
+  molecularWeight: 292.37
+  name: (2R,3S,4S,5R,6R)-2-(hydroxymethyl)-6-octoxyoxane-3,4,5-triol
+  image: https://www.ebi.ac.uk/chembl/api/data/image/CHEMBL446037?dimensions=200
+  alternateName:
+  - β-octylglucoside
+  - 1-octyl-β-D-glucopyranoside
+  - 1-O-n-octyl-β-D-glucopyranoside
+  - β-D-octyl glucoside
+  - 1-O-octyl-β-D-glucopyranoside
+  - octyl-β-D-glucoside
+  - Oct β-Glc

src/fairmd/lipids/data/ToyData/Molecules/membrane/BOG/naming.yaml

@@ -0,0 +1,116 @@
+M_G0C1_M:
+  FRAGMENT: glucose
+  SMILEIDX: 9
+M_G0C1H1_M:
+  FRAGMENT: glucose
+M_G0C1O1_M:
+  FRAGMENT: glucose
+  SMILEIDX: 8
+M_G0C2_M:
+  FRAGMENT: glucose
+  SMILEIDX: 10
+M_G0C2H1_M:
+  FRAGMENT: glucose
+M_G0C2O1_M:
+  FRAGMENT: glucose
+  SMILEIDX: 19
+M_G0C2O1H1_M:
+  FRAGMENT: glucose
+M_G0C3_M:
+  FRAGMENT: glucose
+  SMILEIDX: 11
+M_G0C3H1_M:
+  FRAGMENT: glucose
+M_G0C3O1_M:
+  FRAGMENT: glucose
+  SMILEIDX: 18
+M_G0C3O1H1_M:
+  FRAGMENT: glucose
+M_G0C4_M:
+  FRAGMENT: glucose
+  SMILEIDX: 12
+M_G0C4H1_M:
+  FRAGMENT: glucose
+M_G0C4O1_M:
+  FRAGMENT: glucose
+  SMILEIDX: 17
+M_G0C4O1H1_M:
+  FRAGMENT: glucose
+M_G0C5_M:
+  FRAGMENT: glucose
+  SMILEIDX: 13
+M_G0C5H1_M:
+  FRAGMENT: glucose
+M_G0C5O1_M:
+  FRAGMENT: glucose
+  SMILEIDX: 14
+M_G0C6_M:
+  FRAGMENT: glucose
+  SMILEIDX: 15
+M_G0C6H1_M:
+  FRAGMENT: glucose
+M_G0C6H2_M:
+  FRAGMENT: glucose
+M_G0C6O1_M:
+  FRAGMENT: glucose
+  SMILEIDX: 16
+M_G0C6O1H1_M:
+  FRAGMENT: glucose
+M_C1_M:
+  FRAGMENT: tail
+  SMILEIDX: 7
+M_C1H1_M:
+  FRAGMENT: tail
+M_C1H2_M:
+  FRAGMENT: tail
+M_C2_M:
+  FRAGMENT: tail
+  SMILEIDX: 6
+M_C2H1_M:
+  FRAGMENT: tail
+M_C2H2_M:
+  FRAGMENT: tail
+M_C3_M:
+  FRAGMENT: tail
+  SMILEIDX: 5
+M_C3H1_M:
+  FRAGMENT: tail
+M_C3H2_M:
+  FRAGMENT: tail
+M_C4_M:
+  FRAGMENT: tail
+  SMILEIDX: 4
+M_C4H1_M:
+  FRAGMENT: tail
+M_C4H2_M:
+  FRAGMENT: tail
+M_C5_M:
+  FRAGMENT: tail
+  SMILEIDX: 3
+M_C5H1_M:
+  FRAGMENT: tail
+M_C5H2_M:
+  FRAGMENT: tail
+M_C6_M:
+  FRAGMENT: tail
+  SMILEIDX: 2
+M_C6H1_M:
+  FRAGMENT: tail
+M_C6H2_M:
+  FRAGMENT: tail
+M_C7_M:
+  FRAGMENT: tail
+  SMILEIDX: 1
+M_C7H1_M:
+  FRAGMENT: tail
+M_C7H2_M:
+  FRAGMENT: tail
+M_C8_M:
+  FRAGMENT: tail
+  SMILEIDX: 0
+M_C8H1_M:
+  FRAGMENT: tail
+M_C8H2_M:
+  FRAGMENT: tail
+M_C8H3_M:
+  FRAGMENT: tail

src/fairmd/lipids/molecules.py

@@ -40,10 +40,13 @@ class MoleculeMappingError(MoleculeError):
     def __init__(self, message: str, mol=None) -> None:
         if mol is None:
             msg = message
-        elif not mol.is_mapping_registered:
+        elif not mol._can_load_mapping():
             msg = f"From {mol}: {message}"
         else:
-            disp_name = os.path.relpath(mol._mapping_fpath, FMDL_MOL_PATH)
+            if mol._mapping_fpath is None:
+                disp_name = "naming.yaml"
+            else:
+                disp_name = os.path.relpath(mol._mapping_fpath, FMDL_MOL_PATH)
             msg = f"From {mol}[{disp_name}]: {message}" if mol is not None else message
         super().__init__(msg, mol=mol)
 
@@ -114,18 +117,22 @@ def check_mapping(self, u: mda.Universe, name: str) -> bool:
     @property
     def mapping_dict(self) -> dict:
         """Return mapping dictionary (load on first call)"""
-        if self._mapping_fpath is None:
+        if not self._can_load_mapping():
             msg = "Mapping file is not registered!"
             raise MoleculeError(msg, mol=self)
-        if self._mapping_dict is None:
-            try:
-                with open(self._mapping_fpath) as yaml_file:
-                    self._mapping_dict = yaml.safe_load(yaml_file)  # yaml.load(yaml_file, Loader=yaml.FullLoader)
-            except OSError as e:
-                msg = "Error opening mapping-file!"
-                raise MoleculeError(msg, mol=self) from e
+        if self._mapping_dict is None:  # load on first request
+            self._load_mapping_dict()
         return self._mapping_dict
 
+    def _load_mapping_dict(self) -> None:
+        """Load mapping dictionary from the registered mapping file."""
+        try:
+            with open(self._mapping_fpath) as yaml_file:
+                self._mapping_dict = yaml.safe_load(yaml_file)  # yaml.load(yaml_file, Loader=yaml.FullLoader)
+        except OSError as e:
+            msg = "Error opening mapping-file!"
+            raise MoleculeError(msg, mol=self) from e
+
     def md2uan(self, mdatomname: str, mdresname: str | None = None) -> str:
         """
         Convert MD atom name to the Universal Atom Name.
@@ -134,6 +141,12 @@ def md2uan(self, mdatomname: str, mdresname: str | None = None) -> str:
         :return: Universal Atom Name (str)
         """
         for universal_name, mrecord in self.mapping_dict.items():
+            if "ATOMNAME" not in mrecord:
+                msg = (
+                    f"ATOMNAME field is missing for {universal_name} in mapping dictionary."  # !
+                    " MD mapping is not possible."
+                )
+                raise MoleculeMappingError(msg, mol=self)
             mapping_aname = mrecord["ATOMNAME"]
             # MDAnalysis uses fnmatch patterns for selection language
             # https://userguide.mdanalysis.org/stable/selections.html
@@ -153,6 +166,9 @@ def uan2selection(self, uname: str, resname: str) -> str:
         :raises KeyError: if the universal name is not found in the mapping.
         :return: selection string for MDAnalysis
         """
+        if "ATOMNAME" not in self.mapping_dict[uname]:
+            msg = f"ATOMNAME field is missing for {uname} in mapping dictionary. MD mapping is not possible."
+            raise MoleculeMappingError(msg, mol=self)
         anm = self.mapping_dict[uname]["ATOMNAME"]
         selstr = f"name {anm}"
         if "RESIDUE" in self.mapping_dict[uname]:
@@ -211,8 +227,7 @@ def name(self) -> str:
         """Molecule name"""
         return self._molname
 
-    @property
-    def is_mapping_registered(self) -> bool:
+    def _can_load_mapping(self) -> bool:
         """Is mapping registered for the molecule?"""
         return self._mapping_fpath is not None
 
@@ -253,6 +268,40 @@ def _populate_meta_data(self) -> None:
             msg = f"Metadata file not found for {self.name}."
             raise FileNotFoundError(msg)
 
+    def _can_load_mapping(self) -> bool:
+        super_is_reg = super()._can_load_mapping()
+        naming_path = os.path.join(self._get_path(), "naming.yaml")
+        return super_is_reg or os.path.isfile(naming_path)
+
+    def _load_mapping_dict(self) -> None:
+        if super()._can_load_mapping():
+            super()._load_mapping_dict()
+        else:
+            self._mapping_dict = {}
+        self._load_naming_dict()
+
+    def _load_naming_dict(self) -> None:
+        """
+        Load naming dictionary from `naming.yaml` file if it exists.
+
+        Updates existing mapping dictionary.
+        """
+        naming_path = os.path.join(self._get_path(), "naming.yaml")
+        if os.path.isfile(naming_path):
+            with open(naming_path) as yaml_file:
+                _naming_dict = yaml.load(yaml_file, Loader=yaml.FullLoader)
+            for unm, record in _naming_dict.items():
+                self._mapping_dict.setdefault(unm, {}).update(record)
+
+    @property
+    def fragments(self) -> list[str]:
+        """Return list of fragments for the lipid."""
+        frags = set()
+        for mrecord in self.mapping_dict.values():
+            frag = mrecord.get("FRAGMENT", "total")
+            frags.add(frag)
+        return sorted(frags)
+
     @property
     def metadata(self) -> dict:
         """

tests/test_molecules.py

@@ -7,7 +7,7 @@
 # run only on sim2 mocking data
 pytestmark = [pytest.mark.sim2, pytest.mark.min]
 
-LIPIDS_SET_LENGTH = 5
+LIPIDS_SET_LENGTH = 6
 POPE_MOLECULAR_WEIGHT = 717.5
 
 
@@ -33,15 +33,35 @@ def toy_mols_no_mapping() -> dict:
     from fairmd.lipids import FMDL_MOL_PATH
     from fairmd.lipids.molecules import lipids_set
 
-    mol1 = lipids_set.get("POPE")
-    mol2 = lipids_set.get("POPC")
+    mol1 = lipids_set.get("BOG")
 
-    return {"pope": mol1, "popc": mol2}
+    return {"bog": mol1}
 
 
-def tests_molecule_fragments(toy_mols_no_mapping):
-    pope = toy_mols_no_mapping["pope"]
+def tests_molecule_fragments(toy_mols_no_mapping, toy_mols_w_mapping):
+    bog = toy_mols_no_mapping["bog"]
+    check.is_instance(bog.fragments, list, "Should return a list of fragment names")
+    check.equal(bog.fragments, ["glucose", "tail"], "BOG fragments should be ['head', 'tail']")
+
+    pope = toy_mols_w_mapping["pope/charmm"]
     check.is_instance(pope.fragments, list, "Should return a list of fragment names")
+    check.equal(
+        sorted(pope.fragments),
+        sorted(["headgroup", "glycerol backbone", "sn-1", "sn-2"]),
+        "POPE fragments are improper",
+    )
+
+
+def tests_try_access_nomap(toy_mols_no_mapping):
+    from fairmd.lipids.molecules import Lipid, MoleculeMappingError
+
+    bog: Lipid = toy_mols_no_mapping["bog"]
+    bog.mapping_dict  # should not raise because mapping naming mapping exists
+    with pytest.raises(MoleculeMappingError, match="MD mapping is not possible"):
+        _ = bog.uan2selection("M_G0C1_M", "BOG")
+    # md2uan should also raise
+    with pytest.raises(MoleculeMappingError, match="MD mapping is not possible"):
+        _ = bog.md2uan("C1", "BOG")
 
 
 def test_mapping_dict():

tests/test_molecules_rdkit.py

@@ -30,8 +30,39 @@ def test_lipid_smarts():
     with check.raises(KeyError):
         popc.atoms_by("N", 1)  # 0 is max here
 
+    popc.register_mapping("mappingPOPCcharmm.yaml")
+
     # Test that we can retrieve atoms by SMARTS
     nitrg_atoms = popc.atoms_by("N", 0)  # can call
     check.is_instance(nitrg_atoms, list, "Should return a list of atom universal names")
     check.equal(len(nitrg_atoms), 1, "Should have at least one nitrogen atom")
     check.equal(nitrg_atoms[0], "M_G3N6_M", "Bad atom selected by N from POPC")
+
+
+def test_lipid_nomapping():
+    """Test atoms_by() method with SMARTS queries."""
+    from fairmd.lipids.molecules import Lipid, MoleculeError
+
+    popc = Lipid("BOG")
+    _ = popc.rdkit_object  # rdkit object can be created without mapping
+
+    _save = popc._metadata["bioschema_properties"]["smiles"]
+    del popc._metadata["bioschema_properties"]["smiles"]
+    with check.raises(MoleculeError):
+        _ = popc.rdkit_object  # rdkit object cannot be created without smiles
+    popc._metadata["bioschema_properties"]["smiles"] = _save
+
+    with check.raises(ValueError):
+        popc.atoms_by("~!x[", 0)
+
+    with check.raises(KeyError):
+        popc.atoms_by("N", -1)
+
+    with check.raises(KeyError):
+        popc.atoms_by("N", 1)  # 0 is max here
+
+    # Test that we can retrieve atoms by SMARTS
+    ox_atoms = popc.atoms_by("O", 0)  # can call
+    check.is_instance(ox_atoms, list, "Should return a list of atom universal names")
+    check.equal(len(ox_atoms), 6, "Should have 6 oxygen atoms")
+    check.equal(ox_atoms[0], "M_G0C1O1_M", "Bad atom selected by O from BOG")