diff --git a/CHANGELOG.md b/CHANGELOG.md
index 413456e7c..825727203 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,7 @@
 # Changelog
 
 =======
+- 2026-07-13: Added option to select models/clusters in reverse order - Issue #1620
 - 2026-07-10: Corrected ranking in caprieval for reverse sorting - Issue #1621
 - 2026-07-09: Added `deeprank` scoring module using deeprank-gnn-esm - Issue #569
 - 2026-07-07: Re-add `gdock` as a sampling module
diff --git a/src/haddock/modules/analysis/seletop/__init__.py b/src/haddock/modules/analysis/seletop/__init__.py
index bc60ef91b..dcc934975 100644
--- a/src/haddock/modules/analysis/seletop/__init__.py
+++ b/src/haddock/modules/analysis/seletop/__init__.py
@@ -49,9 +49,16 @@ def _run(self) -> None:
         models_to_select: list[PDBFile] = [
             p for p in self.previous_io.output if p.file_type == Format.PDB
         ]
+        if self.params["sort_ascending"]:
+            reverse = False
+            self.log("Selecting models with the lowest score")
+        else:
+            reverse = True
+        if reverse:
+            self.log("Selecting models with the highest score")
 
         # sort the models based on their score
-        models_to_select.sort(key=lambda x: x.score)
+        models_to_select.sort(key=lambda x: x.score, reverse=reverse)
 
         if len(models_to_select) < self.params["select"]:
             self.log(
diff --git a/src/haddock/modules/analysis/seletop/defaults.yaml b/src/haddock/modules/analysis/seletop/defaults.yaml
index a96700476..8e2a5281f 100644
--- a/src/haddock/modules/analysis/seletop/defaults.yaml
+++ b/src/haddock/modules/analysis/seletop/defaults.yaml
@@ -8,3 +8,11 @@ select:
   long: Number of best-ranked models to select for further steps of the workflow. Typically only 200 models produced at the rigid-body docking stage are selected for further refinement.
   group: analysis
   explevel: easy
+sort_ascending:
+  default: true
+  type: boolean
+  title: Sort in ascending order
+  short: Sort in ascending order.
+  long: Sort in ascending order. “sort_ascending = false” should be used when a higher score corresponds to a better predicted model quality.  
+  group: analysis
+  explevel: easy
diff --git a/src/haddock/modules/analysis/seletopclusts/__init__.py b/src/haddock/modules/analysis/seletopclusts/__init__.py
index eca5d2a02..b9340d5e7 100644
--- a/src/haddock/modules/analysis/seletopclusts/__init__.py
+++ b/src/haddock/modules/analysis/seletopclusts/__init__.py
@@ -21,7 +21,7 @@
 from haddock.modules.analysis.seletopclusts.seletopclusts import (
     select_top_clusts_models,
     write_selected_models,
-    )
+)
 
 
 RECIPE_PATH = Path(__file__).resolve().parent
@@ -33,12 +33,14 @@ class HaddockModule(BaseHaddockModule):
 
     name = RECIPE_PATH.name
 
-    def __init__(self,
-                 order: int,
-                 path: Path,
-                 *ignore: Any,
-                 init_params: FilePath = DEFAULT_CONFIG,
-                 **everything: Any) -> None:
+    def __init__(
+        self,
+        order: int,
+        path: Path,
+        *ignore: Any,
+        init_params: FilePath = DEFAULT_CONFIG,
+        **everything: Any,
+    ) -> None:
         super().__init__(order, path, init_params)
 
     @classmethod
@@ -57,6 +59,15 @@ def _run(self) -> None:
             _msg = "top_clusters must be an integer."
             self.finish_with_error(_msg)
 
+        # The reverse option only applies to score-based ranking
+        reverse = not self.params["sort_ascending"]
+        if reverse and self.params["sortby"] != "size":
+            self.log("Selecting clusters with the highest score first")
+        elif not reverse and self.params["sortby"] != "size":
+            self.log("Selecting clusters with the lowest score first")
+        else:
+            self.log("Selecting clusters with the largest size first")
+            self.log("Selecting clusters with the lowest score first")
         # Retrieve list of previous models
         models_to_select = self.previous_io.retrieve_models()
 
@@ -65,16 +76,18 @@ def _run(self) -> None:
             _msg = (
                 "Impossible to obtain cluster information. Please consider "
                 "running a clustering method prior to this module."
-                )
+            )
             self.finish_with_error(_msg)
 
         # Make model selection
         selected_models, _notes = select_top_clusts_models(
             self.params["sortby"],
+            reverse,
             models_to_select,
             self.params["top_clusters"],
             self.params["top_models"],
-            )
+            self.params["cluster_score_threshold"],
+        )
         # Log notes
         for note in _notes:
             log.info(note)
@@ -84,7 +97,7 @@ def _run(self) -> None:
             "seletopclusts.txt",
             selected_models,
             self.path,
-            )
+        )
 
         # Make these new models the output of this module
         self.output_models = renamed_models
diff --git a/src/haddock/modules/analysis/seletopclusts/defaults.yaml b/src/haddock/modules/analysis/seletopclusts/defaults.yaml
index 6f6ece8b7..1126541ec 100644
--- a/src/haddock/modules/analysis/seletopclusts/defaults.yaml
+++ b/src/haddock/modules/analysis/seletopclusts/defaults.yaml
@@ -31,3 +31,21 @@ sortby:
   long: if the selection is done by 'score' the average score of the top (4) models of each cluster is used to define the cluster rank. When clustering by 'size', a bigger cluster size corresponds to a higher rank. By default, 'score' is selected.
   group: analysis
   explevel: easy
+sort_ascending:
+  default: true
+  type: boolean
+  title: Sort in ascending order
+  short: Sort in ascending order.
+  long: Sort in ascending order.  “sort_ascending = false” should be used when a higher score corresponds to a better predicted model quality and is only used if “sortby” is set to "score".
+  group: analysis
+  explevel: easy
+cluster_score_threshold:
+  default: 4
+  type: integer
+  min: 1
+  max: 99999
+  title: Number of models used to compute a cluster's average score
+  short: Number of best-scoring models per cluster used to compute its average score. By default 4.
+  long: When ranking clusters by 'score', the average score of this many best-scoring models of each cluster is used to define the cluster rank. By default 4.
+  group: analysis
+  explevel: easy
diff --git a/src/haddock/modules/analysis/seletopclusts/seletopclusts.py b/src/haddock/modules/analysis/seletopclusts/seletopclusts.py
index ffde8259e..fee479f65 100644
--- a/src/haddock/modules/analysis/seletopclusts/seletopclusts.py
+++ b/src/haddock/modules/analysis/seletopclusts/seletopclusts.py
@@ -4,15 +4,18 @@
 import math
 from pathlib import Path
 from haddock.core.typing import Union
+from haddock.libs.libclust import rank_clusters
 from haddock.libs.libontology import PDBFile
 
 
 def select_top_clusts_models(
-        sortby: str,
-        models_to_select: list[PDBFile],
-        top_clusters: int,
-        top_models: Union[int, float],
-        ) -> tuple[list[PDBFile], list[str]]:
+    sortby: str,
+    reverse: bool,
+    models_to_select: list[PDBFile],
+    top_clusters: int,
+    top_models: Union[int, float],
+    cluster_score_threshold: int = 4,
+) -> tuple[list[PDBFile], list[str]]:
     """Select best clusters based on structures scores.
 
     Parameters
@@ -25,6 +28,9 @@ def select_top_clusts_models(
         Number of best clusters to take into account.
     top_models : int
         Number of best models in each cluster to take into account.
+    cluster_score_threshold : int
+        Number of best-scoring models per cluster used to compute its
+        average score when ranking by `score`.
 
     Returns
     -------
@@ -40,7 +46,11 @@ def select_top_clusts_models(
     if sortby == "size":
         cluster_rankings = size_clust_order(by_clusters)
     else:
-        cluster_rankings = rank_clust_order(by_clusters)
+        cluster_rankings = rank_clust_order(
+            by_clusters,
+            reverse,
+            cluster_score_threshold,
+        )
 
     # Check if number of clusters >= set of rank
     if top_clusters >= len(cluster_rankings):
@@ -51,23 +61,22 @@ def select_top_clusts_models(
         # select top_clusters clusters
         cluster_rankings = cluster_rankings[:top_clusters]
         cluster_rankins_str = ",".join(map(str, cluster_rankings))
-        notes.append(
-            f"Selecting top {top_clusters} clusters: "
-            f"{cluster_rankins_str}"
-            )
+        notes.append(f"Selecting top {top_clusters} clusters: {cluster_rankins_str}")
 
     # Initiate set of selected models to export
     models_to_export: list[PDBFile] = []
-    # Loop over cluster ranks
-    for clt_rank in cluster_rankings:
+    # Loop over cluster ranks. `new_clt_rank` reflects the ordering computed
+    # above (by average score or size), so the exported model names follow
+    # that ranking rather than the original cluster rank.
+    for new_clt_rank, clt_key in enumerate(cluster_rankings, start=1):
         # Sort models by model rank
-        clt_mdls, note = sort_models(by_clusters[clt_rank])
+        clt_mdls, note = sort_models(by_clusters[clt_key])
         if note:
             notes.append(note)
 
         # Set new ranks to models
         for mdl_rank, pdb in enumerate(clt_mdls, start=1):
-            pdb.clt_rank = clt_rank
+            pdb.clt_rank = new_clt_rank
             pdb.clt_model_rank = mdl_rank
         # In case number of models is not set (nan.)
         if math.isnan(top_models):
@@ -81,9 +90,7 @@ def select_top_clusts_models(
     return models_to_export, notes
 
 
-def sort_models(
-        models: list[PDBFile]
-        ) -> tuple[list[PDBFile], Union[None, str]]:
+def sort_models(models: list[PDBFile]) -> tuple[list[PDBFile], Union[None, str]]:
     """Sort models based on their rank in cluster.
 
     Parameters
@@ -101,16 +108,18 @@ def sort_models(
         sorted_mdls = sorted(
             models,
             key=lambda k: k.clt_model_rank,
-            )
+        )
     except TypeError:
-        note = 'model rank unavailable, falling back to input order'
+        note = "model rank unavailable, falling back to input order"
         sorted_mdls = models
     return sorted_mdls, note
-    
+
 
 def rank_clust_order(
-        by_clusters: dict[int, list[PDBFile]],
-        ) -> list[int]:
+    by_clusters: dict[int, list[PDBFile]],
+    reverse: bool = False,
+    cluster_score_threshold: int = 4,
+) -> list[int]:
     """Select best clusters based on structures scores.
 
     Parameters
@@ -121,6 +130,11 @@ def rank_clust_order(
         Number of best clusters to take into account.
     top_models : int
         Number of best models in each cluster to take into account.
+    reverse : boolean
+        Reverse the ranking (from high to low average score)
+    cluster_score_threshold : int
+        Number of best-scoring models per cluster used to compute its
+        average score.
 
     Returns
     -------
@@ -129,14 +143,20 @@ def rank_clust_order(
     notes : list[str]
         List of notes to be printed.
     """
-    # Generate set of all cluster rank available
-    cluster_rankings = sorted(by_clusters)
+    # Rank clusters by their average score (best/lowest score first)
+    _score_dic, sorted_score_dic = rank_clusters(
+        by_clusters,
+        cluster_score_threshold,
+    )
+    cluster_rankings = [clt_key for clt_key, _score in sorted_score_dic]
+    if reverse:
+        cluster_rankings.reverse()
     return cluster_rankings
 
 
 def size_clust_order(
-        by_clusters: dict[int, list[PDBFile]],
-        ) -> list[int]:
+    by_clusters: dict[int, list[PDBFile]],
+) -> list[int]:
     """Select best clusters based on structures scores.
 
     Parameters
@@ -160,7 +180,7 @@ def size_clust_order(
         by_clusters,
         key=lambda k: len(by_clusters[k]),
         reverse=True,
-        )
+    )
     return cluster_rankings
 
 
@@ -179,9 +199,8 @@ def map_clusters_models(models: list[PDBFile]) -> dict[int, list[PDBFile]]:
     """
     # Preset dictionary keys
     by_clusters: dict[int, list[PDBFile]] = {
-        clrank: []
-        for clrank in list(set([pdb.clt_rank for pdb in models]))
-        }
+        clrank: [] for clrank in list(set([pdb.clt_rank for pdb in models]))
+    }
     # Loop over models
     for pdb in models:
         # Add model to cluster
@@ -190,10 +209,10 @@ def map_clusters_models(models: list[PDBFile]) -> dict[int, list[PDBFile]]:
 
 
 def write_selected_models(
-        output_path: Union[str, Path],
-        models: list[PDBFile],
-        module_path: Union[str, Path],
-        ) -> list[PDBFile]:
+    output_path: Union[str, Path],
+    models: list[PDBFile],
+    module_path: Union[str, Path],
+) -> list[PDBFile]:
     """Dump selected models and new names in a file.
 
     Parameters
@@ -211,26 +230,20 @@ def write_selected_models(
         Updated list of selected models.
     """
     # dump the models to disk and change their attributes
-    with open(output_path, 'w') as fh:
+    with open(output_path, "w") as fh:
         fh.write("rel_path\tori_name\tcluster_name\tmd5" + os.linesep)
         for model in models:
-            name = (
-                f"cluster_{model.clt_rank}_model"
-                f"_{model.clt_model_rank}.pdb"
-                )
+            name = f"cluster_{model.clt_rank}_model_{model.clt_model_rank}.pdb"
             # writing name
             fh.write(
-                f"{model.rel_path}\t"
-                f"{model.ori_name}\t"
-                f"{name}\t"
-                f"{model.md5}" + os.linesep
-                )
+                f"{model.rel_path}\t{model.ori_name}\t{name}\t{model.md5}" + os.linesep
+            )
             # changing attributes
             name_path = Path(name)
             name_path.write_text(model.rel_path.read_text())
             model.ori_name = model.file_name
             model.file_name = name
             model.full_name = name
-            model.rel_path = Path('..', Path(module_path).name, name)
+            model.rel_path = Path("..", Path(module_path).name, name)
             model.path = str(Path(".").resolve())
     return models
diff --git a/tests/test_module_seletop.py b/tests/test_module_seletop.py
new file mode 100644
index 000000000..26fd68b1d
--- /dev/null
+++ b/tests/test_module_seletop.py
@@ -0,0 +1,93 @@
+"""Test related to the [seletop] module."""
+
+import shutil
+import tempfile
+from pathlib import Path
+
+import pytest
+
+from haddock.libs.libontology import PDBFile
+from haddock.modules.analysis.seletop import DEFAULT_CONFIG
+from haddock.modules.analysis.seletop import HaddockModule as SeletopModule
+
+from . import golden_data
+
+
+class MockPreviousIO:
+    """Mocking class to retrieve models."""
+
+    def __init__(self, models):
+        self.models = models
+        self.output = models
+
+    def retrieve_models(self, individualize: bool = False):
+        """Provide a set of models."""
+        return self.models
+
+
+@pytest.fixture(name="scored_models")
+def fixture_scored_models() -> list[PDBFile]:
+    """Set of PDBfiles to be selected."""
+    models: list[PDBFile] = []
+    # Create 5 models with 5 different scores
+    for i in range(5, 0, -1):  # 5, 4, 3, 2, 1
+        with tempfile.NamedTemporaryFile(delete=False, suffix=".pdb") as temp_f:
+            dst = temp_f.name
+            src = Path(golden_data, "protprot_complex_1.pdb")
+            shutil.copy(src, dst)
+            pdbfile = PDBFile(file_name=dst, path=Path(dst).parent)
+            pdbfile.score = i
+            models.append(pdbfile)
+    return models
+
+
+@pytest.fixture(name="seletop")
+def fixture_seletop(monkeypatch):
+    """Test module __init__()."""
+    with tempfile.TemporaryDirectory() as tmpdir:
+        monkeypatch.chdir(tmpdir)
+        yield SeletopModule(
+            order=1,
+            path=Path("."),
+            initial_params=DEFAULT_CONFIG,
+        )
+
+
+def test_confirm_installation(seletop):
+    """Test confirm install."""
+    assert seletop.confirm_installation() is None
+
+
+def test_seletop_default_sort(seletop, mocker, scored_models):
+    """By default the models with the lowest score are selected first."""
+    seletop.params["select"] = 3
+    seletop.previous_io = MockPreviousIO(scored_models)
+    mocker.patch(
+        "haddock.modules.analysis.seletop.HaddockModule.previous_path",
+        return_value="0_rigidbody",
+    )
+    mocker.patch(
+        "haddock.modules.BaseHaddockModule.export_io_models",
+        return_value=None,
+    )
+    seletop.run()
+    # Lowest scores selected: 1, 2, 3
+    assert [m.score for m in seletop.output_models] == [1, 2, 3]
+
+
+def test_seletop_reversed_sort(seletop, mocker, scored_models):
+    """With sort_ascending disabled the highest scores are selected first."""
+    seletop.params["select"] = 3
+    seletop.params["sort_ascending"] = False
+    seletop.previous_io = MockPreviousIO(scored_models)
+    mocker.patch(
+        "haddock.modules.analysis.seletop.HaddockModule.previous_path",
+        return_value="0_rigidbody",
+    )
+    mocker.patch(
+        "haddock.modules.BaseHaddockModule.export_io_models",
+        return_value=None,
+    )
+    seletop.run()
+    # Highest scores selected: 5, 4, 3
+    assert [m.score for m in seletop.output_models] == [5, 4, 3]
diff --git a/tests/test_module_seletopclusts.py b/tests/test_module_seletopclusts.py
index 0cc75f290..73ab38889 100644
--- a/tests/test_module_seletopclusts.py
+++ b/tests/test_module_seletopclusts.py
@@ -4,14 +4,14 @@
 import os
 import shutil
 import tempfile
+from copy import deepcopy
 from pathlib import Path
 
 import pytest
 
 from haddock.libs.libontology import PDBFile
 from haddock.modules.analysis.seletopclusts import DEFAULT_CONFIG
-from haddock.modules.analysis.seletopclusts import \
-    HaddockModule as SeleTopClustModule
+from haddock.modules.analysis.seletopclusts import HaddockModule as SeleTopClustModule
 from haddock.modules.analysis.seletopclusts.seletopclusts import (
     map_clusters_models,
     rank_clust_order,
@@ -19,7 +19,7 @@
     size_clust_order,
     sort_models,
     write_selected_models,
-    )
+)
 
 from . import golden_data
 
@@ -57,6 +57,8 @@ def fixture_clustered_models() -> list[PDBFile]:
                 # Add attributes
                 pdbfile.clt_rank = clt_rank
                 pdbfile.clt_model_rank = mdl_rank
+                # Lower cluster rank -> better (lower) average score
+                pdbfile.score = float(clt_rank * 100 + mdl_rank)
                 models.append(pdbfile)
     return models
 
@@ -226,6 +228,7 @@ def test_select_top_ranked_clusts_models(clustered_models):
     # Test 1 model, 1 cluster
     selected_models, _notes = select_top_clusts_models(
         "score",
+        False,
         clustered_models,
         1,
         1,
@@ -237,6 +240,7 @@ def test_select_top_ranked_clusts_models(clustered_models):
     # Test max models
     selected_models, _notes = select_top_clusts_models(
         "score",
+        False,
         clustered_models,
         2,
         8,
@@ -253,6 +257,7 @@ def test_select_top_ranked_clusts_models(clustered_models):
     # Test max models and max clusters
     selected_models, _notes = select_top_clusts_models(
         "score",
+        False,
         clustered_models,
         NB_CLUSTERS + 2,
         NB_CLUSTERS + MDL_RANK_INCREMENT + 1,
@@ -268,6 +273,7 @@ def test_select_top_ranked_clusts_models(clustered_models):
     # Test unspecified nb models, max clusters
     selected_models, _notes = select_top_clusts_models(
         "score",
+        False,
         clustered_models,
         NB_CLUSTERS + 2,
         math.nan,
@@ -281,23 +287,47 @@ def test_select_top_ranked_clusts_models(clustered_models):
             ind += 1
 
 
+def test_select_top_ranked_clusts_models_reverse(clustered_models):
+    """Naming follows the reversed (high-to-low score) cluster ranking."""
+    selected_models, _notes = select_top_clusts_models(
+        "score",
+        True,
+        clustered_models,
+        NB_CLUSTERS + 2,
+        math.nan,
+    )
+    assert len(selected_models) == len(clustered_models)
+    ind = 0
+    # Reverse score order selects worst-scoring cluster first (orig ranks
+    # 3, 2, 1), renumbered to new ranks 1, 2, 3.
+    for new_rank, orig_rank in enumerate(range(NB_CLUSTERS, 0, -1), start=1):
+        nb_models = orig_rank + MDL_RANK_INCREMENT - 1
+        for mdl_rank in range(1, nb_models + 1):
+            assert selected_models[ind].clt_rank == new_rank
+            assert selected_models[ind].clt_model_rank == mdl_rank
+            ind += 1
+
+
 def test_select_top_sized_clusts_models(clustered_models):
     """Test select_top_clusts_models behavior with size ordered clusters."""
     # Test 1 model, 1 cluster
     selected_models, _notes = select_top_clusts_models(
         "size",
-        clustered_models,
+        False,
+        deepcopy(clustered_models),
         1,
         1,
     )
     assert len(selected_models) == 1
-    assert selected_models[0].clt_rank == 3
+    # Biggest cluster (orig rank 3) is renamed to rank 1
+    assert selected_models[0].clt_rank == 1
     assert selected_models[0].clt_model_rank == 1
 
     # Test max models
     selected_models, _notes = select_top_clusts_models(
         "size",
-        clustered_models,
+        False,
+        deepcopy(clustered_models),
         2,
         8,
     )
@@ -313,30 +343,36 @@ def test_select_top_sized_clusts_models(clustered_models):
     # Test max models and max clusters
     selected_models, _notes = select_top_clusts_models(
         "size",
-        clustered_models,
+        False,
+        deepcopy(clustered_models),
         NB_CLUSTERS + 2,
         NB_CLUSTERS + MDL_RANK_INCREMENT + 1,
     )
     assert len(selected_models) == len(clustered_models)
     ind = 0
-    for clt_rank in range(NB_CLUSTERS, 0, -1):
-        for mdl_rank in range(1, clt_rank + MDL_RANK_INCREMENT):
-            assert selected_models[ind].clt_rank == clt_rank
+    # Size order selects biggest cluster first (orig ranks 3, 2, 1),
+    # renumbered to new ranks 1, 2, 3.
+    for new_rank, orig_rank in enumerate(range(NB_CLUSTERS, 0, -1), start=1):
+        nb_models = orig_rank + MDL_RANK_INCREMENT - 1
+        for mdl_rank in range(1, nb_models + 1):
+            assert selected_models[ind].clt_rank == new_rank
             assert selected_models[ind].clt_model_rank == mdl_rank
             ind += 1
 
     # Test unspecified nb models, max clusters
     selected_models, _notes = select_top_clusts_models(
         "size",
-        clustered_models,
+        False,
+        deepcopy(clustered_models),
         NB_CLUSTERS + 2,
         math.nan,
     )
     assert len(selected_models) == len(clustered_models)
     ind = 0
-    for clt_rank in range(NB_CLUSTERS, 0, -1):
-        for mdl_rank in range(1, clt_rank + MDL_RANK_INCREMENT):
-            assert selected_models[ind].clt_rank == clt_rank
+    for new_rank, orig_rank in enumerate(range(NB_CLUSTERS, 0, -1), start=1):
+        nb_models = orig_rank + MDL_RANK_INCREMENT - 1
+        for mdl_rank in range(1, nb_models + 1):
+            assert selected_models[ind].clt_rank == new_rank
             assert selected_models[ind].clt_model_rank == mdl_rank
             ind += 1