diff --git a/CHANGELOG.txt b/CHANGELOG.txt
index fd3b597bfb594cdedbb695416f8192dbd8b71e47..f5ed1821edbc59b0a5222f77bbee25603b04d223 100644
--- a/CHANGELOG.txt
+++ b/CHANGELOG.txt
@@ -1,7 +1,38 @@
 Changes in Release 2.9.0
 --------------------------------------------------------------------------------
- 
+
  * lDDT-PLI now penalizes added model contacts by default.
+ * Updated unassigned reasons (model and target) to better reflect changes
+   to ligand scoring made in OST 2.8.0.
+ * Added CSV output (--output-format csv) and by model (rather than target)
+   ligand output (--by-model-ligand-output) to compare-ligand-structures
+   action.
+ * Improved logging and output of compare-ligand-structures action.
+ * Residue SetChemType() is now exposed in Python.
+ * Allow reading of BIRD compounds in PRDCC format in the compound library.
+   Compounds in PRD formats cannot be read and are rejected cleanly.
+ * The chemdict_tool executable now exits with an error status if no compound
+   were imported.
+ * The SDF writer populates the program name and time line (2) as per SDF
+   specification.
+ * The mmCIF reader now refuses to read in files with multiple data blocks
+   (except in fault tolerant mode), and warns about mmCIF files containing
+   more than one model (atom_site.pdbx_PDB_model_num).
+ * GDT produces slightly higher scores. It's an optimization problem in the
+   end to find the highest possible fraction of CA atoms that superpose within
+   a specified threshold. OpenStructure GDTTS scores are now on average 0.21
+   points lower (range [0, 100]) than LGA when computed on all CASP15 TS models.
+ * Enable DockQ on protein-nucleotide and nucleotide-nucleotide interfaces.
+ * Enable per-atom lDDT scores.
+ * Optionally compute ICS/IPS scores on trimmed model. Reason for that is
+   scoring with incomplete target structures. The full ICS/IPS scores may
+   consider contacts between residues that are not resolved as false positives,
+   even though there is no experimental evidence. The approach with trimmed
+   models removes any model residue that cannot be mapped to the target
+   structure before scoring.
+ * The SDF reader is now aware of IOProfiles and can read files with invalid
+   bond types (order) from RDKit in fault tolerant mode.
+ * Several bug fixes and improvements.
 
 Changes in Release 2.8.0
 --------------------------------------------------------------------------------
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1ff94f97da8a9eced92f710db9139ef91150ea2e..188ade39134d121af61b5cc7c57b5390a2d6d2a2 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 8)
+set (OST_VERSION_MINOR 9)
 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 77739dd4111e8ad1d8725289966aa3fb882a44f2..3c58acd5951f393ec3142219407045d99b9bba4c 100644
--- a/actions/ost-compare-ligand-structures
+++ b/actions/ost-compare-ligand-structures
@@ -283,6 +283,7 @@ def _ParseArgs():
         "--coverage-delta",
         dest="coverage_delta",
         default=0.2,
+        type=float,
         help=("Coverage delta for partial ligand assignment."))
 
     parser.add_argument(
@@ -314,6 +315,7 @@ def _ParseArgs():
         "--lddt-pli-radius",
         dest="lddt_pli_radius",
         default=6.0,
+        type=float,
         help=("lDDT inclusion radius for lDDT-PLI."))
 
     parser.add_argument(
@@ -343,6 +345,7 @@ def _ParseArgs():
         "--radius",
         dest="radius",
         default=4.0,
+        type=float,
         help=("Inclusion radius to extract reference binding site that is used "
               "for RMSD computation. Any residue with atoms within this "
               "distance of the ligand will be included in the binding site."))
@@ -351,6 +354,7 @@ def _ParseArgs():
         "--lddt-lp-radius",
         dest="lddt_lp_radius",
         default=15.0,
+        type=float,
         help=("lDDT inclusion radius for lDDT-LP."))
 
     parser.add_argument(
@@ -362,6 +366,15 @@ def _ParseArgs():
         help=("Enumerate all potential binding sites in the model when "
               "searching rigid superposition for RMSD computation"))
 
+    parser.add_argument(
+        "-ms",
+        "--max--symmetries",
+        dest="max_symmetries",
+        default=1e4,
+        type=int,
+        help=("If more than that many isomorphisms exist for a target-ligand "
+              "pair, it will be ignored and reported as unassigned."))
+
     args = parser.parse_args()
     if args.output is None:
         args.output = "out.%s" % args.output_format
@@ -512,7 +525,8 @@ def _SetupLDDTPLIScorer(model, model_ligands, reference, reference_ligands, args
                                                 substructure_match = args.substructure_match,
                                                 coverage_delta = args.coverage_delta,
                                                 lddt_pli_radius = args.lddt_pli_radius,
-                                                add_mdl_contacts = args.lddt_pli_add_mdl_contacts)
+                                                add_mdl_contacts = args.lddt_pli_add_mdl_contacts,
+                                                max_symmetries = args.max_symmetries)
 
 def _SetupSCRMSDScorer(model, model_ligands, reference, reference_ligands, args):
     return ligand_scoring_scrmsd.SCRMSDScorer(model, reference,
@@ -524,7 +538,8 @@ def _SetupSCRMSDScorer(model, model_ligands, reference, reference_ligands, args)
                                               coverage_delta = args.coverage_delta,
                                               bs_radius = args.radius,
                                               lddt_lp_radius = args.lddt_lp_radius,
-                                              full_bs_search = args.full_bs_search)
+                                              full_bs_search = args.full_bs_search,
+                                              max_symmetries = args.max_symmetries)
 
 def _Process(model, model_ligands, reference, reference_ligands, args):
 
diff --git a/docker/Dockerfile b/docker/Dockerfile
index 50a5b534eb98837fc5dd2f4d848f70e29a685fba..ec5a420b47cdf07e2c4821220cf84d78133eb73c 100644
--- a/docker/Dockerfile
+++ b/docker/Dockerfile
@@ -2,7 +2,7 @@ FROM ubuntu:22.04
 
 # ARGUMENTS
 ###########
-ARG OPENSTRUCTURE_VERSION="2.8.0"
+ARG OPENSTRUCTURE_VERSION="2.9.0"
 ARG SRC_FOLDER="/usr/local/src"
 ARG CPUS_FOR_MAKE=2
 
diff --git a/modules/doc/actions.rst b/modules/doc/actions.rst
index 89b96d8236bc38afa15739251785d51187e4c995..56afd04dcda6294dfeec0eefd5d64e1bced8fcd8 100644
--- a/modules/doc/actions.rst
+++ b/modules/doc/actions.rst
@@ -37,20 +37,22 @@ Details on the usage (output of ``ost compare-structures --help``):
                                 [-mb MODEL_BIOUNIT] [-rb REFERENCE_BIOUNIT]
                                 [-rna] [-ec] [-d] [-ds DUMP_SUFFIX] [-ft]
                                 [-c CHAIN_MAPPING [CHAIN_MAPPING ...]] [--lddt]
-                                [--local-lddt] [--bb-lddt] [--bb-local-lddt]
-                                [--ilddt] [--cad-score] [--local-cad-score]
-                                [--cad-exec CAD_EXEC]
+                                [--local-lddt] [--aa-local-lddt] [--bb-lddt]
+                                [--bb-local-lddt] [--ilddt] [--cad-score]
+                                [--local-cad-score] [--cad-exec CAD_EXEC]
                                 [--usalign-exec USALIGN_EXEC]
                                 [--override-usalign-mapping] [--qs-score]
                                 [--dockq] [--dockq-capri-peptide] [--ics]
-                                [--ips] [--rigid-scores] [--patch-scores]
-                                [--tm-score] [--lddt-no-stereochecks]
+                                [--ics-trimmed] [--ips] [--ips-trimmed]
+                                [--rigid-scores] [--patch-scores] [--tm-score]
+                                [--lddt-no-stereochecks]
                                 [--n-max-naive N_MAX_NAIVE]
                                 [--dump-aligned-residues] [--dump-pepnuc-alns]
                                 [--dump-pepnuc-aligned-residues]
                                 [--min-pep-length MIN_PEP_LENGTH]
                                 [--min-nuc-length MIN_NUC_LENGTH] [-v VERBOSITY]
                                 [--lddt-add-mdl-contacts]
+                                [--lddt-inclusion-radius LDDT_INCLUSION_RADIUS]
   
   Evaluate model against reference 
   
@@ -113,7 +115,9 @@ Details on the usage (output of ``ost compare-structures --help``):
    * "min_pep_length"
    * "min_nuc_length"
    * "lddt_add_mdl_contacts"
+   * "lddt_inclusion_radius"
    * "dockq_capri_peptide"
+   * "ost_version"
   
   The pairwise sequence alignments are computed with Needleman-Wunsch using
   BLOSUM62 (NUC44 for nucleotides). Many benchmarking scenarios preprocess the
@@ -329,6 +333,29 @@ Details on the usage (output of ``ost compare-structures --help``):
                           metrics are available as keys
                           "per_interface_ics_precision",
                           "per_interface_ics_recall" and "per_interface_ics".
+    --ics-trimmed         Computes interface contact similarity (ICS) related
+                          scores but on a trimmed model. That means that a
+                          mapping between model and reference is performed and
+                          all model residues without reference counterpart are
+                          removed. As a consequence, model contacts for which we
+                          have no experimental evidence do not affect the score.
+                          The effect of these added model contacts without
+                          mapping to target would be decreased precision and
+                          thus lower ics. Recall is not affected. Enabling this
+                          flag adds the following keys: "ics_trimmed",
+                          "ics_precision_trimmed", "ics_recall_trimmed",
+                          "model_contacts_trimmed". The reference contacts and
+                          reference interfaces are the same as for ics and
+                          available as keys: "reference_contacts",
+                          "contact_reference_interfaces". All these measures are
+                          also available on a per-interface basis for each
+                          interface in the reference structure that are defined
+                          as chain pairs with at least one contact (available as
+                          key "contact_reference_interfaces"). The respective
+                          metrics are available as keys
+                          "per_interface_ics_precision_trimmed",
+                          "per_interface_ics_recall_trimmed" and
+                          "per_interface_ics_trimmed".
     --ips                 Computes interface patch similarity (IPS) related
                           scores. They focus on interface residues. They are
                           defined as having at least one contact to a residue
@@ -352,6 +379,7 @@ Details on the usage (output of ``ost compare-structures --help``):
                           metrics are available as keys
                           "per_interface_ips_precision",
                           "per_interface_ips_recall" and "per_interface_ips".
+    --ips-trimmed         The IPS equivalent of ICS on trimmed models.
     --rigid-scores        Computes rigid superposition based scores. They're
                           based on a Kabsch superposition of all mapped CA
                           positions (C3' for nucleotides). Makes the following
@@ -425,7 +453,7 @@ Details on the usage (output of ``ost compare-structures --help``):
                           sequences can be problematic as they may produce high
                           sequence identity alignments by pure chance.
     -v VERBOSITY, --verbosity VERBOSITY
-                          Set verbosity level. Defaults to 3 (Script).
+                          Set verbosity level. Defaults to 2 (Script).
     --lddt-add-mdl-contacts
                           Only using contacts in lDDT thatare within a certain
                           distance threshold in the reference does not penalize
@@ -435,6 +463,9 @@ Details on the usage (output of ``ost compare-structures --help``):
                           necessarily in the reference. No contact will be added
                           if the respective atom pair is not resolved in the
                           reference.
+    --lddt-inclusion-radius LDDT_INCLUSION_RADIUS
+                          Passed to lDDT scorer. Affects all lDDT scores but not
+                          chain mapping.
 
 .. _ost compare ligand structures:
 
@@ -456,14 +487,19 @@ Details on the usage (output of ``ost compare-ligand-structures --help``):
                                        [-rl [REFERENCE_LIGANDS ...]] [-o OUTPUT]
                                        [-mf {pdb,cif,mmcif}]
                                        [-rf {pdb,cif,mmcif}] [-of {json,csv}]
-                                       [-csvm] [-mb MODEL_BIOUNIT]
+                                       [-csvm]
+                                       [--csv-extra-header CSV_EXTRA_HEADER]
+                                       [--csv-extra-data CSV_EXTRA_DATA]
+                                       [-mb MODEL_BIOUNIT]
                                        [-rb REFERENCE_BIOUNIT] [-ft] [-rna]
                                        [-sm] [-cd COVERAGE_DELTA] [-v VERBOSITY]
                                        [--full-results] [--lddt-pli]
                                        [--lddt-pli-radius LDDT_PLI_RADIUS]
-                                       [--lddt-pli-amc] [--rmsd]
+                                       [--lddt-pli-add-mdl-contacts]
+                                       [--no-lddt-pli-add-mdl-contacts] [--rmsd]
                                        [--radius RADIUS]
                                        [--lddt-lp-radius LDDT_LP_RADIUS] [-fbs]
+                                       [-ms MAX_SYMMETRIES]
   
   Evaluate model with non-polymer/small molecule ligands against reference.
   
@@ -500,8 +536,8 @@ Details on the usage (output of ``ost compare-ligand-structures --help``):
   
   Output can be written in two format: JSON (default) or CSV, controlled by the
   --output-format/-of argument.
-
-  Without additional options, the JSON ouput is a dictionary with three keys:
+  
+  Without additional options, the JSON ouput is a dictionary with four keys:
   
    * "model_ligands": A list of ligands in the model. If ligands were provided
      explicitly with --model-ligands, elements of the list will be the paths to
@@ -514,6 +550,7 @@ Details on the usage (output of ``ost compare-ligand-structures --help``):
    * "status": SUCCESS if everything ran through. In case of failure, the only
      content of the JSON output will be "status" set to FAILURE and an
      additional key: "traceback".
+   * "ost_version": The OpenStructure version used for computation.
   
   Each score is opt-in and the respective results are available in three keys:
   
@@ -538,32 +575,32 @@ Details on the usage (output of ``ost compare-ligand-structures --help``):
   items follow the same structure as in "assigned_scores". If no score for a
   specific pair of ligands could be computed, "score" and "coverage" are set to
   null and a key "reason" is added giving an educated guess why this happened.
-
+  
   CSV output is a table of comma-separated values, with one line for each
   reference ligand (or one model ligand if the --by-model-ligand-output flag was
   set).
-
+  
   The following column is always available:
-
+  
    * reference_ligand/model_ligand: If reference ligands were provided explicitly
      with --reference-ligands, elements of the list will be the paths to the
      ligand SDF file(s). Otherwise, they will be the chain name, residue number
      and insertion code of the ligand, separated by a dot. If the
      --by-model-ligand-output flag was set, this will be model ligand instead,
      following the same rules.
-
+  
   If lDDT-PLI was enabled with --lddt-pli, the following columns are added:
-
+  
    * "lddt_pli", "lddt_pli_coverage" and "lddt_pli_(model|reference)_ligand"
      are the lDDT-PLI score result, the corresponding coverage and assigned model
      ligand (or reference ligand if the --by-model-ligand-output flag was set)
      if an assignment was found, respectively, empty otherwise.
    * "lddt_pli_unassigned" is empty if an assignment was found, otherwise it
      lists the short reason this reference ligand was unassigned.
-
+  
   If BiSyRMSD was enabled with --rmsd, the following columns are added:
-
-   * "rmsd", "rmsd_coverage". "rmsd_lddt_lp" "rmsd_bb_rmsd" and
+  
+   * "rmsd", "rmsd_coverage". "lddt_lp" "bb_rmsd" and
      "rmsd_(model|reference)_ligand" are the BiSyRMSD, the corresponding
      coverage, lDDT-LP, backbone RMSD and assigned model ligand (or reference
      ligand if the --by-model-ligand-output flag was set) if an assignment
@@ -582,8 +619,8 @@ Details on the usage (output of ``ost compare-ligand-structures --help``):
     -rl [REFERENCE_LIGANDS ...], --ref-ligands [REFERENCE_LIGANDS ...], --reference-ligands [REFERENCE_LIGANDS ...]
                           Path to reference ligand files.
     -o OUTPUT, --out OUTPUT, --output OUTPUT
-                          Output file name. Default depends on format: out.json or
-                          out.csv
+                          Output file name. Default depends on format: out.json
+                          or out.csv
     -mf {pdb,cif,mmcif}, --mdl-format {pdb,cif,mmcif}, --model-format {pdb,cif,mmcif}
                           Format of model file. pdb reads pdb but also pdb.gz,
                           same applies to cif/mmcif. Inferred from filepath if
@@ -593,11 +630,18 @@ Details on the usage (output of ``ost compare-ligand-structures --help``):
                           pdb.gz, same applies to cif/mmcif. Inferred from
                           filepath if not given.
     -of {json,csv}, --out-format {json,csv}, --output-format {json,csv}
-                          Output format, JSON or CSV, in lowercase. default: json
+                          Output format, JSON or CSV, in lowercase. default:
+                          json
     -csvm, --by-model-ligand, --by-model-ligand-output
                           For CSV output, this flag changes the output so that
                           each line reports one model ligand, instead of a
                           reference ligand. Has no effect with JSON output.
+    --csv-extra-header CSV_EXTRA_HEADER
+                          Extra header prefix for CSV output. This allows adding
+                          additional annotations (such as target ID, group, etc)
+                          to the output
+    --csv-extra-data CSV_EXTRA_DATA
+                          Additional data (columns) for CSV output.
     -mb MODEL_BIOUNIT, --model-biounit MODEL_BIOUNIT
                           Only has an effect if model is in mmcif format. By
                           default, the asymmetric unit (AU) is used for scoring.
@@ -622,13 +666,16 @@ Details on the usage (output of ``ost compare-ligand-structures --help``):
     -cd COVERAGE_DELTA, --coverage-delta COVERAGE_DELTA
                           Coverage delta for partial ligand assignment.
     -v VERBOSITY, --verbosity VERBOSITY
-                          Set verbosity level. Defaults to 3 (INFO).
+                          Set verbosity level. Defaults to 2 (Script).
     --full-results        Outputs scoring results for all model/reference ligand
                           pairs and store as key "full_results"
     --lddt-pli            Compute lDDT-PLI scores and store as key "lddt_pli".
     --lddt-pli-radius LDDT_PLI_RADIUS
                           lDDT inclusion radius for lDDT-PLI.
-    --lddt-pli-amc        Add model contacts (amc) when computing lDDT-PLI.
+    --lddt-pli-add-mdl-contacts
+                          Add model contacts when computing lDDT-PLI.
+    --no-lddt-pli-add-mdl-contacts
+                          DO NOT add model contacts when computing lDDT-PLI.
     --rmsd                Compute RMSD scores and store as key "rmsd".
     --radius RADIUS       Inclusion radius to extract reference binding site
                           that is used for RMSD computation. Any residue with
@@ -640,3 +687,8 @@ Details on the usage (output of ``ost compare-ligand-structures --help``):
                           Enumerate all potential binding sites in the model
                           when searching rigid superposition for RMSD
                           computation
+    -ms MAX_SYMMETRIES, --max--symmetries MAX_SYMMETRIES
+                          If more than that many isomorphisms exist for a
+                          target-ligand pair, it will be ignored and reported as
+                          unassigned.
+
diff --git a/modules/io/doc/io.rst b/modules/io/doc/io.rst
index d4c364c2ce328aedf03b35aa1ce16c115a3bc49e..933f33cc3639d771808b3be71c64a26b5d23f7be 100644
--- a/modules/io/doc/io.rst
+++ b/modules/io/doc/io.rst
@@ -119,22 +119,18 @@ behaviour.
 
   :rtype: :class:`~ost.mol.EntityHandle`.
 
+.. autofunction:: ost.io.LoadSDF
 
-.. function:: LoadSDF(filename)
-
-  Load an SDF file and return an entity.
-
-  :param filename: File to be loaded
-  :type filename: :class:`str`
-
-  :rtype: :class:`~ost.mol.EntityHandle`
-
-.. function:: SDFStrToEntity(sdf_string)
+.. function:: SDFStrToEntity(sdf_string, profile=IOProfile())
 
   Load entity from a string in SDF format.
 
   :param pdb_string: A SDF file as a string.
 
+  :param profile: The IO Profile to read the entity with. For more information
+      on the IO Profiles available, see :doc:`profile`.
+  :type profile: :class:`ost.io.IOProfile`
+
   :rtype: :class:`~ost.mol.EntityHandle`.
 
 .. class:: ost.io.OMF
diff --git a/modules/io/doc/mmcif.rst b/modules/io/doc/mmcif.rst
index 5c18744549199a3094bb4f7d69b6cd7e12369fc4..23bd2a530692788147b3f248578710a3773d5e34 100644
--- a/modules/io/doc/mmcif.rst
+++ b/modules/io/doc/mmcif.rst
@@ -105,7 +105,7 @@ Notes:
   It is a known limitation of the mmCIF format to allow ambiguous identifiers
   for waters (and ligands to some extent) and so we have to require these
   additional identifiers.
-* An mmCIF file can contain several models (``atom.site.pdbx_PDB_model_num``).
+* An mmCIF file can contain several models (``atom_site.pdbx_PDB_model_num``).
   Only the first model occurring in the mmCIF file is read (regardless of the
   actual model number). If extra models are ignored, a warning is logged.
 
diff --git a/modules/io/doc/structure_formats.rst b/modules/io/doc/structure_formats.rst
index 2fe0eaf676f9bc5f146225563a398011e0614a71..6903ed50f015523be82296cd0b61302bada00a86 100644
--- a/modules/io/doc/structure_formats.rst
+++ b/modules/io/doc/structure_formats.rst
@@ -51,12 +51,13 @@ radii.
   
 SDF - Structured Data File
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-Chemical table (Ctab) file format (V2000; read-only V3000 experimental support),
-aka MDL Molfile.
-The SDF format does not support residues, chains or atom names natively.
+Chemical table (Ctab) file format, aka SDF, aka MDL Molfile, originally named
+after Molecular Design Limited, Inc. Follows the "BIOVIA Databases 2020" format
+specification (V2000; read-only V3000 experimental support).
 The SDF importer supports loading gzipped files, which are auto-detected by the
 .gz file extension.
 
+The SDF format does not support residues, chains or atom names natively.
 The reader assigns 1-based atom indices as atom names.
 SDF files containing several molecules are loaded into distinct chains,
 named after the molecule name in the MOLfile header with a numerical prefix.
@@ -66,5 +67,5 @@ Chains are written as separate molecules. If a chain contains more than one
 residue, they will be merged into a single molecule.
 
 *Recognized File Extensions*
-  .sdf, .sdf.gz
+  .sdf, .sdf.gz, .mol, .mol.gz
   
diff --git a/modules/io/pymod/CMakeLists.txt b/modules/io/pymod/CMakeLists.txt
index ead77d54b9d4cc0285845a166a52287774d73f67..02b392754d6b44ec9948b621ef86903ebf3ce060 100644
--- a/modules/io/pymod/CMakeLists.txt
+++ b/modules/io/pymod/CMakeLists.txt
@@ -4,6 +4,7 @@ set(OST_IO_PYMOD_SOURCES
   export_mmcif_io.cc
   export_omf_io.cc
   export_map_io.cc
+  export_sdf_io.cc
 )
 
 set(OST_IO_PYMOD_MODULES
diff --git a/modules/io/pymod/__init__.py b/modules/io/pymod/__init__.py
index ea3be23a85b688a37946a55d298681aa5472979d..dca962c0bcd74611721a5e439e28d6eb1511f2d0 100644
--- a/modules/io/pymod/__init__.py
+++ b/modules/io/pymod/__init__.py
@@ -556,3 +556,46 @@ def _PDBize(biounit, asu, seqres=None, min_polymer_size=None,
   return pdb_bu
 
 MMCifInfoBioUnit.PDBize = _PDBize
+
+
+def LoadSDF(filename, fault_tolerant=None, profile='DEFAULT'):
+  """
+  Load SDF file from disk and return an entity.
+
+  :param filename: File to be loaded
+  :type filename: :class:`str`
+
+  :param fault_tolerant: Enable/disable fault-tolerant import. If set, overrides
+                         the value of :attr:`IOProfile.fault_tolerant`.
+  :type fault_tolerant: :class:`bool`
+
+  :param profile: Aggregation of flags and algorithms to control import and
+                  processing of molecular structures. Can either be a
+                  :class:`str` specifying one of the default profiles
+                  ['DEFAULT', 'SLOPPY', 'CHARMM', 'STRICT'] or an actual object
+                  of type :class:`ost.io.IOProfile`.
+                  See :doc:`profile` for more info.
+  :type profile: :class:`str`/:class:`ost.io.IOProfile`
+
+  :raises: :exc:`~ost.io.IOException` if the import fails due to an erroneous or
+      inexistent file
+  """
+  def _override(val1, val2):
+    if val2 != None:
+      return val2
+    else:
+      return val1
+
+  if isinstance(profile, str):
+    prof = profiles.Get(profile)
+  elif isinstance(profile, IOProfile):
+    prof = profile.Copy()
+  else:
+    raise TypeError('profile must be of type string or IOProfile, ' + \
+                    'instead of %s' % type(profile))
+  prof.fault_tolerant = _override(prof.fault_tolerant, fault_tolerant)
+
+  reader = SDFReader(filename, prof)
+  ent = mol.CreateEntity()
+  reader.Import(ent)
+  return ent
diff --git a/modules/io/pymod/export_sdf_io.cc b/modules/io/pymod/export_sdf_io.cc
new file mode 100644
index 0000000000000000000000000000000000000000..08de790e2850229a568920c790767992491c075c
--- /dev/null
+++ b/modules/io/pymod/export_sdf_io.cc
@@ -0,0 +1,64 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2020 by the OpenStructure authors
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation; either version 3.0 of the License, or (at your option)
+// any later version.
+// This library is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this library; if not, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+//------------------------------------------------------------------------------
+#include <boost/python.hpp>
+#include <boost/shared_ptr.hpp>
+#include <boost/python/suite/indexing/vector_indexing_suite.hpp>
+#include <boost/python/suite/indexing/map_indexing_suite.hpp>
+using namespace boost::python;
+
+#include <ost/export_helper/pair_to_tuple_conv.hh>
+#include <ost/io/mol/entity_io_sdf_handler.hh>
+#include <ost/io/mol/io_profile.hh>
+#include <ost/io/mol/sdf_reader.hh>
+#include <ost/io/mol/sdf_writer.hh>
+#include <ost/io/sdf_str.hh>
+using namespace ost;
+using namespace ost::io;
+using namespace ost::mol;
+
+String (*sdf_str_a)(const mol::EntityHandle&)=&EntityToSDFString;
+String (*sdf_str_b)(const mol::EntityView&)=&EntityToSDFString;
+
+void (*save_sdf_handle)(const mol::EntityHandle& entity, const String& filename)=&SaveSDF;
+void (*save_sdf_view)(const mol::EntityView& entity, const String& filename)=&SaveSDF;
+
+void (SDFWriter::*write_handle)(const mol::EntityHandle&)=&SDFWriter::Write;
+void (SDFWriter::*write_view)(const mol::EntityView&)=&SDFWriter::Write;
+
+void export_sdf_io()
+{
+  class_<SDFReader, boost::noncopyable>("SDFReader", init<String, const IOProfile&>())
+    .def("Import", &SDFReader::Import)
+  ;
+  
+  class_<SDFWriter, boost::noncopyable>("SDFWriter", init<String>())
+    .def("Write", write_handle)
+    .def("Write", write_view)
+  ;
+
+  // def("LoadSDF", &LoadSDF);
+  def("SaveSDF", save_sdf_view);
+  def("SaveSDF", save_sdf_handle);
+
+  def("EntityToSDFStr", sdf_str_a);
+  def("EntityToSDFStr", sdf_str_b);
+
+  def("SDFStrToEntity", &SDFStringToEntity, (arg("SDF_string"),
+                                             arg("profile")=IOProfile()));
+}
diff --git a/modules/io/pymod/wrap_io.cc b/modules/io/pymod/wrap_io.cc
index bae73b1e0fd9168c7673fc39da52ce678d6878d9..b7cade00c33e39eefd2ecc982b7e2c4e2d5a2a6f 100644
--- a/modules/io/pymod/wrap_io.cc
+++ b/modules/io/pymod/wrap_io.cc
@@ -30,9 +30,7 @@ using namespace boost::python;
 #include <ost/io/mol/entity_io_crd_handler.hh>
 #include <ost/io/mol/entity_io_pqr_handler.hh>
 #include <ost/io/mol/entity_io_mae_handler.hh>
-#include <ost/io/mol/entity_io_sdf_handler.hh>
 #include <ost/io/mol/pdb_reader.hh>
-#include <ost/io/mol/sdf_str.hh>
 #include <ost/io/mol/dcd_io.hh>
 #include <ost/io/stereochemical_params_reader.hh>
 using namespace ost;
@@ -66,18 +64,13 @@ BOOST_PYTHON_FUNCTION_OVERLOADS(save_entity_view_ov,
 ost::mol::alg::StereoChemicalProps (*read_props_a)(String filename, bool check) = &ReadStereoChemicalPropsFile;
 ost::mol::alg::StereoChemicalProps (*read_props_b)(bool check) = &ReadStereoChemicalPropsFile;
 
-String (*sdf_str_a)(const mol::EntityHandle&)=&EntityToSDFString;
-String (*sdf_str_b)(const mol::EntityView&)=&EntityToSDFString;
-
-void (*save_sdf_handle)(const mol::EntityHandle& entity, const String& filename)=&SaveSDF;
-void (*save_sdf_view)(const mol::EntityView& entity, const String& filename)=&SaveSDF;
-
 }
 
 void export_pdb_io();
 void export_mmcif_io();
 void export_omf_io();
 void export_map_io();
+void export_sdf_io();
 BOOST_PYTHON_MODULE(_ost_io)
 {
   class_<IOManager, boost::noncopyable>("IOManager", no_init)
@@ -122,14 +115,6 @@ BOOST_PYTHON_MODULE(_ost_io)
       (arg("seq_list"), arg("filename"), arg("format")="auto"));
   def("SaveSequence", &SaveSequence,
       (arg("sequence"), arg("filename"), arg("format")="auto"));
-  def("LoadSDF", &LoadSDF);
-  def("SaveSDF", save_sdf_view);
-  def("SaveSDF", save_sdf_handle);
-
-  def("EntityToSDFStr", sdf_str_a);
-  def("EntityToSDFStr", sdf_str_b);
-
-  def("SDFStrToEntity", &SDFStringToEntity);
 
   def("LoadCRD", &LoadCRD);
   def("LoadCHARMMTraj_", &LoadCHARMMTraj, (arg("ent"), arg("trj_filename"), 
@@ -149,6 +134,7 @@ BOOST_PYTHON_MODULE(_ost_io)
   export_mmcif_io();
   export_omf_io();
   export_map_io();
+  export_sdf_io();
   def("SaveCHARMMTraj", &SaveCHARMMTraj, 
       (arg("traj"), arg("pdb_filename"), arg("dcd_filename"), arg("stride")=1, 
        arg("profile")=IOProfile()));
diff --git a/modules/io/src/mol/entity_io_sdf_handler.cc b/modules/io/src/mol/entity_io_sdf_handler.cc
index 129ce2c86be7a7ec56e1cbee52488a231411067a..fed25f86b2b6c745341a51a4fe6fa770feff500f 100644
--- a/modules/io/src/mol/entity_io_sdf_handler.cc
+++ b/modules/io/src/mol/entity_io_sdf_handler.cc
@@ -21,6 +21,7 @@
  */
 
 #include <ost/log.hh>
+#include <ost/profile.hh>
 #include <ost/io/mol/sdf_writer.hh>
 #include <ost/io/mol/sdf_reader.hh>
 
@@ -37,14 +38,14 @@ bool EntityIOSDFHandler::RequiresProcessor() const
 void EntityIOSDFHandler::Import(mol::EntityHandle& ent,
                                 std::istream& instream)
 {
-  SDFReader reader(instream);
+  SDFReader reader(instream, IOProfileRegistry::Instance().GetDefault());
   reader.Import(ent);
 }
 
 void EntityIOSDFHandler::Import(mol::EntityHandle& ent,
                                 const boost::filesystem::path& loc)
 {
-  SDFReader reader(loc);
+  SDFReader reader(loc, IOProfileRegistry::Instance().GetDefault());
   reader.Import(ent);
 }
 
@@ -69,7 +70,8 @@ bool sdf_handler_is_responsible_for(const boost::filesystem::path& loc,
   if(type=="auto") {
 	String match_suf_string=loc.string();
     std::transform(match_suf_string.begin(),match_suf_string.end(),match_suf_string.begin(),tolower);
-    if(detail::FilenameEndsWith(match_suf_string,".sdf") || detail::FilenameEndsWith(match_suf_string,".sdf.gz")) {
+    if(detail::FilenameEndsWith(match_suf_string,".sdf") || detail::FilenameEndsWith(match_suf_string,".sdf.gz") ||
+       detail::FilenameEndsWith(match_suf_string,".mol") || detail::FilenameEndsWith(match_suf_string,".mol.gz")) {
       return true;
     }
 
diff --git a/modules/io/src/mol/mmcif_reader.cc b/modules/io/src/mol/mmcif_reader.cc
index 077e4c0727f59a929dbd1110607ed14f2bfcf713..03abbad79f9b940e99c556c7ab0bb654d20c648c 100644
--- a/modules/io/src/mol/mmcif_reader.cc
+++ b/modules/io/src/mol/mmcif_reader.cc
@@ -92,7 +92,7 @@ void MMCifReader::SetRestrictChains(const String& restrict_chains)
 
 bool MMCifReader::OnBeginData(const StringRef& data_name) 
 {
-  LOG_DEBUG("MCIFFReader: " << profile_);
+  LOG_DEBUG("MMCifReader: " << profile_);
   Profile profile_import("MMCifReader::OnBeginData");
   if (chain_count_ > 0) {
       std::stringstream ss;
diff --git a/modules/io/src/mol/sdf_reader.cc b/modules/io/src/mol/sdf_reader.cc
index 70ee382b1671388ee9a22aab0ca43020364a76f6..edafc0ff1cb401af7aa63878a8f5074799aa3bf0 100644
--- a/modules/io/src/mol/sdf_reader.cc
+++ b/modules/io/src/mol/sdf_reader.cc
@@ -37,20 +37,20 @@ namespace ost { namespace io {
 
 using boost::format;
 
-SDFReader::SDFReader(const String& filename)
-  : infile_(filename), instream_(infile_)
+SDFReader::SDFReader(const String& filename, const IOProfile& profile)
+  : infile_(filename), instream_(infile_), profile_(profile)
 {
   this->ClearState(boost::filesystem::path(filename));
 }
 
-SDFReader::SDFReader(const boost::filesystem::path& loc)
-  : infile_(loc), instream_(infile_)
+SDFReader::SDFReader(const boost::filesystem::path& loc, const IOProfile& profile)
+  : infile_(loc), instream_(infile_), profile_(profile)
 {
   this->ClearState(loc);
 }
 
-SDFReader::SDFReader(std::istream& instream)
-  : infile_(), instream_(instream)
+SDFReader::SDFReader(std::istream& instream, const IOProfile& profile)
+  : infile_(), instream_(instream), profile_(profile)
 {
   this->ClearState(boost::filesystem::path(""));
 }
@@ -206,7 +206,7 @@ void SDFReader::ParseHeader(const String& line, int line_num,
         version_=version_str;
       }
       else {
-        String msg="Unsupported SDF version: %s.";
+        String msg="Invalid SDF file or unsupported SDF version: %s.";
         throw IOException(str(format(msg) % version_str));
       }
       // Counts will be overridden in V3000
@@ -355,15 +355,40 @@ void SDFReader::AddBond(const bond_data& bond_tuple, int line_num, mol::EntityHa
 
   try {
     type=boost::lexical_cast<int>(boost::trim_copy(s_type));
-    if (type<1 || type>8) {
-      String msg="Bad bond line %d: Bond type number"
-                      " '%s' not within accepted range (1-8).";
-      throw IOException(str(format(msg) % line_num % s_type));
+    // From SDF spec:
+    // bond type      1 = Single, 2 = Double,       [Q] Values 4 through 8 are
+    //                3 = Triple, 4 = Aromatic,     for SSS queries oniy.
+    //                5 = Single or Double,
+    //                6 = Single or Aromatic,
+    //                7 = Double or Aromatic,
+    //                8 = Any
+    if (type < 1 || type > 8) {
+      std::stringstream ss;
+      ss << "Bad bond line " << line_num << ": Bond type number "
+         << std::to_string(type) << " not within accepted range (1-8).";
+      if (profile_.fault_tolerant) {
+        LOG_ERROR(ss.str());
+      } else {
+        throw IOException(ss.str());
+      }
+    } else if (type > 3) {
+      std::stringstream ss;
+      ss << "Bad bond line " << line_num << ": Bond type number "
+         << std::to_string(type) << ": values 4-8 are reserved for queries, "
+         << "should not appear in an SDF file.";
+        LOG_WARNING(ss.str());
     }
   } catch(boost::bad_lexical_cast&) {
-    String msg="Bad bond line %d: Can't convert bond type number"
-                " '%s' to integral constant.";
-    throw IOException(str(format(msg) % line_num % s_type));
+    std::stringstream ss;
+    ss << "Bad bond line " << line_num << ": Can't convert bond type number '"
+       << s_type << "' to integral constant.";
+    if (profile_.fault_tolerant) {
+      ss << " Assuming single bond in fault tolerant mode.";
+      LOG_ERROR(ss.str());
+      type = 1;
+    } else {
+      throw IOException(ss.str());
+    }
   }
 
   mol::AtomHandle first,second;
diff --git a/modules/io/src/mol/sdf_reader.hh b/modules/io/src/mol/sdf_reader.hh
index 59e733a937a3ed0dcb7171deec9ff9a32dd808ae..d524a8bf1561323dcbb2b91c10c70097673b352b 100644
--- a/modules/io/src/mol/sdf_reader.hh
+++ b/modules/io/src/mol/sdf_reader.hh
@@ -28,6 +28,7 @@
 #include <ost/mol/chain_handle.hh>
 #include <ost/mol/residue_handle.hh>
 #include <ost/io/module_config.hh>
+#include <ost/io/mol/io_profile.hh>
 
 namespace ost { namespace io {
 
@@ -36,11 +37,9 @@ namespace ost { namespace io {
 
 class DLLEXPORT_OST_IO SDFReader {
 public:
-  SDFReader(const String& filename);
-  SDFReader(const boost::filesystem::path& loc);
-  SDFReader(std::istream& instream);
-
-  bool HasNext();
+  SDFReader(const String& filename, const IOProfile& profile);
+  SDFReader(const boost::filesystem::path& loc, const IOProfile& profile);
+  SDFReader(std::istream& instream, const IOProfile& profile);
 
   void Import(mol::EntityHandle& ent);
 
@@ -97,6 +96,7 @@ private:
   boost::filesystem::ifstream infile_;
   std::istream& instream_;
   boost::iostreams::filtering_stream<boost::iostreams::input>  in_;
+  IOProfile profile_;
   String version_;
   bool v3000_atom_block_;
   bool v3000_bond_block_;
diff --git a/modules/io/src/mol/sdf_str.cc b/modules/io/src/mol/sdf_str.cc
index a2977c432bf3ec90647de1b05dc8eb5a2afb3f74..0441b672f0edbc3604888cec50abc35f75e0e54b 100644
--- a/modules/io/src/mol/sdf_str.cc
+++ b/modules/io/src/mol/sdf_str.cc
@@ -37,9 +37,9 @@ String EntityToSDFString(const mol::EntityView& ent) {
   return stream.str();
 }
 
-mol::EntityHandle SDFStringToEntity(const String& sdf) {
+mol::EntityHandle SDFStringToEntity(const String& sdf, const IOProfile& profile) {
   std::stringstream stream(sdf);
-  SDFReader reader(stream);
+  SDFReader reader(stream, profile);
   mol::EntityHandle ent = mol::CreateEntity();
   reader.Import(ent);
   return ent;
diff --git a/modules/io/src/mol/sdf_str.hh b/modules/io/src/mol/sdf_str.hh
index 87987679ed7ad28a6e3023f536ad6c6fb716f7f7..8cc6974593a202be792cc8fcf788fae6f9b780df 100644
--- a/modules/io/src/mol/sdf_str.hh
+++ b/modules/io/src/mol/sdf_str.hh
@@ -22,6 +22,7 @@
 #include <ost/io/module_config.hh>
 #include <ost/mol/entity_view.hh>
 #include <ost/mol/entity_handle.hh>
+#include <ost/io/mol/io_profile.hh>
 
 namespace ost { namespace io {
 
@@ -33,7 +34,7 @@ String DLLEXPORT_OST_IO
 EntityToSDFString(const mol::EntityView& ent);
 
 mol::EntityHandle DLLEXPORT_OST_IO
-SDFStringToEntity(const String& pdb);
+SDFStringToEntity(const String& pdb, const IOProfile& profile);
 
 }}
 
diff --git a/modules/io/tests/test_io_sdf.py b/modules/io/tests/test_io_sdf.py
index 7277a399d5a58ca696f50fb19e0775345999313a..40cee207926556fea2916d0a420fcf7f421dee24 100644
--- a/modules/io/tests/test_io_sdf.py
+++ b/modules/io/tests/test_io_sdf.py
@@ -47,7 +47,41 @@ class TestSDF(unittest.TestCase):
     # Charge from atom line is ignored
     o_at = ent.FindAtom("00001_Simple Ligand", 1, "3")
     self.assertEqual(o_at.charge, 0)
-    
+
+  def test_fault_tolerant(self):
+    """This file has a "dative" bond (type = 9).
+    This is a non-standard extension from RDKit which should go through only
+    in fault tolerant mode"""
+
+    with self.assertRaises(Exception):
+      ent = io.LoadSDF('testfiles/sdf/dative_bond.sdf')
+
+    # Directly with fault_tolerant
+    PushVerbosityLevel(-1)  # Expect message at Error level
+    ent = io.LoadSDF('testfiles/sdf/dative_bond.sdf', fault_tolerant=True)
+    PopVerbosityLevel()
+    self.assertEqual(ent.FindAtom("00001_Simple Ligand", 1, "5").bonds[0].bond_order, 9)
+
+    # Sloppy profile
+    PushVerbosityLevel(-1)  # Expect message at Error level
+    ent = io.LoadSDF('testfiles/sdf/dative_bond.sdf', profile="SLOPPY")
+    PopVerbosityLevel()
+    self.assertEqual(ent.FindAtom("00001_Simple Ligand", 1, "5").bonds[0].bond_order, 9)
+
+    # Sloppy profile set as default
+    old_profile = io.profiles['DEFAULT'].Copy()
+    io.profiles['DEFAULT'] = "SLOPPY"
+    PushVerbosityLevel(-1)  # Expect message at Error level
+    ent = io.LoadSDF('testfiles/sdf/dative_bond.sdf')
+    PopVerbosityLevel()
+    self.assertEqual(ent.FindAtom("00001_Simple Ligand", 1, "5").bonds[0].bond_order, 9)
+
+    # Test that a restored default profile has fault_tolerant again
+    io.profiles['DEFAULT'] = old_profile
+    with self.assertRaises(Exception):
+      ent = io.LoadSDF('testfiles/sdf/dative_bond.sdf')
+
+
 if __name__== '__main__':
   from ost import testutils
   testutils.RunTests()
diff --git a/modules/io/tests/testfiles/sdf/dative_bond.sdf b/modules/io/tests/testfiles/sdf/dative_bond.sdf
new file mode 100644
index 0000000000000000000000000000000000000000..30ac15cde8ddd239f13f7c47b33f584df0533591
--- /dev/null
+++ b/modules/io/tests/testfiles/sdf/dative_bond.sdf
@@ -0,0 +1,18 @@
+Simple Ligand
+
+ Teststructure
+  6  6  0  0  1  0            999 V2000
+    0.0000    0.0000    0.0000 N   0  3  0  0  0  0
+    1.0000    0.0000    0.0000 C   0  0  0  0  0  0
+    0.0000    1.0000    0.0000 O   0  0  0  0  0  0
+    1.0000    1.0000    0.0000 S   0  0  0  0  0  0
+    2.0000    2.0000    0.0000 C   0  0  0  0  0  0
+   -1.0000   -1.0000    0.0000 Cl  0  0  0  0  0  0
+  1  2  2  0  0  0
+  1  3  1  0  0  0
+  1  6  1  0  0  0
+  2  4  1  0  0  0
+  3  4  1  0  0  0
+  4  5  9  0  0  0
+M  END
+$$$$
diff --git a/modules/mol/alg/doc/qsscore.rst b/modules/mol/alg/doc/qsscore.rst
index 7b6f1711f1441ce8d7631c1bb28fab471d3568c6..4c0a7c7f9589720a6c2f76e3178e8df0e12dafd5 100644
--- a/modules/mol/alg/doc/qsscore.rst
+++ b/modules/mol/alg/doc/qsscore.rst
@@ -2,6 +2,7 @@
 --------------------------------------------------------------------------------
 
 .. module:: ost.mol.alg.qsscore
+   :noindex:
    :synopsis: New QS score implementation
 
 .. note::
diff --git a/modules/mol/alg/pymod/chain_mapping.py b/modules/mol/alg/pymod/chain_mapping.py
index 74094ac2131f1a82afcd0073141c51b21508812a..d30d2bcec66bd5239de73e0b84c4536c0e154262 100644
--- a/modules/mol/alg/pymod/chain_mapping.py
+++ b/modules/mol/alg/pymod/chain_mapping.py
@@ -190,9 +190,8 @@ class ReprResult:
         self._mdl_bb_pos = None
         self._ref_full_bb_pos = None
         self._mdl_full_bb_pos = None
-        self._transform = None
+        self._superposition = None
         self._superposed_mdl_bb_pos = None
-        self._bb_rmsd = None
         self._ost_query = None
         self._flat_mapping = None
         self._inconsistent_residues = None
@@ -314,24 +313,33 @@ class ReprResult:
             self._mdl_full_bb_pos = self._GetFullBBPos(self.mdl_residues)
         return self._mdl_full_bb_pos
 
+
     @property
-    def transform(self):
-        """ Transformation to superpose mdl residues onto ref residues
+    def superposition(self):
+        """ Superposition of mdl residues onto ref residues
 
         Superposition computed as minimal RMSD superposition on
         :attr:`ref_bb_pos` and :attr:`mdl_bb_pos`. If number of positions is
         smaller 3, the full_bb_pos equivalents are used instead.
 
-        :type: :class:`ost.geom.Mat4`
+        :type: :class:`ost.mol.alg.SuperpositionResult`
         """
-        if self._transform is None:
+        if self._superposition is None:
             if len(self.mdl_bb_pos) < 3:
-                self._transform = _GetTransform(self.mdl_full_bb_pos,
+                self._superposition = _GetSuperposition(self.mdl_full_bb_pos,
                                                 self.ref_full_bb_pos, False)
             else:
-                self._transform = _GetTransform(self.mdl_bb_pos,
+                self._superposition = _GetSuperposition(self.mdl_bb_pos,
                                                 self.ref_bb_pos, False)
-        return self._transform
+        return self._superposition
+
+    @property
+    def transform(self):
+        """ Transformation to superpose mdl residues onto ref residues
+
+        :type: :class:`ost.geom.Mat4`
+        """
+        return self.superposition.transformation
 
     @property
     def superposed_mdl_bb_pos(self):
@@ -346,13 +354,11 @@ class ReprResult:
 
     @property
     def bb_rmsd(self):
-        """ RMSD between :attr:`ref_bb_pos` and :attr:`superposed_mdl_bb_pos`
+        """ RMSD of the binding site backbone atoms after :attr:`superposition`
 
         :type: :class:`float`
         """
-        if self._bb_rmsd is None:
-            self._bb_rmsd = self.ref_bb_pos.GetRMSD(self.superposed_mdl_bb_pos)
-        return self._bb_rmsd
+        return self.superposition.rmsd
 
 
     @property
@@ -1114,7 +1120,8 @@ class ChainMapper:
               for t_pos, t in zip(trg_pos, trg_chains):
                   for m_pos, m in zip(mdl_pos, mdl_chains):
                       if len(t_pos) >= 3 and len(m_pos) >= 3:
-                          transform = _GetTransform(m_pos, t_pos, False)
+                          transform = _GetSuperposition(m_pos, t_pos,
+                                                        False).transformation
                           initial_transforms.append(transform)
                           initial_mappings.append((t,m))
 
@@ -2903,8 +2910,8 @@ def _IterativeRigidRMSD(initial_transforms, initial_mappings, chem_groups,
                 trg_chain_groups[best_sc_group_idx].remove(best_sc_mapping[0])
                 mdl_chain_groups[best_sc_group_idx].remove(best_sc_mapping[1])
 
-                transform = _GetTransform(mapped_mdl_pos, mapped_trg_pos,
-                                          False)
+                transform = _GetSuperposition(mapped_mdl_pos, mapped_trg_pos,
+                                              False).transformation
 
         # compute overall RMSD for current transform
         mapped_mdl_pos.ApplyTransform(transform)
@@ -3219,7 +3226,7 @@ def _ChainMappings(ref_chains, mdl_chains, n_max=None):
     return _ConcatIterators(iterators)
 
 
-def _GetTransform(pos_one, pos_two, iterative):
+def _GetSuperposition(pos_one, pos_two, iterative):
     """ Computes minimal RMSD superposition for pos_one onto pos_two
 
     :param pos_one: Positions that should be superposed onto *pos_two*
@@ -3230,7 +3237,7 @@ def _GetTransform(pos_one, pos_two, iterative):
                 potentially raises, uses standard superposition as fallback.
     :type iterative: :class:`bool`
     :returns: Transformation matrix to superpose *pos_one* onto *pos_two*
-    :rtype: :class:`geom.Mat4`
+    :rtype: :class:`ost.mol.alg.SuperpositionResult`
     """
     res = None
     if iterative:
@@ -3240,7 +3247,7 @@ def _GetTransform(pos_one, pos_two, iterative):
             pass # triggers fallback below
     if res is None:
         res = mol.alg.SuperposeSVD(pos_one, pos_two)
-    return res.transformation
+    return res
 
 # specify public interface
 __all__ = ('ChainMapper', 'ReprResult', 'MappingResult')
diff --git a/modules/mol/alg/pymod/ligand_scoring_base.py b/modules/mol/alg/pymod/ligand_scoring_base.py
index 8710b05d1452f1ee17c1ee68efc2513f77effa23..45376fde5f3d43a7fffb7ce7bbe518b570e8ea9a 100644
--- a/modules/mol/alg/pymod/ligand_scoring_base.py
+++ b/modules/mol/alg/pymod/ligand_scoring_base.py
@@ -43,9 +43,10 @@ class LigandScorer:
 
     * :class:`ost.mol.alg.ligand_scoring_lddtpli.LDDTPLIScorer`
       that assesses the conservation of protein-ligand
-      contacts
+      contacts (lDDT-PLI);
     * :class:`ost.mol.alg.ligand_scoring_scrmsd.SCRMSDScorer`
-      that computes a binding-site superposed, symmetry-corrected RMSD.
+      that computes a binding-site superposed, symmetry-corrected RMSD
+      (BiSyRMSD) and ligand pocket lDDT (lDDT-LP).
 
     All versus all scores are available through the lazily computed
     :attr:`score_matrix`. However, many things can go wrong... be it even
@@ -65,7 +66,7 @@ class LigandScorer:
 
     A common use case is to derive a one-to-one mapping between ligands in
     the model and the target for which :class:`LigandScorer` provides an
-    automated assignment procedure.
+    automated :attr:`assignment` procedure.
     By default, only exact matches between target and model ligands are
     considered. This is a problem when the target only contains a subset
     of the expected atoms (for instance if atoms are missing in an
@@ -85,7 +86,7 @@ class LigandScorer:
 
     :class:`LigandScorer` generally assumes that the
     :attr:`~ost.mol.ResidueHandle.is_ligand` property is properly set on all
-    the ligand atoms, and only ligand atoms. This is typically the case for
+    the ligand residues, and only ligand atoms. This is typically the case for
     entities loaded from mmCIF (tested with mmCIF files from the PDB and
     SWISS-MODEL). Legacy PDB files must contain `HET` headers (which is usually
     the case for files downloaded from the PDB but not elsewhere).
@@ -1225,7 +1226,7 @@ def ComputeSymmetries(model_ligand, target_ligand, substructure_match=False,
             symmetries.append((list(isomorphism.values()),
                                list(isomorphism.keys())))
         assert len(symmetries) > 0
-        LogDebug("Found %s isomorphic mappings (symmetries)" % len(symmetries))
+        LogVerbose("Found %s isomorphic mappings (symmetries)" % len(symmetries))
     elif gm.subgraph_is_isomorphic() and substructure_match:
         if not return_symmetries:
             return True
@@ -1241,14 +1242,14 @@ def ComputeSymmetries(model_ligand, target_ligand, substructure_match=False,
         # Assert that all the atoms in the target are part of the substructure
         assert len(symmetries[0][0]) == len(target_ligand.atoms)
         n_sym = len(symmetries)
-        LogDebug("Found %s subgraph isomorphisms (symmetries)" % n_sym)
+        LogVerbose("Found %s subgraph isomorphisms (symmetries)" % n_sym)
     elif gm.subgraph_is_isomorphic():
-        LogDebug("Found subgraph isomorphisms (symmetries), but"
+        LogVerbose("Found subgraph isomorphisms (symmetries), but"
                  " ignoring because substructure_match=False")
         raise NoIsomorphicSymmetryError("No symmetry between %s and %s" % (
             str(model_ligand), str(target_ligand)))
     else:
-        LogDebug("Found no isomorphic mappings (symmetries)")
+        LogVerbose("Found no isomorphic mappings (symmetries)")
         raise NoSymmetryError("No symmetry between %s and %s" % (
             str(model_ligand), str(target_ligand)))
 
diff --git a/modules/mol/alg/pymod/ligand_scoring_lddtpli.py b/modules/mol/alg/pymod/ligand_scoring_lddtpli.py
index be245c63f22353b344b1345210b8cb44c8d4a08d..467066471a480227480f98f4be4adae065242b22 100644
--- a/modules/mol/alg/pymod/ligand_scoring_lddtpli.py
+++ b/modules/mol/alg/pymod/ligand_scoring_lddtpli.py
@@ -26,20 +26,16 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
     and an lDDT-PLI score is computed. The best possible lDDT-PLI score is
     returned.
 
-    By default, classic lDDT is computed. That means, contacts within
-    *lddt_pli_radius* are identified in the target and checked if they're
-    conserved in the model. Added contacts are not penalized. That means if
-    the ligand is nicely placed in the correct pocket, but that pocket now
-    suddenly interacts with MORE residues in the model, you still get a high
-    score. You can penalize for these added contacts with the
-    *add_mdl_contacts* flag. This additionally considers contacts within
-    *lddt_pli_radius* in the model but only if the involved atoms can
-    be mapped to the target. This is a requirement to 1) extract the respective
-    reference distance from the target 2) avoid usage of contacts for which
-    we have no experimental evidence. One special case are
-    contacts from chains that are NOT mapped to the target binding site. It is
-    very well possible that we have experimental evidence for this chain though
-    its just too far away from the target binding site.
+    The lDDT-PLI score is a variant of lDDT with a custom inclusion radius
+    (`lddt_pli_radius`), no stereochemistry checks, and which penalizes
+    contacts added in the model within `lddt_pli_radius` by default
+    (can be changed with the `add_mdl_contacts` flag) but only if the involved
+    atoms can be mapped to the target. This is a requirement to
+    1) extract the respective reference distance from the target
+    2) avoid usage of contacts for which we have no experimental evidence.
+    One special case are contacts from chains that are not mapped to the target
+    binding site. It is very well possible that we have experimental evidence
+    for this chain though its just too far away from the target binding site.
     We therefore try to map these contacts to the chain in the target with
     equivalent sequence that is closest to the target binding site. If the
     respective atoms can be mapped there, the contact is considered not
@@ -47,7 +43,7 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
 
     Populates :attr:`LigandScorer.aux_data` with following :class:`dict` keys:
 
-    * lddt_pli: The score
+    * lddt_pli: The LDDT-PLI score
     * lddt_pli_n_contacts: Number of contacts considered in lDDT computation
     * target_ligand: The actual target ligand for which the score was computed
     * model_ligand: The actual model ligand for which the score was computed
@@ -102,7 +98,7 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
     def __init__(self, model, target, model_ligands=None, target_ligands=None,
                  resnum_alignments=False, rename_ligand_chain=False,
                  substructure_match=False, coverage_delta=0.2,
-                 max_symmetries=1e5, lddt_pli_radius=6.0,
+                 max_symmetries=1e4, lddt_pli_radius=6.0,
                  add_mdl_contacts=True,
                  lddt_pli_thresholds = [0.5, 1.0, 2.0, 4.0],
                  lddt_pli_binding_site_radius=None):
diff --git a/modules/mol/alg/pymod/ligand_scoring_scrmsd.py b/modules/mol/alg/pymod/ligand_scoring_scrmsd.py
index 9ae0546d6ac62e88237bdb6aba89ad5792683de5..15ef917b48eb4bbcd4c0d97fcf3900ceb8c88d40 100644
--- a/modules/mol/alg/pymod/ligand_scoring_scrmsd.py
+++ b/modules/mol/alg/pymod/ligand_scoring_scrmsd.py
@@ -8,14 +8,14 @@ from ost.mol.alg import ligand_scoring_base
 
 
 class SCRMSDScorer(ligand_scoring_base.LigandScorer):
-    """ :class:`LigandScorer` implementing symmetry corrected RMSD.
+    """ :class:`LigandScorer` implementing symmetry corrected RMSD (BiSyRMSD).
 
     :class:`SCRMSDScorer` computes a score for a specific pair of target/model
     ligands.
 
     The returned RMSD is based on a binding site superposition.
     The binding site of the target structure is defined as all residues with at
-    least one atom within *bs_radius* around the target ligand.
+    least one atom within `bs_radius` around the target ligand.
     It only contains protein and nucleic acid residues from chains that
     pass the criteria for the
     :class:`chain mapping <ost.mol.alg.chain_mapping>`. This means ignoring
@@ -24,28 +24,29 @@ class SCRMSDScorer(ligand_scoring_base.LigandScorer):
     The respective model binding site for superposition is identified by
     naively enumerating all possible mappings of model chains onto their
     chemically equivalent target counterparts from the target binding site.
-    The *binding_sites_topn* with respect to lDDT score are evaluated and 
+    The `binding_sites_topn` with respect to lDDT score are evaluated and 
     an RMSD is computed.
     You can either try to map ALL model chains onto the target binding site by
-    enabling *full_bs_search* or restrict the model chains for a specific
+    enabling `full_bs_search` or restrict the model chains for a specific
     target/model ligand pair to the chains with at least one atom within
     *model_bs_radius* around the model ligand. The latter can be significantly
     faster in case of large complexes.
     Symmetry correction is achieved by simply computing an RMSD value for
     each symmetry, i.e. atom-atom assignments of the ligand as given by
-    :class:`LigandScorer`. The lowest RMSD value is returned
+    :class:`LigandScorer`. The lowest RMSD value is returned.
 
     Populates :attr:`LigandScorer.aux_data` with following :class:`dict` keys:
 
-    * rmsd: The score
-    * lddt_lp: lDDT of the binding site used for superposition
+    * rmsd: The BiSyRMSD score
+    * lddt_lp: lDDT of the binding pocket used for superposition (lDDT-LP)
     * bs_ref_res: :class:`list` of binding site residues in target
     * bs_ref_res_mapped: :class:`list` of target binding site residues that
       are mapped to model
     * bs_mdl_res_mapped: :class:`list` of same length with respective model
       residues
-    * bb_rmsd: Backbone RMSD (CA, C3' for nucleotides) for mapped residues
-      given transform
+    * bb_rmsd: Backbone RMSD (CA, C3' for nucleotides; full backbone for 
+      binding sites with fewer than 3 residues) for mapped binding site
+      residues after superposition
     * target_ligand: The actual target ligand for which the score was computed
     * model_ligand: The actual model ligand for which the score was computed
     * chain_mapping: :class:`dict` with a chain mapping of chains involved in
diff --git a/singularity/Singularity b/singularity/Singularity
index 395c512ca3289341613932f39e24ae9aed0fc264..8ed1bb886168a62cb735a7528204265376e8b580 100644
--- a/singularity/Singularity
+++ b/singularity/Singularity
@@ -1,5 +1,5 @@
 BootStrap: docker
-From: registry.scicore.unibas.ch/schwede/openstructure:2.8.0
+From: registry.scicore.unibas.ch/schwede/openstructure:2.9.0
 %post
 ##############################################################################
 # POST