Frequently Asked Questions

We collected here a list of frequently occurring problems and their solutions. The following topics are currently available:

If your problem falls outside of the topics, please see the Getting support / How to ask for help section.

What about missing atoms?

Missing atoms will be automatically detected (if residues with these missing atoms are in the list of supported particles) and re-generated when running the [topoaa] module. For this reason, it is always used as the first module in a haddock3 workflow configuration file, not only to generate the topology of the input molecules but also to add and reconstruct missing atoms.

What about chain breaks?

In case of missing residues, chain breaks will be introduced. This might cause segments of your molecule to move with respect to each other during the refinement stages. To avoid that, you can define a few specific distance restraints, for example between CA atoms. This can be easily performed by the haddock3-restraints command line interface supporting the restrain_bodies subcommand that allows the detection of such breaks and define distance restraints. Here is the documentation to the haddock3-restraints restrain_bodies subcommand. Those restraints can then be provided to haddock3 as unambiguous restraints for example (using the unambig_fname parameter in CNS modules).

What about point mutations?

To introduce mutations in your input PDB files you can do the following:

  • edit the PDB file and rename the mutated residue to the proper amino acid name
  • keep or rename appropriately the matching side-chain atoms

The extra/missing atoms will be automatically detected and the corrected topology and coordinates will be regenerated by the [topoaa] module. It is important to have at least the backbone atoms and at least the CB atom along the side-chain defined since their average position will be used as a starting point to "grow" the missing atoms. Always check that the sequence of the various PDB files matches!

Note that this approach is only functional for residues supported by HADDOCK.

What about ions?

Some proteins contain ions such as for example calcium. Their inclusion might be important for docking purposes, in particular for proper electrostatics! In principle, they should be recognized when running the [topoaa] module, provided their names in the PDB file matches the ion names in the list of supported ions.

What about disulfide bonds?

Disulfide bonds within a molecule are automatically detected by the [topoaa] module: a bond is created between two cysteines whose SG atoms are closer than disulphide_dist (default: 2.2 Å). There is no option to select specific cysteines, but cysteines that are not close in space will not be bonded. If a disulfide bond is not detected because of a poor geometry of the input structure, you can increase disulphide_dist; setting it to 0 turns the detection off.

[topoaa]
disulphide_dist = 2.5

Since each input molecule is processed separately, no disulfide bond is created between two different molecules. To bring two cysteines from different molecules together, define a few distance restraints between their SG atoms (for example in the unambig_fname file).

Domain definition for docking

In general, it is recommended to remove any part of your system such as flexible linkers that are not involved in the interaction with the partner for docking. Keeping these might cause trouble in the sorting of solutions. For example, such a linker can make contact with the partner molecule, resulting in lower total energy, and, in that way, "bad" solutions could still be kept.

The same applies to AlphaFold2 spaghetti like disordered regions that often surround the domain of interest. Indeed, these regions may induce van der Walls forces due to sterical clashes before the two domains of interest could even interact. Removing regions with low pLDDT (~< 60) can be an appropriate solution so use AlphaFold2 models for docking.

Docking a molecule made of several chains

When one of the partners consists of several chains (e.g. the heavy and light chains of an antibody, a homotrimer or a multimeric receptor), it is best to define it as a single molecule with a single chain. This is also important for the scoring, since the HADDOCK score is calculated over all interfaces, including those between the chains of the same partner.

To do so:

  • shift the residue numbering of the chains so that they do not overlap, using pdb_shiftres;
  • merge them into a single chain, using pdb_chain and pdb_merge, followed by pdb_tidy;
  • keep the chains together during the refinement by defining distance restraints between them, using the haddock3-restraints restrain_bodies subcommand, and provide them as unambiguous restraints (unambig_fname).

Do not use pdb_reres on such a file if the total number of residues exceeds 9999, as PDB files cannot hold larger residue numbers. The antibody-antigen tutorial illustrates this procedure step by step.

Docking nucleic acids

A few settings are recommended when docking DNA or RNA:

  • to keep the base pairing and the conformation of the nucleic acids during the refinement, set dnarest_on = true in the refinement modules (see DNA & RNA restraints);
  • to account for the screening of the charged phosphate backbone, use a constant dielectric with a higher dielectric constant: dielec = "cdie" and epsilon = 78. Note that these settings apply to all interfaces of the system, including protein-protein ones;
  • the state of the 5' end of the nucleic acids is defined per molecule in the [topoaa] module, e.g. 5_phosphate = true under [topoaa.mol1].

See the protein-DNA docking scenario for a complete example.

Note that modified nucleotides are currently not supported: if defined as HETATM, they will not be connected to the rest of the strand.

My glycan falls apart after [topoaa]

Glycosidic linkages are not read from the input file, but detected from the coordinates by the [topoaa] module. A linkage is only created between supported carbohydrates when the anomeric carbon of the donor lies within 2 Å of the acceptor oxygen. If your glycan is split into separate monosaccharides:

  • check that all residue names match the supported carbohydrate codes (glycans built with other tools, e.g. GLYCAM, may use different names);
  • check that the linkage(s) in your glycan are part to the supported glycan linkages and that the geometry of the input is chemically sensible.

Chemically modified sugars are not supported: in that case the whole glycan should be treated as a ligand, providing its own topology and parameters (see Cofactors / Small-ligand docking with HADDOCK).

Clustering issues

When performing RMSD clustering, two modules can be used to compute the RMSD matrix:

  • [rmsdmatrix]: computing the full complex (or single chain) RMSD matrix
  • [ilrmsdmatrix]: computing the interface-ligand RMSD matrix

The [rmsdmatrix] module allows you to define a subset of residues used to perform both the structural alignment and the RMSD computation. For this, you need to specify a list of residues for each chain, using the parameter resdic_*, where * is the chainID. As an example, to perform the selection of residues 12, 13, 14 and 15 from chain A and 1, 2, 3 from chain B, refine the following parameters:

[rmsdmatrix]
resdic_A = [12, 13, 14, 15]
resdic_B = [1, 2, 3]

This will result in the selection of those 7 residues to perform the structural alignment onto the reference and then compute the RMSD.

While for the [ilrmsdmatrix] module, a different approach is taken. Two parameters must be defined

  • receptor_chain: defining the chainID of the receptor. By default "A".
  • ligand_chains: a list of other chain IDs that should represent the "ligands". If not set, all the remaining chains will be considered as ligand.

During the computational workflow, first, all the residue-residue contacts between the receptor and ligand are selected. This selection is then used to perform later structural alignment and RMSD computation.

Those two modules must be followed by the [clustrmsd] module, otherwise, only the pair-wise RMSD matrix will be computed, and clustering not performed. Note that this is not an issue if fractions of common contact (FCC) clustering ([clustfcc] module) is used as the matrix is computed within the clustering module directly (as much faster).

Scoring existing models of complexes

To compute the HADDOCK score of a single complex, obtained with HADDOCK or any other method, use the haddock3-score command line:

haddock3-score complex.pdb

To score many models, merge them into a single ensemble file (the number of input files is limited to 20, but an ensemble can contain any number of models):

pdb_mkensemble model_*.pdb | pdb_tidy > ensemble.pdb

and score them with a short workflow:

run_dir = "scoring"
molecules = ["ensemble.pdb"]

[topoaa]
[emscoring]
[caprieval]

A few points to keep in mind:

  • the scoring modules expect each complex as one PDB file, containing all chains, and not one file per chain;
  • the models are energy minimised before the scoring. If a model contains severe clashes, the minimisation can fail and no score is reported (None);
  • the HADDOCK score covers all interfaces of the complex. Set per_interface_scoring = true in [emscoring] to also obtain the HADDOCK score of each individual interface.

More examples are available in the scoring workflows section.

Comparing or clustering models from several runs

Each run is analysed separately. To compare or cluster models coming from several runs (e.g. docking runs with different conformations of the partners), merge them into a single ensemble and run a scoring workflow on it:

run_dir = "rescore-and-cluster"
molecules = ["ensemble-of-all-runs.pdb"]

[topoaa]
[emscoring]
[clustfcc]
[seletopclusts]
[caprieval]

Note that scores from individual runs can only be compared between runs using the same restraints, since the restraint energy is part of the HADDOCK score. No restraints are used when running a scoring scenario like the one above.

What does the HADDOCK score mean?

The HADDOCK score is a weighted sum of energy terms and restraint violations (see the scoring modules). It is expressed in arbitrary units: it is meant to rank models of the same complex, and is neither a binding free energy nor a dissociation constant. For this reason, scores should not be used to compare the affinity of different complexes (see the Dos and Don'ts of the analysis). To predict binding affinities, use the [prodigyprotein] and [prodigyligand] modules.

A positive HADDOCK score usually means that many restraints are not satisfied, resulting in a large restraint energy. This can be the case when many restraints were defined that cannot be satisfied at the same time, or when unambiguous restraints have incorrect distances.

Note that the Z-score reported by HADDOCK2.4 is not provided by haddock3: it is the number of standard deviations of the cluster score from the mean, and can easily be calculated from the cluster statistics.

Ranking models based on restraints only

To rank models based only on how well they satisfy a set of restraints (e.g. to filter the models of an ab-initio run with uncertain data), use a scoring module in which all weights are set to 0, except the one of the restraints:

[emscoring]
ambig_fname = "restraints.tbl"
w_vdw = 0
w_elec = 0
w_desolv = 0
w_bsa = 0
w_air = 1

Using [caprieval] without a reference structure

The [caprieval] module can be used without a reference structure: simply do not define reference_fname. In that case, the best scoring model of the previous step is used as reference.

[caprieval]

Even without a reference, it is useful to add it to the workflow, as it reports the scores and their components, and generates the plots of the analysis.

Nothing seems to be docked

If the molecules remain apart, or are not brought into contact, check that the docking is driven by some information. Without restraints (ambig_fname or unambig_fname) and without an ab-initio option (cmrest, surfrest or ranair), nothing brings the molecules together. See the ab-initio docking section if you do not have information about the interface.

Note that restraint files must be provided to each module using them (e.g. [rigidbody], [flexref] and [emref]), and not only to the first one.

How can I make my runs faster?

The computational cost of a run mostly depends on the number of models generated and refined and the number of restraints defined and their ambiguity level (the more ambiguous the more costly). To reduce it:

  • increase the number of cores used with the ncores parameter;
  • provide meaningful restraints, which allows to safely reduce the sampling in [rigidbody] (default: 1000 models);
  • reduce the number of models selected for refinement (select parameter in [seletop], default: 200);
  • cluster the models ([clustfcc]) and refine only the best models of each cluster ([seletopclusts]);
  • in ab-initio docking, do not combine cmrest and surfrest, as this increases the computational time without benefit.

Running HADDOCK on a cluster using a queuing system (e.g. Torque or Slurm)

In order to submit to the queuing system we typically use a wrapper script that will add some directives to the job files.

First, we must define a haddock3 workflow (e.g.: haddock_run.cfg):

#################################
# GLOBAL PARAMETERS
#################################
run_dir = "amazing_docking_experiment"
molecules ["protein1.pdb", "protein2.pdb"]
# Here we define the maximum number of available cores to use
ncores = 40

#################################
# WORKFLOW MODULES PARAMETERS
#################################
[topoaa]
[rigidbody]
[seletop]
[flexref]
[emref]
[clustfcc]
[seletopclusts]
[contactmap]
[caprieval]

Here is one example of such a wrapper script (named haddock3_run.job) that would submit to the slurm queue:

#!/bin/bash
#SBATCH --nodes=1
#SBATCH --tasks-per-node=40
#SBATCH -J haddock3-run
#SBATCH -p short

# 1. Active the haddock3 virtual environment
# From venv
source /path/to/haddock3/install/dir/.haddock3-env/bin/activate
# Or using conda
source /path/to/conda/install/dir/bin/activate
conda activate haddock3-env

# Go to the base directory where the `haddock_run.cfg` workflow is written
cd /path/to/workflow/

# Execute haddock3 on the `haddock_run.cfg` workflow
haddock3 haddock_run.cfg

Note here that we set up the number of cores for both the haddock3 run (in haddock_run.cfg) and slurm job (in haddock3_run.job)to 40

Cofactors / Small-ligand docking with HADDOCK

It's possible to dock small ligands or cofactor using haddock3, but for that topology and parameter files for the ligand are required in CNS format. Residues for which no topology is available are not kept in the models, which is why a ligand present in an input file can disappear from the output models.

The simplest option is to let haddock3 generate them automatically with PRODRG, by setting autotoppar = true in the [topoaa] module. This also works for ligands bound to a protein, and for several different ligands in the same system (see Dealing with non-standard molecules).

Note that the ligand will not be included if coarse-graining is used. They will however be present after the back conversion from CG to AA models (the [cgtoaa] module).

Alternatively, the topology and parameter files can be generated with other tools. Several sources exist to find such files:

  • ccp4-prodrg: ccp4-prodrg.

  • the PRODRG server was maintained by Daan van Aalten at Dundee University. This server allows you to draw your molecule or paste coordinates and will return topologies and parameter files in various formats, including CNS. You should turn on the electrostatic to obtain partial charges. Save the resulting PDB file and the corresponding CNS parameter and topology files to use in HADDOCK.

    Important: The generated parameter file contains a CNS NBONds statement which should be removed prior to their use in HADDOCK. Look in the parameter file for:

     NBONds

CUTNB=7.0 WMIN=1.5 REPEL=1.0 REXPONENT=4 IREXPONENT=1 RCONST=16.0 TOLERANCE=0.5 NBXMOD=5 CTONNB=5.5 CTOFNB=6.0 END

and remove or comment it out (by adding ! before each line).

  • the Automated Topology Builder (ATB) and Repository developed in Prof. Alan Mark's group at the University of Queensland in Brisbane: https://atb.uq.edu.au/
    Note: we have not yet tested those parameters in HADDOCK.

For docking small ligands with haddock3 using custom-made topology and parameter files, you should:

  • Define the path to the files in CNS modules ([topoaa], [flexref], [emref], [mdref], [emscoring], [mdscoring])
    • Input the topology file using ligand_top_fname parameter.
    • Input the parameter file using ligand_param_fname parameter.
  • When using several ligands, combine all their topologies in a single file, and all their parameters in another one. Their residue names and atom types must be unique, and must not overlap with the ones already defined in HADDOCK.

Note for D-amino acids: they do not need specific topologies. Name them as the corresponding L-amino acid (e.g. SER for D-serine): the chirality is taken from the input structure.

Also, we recommend setting the number of MD steps for the first two parts (rigid-body high temperature dynamic and slow cooling annealing) of the [flexref] module to 0.
This is performed by tuning the mdsteps_rigid and mdsteps_cool1 parameters and setting their values to 0.

Haddock3 comes with an example for protein-ligand docking. Check the setting in that example.

Important: When starting a run, always check for error messages in the 0_topoaa directory in the various generated .out files, especially for your ligand.

Beads dummy atoms docking with haddock3

Dummy atoms can be used in haddock3 and can be useful as distance restraints can be built towards them. This is used, for example:

Dummy atoms (also called shape beads) must be defined in a separate PDB file and have the following naming convention:

  • Start with the ATOM
  • using SHA for both atom and residue name
  • defined as chain S
  • have the same atom and residue index
ATOM      1  SHA SHA S   1      24.222  -6.426 -14.545  1.00  1.00
ATOM      2  SHA SHA S   2      23.059  -6.675 -14.930  1.00  1.00
...

Because such types of dummy atoms neither have topology nor force-field parameters, they must be explicitly defined as shapes. To do so, two parameters must be set in your configuration file:

  • mol_shape_X = true: allows to tell haddock3 that molecule X is a shape.
  • mol_fix_origin_X = true: allows telling haddock3 not to move molecule X, and keep the original coordinates.

Where X is a number that corresponds to the molecule position in the input list of molecules in the configuration file.

Here is an example, where the shapes will be input at the second position in the molecules:

run_dir = "test_shape"
molecules = ["protein.pdb", "shapes.pdb"]
[topoaa]
[rigidbody]
  # `_2` as the shape is placed second in the input molecules
mol_shape_2 = true  # Defines second input molecule as `shape`
mol_fix_origin_2 = true  # Fix origin/input coordinates of the second input molecule

Typical haddock3 error messages

In some cases, the haddock3 execution can stop for a given reason. While we are already trying the handle possible errors, some of them will lead to critical failure, terminating the workflow. If such an error occurs, we report it in the log file. The log file can be found at two locations:

  • printed on your screen as standard output.
  • written in a file named log located in the workflow run directory.

We often try to provide a meaningful error message that can help you figure out what could be the issue related to it. If not, please refer to the Getting support / How to ask for help section to get assistance.

Here is a list of the most common errors:

Tolerance issue

Here is a typical tolerance issue log error message:

[2024-09-09 20:04:33,709 libutil ERROR] 100.00% of output was not generated for this module and tolerance was set to 5.00%.
Traceback (most recent call last):
  File "/data/haddock3/src/haddock/libs/libutil.py", line 335, in log_error_and_exit
    yield
  File "/data/haddock3/src/haddock/clis/cli.py", line 192, in main
    workflow.run()
 File "/data/haddock3/src/haddock/libs/libworkflow.py", line 43, in run
    step.execute()
 File "/data/haddock3/src/haddock/libs/libworkflow.py", line 162, in execute
    self.module.run()  # type: ignore
 File "/data/haddock3/src/haddock/modules/base_cns_module.py", line 61, in run
    self._run()
 File "/data/haddock3/src/haddock/modules/sampling/rigidbody/__init__.py", line 246, in _run
 self.export_io_models(faulty_tolerance=self.params["tolerance"])
  File "/data/haddock3/src/haddock/modules/__init__.py", line 300, in export_io_models
    self.finish_with_error(_msg)
  File "/data/haddock3/src/haddock/modules/__init__.py", line 308, in finish_with_error
    raise RuntimeError(reason)
RuntimeError: 100.00% of output was not generated for this module and tolerance was set to 5.00%.

This means that models that should have been generated by a module are missing from the file system, and therefore were not written. This can come from multiple reasons:

  • There is an issue with the parameters/topology generation.
  • The model contained clashes that led to extremely high energetics that blew up the system.
  • The tolerance threshold was set too low.

If 100.00% of output was not generated, there is probably an issue with the input molecules:

  • unrecognized amino acids or missing parameters/topology for residues outside of the supported particles list
  • huge sterical clash

If the value is inferior to 100.00%, it can come from either one of the input conformers or a random error from the molecular dynamic simulation. In this case, you could increase the tolerance threshold to a higher value (e.g.: tolerance = 10), for the workflow to continue, as missing a few models can be acceptable. Note that you could also restart the workflow from the next module (e.g.: haddock3 workflow.cfg --restart 5 to restart from module 5 if the module did not meet the tolerance threshold), allowing you to save computational resources not having to recompute data from module 0 to 4.

To find the origin of the problem:

  • set debug = true in the global parameters of the configuration file, so that the input and output files of the CNS jobs are kept;
  • check the .out files in the directory of the failing module (and in 0_topoaa) for error messages;
  • run one of the examples provided with haddock3, to make sure that the installation is working.

Configuration file syntax errors

The configuration file follows the TOML syntax, and syntax errors are reported with their position in the file. Common causes are:

  • a missing comma between the elements of a list, e.g. in molecules = ["mol1.pdb", "mol2.pdb"] (Found tokens after a closed string);
  • special characters or "smart" quotes introduced by a text editor or word processor (Reserved escape sequence used, or quotes not recognised). Edit configuration files with a plain-text editor.