Ligand preparation. Ligand structures identified from PDB were built in their native form restoring the active warheads using the primary literature citations, (e.g.: reconstruction of the β-lactam ring, addition of leaving groups, restoring unsaturated bonds, etc.). Initial 3D coordinates were generated using OpenBabel modeling protonation states at pH 7.4, then minimized (MMFF94s; 300 steps Steepest Descent; 300 Conjugated Gradient). Partial charges, torsions, standard and reactive atom type parameters were assigned according to the AutoDock protocol using Meeko () with SMARTS patterns to define the warhead atoms and assign the reactive docking force field parameters (Fig. ). If present, macrocyclic structures were modeled as flexible by default during docking . Target preparation. Target structures were retrieved from the Protein Data Bank and hydrogens were added with Reduce . For oligomeric structures, dockings were performed only on the first chain. Partial charges, torsions, and standard and reactive atom type parameters were assigned using Meeko according to the AutoDock protocol . For each site, a cubic docking box of 30 Å side (80 points in the AutoGrid parameter file ) was defined and centered on the Cα of the target residue. Reactive docking. By default, the residues to be evaluated by reactive docking are automatically calculated on each target structure using MSMS to identify all solventaccessible residues (default probe radius: 1.5 Å) of a given type (i.e.: cysteines) including buried cavities (Fig. ). Optionally, a user-defined list can be provided. Then, each ligand is docked against the individual residues to be evaluated. During docking, ligands are modeled in their unmodified form with reactive warheads in place using a conventional (untethered) docking method, while the side chain of the target residue is modeled as flexible . Dockings are performed using AutoDock-GPU , generating 50 poses for each ligand using the default Lamarckian Genetic Algorithm (LGA) parameters . A ligand and a residue are considered reacting if the distance between their reactive atoms in the lowest energy pose is ≤2.0 Å. Residues are then ranked by their likelihood of reacting based on the best energy of the ligand(s) predicted to react with it. Reactive docking parameters do not affect the length nor the complexity of the calculation, resulting in nearly identical docking times as conventional dockings. Datasets. In order to calibrate and validate the reactive docking parameters, covalent complexes were collected from the PDB , containing a very diverse pool of residues (Fig. ) and warheads (Fig. ). Complexes were visually inspected to discard highly distorted or problematic structures selecting a dataset of 431 structurally diverse ligands and chemical reactions. The dataset was then subdivided into specialized sets, with no overlap between them. The list of all PDBs and warheads in the dataset is available in Supplementary Material. Training and test sets. First, we selected a training set of 80 complexes using the most represented residues (cysteine, serine) for which the largest number of complexes with the most diverse ligands and protein families is available. The training set was built by randomly selecting 20 complexes representative of each of the four most abundant warheads (Fig. ): acrylamides and chloroacetamides for cysteine; boronic acids and βlactams for serine. The training set was used to calibrate the reactive docking force field parameters. The remaining set of complexes was used to build the test set of 351 complexes, which includes 8 residue types (cystine, serine, threonine, lysine, tyrosine, glutamic acid, histidine, and methionine, Fig. ) and 37 warheads (Fig. ). The test set was used to assess the performance and transferability of the parameters obtained with the calibration process on the training set. Virtual screening set. A virtual screening set was built based on the experimental data published by Resnik et al. which used a library of 993 reactive fragments functionalized with mild electrophiles including chloroacetamides (n = 752, 76%) and acrylamides (n=241, 24%) targeting cysteines. This virtual screening set contains three targets for which structure coordinates were available and potent and selective fragments were identified: two enzymes, deubiquitinase OTUB2 and the pyrophosphatase NUDT7 , and K-Ras . Parameter selection. The training set was used to calibrate the reactive docking parameters to obtain a single set of weights to be used to maximize the predictive accuracy. The req value of 1.7 Å is fixed for all reactions and was chosen for being between the lengths of the shortest bond in our dataset (C-C bond, ~1.5 Å) and the longest (C-S bond, ~1.8 Å). After preliminary tests, values were obtained by performing an exhaustive search with respect to ε value (from 2.0 to 10.0 kcal/mol, 0.5 intervals), and the scaling factors W1, 3, and W1, 4 (0.0 to 0.9, 0.1 interval). During the calibration, each parameter set was tested with the entire training set by docking each ligand on its cognate target including all the solvent-accessible residues that would be modified by its warhead (i.e.: all solvent-accessible cysteines for acrylamides). Virtual screening protocol. Firstly, we assessed the ability of the reactive docking to find the correct modified cysteine from the pool of all solvent-accessible cysteines in each structure. We calculated the average, median, and standard deviation of the docking score of the whole library for each (Table ). Then the average docking score of the whole library docking against each target calculation was used to rank the most likely cysteine to be alkylated by the fragments. True positive success rates for each target were then calculated by ranking ligands predicted to react with the target residue by their best docking score, and true hit rates of known binders were calculated for the top 0.5%, 1%, and 10% of the docking results.