diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6f57cff01..7e22505b6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,6 @@
 # Changelog
 
+- 2026-08-10: Fixed nucleobase selection for DNA default caprieval - Issue #1658
 - 2026-08-05: Harmonised running mode of test config files - Issue #1655
 - 2026-08-10: Added missing improper for HYP - Issue #1662
 - 2026-08-05: Removed identical duplicate `ANGLe` statements from `protein-allhdg5-4.param` and `protein-CG-Martini-2-2.param` - Issue #1589
diff --git a/src/haddock/libs/libalign.py b/src/haddock/libs/libalign.py
index 0af3740c9..61ebd7a40 100644
--- a/src/haddock/libs/libalign.py
+++ b/src/haddock/libs/libalign.py
@@ -42,11 +42,97 @@
     slc_y,
     slc_z,
     split_by_chain,
-    )
+)
 
 
 RES_TO_BE_IGNORED = ["SHA", "WAT"]
 
+ResCode = Literal[
+    "C",
+    "D",
+    "S",
+    "Q",
+    "K",
+    "I",
+    "P",
+    "T",
+    "F",
+    "N",
+    "G",
+    "H",
+    "L",
+    "R",
+    "W",
+    "A",
+    "V",
+    "E",
+    "Y",
+    "M",
+    "A",
+    "G",
+    "C",
+    "T",
+    "X",
+]
+"""
+The single letter code of a residue.
+
+Unrecognized residues' code is `X`.
+"""
+
+RES_3L1L = dict(
+    [
+        ("CYS", "C"),
+        ("ASP", "D"),
+        ("SER", "S"),
+        ("GLN", "Q"),
+        ("LYS", "K"),
+        ("ILE", "I"),
+        ("PRO", "P"),
+        ("THR", "T"),
+        ("PHE", "F"),
+        ("ASN", "N"),
+        ("GLY", "G"),
+        ("HIS", "H"),
+        ("LEU", "L"),
+        ("ARG", "R"),
+        ("TRP", "W"),
+        ("ALA", "A"),
+        ("VAL", "V"),
+        ("GLU", "E"),
+        ("TYR", "Y"),
+        ("MET", "M"),
+        ("DA", "A"),
+        ("DG", "G"),
+        ("DC", "C"),
+        ("DT", "T"),
+        # 9/8/2023: adding non-standard amino-acids (src/haddock/cns/toppar/protein-allhdg5-4.top) # noqa: E501
+        ("ALY", "K"),
+        ("ASH", "D"),
+        ("CFE", "C"),
+        ("CSP", "C"),
+        ("CYC", "C"),
+        ("CYF", "C"),
+        ("CYM", "C"),
+        ("DDZ", "A"),
+        ("GLH", "E"),
+        ("HLY", "P"),
+        ("HY3", "P"),
+        ("HYP", "P"),
+        ("M3L", "K"),
+        ("MLY", "K"),
+        ("MLZ", "K"),
+        ("MSE", "M"),
+        ("NEP", "H"),
+        ("PNS", "S"),
+        ("PTR", "Y"),
+        ("SEP", "S"),
+        ("TOP", "T"),
+        ("TYP", "Y"),
+        ("TYS", "Y"),
+    ]
+)
+
 PROT_RES = [
     "ALA",
     "ARG",
@@ -75,9 +161,9 @@
 PROT_ATOMS = ["C", "N", "CA", "O"]
 PROT_ATOMS_MARTINI2 = ["BB"]
 PROT_ATOMS_all = {
-    "aa" : PROT_ATOMS,
-    "martini2" : PROT_ATOMS_MARTINI2,
-    "martini3" : PROT_ATOMS_MARTINI2 # The backbone particle is identical between MARTINI 2 and 3
+    "aa": PROT_ATOMS,
+    "martini2": PROT_ATOMS_MARTINI2,
+    "martini3": PROT_ATOMS_MARTINI2,  # The backbone particle is identical between MARTINI 2 and 3
 }
 
 # Side chains
@@ -133,78 +219,78 @@
 }
 
 PROT_SIDE_CHAINS_DICT_MARTINI2 = {
-    "GLY": ['BB'],
-    "ALA": ['BB'],
-    "CYS": ['BB', 'SC1'],
-    "VAL": ['BB', 'SC1'],
-    "LEU": ['BB', 'SC1'],
-    "ILE": ['BB', 'SC1'],
-    "MET": ['BB', 'SC1'],
-    "PRO": ['BB', 'SC1'],
-    "ASN": ['BB', 'SC1'],
-    "GLN": ['BB', 'SC1'],
-    "ASP": ['BB', 'SC1'],
-    "ASP0": ['BB', 'SC1'],
-    "GLU": ['BB', 'SC1'],
-    "GLU0": ['BB', 'SC1'],
-    "THR": ['BB', 'SC1'],
-    "SER": ['BB', 'SC1'],
-    "LYS": ['BB', 'SC1', 'SC2'],
-    "LSN": ['BB', 'SC1', 'SC2'],
-    "ARG": ['BB', 'SC1', 'SC2'],
-    "ARG0": ['BB', 'SC1', 'SC2'],
-    "HIS": ['BB', 'SC1', 'SC2', 'SC3'],
-    "HSE": ['BB', 'SC1', 'SC2', 'SC3'],
-    "HSD": ['BB', 'SC1', 'SC2', 'SC3'],
-    "HSP": ['BB', 'SC1', 'SC2', 'SC3'],
-    "PHE": ['BB', 'SC1', 'SC2', 'SC3'],
-    "TYR": ['BB', 'SC1', 'SC2', 'SC3'],
-    "TRP": ['BB', 'SC1', 'SC2', 'SC3', 'SC4'],
+    "GLY": ["BB"],
+    "ALA": ["BB"],
+    "CYS": ["BB", "SC1"],
+    "VAL": ["BB", "SC1"],
+    "LEU": ["BB", "SC1"],
+    "ILE": ["BB", "SC1"],
+    "MET": ["BB", "SC1"],
+    "PRO": ["BB", "SC1"],
+    "ASN": ["BB", "SC1"],
+    "GLN": ["BB", "SC1"],
+    "ASP": ["BB", "SC1"],
+    "ASP0": ["BB", "SC1"],
+    "GLU": ["BB", "SC1"],
+    "GLU0": ["BB", "SC1"],
+    "THR": ["BB", "SC1"],
+    "SER": ["BB", "SC1"],
+    "LYS": ["BB", "SC1", "SC2"],
+    "LSN": ["BB", "SC1", "SC2"],
+    "ARG": ["BB", "SC1", "SC2"],
+    "ARG0": ["BB", "SC1", "SC2"],
+    "HIS": ["BB", "SC1", "SC2", "SC3"],
+    "HSE": ["BB", "SC1", "SC2", "SC3"],
+    "HSD": ["BB", "SC1", "SC2", "SC3"],
+    "HSP": ["BB", "SC1", "SC2", "SC3"],
+    "PHE": ["BB", "SC1", "SC2", "SC3"],
+    "TYR": ["BB", "SC1", "SC2", "SC3"],
+    "TRP": ["BB", "SC1", "SC2", "SC3", "SC4"],
 }
 
 PROT_SIDE_CHAINS_DICT_MARTINI3 = {
-    "GLY": ['BB'],
-    "ALA": ['BB', 'SC1'],
-    "CYS": ['BB', 'SC1'],
-    "VAL": ['BB', 'SC1'],
-    "LEU": ['BB', 'SC1'],
-    "ILE": ['BB', 'SC1'],
-    "MET": ['BB', 'SC1'],
-    "PRO": ['BB', 'SC1'],
-    "HYP": ['BB', 'SC1'],
-    "ASN": ['BB', 'SC1'],
-    "GLN": ['BB', 'SC1'],
-    "ASP": ['BB', 'SC1'],
-    "ASPP": ['BB', 'SC1'],
-    "ASH": ['BB', 'SC1'],
-    "GLU": ['BB', 'SC1'],
-    "GLUP": ['BB', 'SC1'],
-    "GLH": ['BB', 'SC1'],
-    "LYS": ['BB', 'SC1', 'SC2'],
-    "LSN": ['BB', 'SC1', 'SC2'],
-    "LYN": ['BB', 'SC1', 'SC2'],
-    "THR": ['BB', 'SC1'],
-    "SER": ['BB', 'SC1'],
-    "ARG": ['BB', 'SC1', 'SC2'],
-    "HIS": ['BB', 'SC1', 'SC2', 'SC3'],
-    "HIE": ['BB', 'SC1', 'SC2', 'SC3'],
-    "HSE": ['BB', 'SC1', 'SC2', 'SC3'],
-    "HSD": ['BB', 'SC1', 'SC2', 'SC3'],
-    "HID": ['BB', 'SC1', 'SC2', 'SC3'],
-    "HSP": ['BB', 'SC1', 'SC2', 'SC3'],
-    "HIP": ['BB', 'SC1', 'SC2', 'SC3'],
-    "PHE": ['BB', 'SC1', 'SC2', 'SC3'],
-    "TYR": ['BB', 'SC1', 'SC2', 'SC3', 'SC4'],
-    "TRP": ['BB', 'SC1', 'SC2', 'SC3', 'SC4', 'SC5'],    
+    "GLY": ["BB"],
+    "ALA": ["BB", "SC1"],
+    "CYS": ["BB", "SC1"],
+    "VAL": ["BB", "SC1"],
+    "LEU": ["BB", "SC1"],
+    "ILE": ["BB", "SC1"],
+    "MET": ["BB", "SC1"],
+    "PRO": ["BB", "SC1"],
+    "HYP": ["BB", "SC1"],
+    "ASN": ["BB", "SC1"],
+    "GLN": ["BB", "SC1"],
+    "ASP": ["BB", "SC1"],
+    "ASPP": ["BB", "SC1"],
+    "ASH": ["BB", "SC1"],
+    "GLU": ["BB", "SC1"],
+    "GLUP": ["BB", "SC1"],
+    "GLH": ["BB", "SC1"],
+    "LYS": ["BB", "SC1", "SC2"],
+    "LSN": ["BB", "SC1", "SC2"],
+    "LYN": ["BB", "SC1", "SC2"],
+    "THR": ["BB", "SC1"],
+    "SER": ["BB", "SC1"],
+    "ARG": ["BB", "SC1", "SC2"],
+    "HIS": ["BB", "SC1", "SC2", "SC3"],
+    "HIE": ["BB", "SC1", "SC2", "SC3"],
+    "HSE": ["BB", "SC1", "SC2", "SC3"],
+    "HSD": ["BB", "SC1", "SC2", "SC3"],
+    "HID": ["BB", "SC1", "SC2", "SC3"],
+    "HSP": ["BB", "SC1", "SC2", "SC3"],
+    "HIP": ["BB", "SC1", "SC2", "SC3"],
+    "PHE": ["BB", "SC1", "SC2", "SC3"],
+    "TYR": ["BB", "SC1", "SC2", "SC3", "SC4"],
+    "TRP": ["BB", "SC1", "SC2", "SC3", "SC4", "SC5"],
 }
 
 PROT_SIDE_CHAINS_DICT_all = {
-    "aa" : PROT_SIDE_CHAINS_DICT,
-    "martini2" : PROT_SIDE_CHAINS_DICT_MARTINI2,
-    "martini3" : PROT_SIDE_CHAINS_DICT_MARTINI3
+    "aa": PROT_SIDE_CHAINS_DICT,
+    "martini2": PROT_SIDE_CHAINS_DICT_MARTINI2,
+    "martini3": PROT_SIDE_CHAINS_DICT_MARTINI3,
 }
 
-# Bases
+# For DNA, we use nucleo-bases as default atom selection
 DNA_ATOMS = [
     "C5",
     "N9",
@@ -223,19 +309,20 @@
     "C4",
     "O6",
 ]
-DNA_ATOMS_MARTINI2 = ["BB1", "BB2", "BB3"]
-
+DNA_ATOMS_MARTINI2 = ["SC1", "SC2", "SC3", "SC4"]
+DNA_ATOMS_MARTINI3 = ["SC1", "SC2", "SC3", "SC4"]
 DNA_ATOMS_all = {
-    "aa" : DNA_ATOMS,
-    "martini2" : DNA_ATOMS_MARTINI2,
-    "martini3" : DNA_ATOMS_MARTINI2
+    "aa": DNA_ATOMS,
+    "martini2": DNA_ATOMS_MARTINI2,
+    "martini3": DNA_ATOMS_MARTINI3,
 }
 
+# Definition of all atoms in DNA
 DNA_FULL_DICT = {
     "DA": [
         "P",
-        "O1P",
-        "O2P",
+        "OP1",
+        "OP2",
         "O5'",
         "C5'",
         "C4'",
@@ -263,8 +350,8 @@
     ],
     "DG": [
         "P",
-        "O1P",
-        "O2P",
+        "OP1",
+        "OP2",
         "O5'",
         "C5'",
         "C4'",
@@ -292,8 +379,8 @@
     ],
     "DC": [
         "P",
-        "O1P",
-        "O2P",
+        "OP1",
+        "OP2",
         "O5'",
         "C5'",
         "C4'",
@@ -321,8 +408,8 @@
     ],
     "DT": [
         "P",
-        "O1P",
-        "O2P",
+        "OP1",
+        "OP2",
         "O5'",
         "C5'",
         "C4'",
@@ -351,28 +438,31 @@
 }
 
 DNA_FULL_DICT_MARTINI2 = {
-    "DA" : ["BB1", "BB2", "BB3", "SC1", "SC2", "SC3", "SC4"],
-    "DC" : ["BB1", "BB2", "BB3", "SC1", "SC2", "SC3"],
-    "DG" : ["BB1", "BB2", "BB3", "SC1", "SC2", "SC3", "SC4"],
-    "DT" : ["BB1", "BB2", "BB3", "SC1", "SC2", "SC3"]
+    "DA": ["BB1", "BB2", "BB3", "SC1", "SC2", "SC3", "SC4"],
+    "DC": ["BB1", "BB2", "BB3", "SC1", "SC2", "SC3"],
+    "DG": ["BB1", "BB2", "BB3", "SC1", "SC2", "SC3", "SC4"],
+    "DT": ["BB1", "BB2", "BB3", "SC1", "SC2", "SC3"],
 }
-
+DNA_FULL_DICT_MARTINI3 = DNA_FULL_DICT_MARTINI2
 DNA_FULL_DICT_all = {
-    "aa" : DNA_FULL_DICT,
-    "martini2" : DNA_FULL_DICT_MARTINI2,
-    "martini3" : DNA_FULL_DICT_MARTINI2
+    "aa": DNA_FULL_DICT,
+    "martini2": DNA_FULL_DICT_MARTINI2,
+    "martini3": DNA_FULL_DICT_MARTINI3,
 }
 
+# RNA definitions
+# Backbone
 RNA_RES = ["A", "G", "C", "U"]
 RNA_ATOMS = ["P", "O5'", "C5'", "C4'", "C3'", "O3'"]
 RNA_ATOMS_MARTINI2 = ["BB1", "BB2", "BB3"]
-
+RNA_ATOMS_MARTINI3 = RNA_ATOMS_MARTINI2
 RNA_ATOMS_all = {
-    "aa" : RNA_ATOMS,
-    "martini2" : RNA_ATOMS_MARTINI2,
-    "martini3" : RNA_ATOMS_MARTINI2
+    "aa": RNA_ATOMS,
+    "martini2": RNA_ATOMS_MARTINI2,
+    "martini3": RNA_ATOMS_MARTINI3,
 }
 
+# All atoms in RNA
 RNA_FULL_DICT = {
     "A": [
         "P",
@@ -436,13 +526,11 @@
         "C2'",
         "O2'",
         "C1'",
-        "N9",
         "C5",
-        "C6",
         "O6",
         "N1",
         "C2",
-        "N2",
+        "O2",
         "N3",
         "C4",
         "N4",
@@ -472,16 +560,16 @@
 }
 
 RNA_FULL_DICT_MARTINI2 = {
-    "A" : ["BB1", "BB2", "BB3", "SC1", "SC2", "SC3", "SC4"],
-    "C" : ["BB1", "BB2", "BB3", "SC1", "SC2", "SC3"],
-    "G" : ["BB1", "BB2", "BB3", "SC1", "SC2", "SC3", "SC4"],
-    "U" : ["BB1", "BB2", "BB3", "SC1", "SC2", "SC3"]
+    "A": ["BB1", "BB2", "BB3", "SC1", "SC2", "SC3", "SC4"],
+    "C": ["BB1", "BB2", "BB3", "SC1", "SC2", "SC3"],
+    "G": ["BB1", "BB2", "BB3", "SC1", "SC2", "SC3", "SC4"],
+    "U": ["BB1", "BB2", "BB3", "SC1", "SC2", "SC3"],
 }
-
+RNA_FULL_DICT_MARTINI3 = RNA_FULL_DICT_MARTINI2
 RNA_FULL_DICT_all = {
-    "aa" : RNA_FULL_DICT,
-    "martini2" : RNA_FULL_DICT_MARTINI2,
-    "martini3" : RNA_FULL_DICT_MARTINI2
+    "aa": RNA_FULL_DICT,
+    "martini2": RNA_FULL_DICT_MARTINI2,
+    "martini3": RNA_FULL_DICT_MARTINI3,
 }
 
 
@@ -583,7 +671,7 @@ def load_coords(
 
     add_resname : bool
         use the residue name in the identifier
-    
+
     keep_hetatm : bool
         Should HETATM lines be considered ?
 
@@ -600,9 +688,12 @@ def load_coords(
     idx: int = 0
     # Set types of coordinates lines to extract
     if keep_hetatm:
-        coordinates_line_to_extract = ("ATOM", "HETATM", )
+        coordinates_line_to_extract = (
+            "ATOM",
+            "HETATM",
+        )
     else:
-        coordinates_line_to_extract = ("ATOM", )
+        coordinates_line_to_extract = ("ATOM",)
     # Check filetype
     if isinstance(pdb_f, PDBFile):
         pdb_f = pdb_f.rel_path
@@ -700,9 +791,7 @@ def load_coords(
     return coord_dic, chain_ranges
 
 
-def get_atoms(pdb: PDBPath, 
-              full: bool = False,
-              ff: str = "aa") -> AtomsDict:
+def get_atoms(pdb: PDBPath, full: bool = False, ff: str = "aa") -> AtomsDict:
     """Identify what is the molecule type of each PDB.
 
     Parameters
@@ -767,41 +856,10 @@ def get_atoms(pdb: PDBPath,
     return atom_dic
 
 
-ResCode = Literal[
-    "C",
-    "D",
-    "S",
-    "Q",
-    "K",
-    "I",
-    "P",
-    "T",
-    "F",
-    "N",
-    "G",
-    "H",
-    "L",
-    "R",
-    "W",
-    "A",
-    "V",
-    "E",
-    "Y",
-    "M",
-    "A",
-    "G",
-    "C",
-    "T",
-    "X",
-]
-"""
-The single letter code of a residue.
-
-Unrecognized residues' code is `X`.
-"""
-
-
-def pdb2fastadic(pdb_f: PDBPath, keep_hetatm: bool = False) -> dict[str, dict[int, str]]:
+def pdb2fastadic(
+    pdb_f: PDBPath,
+    keep_hetatm: bool = False,
+) -> dict[str, dict[int, str]]:
     """
     Write the sequence as a fasta.
 
@@ -817,68 +875,22 @@ def pdb2fastadic(pdb_f: PDBPath, keep_hetatm: bool = False) -> dict[str, dict[in
     seq_dic : dict
         dict of fasta sequences (one per chain)
     """
-    res_codes = dict(
-        [
-            ("CYS", "C"),
-            ("ASP", "D"),
-            ("SER", "S"),
-            ("GLN", "Q"),
-            ("LYS", "K"),
-            ("ILE", "I"),
-            ("PRO", "P"),
-            ("THR", "T"),
-            ("PHE", "F"),
-            ("ASN", "N"),
-            ("GLY", "G"),
-            ("HIS", "H"),
-            ("LEU", "L"),
-            ("ARG", "R"),
-            ("TRP", "W"),
-            ("ALA", "A"),
-            ("VAL", "V"),
-            ("GLU", "E"),
-            ("TYR", "Y"),
-            ("MET", "M"),
-            ("DA", "A"),
-            ("DG", "G"),
-            ("DC", "C"),
-            (
-                "DT",
-                "T",
-            ),  # 9/8/2023: adding non-standard amino-acids (src/haddock/cns/toppar/protein-allhdg5-4.top) # noqa: E501
-            ("ALY", "K"),
-            ("ASH", "D"),
-            ("CFE", "C"),
-            ("CSP", "C"),
-            ("CYC", "C"),
-            ("CYF", "C"),
-            ("CYM", "C"),
-            ("DDZ", "A"),
-            ("GLH", "E"),
-            ("HLY", "P"),
-            ("HY3", "P"),
-            ("HYP", "P"),
-            ("M3L", "K"),
-            ("MLY", "K"),
-            ("MLZ", "K"),
-            ("MSE", "M"),
-            ("NEP", "H"),
-            ("PNS", "S"),
-            ("PTR", "Y"),
-            ("SEP", "S"),
-            ("TOP", "T"),
-            ("TYP", "Y"),
-            ("TYS", "Y"),
-        ]
-    )
 
     seq_dic: dict[str, dict[int, str]] = {}
 
     if isinstance(pdb_f, PDBFile):
         pdb_f = pdb_f.rel_path
 
-    coord_req = ("ATOM", "HETATM", ) if keep_hetatm else ("ATOM", )
-    with open(pdb_f) as fh:
+    coord_req = (
+        (
+            "ATOM",
+            "HETATM",
+        )
+        if keep_hetatm
+        else ("ATOM",)
+    )
+
+    with open(pdb_f, "r") as fh:
         for line in fh.readlines():
             if line.startswith(coord_req):
                 res_num = int(line[slc_resseq])
@@ -887,7 +899,7 @@ def pdb2fastadic(pdb_f: PDBPath, keep_hetatm: bool = False) -> dict[str, dict[in
                 if res_name in RES_TO_BE_IGNORED:
                     continue
                 try:
-                    one_letter = res_codes[res_name]
+                    one_letter = RES_3L1L[res_name]
                 except KeyError:
                     one_letter = "X"
                 if chain not in seq_dic:
@@ -897,7 +909,9 @@ def pdb2fastadic(pdb_f: PDBPath, keep_hetatm: bool = False) -> dict[str, dict[in
 
 
 def get_align(
-    method: str, lovoalign_exec: FilePath, keep_hetatm: bool = False,
+    method: str,
+    lovoalign_exec: FilePath,
+    keep_hetatm: bool = False,
 ) -> partial[dict[str, dict[int, int]]]:
     """
     Get the alignment function.
@@ -958,12 +972,12 @@ def align_strct(
     """
     if lovoalign_exec is None:
         log.error(
-            "Structural alignment needs LovoAlign " "get it at github.com/m3g/lovoalign"
+            "Structural alignment needs LovoAlign get it at github.com/m3g/lovoalign"
         )
         raise ALIGNError("Path to LovoAlign executable required.")
 
     if not lovoalign_exec:
-        raise ALIGNError("lovoalign_exec parameter not defined ")
+        raise ALIGNError("lovoalign_exec parameter not defined.")
 
     if not os.access(lovoalign_exec, os.X_OK):
         raise ALIGNError(f"{lovoalign_exec!r} for LovoAlign is not executable")
@@ -1009,7 +1023,7 @@ def align_strct(
                 alignment_end_index = i - 2
             elif "ERROR" in line:
                 failed_pdb = line.split()[-1]
-                _msg = f"LovoAlign could not read {failed_pdb} " "is it a ligand?"
+                _msg = f"LovoAlign could not read {failed_pdb} is it a ligand?"
                 log.warning(_msg)
                 alignment_pass = False
 
@@ -1019,7 +1033,7 @@ def align_strct(
         if not alignment_pass:
             # This alignment failed, move on to the next
             log.warning(
-                f"Skipping alignment of chain {chain}, " "used sequential matching"
+                f"Skipping alignment of chain {chain}, used sequential matching"
             )
             continue
 
@@ -1377,9 +1391,7 @@ def dump_as_izone(fname, numbering_dic, model2ref_chain_dict=None):
                     unb_chain = model2ref_chain_dict[chain]
                 #
                 izone_str = (
-                    "ZONE "
-                    f"{chain}{bound_res}:{unb_chain}{unbound_res}"
-                    f"{os.linesep}"
+                    f"ZONE {chain}{bound_res}:{unb_chain}{unbound_res}{os.linesep}"
                 )
                 fh.write(izone_str)