Make the coarse-graining deterministic and fix the dummy-bead geometry (#1657)

CHANGELOG.md

@@ -1,5 +1,6 @@
 # Changelog
 
+- 2026-08-07: Fixed the martinise library (libaa2cg) to make it deterministic - Issue #1657
 - 2026-08-05: Removed identical duplicate `ANGLe` statements from `protein-allhdg5-4.param` and `protein-CG-Martini-2-2.param` - Issue #1589
 - 2026-08-04: Added workflow module ordering validation - related to Issue #1530
 - 2026-08-02: Fixed logging/warning leaks - Issue #1647

examples/refine-complex/refine-complex-CG-test.cfg

@@ -0,0 +1,38 @@
+# ====================================================================
+# Refinment of a complex example
+
+# directory in which the scoring will be done
+run_dir = "run1-CG-refinement"
+
+# execution mode
+mode = "local"
+ncores = 10
+
+# molecules to be refined (given as a single complex PDBs)
+molecules =  [
+    "data/e2a-hpr_1GGR_A.pdb",
+    "data/e2a-hpr_1GGR_B.pdb"
+    ]
+
+# ====================================================================
+# Parameters for each stage are defined below, prefer full paths
+# ====================================================================
+[topoaa]
+
+[topocg]
+
+# required to create the complex for cgtoaa
+[emref]
+nemsteps = 0
+
+[cgtoaa]
+# generate five model for each input model
+sampling_factor = 5
+
+[emref]
+
+[caprieval]
+reference_fname = "data/e2a-hpr_1GGR.pdb"
+
+# ====================================================================
+

src/haddock/libs/libaa2cg.py

@@ -56,6 +56,11 @@
 
 CRYST_LINE = "CRYST1 " + os.linesep
 
+# Fallback seed for the ``SCD*`` dummy beads, used by callers that have no
+# ``iniseed`` parameter of their own. Kept equal to the ``topocg`` default; the
+# test suite asserts that.
+DEFAULT_SEED = 917
+
 
 def norm(a):
     """
@@ -373,47 +378,44 @@ def ss_classification(ss, program="dssp"):
 charged = ["ARG", "LYS", "ASP", "GLU"]
 
 
-def add_dummy(bead_list, dist=0.11, n=2):
-    """
+def add_dummy(bead_list, rng, dist, n):
+    """Place the ``SCD*`` dummy beads around the last bead of a residue.
+
+    The beads are laid out along a randomly oriented axis through the parent
+    bead: the first at ``+dist``, the second (if any) at ``-dist``, so that a
+    pair straddles the parent symmetrically. Only the *orientation* is random.
 
     Args:
-        bead_list:
-        dist:
-        n:
+        bead_list: ``(name, coordinates)`` pairs of the residue's real beads.
+            The dummy beads are attached to the last one.
+        rng: ``random.Random`` instance used to orient the dummy beads.
+        dist: Distance from the parent bead, in angstrom (see ``map_cg``).
+        n: Number of dummy beads, 1 or 2.
 
     Returns:
-
+        Mapping of bead name (``SCD1``, ``SCD2``) to its coordinates.
     """
-    new_bead_dic = {}
-
-    # Generate a random vector in a sphere of -1 to +1, to add to the bead position
-    v = [
-        random.random() * 2.0 - 1,
-        random.random() * 2.0 - 1,
-        random.random() * 2.0 - 1,
-    ]
-
-    # Calculated the length of the vector and divide by the final distance of the dummy bead
-    norm_v = norm(v) / dist
-
-    # Resize the vector
-    vn = [i / norm_v for i in v]
-
-    # m sets the direction of the added vector, currently only works when adding one or two beads.
-    m = 1
-    for j in range(n):  # create two new beads
-        bead_s = str(j + 1)
-        new_name = f"SCD{bead_s}"  # set the name of the new bead
-        new_bead_dic[new_name] = [i + (m * j) for i, j in zip(bead_list[-1][1], vn)]
-        m *= -2
-    return new_bead_dic
+    # Random direction, uniform on the unit sphere. Normalising a vector drawn
+    # from a cube would bias the direction towards the cube's corners; a
+    # Gaussian vector is isotropic, so normalising it is unbiased.
+    v = [rng.gauss(0.0, 1.0) for _ in range(3)]
+    scale = dist / norm(v)
+    vn = [i * scale for i in v]
+
+    parent_coord = bead_list[-1][1]
+    # The signs straddle the parent bead, and cap the layout at two beads.
+    return {
+        f"SCD{idx}": [coord + sign * offset for coord, offset in zip(parent_coord, vn)]
+        for idx, sign in enumerate((1, -1)[:n], start=1)
+    }
 
 
-def map_cg(chain):
+def map_cg(chain, rng):
     """
 
     Args:
         chain:
+        rng: ``random.Random`` instance used to orient the dummy beads.
 
     Returns:
 
@@ -500,13 +502,17 @@ def map_cg(chain):
             m_dic[aares][bead_name] = bead_coord, code, restrain
 
     # add dummy beads whenever its needed
+    # Distances are in angstrom, matching the SCd bond lengths of
+    # cns/toppar/protein-CG-Martini-2-2.param: a pair placed at +/-1.4 A is
+    # 2.8 A apart (BOND SCd SCd), and a lone bead sits at the 1.1 A
+    # parent-SCd bond length.
     for r in m_dic:
         if r.resname in polar:
-            d = 0.14  # distance
+            d = 1.4  # distance
             n = 2  # number of dummy beads to be placed
 
         elif r.resname in charged:
-            d = 0.11  # distance
+            d = 1.1  # distance
             n = 1  # number of dummy beads to be placed
 
         else:
@@ -515,7 +521,7 @@ def map_cg(chain):
         # add to data structure
         # this special beads have no HADDOCK code
         bead_list = [(b, m_dic[r][b][0]) for b in m_dic[r]]
-        dummy_bead_dic = add_dummy(bead_list, dist=d, n=n)
+        dummy_bead_dic = add_dummy(bead_list, rng, dist=d, n=n)
         for db in dummy_bead_dic:
             db_coords = dummy_bead_dic[db]
             # code should be the same as the residue
@@ -733,66 +739,6 @@ def identify_pairing(ra, rb):
     return pair
 
 
-def output_cg_restraints(pair_list):
-    """
-
-    Args:
-        pair_list:
-
-    Returns:
-
-    """
-    out = open("dna_restraints.def", "w")
-    for i, e in enumerate(pair_list):
-        idx = i + 1
-        res_a = e[0][0]
-        segid_a = e[0][1]
-        res_b = e[1][0]
-        segid_b = e[1][1]
-        out.write(
-            f"{{===>}} base_a_{idx}=(resid {res_a} and segid {segid_a});\n"
-            f"{{===>}} base_b_{idx}=(resid {res_b} and segid {segid_b});\n\n"
-        )
-    out.close()
-
-
-def extract_groups(pair_list):
-    """
-
-    Args:
-        pair_list:
-
-    Returns:
-
-    """
-    # this will be used to define AA restraints
-    out = open("dna-aa_groups.dat", "w")
-    # extract groups
-    group_a = [a[0][0] for a in pair_list]
-    segid_a = list(set([a[0][1] for a in pair_list]))
-
-    group_b = [a[1][0] for a in pair_list]
-    segid_b = list(set([a[0][1] for a in pair_list]))
-
-    if len(segid_a) != 1:
-        emsg = "Something is wrong with SEGID A"
-        raise ModuleError(emsg)
-
-    if len(segid_b) != 1:
-        emsg = "Something is wrong with SEGID B"
-        raise ModuleError(emsg)
-
-    segid_a = segid_a[0]
-    segid_b = segid_b[0]
-
-    group_a.sort()
-    group_b.sort()
-    out.write(
-        f"{group_a[0]}:{group_a[-1]}\n{segid_a}\n{group_b[0]}:{group_b[-1]}\n{segid_b}"
-    )
-    out.close()
-
-
 def create_file_with_cryst(pdb_file: str) -> None:
     """
     This function creates a new pdb because the CRYST line is missing from the pdf file.
@@ -930,6 +876,7 @@ def martinize(
     input_pdb: str,
     output_path: str,
     skipss: bool,
+    seed: int = DEFAULT_SEED,
 ) -> tuple[str, bool]:
     """
     Converts an all-atom (AA) PDB structure into a coarse-grained (CG) model
@@ -947,6 +894,8 @@ def martinize(
             If True, skips secondary structure assignment (DSSP step).
             If False, assigns secondary structure and encodes it
             into HADDOCK-compatible B-factors.
+        seed (int):
+            Pseudo-random seed used to orient the ``SCD*`` dummy beads.
 
     Returns:
         tuple[str, bool]:
@@ -959,6 +908,11 @@ def martinize(
         emsg = "No input file detected"
         raise ModuleError(emsg)
 
+    # Dedicated generator, built fresh per call: never touch the global
+    # `random` state, which would make the result depend on whatever else ran
+    # first in this interpreter.
+    rng = random.Random(seed)
+
     p = PDBParser()
     io = PDBIO()
 
@@ -984,10 +938,12 @@ def martinize(
     # WARNING, THIS ASSUMES THAT INPUT DNA/RNA IS 3-LETTER CODE
     rename_nucbases(aa_model)
 
-    # Assign HADDOCK code for hydrogen bonding capable nucleotides (0-1)
-    pair_list = determine_hbonds(aa_model)
-    if pair_list:
-        output_cg_restraints(pair_list)
+    # Assign HADDOCK code for hydrogen bonding capable nucleotides (0-1).
+    # The returned pair list is not used here: base-pair restraints for CG
+    # models are derived in CNS by `dna-rna_restraints.cns`. What matters is
+    # the side effect, marking the paired bases with bfactor 1 so that
+    # `patch-types-cg-hbond-dna-rna.cns` can patch their bead types.
+    determine_hbonds(aa_model)
 
     # Map CG beads to AA structure
     structure_builder = StructureBuilder()
@@ -1003,7 +959,7 @@ def martinize(
             structure_builder.init_chain(chain.id)
             structure_builder.init_seg(chain.id)
 
-            mapping_dic = map_cg(chain)
+            mapping_dic = map_cg(chain, rng)
 
             for residue in mapping_dic:
                 if residue.id[0] != " ":  # filter HETATMS

src/haddock/modules/refinement/cgtoaa/__init__.py

@@ -27,7 +27,7 @@
 
 
 class HaddockModule(BaseCNSModule):
-    """HADDOCK3 module energy minimization refinement."""
+    """HADDOCK3 module for CG to AA conversion."""
 
     name = RECIPE_PATH.name
 
@@ -45,11 +45,15 @@ def confirm_installation(cls) -> None:
 
     def _run(self) -> None:
         """Execute module."""
+        # Pool of jobs to be executed by the CNS engine
+        jobs: list[CNSJob] = []
+
         # Get the models generated in previous step
         try:
             models_to_refine = self.previous_io.retrieve_models(individualize=True)
         except Exception as e:
             self.finish_with_error(e)
+
         self.output_models = []
         sampling_factor = self.params["sampling_factor"]
         if sampling_factor == 0:
@@ -67,8 +71,6 @@ def _run(self) -> None:
                 " decrease the sampling_factor."
             )
 
-        # Pool of jobs to be executed by the CNS engine
-        jobs: list[CNSJob] = []
         idx = 1
         for model in models_to_refine:
             if isinstance(model, PDBFile):

src/haddock/modules/refinement/cgtoaa/cns/cgtoaa.cns

@@ -25,7 +25,6 @@ evaluate ($ini_count    =1)
 
 evaluate ($data.ncomponents=$ncomponents)
 
-
 ! non-bonded parameter set to use
 evaluate ($toppar.par_nonbonded = "OPLSX" )
 
@@ -90,6 +89,8 @@ else
     set message=off echo=off end
 end if
 
+! initialize random number generator
+set seed $seed end
 
 {* Change segid of CG model and read AA pdb, psf files =========== *}
 

src/haddock/modules/topology/topocg/__init__.py

@@ -56,8 +56,13 @@ def generate_topology(
     write_to_disk: Optional[bool] = True,
     force_field: str = "martini2",
     shape: bool = False,
+    seed: int = 916,
 ) -> Union[Path, str]:
-    """Generate a HADDOCK topology file from input_pdb."""
+    """Generate a HADDOCK topology file from input_pdb.
+
+    ``seed`` is the ``iniseed`` parameter; it seeds the placement of the CG
+    dummy beads so that the coarse-graining is reproducible across runs.
+    """
     # generate params headers
     general_param = load_workflow_params(**defaults)
     input_mols_params = load_workflow_params(param_header="", **mol_params)
@@ -70,7 +75,7 @@ def generate_topology(
 
     if not shape:
         # AA to CG
-        cg_pdb_name = martinize(input_pdb, output_path, False)
+        cg_pdb_name = martinize(input_pdb, output_path, False, seed=seed)
         output = prepare_output(
             output_pdb_filename=f"{Path(cg_pdb_name).stem}_{force_field}{input_pdb.suffix}",
             output_psf_filename=f"{Path(cg_pdb_name).stem}_{force_field}.{Format.TOPOLOGY}",
@@ -262,6 +267,7 @@ def _run(self) -> None:
                     write_to_disk=self.params["debug"],
                     force_field=force_field,
                     shape=shape_dic[i],
+                    seed=self.params["iniseed"],
                 )
                 self.log("Topology CNS input created")