diff --git a/CHANGELOG.txt b/CHANGELOG.txt
index 1fa74533fa8039047f11f5d94c60b19de4bfd225..01e85e9cf5e605b142e1f032f30eb99e1776bf61 100644
--- a/CHANGELOG.txt
+++ b/CHANGELOG.txt
@@ -15,6 +15,7 @@ Changes in Release 2.5.0
    flag to disable stereochemistry checks for lDDT, TMscore (including
    associated chain mapping) computed by USalign, make residues/atoms uniquely
    identifiable (also considering residue number insertion codes).
+ * Read atom charges if present from PDB/mmCIF/SDF files
  * Several bug fixes and improvements.
 
 Changes in Release 2.4.0
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 83aa5d857cd44b803c1de242cfeb5cdbe018624b..aa7daa5413e1039eb3ec3181d824caa050686f3a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -14,7 +14,7 @@ cmake_policy(SET CMP0060 NEW)
 project(OpenStructure CXX C)
 set (CMAKE_EXPORT_COMPILE_COMMANDS 1)
 set (OST_VERSION_MAJOR 2)
-set (OST_VERSION_MINOR 4)
+set (OST_VERSION_MINOR 5)
 set (OST_VERSION_PATCH 0)
 set (OST_VERSION_STRING ${OST_VERSION_MAJOR}.${OST_VERSION_MINOR}.${OST_VERSION_PATCH} )
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${PROJECT_SOURCE_DIR}/cmake_support)
diff --git a/actions/ost-compare-ligand-structures b/actions/ost-compare-ligand-structures
index 3c33dceb50fa6378532e6066067283de524530ac..176cd5e1f6028c096d9b1bdc3dbd4b6f3fdce673 100644
--- a/actions/ost-compare-ligand-structures
+++ b/actions/ost-compare-ligand-structures
@@ -57,7 +57,6 @@ import os
 import traceback
 
 import ost
-from ost import conop
 from ost import io
 from ost.mol.alg import ligand_scoring
 
@@ -366,11 +365,24 @@ def _Process(model, model_ligands, reference, reference_ligands, args):
 
     if model_ligands is not None:
         # Replace model ligand by path
-        assert len(model_ligands) == len(scorer.model_ligands)
-        # Map ligand => path
-        model_ligands_map = {k: v for k, v in zip(scorer.model_ligands,
-                                                  args.model_ligands)}
-        out["model_ligands"] = args.model_ligands
+        if len(model_ligands) == len(scorer.model_ligands):
+            # Map ligand => path
+            model_ligands_map = {k: v for k, v in zip(scorer.model_ligands,
+                                                      args.model_ligands)}
+            out["model_ligands"] = args.model_ligands
+        elif len(model_ligands) < len(scorer.model_ligands):
+            # Multi-ligand SDF files were given
+            # Map ligand => path:idx
+            out["model_ligands"] = []
+            for ligand, filename in zip(model_ligands, args.model_ligands):
+                assert isinstance(ligand, ost.mol.EntityHandle)
+                for i, residue in enumerate(ligand.residues):
+                    out["model_ligands"].append(f"{filename}:{i}")
+        else:
+            # This should never happen and would be a bug
+            raise RuntimeError("Fewer ligands in the model scorer "
+                               "(%d) than given (%d)" % (
+                len(scorer.model_ligands), len(model_ligands)))
     else:
         model_ligands_map = {l: _QualifiedResidueNotation(l)
                              for l in scorer.model_ligands}
@@ -378,11 +390,25 @@ def _Process(model, model_ligands, reference, reference_ligands, args):
 
     if reference_ligands is not None:
         # Replace reference ligand by path
-        assert len(reference_ligands) == len(scorer.target_ligands)
-        # Map ligand => path
-        reference_ligands_map = {k: v for k, v in zip(scorer.target_ligands,
-                                                  args.reference_ligands)}
-        out["reference_ligands"] = args.reference_ligands
+        if len(reference_ligands) == len(scorer.target_ligands):
+            # Map ligand => path
+            reference_ligands_map = {k: v for k, v in zip(scorer.target_ligands,
+                                                      args.reference_ligands)}
+            out["reference_ligands"] = args.reference_ligands
+        elif len(reference_ligands) < len(scorer.target_ligands):
+            # Multi-ligand SDF files were given
+            # Map ligand => path:idx
+            out["reference_ligands"] = []
+            for ligand, filename in zip(reference_ligands, args.reference_ligands):
+                assert isinstance(ligand, ost.mol.EntityHandle)
+                for i, residue in enumerate(ligand.residues):
+                    out["reference_ligands"].append(f"{filename}:{i}")
+        else:
+            # This should never happen and would be a bug
+            raise RuntimeError("Fewer ligands in the reference scorer "
+                               "(%d) than given (%d)" % (
+                len(scorer.target_ligands), len(reference_ligands)))
+
     else:
         reference_ligands_map = {l: _QualifiedResidueNotation(l)
                              for l in scorer.target_ligands}
diff --git a/docker/Dockerfile b/docker/Dockerfile
index 1975e401be9e21d7bcbee708b8ab5b5d068ebc5b..e019bcef45d8742129ea99819072b09b9baea6eb 100644
--- a/docker/Dockerfile
+++ b/docker/Dockerfile
@@ -2,7 +2,7 @@ FROM ubuntu:22.04
 
 # ARGUMENTS
 ###########
-ARG OPENSTRUCTURE_VERSION="2.4.0"
+ARG OPENSTRUCTURE_VERSION="2.5.0"
 ARG SRC_FOLDER="/usr/local/src"
 ARG CPUS_FOR_MAKE=2
 ARG OPENMM_VERSION="7.7.0"
diff --git a/modules/doc/actions.rst b/modules/doc/actions.rst
index a71c8db2ec2cbc978da9de73001857e8718da23d..8fd2dbab5d4be5f12cc6ae9c0443818fa20e8e44 100644
--- a/modules/doc/actions.rst
+++ b/modules/doc/actions.rst
@@ -23,7 +23,6 @@ You can compare two structures from the command line with the
 interface to :class:`ost.mol.alg.scoring.Scorer`
 
 .. warning::
-
   ``compare-structures`` underwent a complete rewrite in OpenStructure
   release 2.4.0. The old version is still available as
   ``compare-structures-legacy`` with documentation available
@@ -412,4 +411,4 @@ Details on the usage (output of ``ost compare-ligand-structures --help``):
 
 Additional information about the scores and output values is available in
 :meth:`rmsd_details <ost.mol.alg.ligand_scoring.LigandScorer.rmsd_details>` and
-:meth:`lddt_pli_details <ost.mol.alg.ligand_scoring.LigandScorer.lddt_pli_details>`.
\ No newline at end of file
+:meth:`lddt_pli_details <ost.mol.alg.ligand_scoring.LigandScorer.lddt_pli_details>`.
diff --git a/modules/io/src/mol/mmcif_info.cc b/modules/io/src/mol/mmcif_info.cc
index 576e05912bb0073c9d93215d12ced3747353a2c8..ef44d4813fef36f2c709888b2acaef4d007342ee 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/src/mol/mmcif_reader.cc b/modules/io/src/mol/mmcif_reader.cc
index 7ddde3b2cf53f36807b659cb943f13207ed0a203..e550d96f83dfc160819f860055561534606b48fd 100644
--- a/modules/io/src/mol/mmcif_reader.cc
+++ b/modules/io/src/mol/mmcif_reader.cc
@@ -141,6 +141,7 @@ bool MMCifReader::OnBeginLoop(const StarLoopDesc& header)
     indices_[AUTH_SEQ_ID]        = header.GetIndex("auth_seq_id");
     indices_[PDBX_PDB_INS_CODE]  = header.GetIndex("pdbx_PDB_ins_code");
     indices_[PDBX_PDB_MODEL_NUM] = header.GetIndex("pdbx_PDB_model_num");
+    indices_[FORMAL_CHARGE]     = header.GetIndex("pdbx_formal_charge");
 
     // post processing
     if (category_counts_[category_] > 0) {
@@ -482,6 +483,7 @@ void MMCifReader::ParseAndAddAtom(const std::vector<StringRef>& columns)
     return;                            
   }
   Real occ = 1.00f, temp = 0;
+  int charge = 0;
   geom::Vec3 apos;
   
   for (int i = CARTN_X; i <= CARTN_Z; ++i) {
@@ -505,6 +507,13 @@ void MMCifReader::ParseAndAddAtom(const std::vector<StringRef>& columns)
                               "atom_site.B_iso_or_equiv");
     }
   }
+  if (indices_[FORMAL_CHARGE] != -1) { // unit test
+    String charge_s = columns[indices_[FORMAL_CHARGE]].str();
+    if (charge_s != "?" && charge_s != ".") {
+      charge = this->TryGetInt(columns[indices_[FORMAL_CHARGE]],
+                               "atom_site.pdbx_formal_charge");
+    }
+  }
 
   // determine element
   String s_ele(columns[indices_[TYPE_SYMBOL]].str());
@@ -665,6 +674,8 @@ void MMCifReader::ParseAndAddAtom(const std::vector<StringRef>& columns)
 
   ah.SetOccupancy(occ);
 
+  ah.SetCharge(charge);
+
   // record type
   ah.SetHetAtom(indices_[GROUP_PDB] == -1 ? false :  
                 columns[indices_[GROUP_PDB]][0]=='H');
diff --git a/modules/io/src/mol/mmcif_reader.hh b/modules/io/src/mol/mmcif_reader.hh
index 1f21c3928cac354f5f69197e763443d37a2dc9f6..f2f2d2dbd78069a778c6cd3ccab9d49abb7b0cb1 100644
--- a/modules/io/src/mol/mmcif_reader.hh
+++ b/modules/io/src/mol/mmcif_reader.hh
@@ -359,7 +359,7 @@ protected:
 private:
   /// \enum magic numbers of this class
   typedef enum {
-    MAX_ITEMS_IN_ROW=18 ///< count for possible items in a loop row
+    MAX_ITEMS_IN_ROW=19 ///< count for possible items in a loop row
   } MMCifMagicNos;
 
   /// \enum items of the atom_site category
@@ -381,7 +381,8 @@ private:
     B_ISO_OR_EQUIV,
     PDBX_PDB_INS_CODE,
     GROUP_PDB,         ///< record name
-    PDBX_PDB_MODEL_NUM ///< model no. (especially NMR structures)
+    PDBX_PDB_MODEL_NUM,///< model no. (especially NMR structures)
+    FORMAL_CHARGE
   } AtomSiteItems;
 
   /// \enum items of the entity category
diff --git a/modules/io/src/mol/pdb_reader.cc b/modules/io/src/mol/pdb_reader.cc
index 3b69915de24795389f6566e5ce43e0cc487d9c8d..930ffbcf453f0563768d8aae0351f65fe6ac07bf 100644
--- a/modules/io/src/mol/pdb_reader.cc
+++ b/modules/io/src/mol/pdb_reader.cc
@@ -697,8 +697,11 @@ void PDBReader::ParseAndAddAtom(const StringRef& line, int line_num,
                             atom_name, alt_loc, record_type, serial)) {
     return;                            
   }
-  std::pair<bool, Real> charge, radius;
-  std::pair<bool, Real> occ, temp;
+
+  std::pair<bool, Real> charge = std::make_pair(false, 0.0);
+  std::pair<bool, Real> radius = std::make_pair(false, 0.0);
+  std::pair<bool, Real> occ = std::make_pair(false, 0.0);
+  std::pair<bool, Real> temp = std::make_pair(false, 0.0);
   geom::Vec3 apos;
 
   for (int i=0;i<3; ++i) {
@@ -728,6 +731,20 @@ void PDBReader::ParseAndAddAtom(const StringRef& line, int line_num,
     if (line.length()>=66) {
       temp=line.substr(60, 6).ltrim().to_float();      
     }
+    if (line.length()>=80) {
+      // example charge formatting: 1+, 2- etc.
+      charge = line.substr(78,1).trim().to_float();
+      if(charge.first) {
+        if(line[79] != '-' && line[79] != '+') {
+          std::stringstream ss;
+          ss << "error on line " << line_num << ": "
+             << "expect charge in format 1+, 2-, etc. got: "
+             << line.substr(78, 2);
+          throw IOException(ss.str());      
+        }
+        if(line[79] == '-') charge.second *= (-1);
+      }
+    }
   }
   LOG_TRACE( "line: [" << line << "]" );
   String s_ele;
@@ -838,6 +855,7 @@ void PDBReader::ParseAndAddAtom(const StringRef& line, int line_num,
   }
   Real b=temp.first ? temp.second : 0.0;
   Real o=occ.first ? occ.second : 1.0;
+
   if (!profile_.quack_mode && alt_loc!=' ') {
     // Check if there is already a atom with the same name.
     mol::AtomHandle me=curr_residue_.FindAtom(aname);
@@ -869,16 +887,16 @@ void PDBReader::ParseAndAddAtom(const StringRef& line, int line_num,
     ah=editor.InsertAtom(curr_residue_, aname, apos, s_ele);
     ++atom_count_;
   }
-  if(is_pqr_) {
-    if (radius.first) {
-      ah.SetRadius(radius.second);
-    }
+
+  if (radius.first) {
+    ah.SetRadius(radius.second);
   }
-  ah.SetBFactor(b);
-  ah.SetOccupancy(o);
   if (charge.first) {
     ah.SetCharge(charge.second);
   }
+  ah.SetBFactor(b);
+  ah.SetOccupancy(o);
+
   ah.SetHetAtom(record_type[0]=='H');
   if(profile_.read_conect && serial != -1) {
     this->amap_[serial] = ah;
diff --git a/modules/io/src/mol/sdf_reader.cc b/modules/io/src/mol/sdf_reader.cc
index 2db564c66de7d6c5990a9b78897d3420b5abfc12..3c28e279cb71a3f89575d5fe3bca70ce621c6c6f 100644
--- a/modules/io/src/mol/sdf_reader.cc
+++ b/modules/io/src/mol/sdf_reader.cc
@@ -159,6 +159,11 @@ void SDFReader::ParseAndAddHeader(const String& line, int line_num,
       break;
     case 4:  // counts line
     {
+      String version_str=line.substr(34, 5);
+      if (version_str != "V2000") {
+        String msg="Unsupported SDF version: %s.";
+        throw IOException(str(format(msg) % version_str));
+      }
       String s_anum=line.substr(0,3);
       try {
         atom_count_=boost::lexical_cast<int>(boost::trim_copy(s_anum));
@@ -188,9 +193,15 @@ void SDFReader::ParseAndAddAtom(const String& line, int line_num,
   LOG_TRACE( "line: [" << line << "]" );
 
   if(line.length()<48 || line.length()>69) {
-    String msg="Bad atom line %d: Not correct number of characters on the"
-               " line: %i (should be between 48 and 69)";
-    throw IOException(str(format(msg) % line_num % line.length()));
+    // Handle the case where we have trailing space characters
+    if (line.length()>69 && boost::trim_copy(line.substr(69)) == "") {
+      LOG_DEBUG( "Ignoring trailing space" );
+    }
+    else {
+      String msg="Bad atom line %d: Not correct number of characters on the"
+                 " line: %i (should be between 48 and 69)";
+      throw IOException(str(format(msg) % line_num % line.length()));
+    }
   }
   int anum = line_num-4;  // start at 1 on fifth line since first four lines are header
   String s_posx=line.substr(0,10);
@@ -253,9 +264,15 @@ void SDFReader::ParseAndAddBond(const String& line, int line_num,
   LOG_TRACE( "line: [" << line << "]" );
 
   if(line.length()<9 || line.length()>21) {
-    String msg="Bad bond line %d: Not correct number of characters on the"
-               " line: %i (should be between 9 and 21)";
-    throw IOException(str(format(msg) % line_num % line.length()));
+    // Handle the case where we have trailing space characters
+    if (line.length()>21 && boost::trim_copy(line.substr(21)) == "") {
+      LOG_DEBUG( "Ignoring trailing space" );
+    }
+    else {
+      String msg="Bad bond line %d: Not correct number of characters on the"
+                 " line: %i (should be between 9 and 21)";
+      throw IOException(str(format(msg) % line_num % line.length()));
+    }
   }
 
   String s_first_name=line.substr(0,3);
diff --git a/modules/io/src/mol/sdf_writer.cc b/modules/io/src/mol/sdf_writer.cc
index 8deaf9a7819f0fff6c069b721af7b55e446a0245..e39444d0cdd394b568f82568622c6b4fbc006545 100644
--- a/modules/io/src/mol/sdf_writer.cc
+++ b/modules/io/src/mol/sdf_writer.cc
@@ -50,7 +50,12 @@ namespace {
               << format("%10.4f") % atom.GetPos()[1]
               << format("%10.4f ") % atom.GetPos()[2]
               << format("%-3s") % SDFAtomWriter::FormatEle(atom.GetElement())
-              << " 0  0  0  0  0  0"
+              << " 0" // Mass difference
+              << format("%-3s") % SDFAtomWriter::FormatCharge(atom.GetCharge()) // Charge
+              << "  0" // Atom stereo parity
+              << "  0" // Hydrogen count + 1
+              << "  0" // Stereo care box
+              << "  0" // Valence
               << std::endl;
         return true;
       }
@@ -66,6 +71,24 @@ namespace {
         }
         return return_ele;
       }
+
+      static String FormatCharge(const Real& chg) {
+        // Format charge according to https://doi.org/10.1021/ci00007a012
+        // 0 = uncharged or value other than these, 1 = +3, 2 = +2, 3 = +1,
+        // 4 doublet (A), 5 = -1, 6 = -2, 7 = -3
+        // Doublet means radical. This function would never return 4.
+        if (chg == 0) {
+          return "  0";
+        }
+        else if (abs(chg) > 3) {
+          String msg = "SDF format only supports charges from -3 to +3, not %g";
+          throw IOException(str(format(msg) % chg));
+        }
+        else {
+          Real chg_sdf = 4 - chg;
+          return str(format("%3.0f") % chg_sdf);
+        }
+      }
     private:
       std::ostream&      ostr_;
       std::map<long, int>& atom_indices_;
diff --git a/modules/io/tests/test_io_sdf.py b/modules/io/tests/test_io_sdf.py
index 718ac0691192f24ddbc04c71d2ef71d15fa631cf..735158fc3c38f8ed2345415ae913143f8b7d1bf1 100644
--- a/modules/io/tests/test_io_sdf.py
+++ b/modules/io/tests/test_io_sdf.py
@@ -16,6 +16,28 @@ class TestSDF(unittest.TestCase):
     ent = io.LoadSDF('testfiles/sdf/6d5w_rank1_crlf.sdf.gz')
     self.assertEqual(len(ent.atoms), 21)
     self.assertEqual(len(ent.bonds), 24)
+
+  def test_Charge(self):
+    ent = io.LoadSDF('testfiles/sdf/simple.sdf')
+    self.assertEqual(ent.FindAtom("00001_Simple Ligand", 1, "6").charge,  0)
+
+    # Write and read charges properly
+    for chg in range(-3, 4):
+      ent.FindAtom("00001_Simple Ligand", 1, "6").charge = chg
+      sdf_str = io.EntityToSDFStr(ent)
+      ent = io.SDFStrToEntity(sdf_str)
+      self.assertEqual(ent.FindAtom("00001_Simple Ligand", 1, "6").charge,  chg)
+
+    # Only -3 to +3 is supported
+    # If M CHG is implemented the following tests can be removed
+    with self.assertRaises(Exception):
+      ent.FindAtom("00001_Simple Ligand", 1, "6").charge = 4
+      io.EntityToSDFStr(ent)
+
+    with self.assertRaises(Exception):
+      ent.FindAtom("00001_Simple Ligand", 1, "6").charge = -4
+      io.EntityToSDFStr(ent)
+
     
 if __name__== '__main__':
   from ost import testutils
diff --git a/modules/io/tests/test_mmcif_info.cc b/modules/io/tests/test_mmcif_info.cc
index 5370d2ce1eeb733bbba23b2ff52fadaf2cd22274..faf5cc10379e88b73db25e58be05bbf64a34ce6b 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);
diff --git a/modules/io/tests/test_mmcif_reader.cc b/modules/io/tests/test_mmcif_reader.cc
index e03eb0073be2133e959420e390e3c6f07ec07fde..6d81b46033ab1143be596f7464020706a4645abd 100644
--- a/modules/io/tests/test_mmcif_reader.cc
+++ b/modules/io/tests/test_mmcif_reader.cc
@@ -1616,5 +1616,22 @@ BOOST_AUTO_TEST_CASE(mmcif_atom_site_B_iso_or_equiv_tests)
   BOOST_TEST_MESSAGE("  done.");
 }
 
+BOOST_AUTO_TEST_CASE(mmcif_formal_charge)
+{
+  mol::EntityHandle eh = mol::CreateEntity();
+  std::ifstream s("testfiles/mmcif/4C79_charged.cif");
+  IOProfile profile;
+  MMCifReader mmcif_p(s, eh, profile);
+  mmcif_p.Parse();
+
+  BOOST_CHECK_EQUAL(eh.FindAtom("A", 49, "OE2").GetCharge(), -1);
+  BOOST_CHECK_EQUAL(eh.FindAtom("A", 49, "OE1").GetCharge(), 0); // '?'
+  BOOST_CHECK_EQUAL(eh.FindAtom("A", 49, "CA").GetCharge(), 0);  // Explicit 0
+  BOOST_CHECK_EQUAL(eh.FindAtom("A", 49, "CB").GetCharge(), 0);  // '.'
+  BOOST_CHECK_EQUAL(eh.FindAtom("C", 1, "ZN").GetCharge(), 2);
+  BOOST_CHECK_EQUAL(eh.FindAtom("D", 1, "NA").GetCharge(), 1);
+
+}
+
 
 BOOST_AUTO_TEST_SUITE_END();
diff --git a/modules/io/tests/testfiles/mmcif/4C79_charged.cif b/modules/io/tests/testfiles/mmcif/4C79_charged.cif
new file mode 100644
index 0000000000000000000000000000000000000000..c5089608ff706068345942fb2240a2f7394c5a52
--- /dev/null
+++ b/modules/io/tests/testfiles/mmcif/4C79_charged.cif
@@ -0,0 +1,289 @@
+data_4C79
+# taken and modified from 4C79.cif
+_entry.id   4C79
+# 
+loop_
+_entity.id 
+_entity.type 
+_entity.src_method 
+_entity.pdbx_description 
+_entity.formula_weight 
+_entity.pdbx_number_of_molecules 
+_entity.pdbx_ec 
+_entity.pdbx_mutation 
+_entity.pdbx_fragment 
+_entity.details 
+1 polymer     man SMOOTHENED   22144.174 2  ? ? 'CYSTEINE-RICH DOMAIN (CRD), RESIDUES 28-210' ? 
+2 non-polymer syn 'ZINC ION'   65.409    1  ? ? ?                                             ? 
+3 non-polymer syn 'SODIUM ION' 22.990    2  ? ? ?                                             ? 
+4 water       nat water        18.015    48 ? ? ?                                             ? 
+# 
+_entity_poly.entity_id                      1 
+_entity_poly.type                           'polypeptide(L)' 
+_entity_poly.nstd_linkage                   no 
+_entity_poly.nstd_monomer                   no 
+_entity_poly.pdbx_seq_one_letter_code       
+;MAVILHPNETIFNDFCKKSTTCEVLKYNTCLGSPLPYTHTSLILAEDSETQEEAFEKLAMWSGLRNAPRCWAVIQPLLCA
+VYMPKCENGKVELPSQHLCQATRNPCSIVERERGWPNFLKCENKEQFPKGCQNEVQKLKFNTSGQCEAPLVKTDIQASWY
+KDVEGCGIQCDNPLFTEDEHSDMHKLEHHHHHH
+;
+_entity_poly.pdbx_seq_one_letter_code_can   
+;MAVILHPNETIFNDFCKKSTTCEVLKYNTCLGSPLPYTHTSLILAEDSETQEEAFEKLAMWSGLRNAPRCWAVIQPLLCA
+VYMPKCENGKVELPSQHLCQATRNPCSIVERERGWPNFLKCENKEQFPKGCQNEVQKLKFNTSGQCEAPLVKTDIQASWY
+KDVEGCGIQCDNPLFTEDEHSDMHKLEHHHHHH
+;
+_entity_poly.pdbx_strand_id                 A
+_entity_poly.pdbx_target_identifier         ? 
+# 
+loop_
+_entity_poly_seq.entity_id 
+_entity_poly_seq.num 
+_entity_poly_seq.mon_id 
+_entity_poly_seq.hetero 
+1 1   MET n 
+1 2   ALA n 
+1 3   VAL n 
+1 4   ILE n 
+1 5   LEU n 
+1 6   HIS n 
+1 7   PRO n 
+1 8   ASN n 
+1 9   GLU n 
+1 10  THR n 
+1 11  ILE n 
+1 12  PHE n 
+1 13  ASN n 
+1 14  ASP n 
+1 15  PHE n 
+1 16  CYS n 
+1 17  LYS n 
+1 18  LYS n 
+1 19  SER n 
+1 20  THR n 
+1 21  THR n 
+1 22  CYS n 
+1 23  GLU n 
+1 24  VAL n 
+1 25  LEU n 
+1 26  LYS n 
+1 27  TYR n 
+1 28  ASN n 
+1 29  THR n 
+1 30  CYS n 
+1 31  LEU n 
+1 32  GLY n 
+1 33  SER n 
+1 34  PRO n 
+1 35  LEU n 
+1 36  PRO n 
+1 37  TYR n 
+1 38  THR n 
+1 39  HIS n 
+1 40  THR n 
+1 41  SER n 
+1 42  LEU n 
+1 43  ILE n 
+1 44  LEU n 
+1 45  ALA n 
+1 46  GLU n 
+1 47  ASP n 
+1 48  SER n 
+1 49  GLU n 
+1 50  THR n 
+1 51  GLN n 
+1 52  GLU n 
+1 53  GLU n 
+1 54  ALA n 
+1 55  PHE n 
+1 56  GLU n 
+1 57  LYS n 
+1 58  LEU n 
+1 59  ALA n 
+1 60  MET n 
+1 61  TRP n 
+1 62  SER n 
+1 63  GLY n 
+1 64  LEU n 
+1 65  ARG n 
+1 66  ASN n 
+1 67  ALA n 
+1 68  PRO n 
+1 69  ARG n 
+1 70  CYS n 
+1 71  TRP n 
+1 72  ALA n 
+1 73  VAL n 
+1 74  ILE n 
+1 75  GLN n 
+1 76  PRO n 
+1 77  LEU n 
+1 78  LEU n 
+1 79  CYS n 
+1 80  ALA n 
+1 81  VAL n 
+1 82  TYR n 
+1 83  MET n 
+1 84  PRO n 
+1 85  LYS n 
+1 86  CYS n 
+1 87  GLU n 
+1 88  ASN n 
+1 89  GLY n 
+1 90  LYS n 
+1 91  VAL n 
+1 92  GLU n 
+1 93  LEU n 
+1 94  PRO n 
+1 95  SER n 
+1 96  GLN n 
+1 97  HIS n 
+1 98  LEU n 
+1 99  CYS n 
+1 100 GLN n 
+1 101 ALA n 
+1 102 THR n 
+1 103 ARG n 
+1 104 ASN n 
+1 105 PRO n 
+1 106 CYS n 
+1 107 SER n 
+1 108 ILE n 
+1 109 VAL n 
+1 110 GLU n 
+1 111 ARG n 
+1 112 GLU n 
+1 113 ARG n 
+1 114 GLY n 
+1 115 TRP n 
+1 116 PRO n 
+1 117 ASN n 
+1 118 PHE n 
+1 119 LEU n 
+1 120 LYS n 
+1 121 CYS n 
+1 122 GLU n 
+1 123 ASN n 
+1 124 LYS n 
+1 125 GLU n 
+1 126 GLN n 
+1 127 PHE n 
+1 128 PRO n 
+1 129 LYS n 
+1 130 GLY n 
+1 131 CYS n 
+1 132 GLN n 
+1 133 ASN n 
+1 134 GLU n 
+1 135 VAL n 
+1 136 GLN n 
+1 137 LYS n 
+1 138 LEU n 
+1 139 LYS n 
+1 140 PHE n 
+1 141 ASN n 
+1 142 THR n 
+1 143 SER n 
+1 144 GLY n 
+1 145 GLN n 
+1 146 CYS n 
+1 147 GLU n 
+1 148 ALA n 
+1 149 PRO n 
+1 150 LEU n 
+1 151 VAL n 
+1 152 LYS n 
+1 153 THR n 
+1 154 ASP n 
+1 155 ILE n 
+1 156 GLN n 
+1 157 ALA n 
+1 158 SER n 
+1 159 TRP n 
+1 160 TYR n 
+1 161 LYS n 
+1 162 ASP n 
+1 163 VAL n 
+1 164 GLU n 
+1 165 GLY n 
+1 166 CYS n 
+1 167 GLY n 
+1 168 ILE n 
+1 169 GLN n 
+1 170 CYS n 
+1 171 ASP n 
+1 172 ASN n 
+1 173 PRO n 
+1 174 LEU n 
+1 175 PHE n 
+1 176 THR n 
+1 177 GLU n 
+1 178 ASP n 
+1 179 GLU n 
+1 180 HIS n 
+1 181 SER n 
+1 182 ASP n 
+1 183 MET n 
+1 184 HIS n 
+1 185 LYS n 
+1 186 LEU n 
+1 187 GLU n 
+1 188 HIS n 
+1 189 HIS n 
+1 190 HIS n 
+1 191 HIS n 
+1 192 HIS n 
+1 193 HIS n
+# 
+_struct.entry_id                  4C79 
+_struct.title                     'Crystal structure of the Smoothened CRD, native' 
+_struct.pdbx_descriptor           SMOOTHENED 
+_struct.pdbx_model_details        ? 
+_struct.pdbx_CASP_flag            ? 
+_struct.pdbx_model_type_details   ? 
+# 
+_struct_keywords.entry_id        4C79 
+_struct_keywords.pdbx_keywords   'SIGNALING PROTEIN' 
+_struct_keywords.text            'SIGNALING PROTEIN'
+# 
+loop_
+_atom_site.group_PDB 
+_atom_site.id 
+_atom_site.type_symbol 
+_atom_site.label_atom_id 
+_atom_site.label_alt_id 
+_atom_site.label_comp_id 
+_atom_site.label_asym_id 
+_atom_site.label_entity_id 
+_atom_site.label_seq_id 
+_atom_site.pdbx_PDB_ins_code 
+_atom_site.Cartn_x 
+_atom_site.Cartn_y 
+_atom_site.Cartn_z 
+_atom_site.occupancy 
+_atom_site.B_iso_or_equiv 
+_atom_site.pdbx_formal_charge 
+_atom_site.auth_seq_id 
+_atom_site.auth_comp_id 
+_atom_site.auth_asym_id 
+_atom_site.auth_atom_id 
+_atom_site.pdbx_PDB_model_num
+ATOM   262  N  N   . GLU A 1 49  ? -25.812 5.207   -4.954  1.00 35.58  ? 75   GLU A N   1
+ATOM   263  C  CA  . GLU A 1 49  ? -26.815 4.149   -4.979  1.00 36.58  0 75   GLU A CA  1
+ATOM   264  C  C   . GLU A 1 49  ? -27.600 4.155   -6.288  1.00 35.91  ? 75   GLU A C   1
+ATOM   265  O  O   . GLU A 1 49  ? -28.206 3.150   -6.663  1.00 35.94  ? 75   GLU A O   1
+ATOM   266  C  CB  . GLU A 1 49  ? -27.774 4.284   -3.794  1.00 38.91  . 75   GLU A CB  1
+ATOM   267  C  CG  . GLU A 1 49  ? -27.114 4.163   -2.426  1.00 40.64  ? 75   GLU A CG  1
+ATOM   268  C  CD  . GLU A 1 49  ? -26.931 2.725   -1.965  1.00 42.04  ? 75   GLU A CD  1
+ATOM   269  O  OE1 . GLU A 1 49  ? -26.905 1.807   -2.812  1.00 42.07  ? 75   GLU A OE1 1
+ATOM   270  O  OE2 . GLU A 1 49  ? -26.811 2.516   -0.740  1.00 42.97 -1 75   GLU A OE2 1
+HETATM 1881 ZN ZN  . ZN  C 2 .   ? -29.714 -4.622  -22.478 1.00 31.76  2 1159 ZN  A ZN  1 
+HETATM 1882 NA NA  . NA  D 3 .   ? -23.050 1.721   -4.584  1.00 26.51  1 1160 NA  A NA  1 
+#
+loop_
+_pdbx_entity_nonpoly.entity_id 
+_pdbx_entity_nonpoly.name 
+_pdbx_entity_nonpoly.comp_id 
+2 'ZINC ION'   ZN  
+3 'SODIUM ION' NA  
+4 water        HOH 
+# 
diff --git a/modules/mol/alg/doc/lddt_deprecated.rst b/modules/mol/alg/doc/lddt_deprecated.rst
index 4293046fd7e9359808d394a815366b1696277453..173bcb5ed5b3d8adbdc5abec02df9f2bc580b765 100644
--- a/modules/mol/alg/doc/lddt_deprecated.rst
+++ b/modules/mol/alg/doc/lddt_deprecated.rst
@@ -3,6 +3,11 @@
 lDDT (deprecated)
 ================================================================================
 
+.. warning::
+
+  These functions in `ost.mol.alg` are deprecated. Consider using the newer
+  implementation in :class:`ost.mol.alg.lDDTScorer` instead.
+
 .. function:: LocalDistDiffTest(model, distance_list, tolerance_list, \
                                 sequence_separation=0, \
                                 local_lddt_property_string="")
diff --git a/modules/mol/alg/doc/molalg.rst b/modules/mol/alg/doc/molalg.rst
index bf451df52b20b045f3a0a1e5e82d793cd537847b..3a43ea5de90b8b5bf14393c15754b09b6e2d37d4 100644
--- a/modules/mol/alg/doc/molalg.rst
+++ b/modules/mol/alg/doc/molalg.rst
@@ -7,14 +7,14 @@
 Local Distance Test scores (lDDT, DRMSD)
 --------------------------------------------------------------------------------
 
-.. warning::
+.. note::
 
-  The code that comes with
+  This is a new implementation of lDDT, introduced in OpenStructure 2.4 with
+  focus on supporting quaternary structure and compounds beyond the 20
+  standard proteinogenic amino acids.
+  The :doc:`previous lDDT code <lddt_deprecated>` that comes with
   `Mariani et al. <https://dx.doi.org/10.1093/bioinformatics/btt473>`_ is
-  considered deprecated. lDDT has been re-implemented with a focus on
-  supporting quaternary structure and compounds beyond the 20 standard
-  proteinogenic amino acids. The old code is still available and documented
-  :doc:`here <lddt_deprecated>`.
+  considered deprecated.
 
 .. note::
 
@@ -135,16 +135,16 @@ Local Distance Test scores (lDDT, DRMSD)
 .. currentmodule:: ost.mol.alg
 
 
-:mod:`qsscoring <ost.mol.alg.qsscore>` -- QS score implementation
+:mod:`qsscore <ost.mol.alg.qsscore>` -- New QS score implementation
 --------------------------------------------------------------------------------
 
-.. warning::
+.. note::
 
-  The code that comes with
+  This is a new implementation of QS Score, introduced in OpenStructure 2.4 and
+  tightly integrated with the chain mapping algorithms.
+  The :doc:`previous qsscoring code <qsscoring_deprecated>` that comes with
   `Bertoni et al. <https://www.nature.com/articles/s41598-017-09654-8>`_ is
-  considered deprecated. QS score has been re-implemented to be tightly
-  integrated with the chain mapping algorithms. The old code is still available
-  and documented :doc:`here <qsscoring_deprecated>`.
+  considered deprecated.
 
 .. automodule:: ost.mol.alg.qsscore
    :members:
diff --git a/modules/mol/alg/pymod/qsscoring.py b/modules/mol/alg/pymod/qsscoring.py
index b47aa4e76e1544ba28caa18482253833b57106c9..8928646b98283584f2f5e5b82a0b4b8e0b49715d 100644
--- a/modules/mol/alg/pymod/qsscoring.py
+++ b/modules/mol/alg/pymod/qsscoring.py
@@ -2,6 +2,11 @@
 Scoring of quaternary structures (QS). The QS scoring is according to the paper
 by `Bertoni et al. <https://dx.doi.org/10.1038/s41598-017-09654-8>`_.
 
+.. warning::
+
+  The `qsscoring` module is deprecated. Consider using the newer implementation
+  in :mod:`~ost.mol.alg.qsscore` instead.
+
 .. note ::
 
   Requirements for use:
diff --git a/singularity/Singularity b/singularity/Singularity
index 1d1a0da33345fc2415fca5d32643f7a15b1bad5a..71f2c00fcad8ca8a2be1939865a2e7254ef0b486 100644
--- a/singularity/Singularity
+++ b/singularity/Singularity
@@ -1,5 +1,5 @@
 BootStrap: docker
-From: registry.scicore.unibas.ch/schwede/openstructure:2.4.0-jammy
+From: registry.scicore.unibas.ch/schwede/openstructure:2.5.0-jammy
 %post
 ##############################################################################
 # POST