From 4097891f3178cc0a261654095e9af10b281f293d Mon Sep 17 00:00:00 2001
From: Stefan Bienert <stefan.bienert@unibas.ch>
Date: Mon, 19 Jun 2023 09:43:46 +0200
Subject: [PATCH] Ignore invalid atoms in brnach links of the mmCIF foemat.

---
 modules/io/src/mol/mmcif_info.cc    |  6 ++++++
 modules/io/tests/test_mmcif_info.cc | 15 +++++++++++----
 2 files changed, 17 insertions(+), 4 deletions(-)

diff --git a/modules/io/src/mol/mmcif_info.cc b/modules/io/src/mol/mmcif_info.cc
index 576e05912..ef44d4813 100644
--- a/modules/io/src/mol/mmcif_info.cc
+++ b/modules/io/src/mol/mmcif_info.cc
@@ -207,6 +207,12 @@ void MMCifInfo::AddEntityBranchLink(String chain_name,
                                     mol::AtomHandle atom2,
                                     unsigned char bond_order)
 {
+  if (!atom1.IsValid() || !atom2.IsValid()) {
+    /* Would love to give details about the atoms... but atom names are not
+       available at this point. */
+    LOG_WARNING("Invalid branch link found in chain '"+chain_name+"'.");
+    return;
+  }
   // check if element already exists
   MMCifInfoEntityBranchLinkMap::iterator blm_it =
                                                entity_branches_.find(chain_name);
diff --git a/modules/io/tests/test_mmcif_info.cc b/modules/io/tests/test_mmcif_info.cc
index 5370d2ce1..faf5cc103 100644
--- a/modules/io/tests/test_mmcif_info.cc
+++ b/modules/io/tests/test_mmcif_info.cc
@@ -341,16 +341,23 @@ BOOST_AUTO_TEST_CASE(mmcif_info)
   mol::ResidueHandle res11 = editor.AppendResidue(ch1, "NAG");
   mol::ResidueHandle res12 = editor.AppendResidue(ch1, "NAG");
   // create AtomHandles for testing
-  mol::AtomHandle atom11 = editor.InsertAtom(res12, "C1",geom::Vec3());
-  mol::AtomHandle atom12 = editor.InsertAtom(res11, "O4",geom::Vec3());
+  mol::AtomHandle atom11 = editor.InsertAtom(res12, "C1", geom::Vec3());
+  mol::AtomHandle atom12 = editor.InsertAtom(res11, "O4", geom::Vec3());
   mol::ChainHandle ch2 = editor.InsertChain("B");
   mol::ResidueHandle res21 = editor.AppendResidue(ch2, "BMA");
   mol::ResidueHandle res22 = editor.AppendResidue(ch2, "MAN");
   // create AtomHandles for testing
-  mol::AtomHandle atom21 = editor.InsertAtom(res22, "C1",geom::Vec3());
-  mol::AtomHandle atom22 = editor.InsertAtom(res21, "O3",geom::Vec3());
+  mol::AtomHandle atom21 = editor.InsertAtom(res22, "C1", geom::Vec3());
+  mol::AtomHandle atom22 = editor.InsertAtom(res21, "O3", geom::Vec3());
+  // create invalid AtomHandle pairs for testing
+  mol::AtomHandle atom_invalid;
   info.AddEntityBranchLink(ch1.GetName(), atom11, atom12, 1);
   info.AddEntityBranchLink(ch2.GetName(), atom21, atom22, 1);
+  /* Sometimes branched PDB entries link two atoms which are available in the
+     compound's definition but not resolved (missing) in the coordinates, e.g.
+     RCSB entry 7zim. Check that in case of invalid atom, no link is created. */
+  info.AddEntityBranchLink(ch2.GetName(), atom11, atom_invalid, 1);
+  info.AddEntityBranchLink(ch2.GetName(), atom_invalid, atom12, 1);
   std::vector<MMCifInfoEntityBranchLink> blinks = info.GetEntityBranchLinks();
 
   BOOST_CHECK(blinks.size() == 2);
-- 
GitLab