BindFlow: a free, user-friendly pipeline for absolute binding free energy calculations using free energy perturbation or MM(PB/GB)SA

BindFlow: a free, user-friendly pipeline for absolute binding free energy calculations using free energy perturbation or MM(PB/GB)SA | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results BindFlow: a free, user-friendly pipeline for absolute binding free energy calculations using free energy perturbation or MM(PB/GB)SA View ORCID Profile Alejandro Martínez León , View ORCID Profile Lucas Andersen , View ORCID Profile Jochen S. Hub doi: https://doi.org/10.1101/2025.09.25.678545 Alejandro Martínez León 1 Theoretical Physics and Center for Biophysics, Universität des Saarlandes , 66123 Saarbrücken, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Alejandro Martínez León Lucas Andersen 1 Theoretical Physics and Center for Biophysics, Universität des Saarlandes , 66123 Saarbrücken, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Lucas Andersen Jochen S. Hub 1 Theoretical Physics and Center for Biophysics, Universität des Saarlandes , 66123 Saarbrücken, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Jochen S. Hub Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF Abstract We present BindFlow, a Python-based software for automated absolute binding free energy (ABFE) calculations at the free energy perturbation (FEP) or at the molecular mechanics Poisson-Boltzmann/generalized Born surface area [MM(PB/GB)SA] level of theory. BindFlow is free, open-source, user-friendly, easily customizable, runs on work-stations or distributed computing platforms, and provides extensive documentation and tutorials. BindFlow uses GROMACS as molecular dynamics engine and provides built-in support for the small-molecule force fields GAFF, OpenFF, and Espaloma. We test BindFlow by computing affinities for 139 ligand/target pairs, involving eight different targets including six soluble proteins, one membrane protein and one non-protein host– guest system. Quantified by Pearson, Kendall, and Spearman correlations coefficients, we find that the agreement of BindFlow predictions with experiments are overall similar to gold standards in the field. Interestingly, we find that MM(PB/GB)SA achieves correlations that, for some systems and force fields, approach those obtained with FEP, while requiring only a fraction of the computational cost. This study establishes BindFlow as a validated and accessible tool for ABFE calculations. Introduction During the early stages of drug discovery, a large number of chemical compounds are evaluated with the aim of finding potent binders, a process that is both time-consuming and resource-intensive. 1 , 2 Thus, computational models are widely used to estimate the binding free energy of small molecules to biological targets in order to prioritize compounds for followup synthesis and experimental evaluation. 3 – 6 In the 1950s, methods collectively known as alchemical free energy perturbation were introduced. 7 The term “free energy perturbation” (FEP) referred originally to a specific class of alchemical methods 8 but has since been more generally applied to alchemical binding free energy methods, a notion followed in this study. FEP enabled the estimation of free energy differences at a fraction of the computational cost compared to conventional MD simulations. Nevertheless, FEP remained computationally demanding for high-throughput studies, rationalizing the development of computationally more efficient yet more approximate endpoint free energy techniques such as the molecular mechanics Poisson-Boltzmann/generalized Born surface area [MM(PB/GB)SA] methods, 9 which achieved varying degrees of success. 10 Recent advances in biomolecular simulation methods 6 , 11 combined with the growth of computational power have enabled increasingly accurate affinity predictions, 6 , 12 making both FEP and MM(PB/GB)SA calculations more routine. Binding free energy estimation have been categorized into relative (RBFE) and absolute binding free energy (ABFE) calculations. 11 RBFE is typically preferred for ranking congeneric molecular series 13 or for evaluating the effects of mutations on ligand binding. 14 Alchemical RBFE methodologies have seen substantial advancements, resulting in numerous implementations. 15 – 29 These developments have enabled alchemical RBFE methods to achieve remarkable accuracy, often with a root mean-squared error (RMSE) of 1–2 kcal mol −1 relative to experimental values, 30 making them highly effective for drug discovery applications. 13 – 15 , 31 – 37 ABFE calculations yield the binding free energy relative to a standard state of the com-pound in solution. 38 ABFE calculations are particularly useful for studying highly diverse molecular sets, 4 , 39 – 41 for binding pose validation, 42 , 43 or for multi-target selectivity prediction. 44 , 45 Alchemical ABFE methods have demonstrated their power in addressing complex challenges across various projects. 4 , 34 , 39 , 40 , 44 – 48 Several tools for alchemical ABFE calculations have been proposed recently including ABFE workflow, 49 A3FE, 50 BAT2 51 (successor to BAT.py 23 ), SAFEP, 52 FEP+, 46 BFEE2 53 , 54 (successor to BFEE 55 ), HTBAC, 56 XFEP, 17 CHARMM-GUI, 29 and FEPSetup. 57 Since such tools automate significant portions of the binding free energy calculation pipeline, large-scale comparisons with experimental data have become accessible, reveling that alchemical ABFE predictions are approaching experimental accuracy. 58 However, the preparation, execution, and analysis of ABFE simulations remain at times tedious, system-dependent, often reliant on expert intervention, and, in some cases, associated with substantial costs due to software licensing. Fully automating the binding free energy pipeline has the potential to overcome these challenges, enabling high-throughput calculations for modern drug design campaigns while enhancing reproducibility and accessibility. We aim to close this gap by introducing BindFlow, a free open-source software that automates the entire ABFE pipeline. BindFlow implements two ABFE methods: (i) the endpoint free energy method MM(PB/GB)SA 9 based on a simulation of the receptor–ligand complex only, also referred to as single trajectory approach, and (ii) a double-decoupling alchemical free energy method that uses thermodynamic integration (TI) 59 or the multistate Bennett acceptance ratio (MBAR) 60 frameworks for free energy estimation. BindFlow was inspired by Alibay et al. 5 and Ries et al. 49 . BindFlow is highly user-friendly, provides extensive tutorials and documentation, natively supports the small-molecule force fields GAFF, OpenFF, and Espaloma. It schedules tasks either on a local desktop or to SLURM-based distributed computing environments, while being resilient against rare hardware failures or simulation instabilities. For advanced users, it allows full controls over the pipeline and MD settings, and it provides various options of further customization. Below, we first present implementation principles of BindFlow and concepts of the Bind-Flow API (application programming interface). Next, we validate BindFlow by computing ABFEs for 139 diverse ligands binding to eight different receptors including six soluble proteins, one membrane protein, and one host–guest system ( Fig. 1b ). The calculations include, both, FEP and MM(PB/GB)SA methods using three widely used open-source small-molecule force fields: GAFF-2.11, 61 OpenFF-2.0.0, 62 and Espaloma-0.3.1, 63 referred to GAFF, OpenFF, and Espaloma in the following. Download figure Open in new tab Figure 1. ( a ) Schematic representation of the BindFlow workflow and ( b ) the systems used in this study for validating ABFE calculations. Overall, our results show that BindFlow delivers predictive performance on par with established gold standards in the field. Strikingly, we observe that for certain targets and force fields, the computationally efficient MM(PB/GB)SA method achieves correlations with experiment that approach those of the more rigorous FEP calculations—yet at only a fraction of the computational cost. By fully automating the binding free energy workflow and minimizing user intervention, BindFlow enhances not only efficiency but also reproducibility. Together, these features make BindFlow a practical, accessible, and reliable platform for large-scale binding affinity prediction, with strong potential for accelerating modern drug discovery campaigns. Theory The theory of FEP and MM(PB/GB)SA have been reviewed frequently, 9 , 10 , 12 , 60 , 64 – 66 hence we provide only a brief summary of the underlying concepts. MM(PB/GB)SA Following the thermodynamic cycle in Fig. 2a , MM(PB/GB)SA estimates the binding free energy Δ G bind as follows: 9 , 10 , 65 , 66 Download figure Open in new tab Figure 2. Thermodynamic cycles for ( a ) MM(PB/GB)SA and ( b ) FEP. denotes the solvation free energy for species X, where superscripts R, L, and RL denote the receptor, ligand, and receptor–ligand complex, respectively. is the binding free energy in vacuum. Superscripts “solv” and “rcp” indicate states with the ligand in solvent or in the receptor, respectively. indicates. Superscripts “dcpl” and “cpl” indicate decoupling or coupling processes, while “coul” and “vdw” specify transitions of Coulomb or Lennard-Jones interactions. Superscripts “on/off” specify the activation or deactivation of Boresch restraints (subscript “rest”). The light pink, gray, and light blue backgrounds represent the implicit solvation model, vacuum, and explicit water model, respectively. Here, denotes the solvation free energy for species X (X = R, L, or RL), where superscripts R, L, and RL denote receptor, ligand, and receptor–ligand complex, respectively. The solvation free energy for species X is decomposed into polar and non-polar terms, is obtained by solving the Poisson-Boltzmann (PB) equation or using the generalized Born (GB) model, 10 while is derived from the solvent-accessible surface area (SASA). 67 denotes the binding free energy in vacuum, computed from the free energies of receptor, ligand, or receptor–ligand complex. The free energies are estimated via where is the molecular mechanics (MM) potential energy of species X comprising bonded and non-bonded terms, as defined by the MM force field. T is the temperature, and S X is the entropy of species X. By defining the free energy of species X as Δ G bind may be written as the free energy difference between products and reactants, thereby rationalizing by the subscripts of Eq. 3 the name MM(PB/GB)SA. 9 , 10 , 65 , 66 BindFlow computes the energy contributions to Δ G bind by averaging over multiple MD simulation frames using the package gmx_MMPBSA. 65 The entropic term may be estimated using normal mode analysis, interaction entropy (IE), 68 or a cumulant approximation to the second order of the exponential average (C2). 69 These entropy estimations have been implemented by gmx_MMPBSA 65 and may therefore be used by BindFlow. BindFlow employs a single-trajectory protocol, that is, all contributions from ligand, receptor, or receptor–ligand complex are computed from MD simulations of the complex. This approach reduces the computational cost and facilitates the cancellation of bonded terms in . However, the approach neglects contributions from structural changes of the ligand or receptor upon binding. 10 , 65 Free energy perturbation (FEP) Alchemical FEP estimates binding free energies by constructing a thermodynamic cycle. Such cycles involve decoupling the ligand from solvent ( Fig. 2b , Steps 1–3) and coupling it back in the binding site ( Fig. 2b , Steps 5–7) in the presence of restraints. The restraints maintain the ligand in the binding pocket, thereby improving sampling and unambiguously defining the reference standard state. BindFlow employs Boresch restraints. 38 The free energy for activating Boresch restraints in the solvent are calculated analytically ( Fig. 2b , Steps 3–4) and the free energy for removing the restraints in the receptor environment are derived numerically ( Fig. 2b , Steps 7–8). Free energy differences Δ G between two states (e.g., between states 2 and 1 in Fig. 2b ) are estimated from simulations across a series of intermediate alchemical states, specified by a parameter λ . Bindflow calculates Δ G values using either thermodynamic integration (TI) 64 or the multistate Bennett acceptance ratio (MBAR) 60 method. TI computes the Δ G values by integrating the mean work along the λ parameter: 64 where ⟨… ⟩ λ denotes the ensemble average at a given value λ and U ( λ ) is the potential energy. In practice, the integral is solved numerically using a finite number of λ points. MBAR estimates Δ G values using information from all alchemical states. MBAR requires solving the following set of equations self-consistently to obtain the reduced free energies f i : 60 Here, M is the number of λ -states, N j is the number of samples from the λ -state j , and u i ( x j,n ) is the reduced potential energy of the n th sample x j,n from the simulation of λ -state j evaluated with the energy function of λ -state i . Implementation General concept of BindFlow As illustrated in Fig. 3 , BindFlow streams the entire binding free energy pipeline, including building of the simulations systems, definition of force fields and MD parameters, careful multi-step equilibration routines, launching of productions simulations, and analysis. For FEP calculations, specifically, production simulations include the definition of Boresch restraints as well as the setup, equilibration, and launching of all λ -windows. The pipeline is organized into tasks with well-defined dependencies, which are deployed to the computing environment. The environment may be a desktop computer, a parallelized high-performance computer, or another distributed computing environment. Task scheduling is carried out by Snakemake, 70 a robust task manager that has been widely used for defining complex bioinformatics analyses, but has to our knowledge been less used in the MD community. Snakemake constructs a direct acyclic graph (DAG), which specifies the dependencies among Bindflow tasks to be executed asynchronously, thereby optimally using the available hardware. Bind-Flow uses GROMACS 71 as molecular dynamics engine. Download figure Open in new tab Figure 3. Task dependency graph for BindFlow. Each rectangle represents a task (or group of tasks), while arrows indicate dependencies between them. Some tasks are shared between methodologies. Purple boxes correspond to ligand-in-solvent FEP simulations, and green boxes to protein–ligand complex FEP simulations. The darkest blue boxes denote MM(PB/GB)SA simulations, which reuse some tasks from the protein–ligand complex FEP workflow. All three approaches share the initial system-building step (light blue) and converge at the final result-gathering stage to report the binding free energy, Δ G bind . BindFlow may compute Δ G bind from a protein PDB file together with a set of ligands provided as MOL files. However, for more complex systems including cofactors or membranes, additional input files may be provided. The pipeline is implemented into a single function calculate within the module runners . The function calculate takes general parameters such as the protein and ligand structures and, through the keyword argument global_config , more specific definitions such as the computational environment, or parameters for simulations analysis. Listing 1 presents a minimalistic Python script for running an FEP pipeline in an HPC environment, while Listing S1 presents a more extended example for running an MM(PB/GB)SA pipeline on a desktop computer. Parameters passed with global_config may be specified within Python ( Listing 1 ) or, more conveniently, using JSON or YAML files (Listing S1). While BindFlow implements default workflows, the user has full controls over equilibration and production settings, GROMACS parameters, and task deployment (Listing S7). Full documentation is provided online. Two classes for deployment are provided, namely for a desktop computer and for the SLURM queuing system, 72 however, alternative deployment managers may be added by the user via the abstract base class (ABC) Scheduler , for instance for using cloud computing services. Thus, BindFlow is highly user-friendly, yet allows extensive customization by the user. Download figure Open in new tab Listing 1: Minimal Python code for running an FEP campaign with BindFlow, as specified by calculation type . The protein specifies the path to the protein PDB file, and the ligands the list of ligand MOL files. The dictionary global_config specifies options for the computing environment. BindFlow has been forked from ABFE Workflow, 49 yet has been largely rewritten with the following aims in mind: options for extensive customization, efficient resource utilization, built-in support for several small-molecule force fields, flexibility with respect to user-provided force fields, comprehensive documentation, implementation of MM(PB/GB)SA, handling of cofactors, as well as support for membrane proteins and non-protein receptors. BindFlow is released as Open Source under the GPL-3.0 license. Input files and force fields For a typical workflow, the protein structure is provided as PDB file and a set of ligands as MOL files, where the ligand coordinates are aligned with the protein binding pocket. Optionally, a cofactor as MOL file or a membrane as PDB file may be provided. Alternatively, to enable the use of custom force fields, the user may provide GROMACS topology (TOP) and structure files (GRO). Running membrane protein systems requires additional care, as described in the online documentation. Examples for system component definitions are given in Listings S2, S3, S4, S5, and S6. By default, proteins are described with Amber99sb-ildn, 73 membranes with SLipids2020, 74 cofactors and ligands with OpenFF-2.0.0, 62 and water molecules with the TIP3P model. 75 Other protein force fields may be selected from the GROMACS distribution or from a user-specified path. Alternative water models may be chosen. The small-molecule force fields OpenFF, 62 GAFF, 61 and Espaloma 63 are natively supported by BindFlow via the TOFF 76 package. Results Validation set involving 139 receptor–ligand pairs and three small-molecule force fields We validated BindFlow by computing the binding affinity of 139 ligand/receptor pairs, for which high-quality experimental affinity data is available ( Fig. 1b ). The receptors included six soluble proteins (Cyclophilin D, MCL1, Thrombin, TYK2, PTP1B, P38), the transmembrane GPCR protein A2A, and the non-protein host–guest system SAMPL6-OA that has been used for a binding affinity prediction challenge. 78 Affinity data for P38, PTP1B, TYK2, Thrombin, and MCL1 have been widely utilized for validating RBFE 15 – 18 and ABFE calculations, 17 , 46 , 79 allowing us to compare BindFlow results with the literature (see below). These systems involve various challenges for binding affinity calculations. (i) Some ligands binding to MCL1, PTP1B, Thrombin, or SAMPL6-OA are charged (Figs. S25, S22, S24, S27). (ii) Ligand binding to PTP1B or P38 may proceed via slow induced-fit conformational transitions, suggesting that simulations may not sufficiency sample the apo state and, thereby, partly miss the free energy cost for the induced-fit transition. 46 (iii) Binding affinities may be influenced by cofactors: for P38, PTP1B, TYK2, and Thrombin, crystallographic water molecules may interact with the ligand ( Table 1 , 6 th column). In A2A, a sodium ion interacts with the ligand, suggesting that imprecise modeling of a salt bridge by the force field may add considerable uncertainty. (iv) Among the ten ligands binding to Cyclophilin D (Fig. S26), four comprise more than 30 heavy atoms, which could potentially increase sampling challenges, while (v) eight ligands contain exotic ring structures, posing challenges for the accurate force field representation. 5 , 49 , 50 Thus, our test sets involve highly diverse lig- and/receptor pairs that cover various common challenges during binding affinity calculations. All ligands are shown in Figs. S20–S27. View this table: View inline View popup Download powerpoint Table 1: Simulation systems used for validating BindFlow. We computed binding for these 139 ligand/receptor pairs in triplicates using FEP or MM(PB/GB)SA, except for the Cyclophilin D set, where five FEP replicates were used instead. For all protein targets, we computed binding affinities using three popular small-molecule force fields: GAFF-2.11, 61 OpenFF-2.0.0, 62 and Espaloma-0.3.1. 63 For the SAMPL6-OA system, only Espaloma-0.3.1 was used. 63 Cyclophilin D was described with the Amber99sbildn force field to allow comparison with a previous study. 5 , 49 , 50 Among our MM(PB/GB)SA calculations, we compared results from the Poisson-Boltzmann with results from the generalized Born model, and, for each solvation model, evaluated the effects from using either no entropy correction or using the IE or C2 entropy correction. Together, our validation sets comprise 1215 FEP calculations and 7254 MM(PB/GB)SA calculations. We quantified the agreement between calculated Δ G calc and experimental Δ G exp affinities using the Pearson ρ , Kendall τ , and Spearman r S correlation coefficient. Here, Pearson ρ quantifies the linear correlation between calculated and experimental values, accounting for magnitudes rather than just ranks, though it remains sensitive to outliers. In contrast, Kendall τ and Spearman r S quantify agreement in ranking and are far less sensitive to outliers. Whereas Spearman r S quantifies whether the ligand ranking according to Δ G exp agrees with the ligand ranking according to Δ G calc , Kendall τ quantifies the concordance between pairs of ligands. Upon comparing correlation coefficients from different receptors, it is critical to note that Pearson ρ is highly sensitive to the dynamic range among the ligands, whereas Kendall τ and Spearman r S are less sensitive to the dynamic range ( Table 1 , 7 th column). In this study, we strongly rely on Kendall τ as key quality measure for Δ G bind calculations since it has been shown to be a robust estimator, 80 while reporting Pearson ρ and Spearman r S as complementary measures. In addition to the correlation coefficients, we computed the root mean-squared error (RMSE), mean signed error (MSE), and mean unsigned error (MUE) of Δ G calc relative to Δ G exp (see SI methods). To focus on the global correlation between Δ G calc and Δ G exp across all ligand-receptor pairs, we furthermore subtracted from each set its corresponding MSE from Δ G calc to yield offset-corrected values and, respectively, offset-corrected RMSE (ocRMSE) and offset-corrected MUE (ocMUE, see Methods). Thereby, ocRMSE and ocMUE ignore systematic offsets between Δ G calc and Δ G exp , for instance owing to insufficiently sampling of the apo state or owing to a systematic bias in protein–ligand interactions as modeled by the force field. ocRMSE and ocMUE values may be compared with results from RBFE calculations since the latter are blind to such systematic offsets. Uncertainties of the statistical measures were derived by bootstrapping among the ligands for each set. Thus, critically, the errors bars for statistical measures are not only caused by limited sampling but furthermore stem from the limited number of 7 to 29 ligands per receptor. Comparison of free energy perturbation (FEP) results with experiments Figure 4 presents correlation plots between offset-corrected calculated affinities ob-tained with GAFF and experimental affinities Δ G exp . Values are either listed for seven protein targets individually (small panels) or accumulated from these seven targets (large panel). Insets report Pearson ρ , Kendall τ , ocMUE, and ocRMSE for each data set. Correlations plots prior to the offset correction are shown in Figs. S2 and S3. Correlation plots obtained with OpenFF or Espaloma are shown in Figs. S6 and S7. Notably, we removed a single extreme outlier given by the ligand lig 4 binding to Thrombin (Fig. S24) from the statistical quality measures to avoid a bias from a single receptor–ligand pair (see Fig. S3). Lig 4 comprises a cationic amidinium moiety forming a salt bridge to Asp-189 of Thrombin (Fig. S28). We speculate that the greatly overestimated affinity of lig 4 may be a consequence of overly stable amidinium–carboxylate bridge modeled by the three force fields. For most receptors and force fields, the affinities were overestimated, as shown by mostly negative MSE values between typically −4 and 0 kcal mol − 1 (Figs. S2, S3, see also Fig. 5d ). Such negative MSE values have been rationalized by insufficient sampling of the apo state, thereby partly missing the free energy cost of the induced-fit conformations transition upon ligand binding. 46 An extreme case was PTP1B with MSE values between −14 and −10 kcal mol − 1 , pointing to a major free energy cost for the induced-fit conformational transition that was not captured during FEP simulation, in line with the value of approximately −10 kcal mol − 1 reported by Chen et al. 46 . Download figure Open in new tab Figure 4. Offset-corrected calculated affinities versus experimental affinities Δ G exp from FEP with GAFF-2.11. Results shown for individual set (small panels, see labels for protein name and number of ligands N ) and collected from all sets (large panel). Insets show Pearson ρ , Kendall τ , osMUE, and ocRMSE (last two in kcal mol − 1 ) for each data set with its corresponding 68 % confident interval. Colors of dots indicate the absolute deviation between Δ G calc and Δ G exp (color bar). Dark and light gray diagonal regions indicate 1 or 2 kcal mol − 1 deviations, respectively. Dashed lines are linear fits shown to guide the eye. A single outlier (lig 4, Thrombin) has been removed. Error bars show uncertainties obtained via three independent replicates. Figure style following Ref. 79. Download figure Open in new tab Figure 5. ( a ) Pearson ρ , ( b ) Kendall τ , ( c ) Spearman r S , and ( d ) deviations between calculated and experimental binding free energies Δ G calc − Δ G exp from FEP using Espaloma-0.3.1 (green), GAFF-2.11 (orange), and OpenFF-2.0.0 (blue). Results are shown for seven different data sets (labels at abscissa). Error bars represent the 68 % confident interval obtained from bootstrapping among the ligands from each set. For the Δ G calc − Δ G exp , dark and light gray horizontal regions indicate 1 or 2 kcal mol − 1 deviations, respectively. Box plots present the median (50 th percentile) as a line within the box, and the lower and upper box edges correspond to the first (25 th percentile) and third (75 th percentile) quartiles. Whiskers extend to the smallest and largest data points within 1.5 × interquartile range (IQR = Q3 − Q1), while outliers beyond this range are shown as circles. After removing the MSE form Δ G calc to obtain , statistical quality measures for the agreement with experiments still greatly depended on the protein target, with ocMUE values far below 1 kcal mol − 1 for Cyclophilin D up to approximately 2 kcal mol − 1 for PTP1B. Excellent agreement was found for SAMPL6-OA (Fig. S7). These variations were likewise reflected by Kendall τ which spanned 0.82 for Cyclophilin D with Espaloma down to poor τ for Thrombin. Pearson ρ and ocRMSE values reported similar trends. The wide range of correlations obtained with different receptor and different force fields are summarized by Fig. 5 . Comparison of FEP results from BindFlow with the literature Ross et al. 81 evaluated the reproducibility of experimental RBFE and reported an overall RMSE (weighted by number of ligands on each set) of 0.91 kcal mol − 1 (95% CI: [0.83, 1.11]) and a weighted average Kendall τ = 0.71 (95% CI: [0.65, 0.74]). Thus, upon comparing Δ G calc with Δ G exp , the considerable intrinsic uncertainty of the Δ G exp must be kept in mind. Accordingly, an RMSE of ∼ 1 kcal mol − 1 and Kendall τ of ∼ 0.7 between Δ G calc and Δ G exp would indicate excellent agreement. In addition, Ross et al. 81 validated the FEP+ software for RBFE calculations in combination with the OPLS4 force field 82 against protein–ligand sets different from the ones used here. Beyond the differences in validation sets, their study employed RBFEs, whereas we performed ABFEs, which are generally more challenging. For these reasons, a direct, quantitative comparison of RMSE and Kendall τ is not meaningful. Nevertheless, placing the results side by side helps to contextualize the current state of the art in FEP calculations. FEP+ achieved an RMSE of 1.25 (95% CI: [1.17, 1.33]) and a Kendall τ of 0.51 (95% CI: [0.48, 0.55]). With Espaloma, GAFF, and OpenFF, we obtained slightly higher ocRMSE values, which may indicate that the closed-source OPLS4 force field 82 performs somewhat better than the currently freely available force fields. However, the Kendall τ values obtained here were even slightly higher than those from FEP+, suggesting that our ABFE calculations are competitive in terms of ligand ranking. Next, we systematically compared ABFE results from BindFlow with FEP+ by restricting the analysis to ligands that were present in both the dataset of Chen et al. 46 and our BindFlow simulations, resulting in 89 ligands across the P38, PTP1B, TYK2, and MCL1 sets. Both BindFlow and FEP+ systematically overestimated affinities for PTP1B (Fig. S11, top row). Thus, to enable a quantitative comparison of ABFE values, we removed PTP1B from the datasets and focused the comparison on the remaining 67 ligands. FEP+ yielded systematically stronger binding affinities on these three sets compared to BindFlow, as quantified by the MSE of for FEP+, versus , and for Espaloma, GAFF, and OpenFF, respectively (Fig. S11, middle row). Keeping the PTP1B set and removing the constant offset for each set by subtracting the set-specific MSE from Δ G calc , FEP+ yielded a considerably lower ocRMSE and ocMUE values compared to BindFlow when used with any of the three force fields. FEP+ also yielded higher correlation coefficients (Fig. S11, bottom row). Thus, in the light of these data, FEP+ may achieve superior ranking at the cost of overestimating the absolute binding affinity. We speculate that the superior ranking by FEP+ is primarily a consequence of using the OPLS4 force field, although alternative sampling algorithms may also play a role. Within this subset of ligands, GAFF achieved better agreement than Espaloma and OpenFF with experimental data, with , and (Fig. S11, bottom row). Deflorian et al. 13 reported RBFE calculations for the A2A system using the FEP+ software. We converted their RBFE values into ABFEs by referencing the experimental binding affinity of ligand 4k and correcting for the MSE. From their data, we ob-tained and . Evidently, since all ligands were similarly composed of three aromatic rings (Fig. S21), RBFE required relatively small perturbations, rationalizing the high τ and exceptionally low ocRMSE. Although our ABFE calculations involved by far larger perturbations, we obtained decent agreement (Fig. S14) with GAFF ( , Fig. 4 ) and with OpenFF ( , Fig. S6). Only with Espaloma, we obtained a rather poor τ of only (Fig. S7), which might indicate inaccurate modeling of the ligand interac-tions with the Na + cofactor inside the A2A binding pocket. The Cyclophilin D set was first introduced by Alibay et al. 5 for ABFE benchmarking and was recently employed to validate ABFE workflow 49 and A3FE. 50 Excellent agreement with experiments was found across the three pipelines —ABFE workflow, A3FE, and BindFlow—, underscoring that the Cyclophilin D set was the least challenging receptor–ligand set considered in this study (Figs. S12, 5). The P38 and TYK2 systems have been used by several authors for validating binding free energy calculations using different methodologies. 17 , 46 , 48 , 79 Figure S13 compares results from BindFlow with previous studies 17 , 46 , 48 , 79 . Here, Lin et al. 17 and Chen et al. 46 achieved the best ranking of ligands, however, at the cost of overestimating the absolute binding free energies as shown by large negative MSE. Among our results, GAFF achieved comparable ranking at a less negative MSE. Taken together, the reasonable agreement with previous data validates BindFlow’s FEP pipeline and demonstrates competitive agreement with experimental data for a wide range of ligands and receptors. Comparison of MM(PB/GB)SA with experiments MM(PB/GB)SA provides a computationally efficient alternative to FEP. Specifically, the MM(PB/GB)SA pipeline used here employed only ∼3 ns of MD simulations; corresponding to a 180-fold lower computing cost compared to the FEP pipeline used here. However, whereas the MMGBSA calculation from the MD frames is highly efficient, the MMPBSA calculation may take considerable additional computing time owing to the cost of Poisson-Boltzmann calculations. As expected from previous studies, 10 absolute binding free energies obtained with MMGBSA (Figs. S4, S5, 6d, S15d, S16d) or MMPBSA (Figs. S17d, S18d, and S19d) revealed poor agreement with experimental values. For instance, MMGBSA without entropy correction overestimated the absolute binding affinities to protein targets by approximately 25 to 75 kcal mol − 1 ( Fig. 6d ). Even upon correcting the binding affinity with a target-specific offset, ocRMSE values remained large (Figs. S8, S9, S10). Download figure Open in new tab Figure 6. Pearson ρ , Kendall τ , Spearman r S , and deviations between calculated and experimental binding free energies (Δ G calc − Δ G exp ) from MMGBSA results without entropy correction. Presentation style according to Fig. 5 . However, MMGBSA without entropy correction achieved overall excellent Pearson, Kendall, and Spearman correlation coefficients for A2A, Cyclophilin D, and SAMPL6-OA, while reaching still good correlation coefficients for PTP1B, TYK2, and MCL1 ( Fig. 6a–c ). In contrast, correlation coefficients for P38 or Thrombin were rather poor. The three force fields GAFF, OpenFF, and Espaloma achieved similar ranking (Figs. S10, S8, S9 and also 6a–c). Thus, while MMGBSA is not suitable for obtaining absolute binding free energies ( Fig. 6d ), it may provide for certain receptors reasonable ranking among ligands and, thereby, serve as a useful tool for screening large datasets of potential binders. Using the Poisson-Boltzmann instead of the generalized Born solvation model had only moderate effects on the correlation coefficients (compare Fig. S17 with Fig. 6 ), suggesting that the computationally more efficient generalized Born model yields a good starting point for computational studies. Using the IE entropy correction had only a small effects on the correlation coefficients (compare Fig. S19 with S16), whereas the C2 entropy correction decreased the correlation coefficients (compare Fig. S18 with S15). Hence, in line with previous findings, 83 entropy corrections should be used with care or tested for the receptor of interest, because entropy corrections may deteriorate the ranking by MM(PB/GB)SA. Ligand ranking by FEP versus MM(PB/GB)SA and by different force fields Figure 7 summarizes the overall performance of FEP versus MM(PB/GB)SA among the three force fields, as collected from all ligands–receptor pairs, thus including easy receptors such as Cyclophilin D and challenging receptors such as Thrombin (see also Table S1). Here, the matrices in Figs. 7b/d visualize whether one combination of method/force field, such FEP/Espaloma, outperforms another combination, such as MMGBSA/GAFF (see Methods for the applied significance test). The following conclusions are drawn from the analysis: (i) In terms of ocMUE, FEP by far outperformed MMGBSA ( Figs. 7d ). These findings suggest that errors caused by major approximations underlying MMGBSA remain after correcting Δ G calc by the MSE. (ii) In terms of ranking by Kendall τ , FEP overall outperformed MMG-BSA ( Figs. 7a/b ). An exception was given by FEP/Espaloma that showed only insignificant differences relative to MMGBSA ( Figs. 7b ). (iii) For a given method (FEP or MMGBSA), the three force field showed no statistically significant differences ( Figs. 7b/d ), except among the poor ocMUE by MMGBSA. Thus, simulations with additional ligands will be required to test whether the slightly lower τ and slightly larger ocMUE obtained with Espaloma are significant or whether they were caused by the limited set of receptor–ligand pairs considered in this study. Download figure Open in new tab Figure 7. Summary of agreement between FEP or MMGBSA calculations and experimental data quantified by ( a / b ) Kendall τ and ( b / c ) offset-corrected mean unsigned error (ocMUE) in kcal mol − 1 , summarized from all receptors. Error bars show 68 % confident intervals obtained from bootstrapping among the ligands from all sets. Pairwise significance difference matrices are for ( b ) τ and ( d ) ocMUE. An arrow indicates a statistically significant difference pointing to the combination of method/force field with the better statistical metric. Equal symbols indicate an insignificant difference. The set SAMPL6-OA is excluded to ensure consistence in the analysis across force fields. The extreme outlier lig 4 from the Thrombin set was excluded. Discussion We introduced BindFlow, a free and user-friendly software for absolute binding free energy (ABFE) calculations. BindFlow offers fully automated pipelines for FEP and MM(GB/PB)SA, requiring minimal user intervention while providing extensive configuration options for advanced users. Comprehensive online tutorials and documentation support accessibility. The software is developed according to best coding practices, ensuring transparency, reproducibility, and a solid foundation for future community-driven development. Beyond the benchmark simulations presented here, BindFlow has recently been applied successfully for combining ∼60.000 MM(GB/PB)SA calculations with Bayesian active learning 41 as well as for interpreting X-ray crystallographic data from the antiviral drug tecoviri-mat binding to the phospholipase F13 from monkeypox virus. 43 These use cases demonstrate that BindFlow is ready for production use. We validated BindFlow by computing Δ G bind values for 139 ligands binding to a one of eight different receptors using FEP together with MBAR or using MM(GB/PB)SA. Overall, we found that BindFlow yields competitive Δ G bind compared to the standards in the field. As expected from previous studies, the accuracy of Δ G bind values and ligand ranking greatly depend on the target and on the chemical moieties of the specific ligand. For the 139 receptor–ligand pairs, in terms of accuracy of Δ G bind predictions quantified by ocMUE or ocRMSE, FEP by far outperformed MM(GB/PB)SA. In terms of ligand ranking, FEP still outperformed MMGBSA significantly; however, for several receptors studied here, ranking by MMGBSA was remarkably good. Only for P38 and Thrombin, MMGBSA achieved poor ranking. Thus, our study underscores that MMGBSA may provide a computationally inexpensive method for pre-screening large sets of ligands prior to validation by the more expensive FEP. However, future studies will need to establish whether favorable correlations between MM(GB/PB)SA and experiments, as observed here, also hold for non-congeneric sets of ligands. BindFlow will be useful to systematically test and improve small molecule force fields. For instance, we spotted that the affinity of lig 4 binding to Thrombin was greatly overestimated by 7 to 13 kcal mol − 1 , possibly caused by an overly strong amidinium–aspartate salt bridge (Fig. S28). This interpretation is supported by the fact that larger partial charges used by Espaloma correlated with an even stronger computed affinity. Since amidinium moieties are quite common in drugs, these findings point towards options for further force field refinements. Future implementations may further enhance the capabilities of BindFlow. (i) Enhanced sampling techniques, such as Hamiltonian replica exchange, 84 – 86 may accelerate convergence and, thereby, reduce the computational cost. Enhanced sampling may, specifically, improve the sampling of the apo state and lead to smaller MSE values. (ii) The adaptive allocation of computational resources may focus sampling to λ windows with increased sampling challenges. 50 (iii) Charged ligand systems (e.g., MCL1, Thrombin, and PTP1B) showed satisfactory accuracy in our validations. However, finite-size effects in charge-imbalanced systems, primarily due to the use of particle-mesh Ewald (PME) method, 87 , 88 remain a concern. 89 , 90 Charge corrections 91 , 92 such as the co-alchemical ion approach 12 represent viable solutions to address this issue. (iv) Water exchange between the binding site and bulk solvent represents a slow sampling process that is difficult to capture during the short MD simulations performed in FEP. 12 Monte Carlo-based water swap methods in a grand canonical ensemble, 93 , 94 methods based on inhomogeneous fluid solvation theory, 95 or specialized hydration-shell generators such as SOLVATE 96 , 97 have shown promise in treating the explicit water. Future versions of BindFlow may include such methodologies. Conclusions We presented BindFlow, an open-source and user-friendly pipeline for absolute binding free energy (ABFE) calculations using FEP and MM(PB/GB)SA. Validation on 139 receptor– ligand pairs showed that BindFlow achieves predictive performance comparable to standards in the field, while offering full automation and extensive customization. FEP achieved more accurate Δ G bind predictions and ligand ranking compared to MM(PB/GB)SA. However, MM(PB/GB)SA provided remarkable ligand rankings for several systems at a fraction of the computational cost, highlighting its value for large-scale computational pre-screening campaigns. Thus, we anticipated that the combination of MM(PB/GB)SA and FEP will be powerful for balancing efficiency and accuracy in future BindFlow applications. By combining automation, flexibility, accessibility, and documentation, BindFlow lowers the entry barrier for routine ABFE calculations and provides a platform for modern drug discovery or for systematically improving small-molecule force fields. Methods Building simulation systems Simulations followed largely the current default settings of BindFlow. The setup protocols are fully automated, yet may be easily adapted by the user via the global input keyword, which is a Python dictionary. This dictionary may be conveniently defined with a YAML file (Listing S7). The simulations presented here were set up as follows. The task Build Simulation System ( Fig. 3 ) assembles the simulation systems and defines force fields. Starting structures were taken from previous studies. 5 , 13 , 30 , 49 , 77 , 78 Soluble proteins, SAMPLE6-OA, and ligands were placed in a dodecahedral simulation box, keeping a distance of 1.5 nm between solute and box boundary. The systems were solvated by explicit water with the GROMACS solvate module and approximately 150 mM of NaCl were added, thereby neutralizing the system. A minimum distance of 1 nm between ions and nonsolvent molecules was used. The structure of the membrane protein A2A was taken from Ref. 13. BindFlow does not build simulation systems of membrane proteins. Thus, the A2A system was built using CHARMM-GUI 98 by placing the protein in a simulation box of a hexagonal prism and by embedding it in a membrane of 172 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) lipids. Ligand interactions were described with Espaloma-0.3.1, 63 GAFF-2.11, 61 or OpenFF-2.0.0. 62 Ligand topologies were generated within BindFlow using the TOFF Python library. 76 Cyclophilin D was described with Amber99sb-ildn 73 as used by previous studies. 5 , 49 , 50 All other proteins were described with Amber14sb, 99 using OpenMM 100 and ParmEd 101 for generating force field topologies. The host molecule octa acid of the SAMPL6-OA system was described with Espaloma-0.3.1 63 and the topologies were generated with TOFF. 76 Water was modeled with the TIP3P model 75 and ion parameters were taken from the definitions provided with Amber99sb-ildn. 73 Simulation parameters The tasks Setup Complex Equilibration, Setup Ligand Equilibration, Setup Complex FEP, Setup Ligand FEP , and Setup Complex MM(PB/GB)SA ( Fig. 3 ) define MD parameters and set up directory structures for the initial equilibration, FEP or MM(GB/PB)SA simulations. All simulations were carried out with GROMACS, version 2022.4. 71 The geometry of water acting as cofactor was constrained with LINCS, 102 while the geometry of all other water molecules were constrained with SETTLE. 103 All other bonds were constrained with LINCS 102 if not stated otherwise. Hydrogen mass repartitioning was used with a mass repartition factor of 2.5. During production simulations a time step of 4 fs was used. Dispersive interactions and short-range repulsion were modeled using a Lennard-Jones potential with a 1 nm cutoff. Electrostatic interactions were calculated with the particle-mesh Ewald 87 , 88 (PME) method, applying a real-space cutoff of 1 nm. The temperature was maintained at 298.15 K. For the initial equilibration of membrane protein–ligand complexes, the temperature is maintained using velocity rescaling ( τ t = 1 ps). Here, the protein and the Na + cofactor were coupled to the same thermostat. During all other simulations, the temperature was maintained using Langevin dynamics ( τ t = 2 ps). The pressure was maintained at 1 bar for membrane protein–ligand complexes using a semi-isotropic stochastic cell rescaling 104 with a time constant τ p = 5 ps during the initial equilibration. The same barostat was used during the perturbation phase for membrane protein–ligand complexes, but the pressure was maintained at 1 atm (1.01325 bar) and τ p = 2 ps was used. For all other simulations the pressure was controlled at 1 atm with the isotropic Berendsen barostat during equilibration ( τ t = 1 ps) 105 and the Parrinello-Rahman barostat during production ( τ p = 2 ps). 106 Initial multi-step equilibration protocol After minimizing the energy with the steepest-descent algorithm, each system was equilibrated with a multi-step protocol (see SI methods). FEP calculations Δ G bind from FEP was computed according the thermodynamic cycle shown in Fig. 2b . Coulomb interactions of the ligand in solvent were decoupled over 11 λ points ( Fig. 2b , transition 1→2), followed by decoupling of Lennard-Jones interactions over 11 λ points (transition 2→3). Both, inter- and intramolecular interactions were decoupled. A 10 ns equilibrium simulation together with the software MDRestraintsGenerator 107 were used to obtain the optimal Boresch restraints, 38 the free energy cost for introducing restraints to the ligand in solvent (transition 3→4), as well as the initial frame for the complex decoupling simulations. Inter- and intramolecular Lennard-Jones interactions of the ligand in the receptor were activated over 21 λ points (transition 5→6), followed by 11 λ points for activating the Coulomb interactions (transition 6→7). Finally, Boresch restraints were removed for the ligand in the receptor over 11 λ points (transition 7→8). For each λ window, the system was energy minimized and equilibrated with a three-step protocol (see SI methods). Each window was simulated for 10 ns using a stochastic dynamics integrator. During simulations that decouple Lennard-Jones interactions, a soft-core potential ( α = 0.5, σ = 0.3, power = 1) was used. Free energy differences were computed with MBAR. 60 As a control, free energies were additionally computed using thermodynamic integration (TI). 64 BindFlow reports a warning if results from MBAR and TI differ by more than 0.5 kcal mol − 1 , as such differences may indicate poor convergence or simulation instabilities. Alchemlyb-2.0.0 was used for MBAR and TI evaluations. 108 The binding free energy is computed via ( Fig. 2B ): Here, and denote free energies for decoupling (dcpl) Coulomb and Lennard-Jones interactions of the ligand in solvent (solv), respectively. denotes the free energy for turning on Boresh restraints (rest) for the ligand in solvent, which is computed nalytically. and denote the free energy for activating (coupling, cpl) Lennard-Jones and Coulomb interactions for the ligand in the receptor (rcp), respectively. denotes the free energy cost for turning off the Boresch restraints for the ligand in the receptor. BindFlow reports the statistical error obtained from MBAR 60 together with Gaussian error propagation. However, owing to autocorrelations, this value greatly underestimates the true uncertainty of Δ G bind . Thus, we carried out the whole pipeline including setup, equilibration, and FEP in three (five for Cyclophilin D 5 ) independent replicates and report the respective standard error of the mean as error bars in the correlation plots. MM(PB/GB)SA calculations The initial multi-step equilibration for MM(PB/GB)SA calculation was carried out as described in the SI methods. Following Su et al. 83 , MM(PB/GB)SA values were computed from multiple short simulations, rather than from a single long simulation. Accordingly, twenty starting frames were taken from a 950 ps equilibrium simulation, using one frame every 50 ps. From each frame, a 100 ps simulation was carried out writing a frame every 5 ps. Here, we follow Genheden and Ryde 10 , who recommended using an output frequency of 1 to 10 ps. MM(PB/GB)SA values were computed with the gmx_MMPBSA software, 65 yielding 20 samples of the binding affinity. Cofactors were defined as part of the receptor for the MM(PB/GB)SA calculation. For each complex, the pipeline was run in three independent replicates. Binding affinities reported here were computed by averaging over 60 samples from the three replicates with 20 samples per replica, and the error was derived as SEM and reported as error bars in correlation plots. PB and GB models were used, each without entropy correction or using the IE or C2 entropy corrections. The main text reports results from the GB model without entropy correction, while all other results are provided in the SI. Code and data availability statement BindFlow is free software published under the GPL-3.0 license. The code is currently hosted at https://github.com/ale94mleon/BindFlow . Documentation is available at https://bindflow.readthedocs.io . Scripts and input files required for reproducing our results and analysis are accessible at: https://github.com/ale94mleon/bindflow-api-paper-si . Supporting Information Available Supporting Information provides Supporting Methods, Listings S1–S7, Figures S1–S28 and Table S1. TOC Graphic Download figure Open in new tab Acknowledgement We thank the authors of ABFE Workflow for publishing the software under a free license, and in particular Benjamin Ries and Aniket Magarkar for helpful discussions. We thank Mario Sergio Valdés-Trasanco for support with the gmx_MMPBSA software. This research was funded by the European Union’s Horizon 2020 research and innovation program under Marie Skłodowska Curie Grant 860592. Funder Information Declared European Union, https://ror.org/019w4f821 , Marie Sklodowska Curie Grant 860592 Footnotes ↵ * E-mail: alejandro.martinezleon{at}uni-saarland.de ; jochen.hub{at}uni-saarland.de https://github.com/ale94mleon/BindFlow https://bindflow.readthedocs.io References (1). ↵ Paul , S. M. ; Mytelka , D. S. ; Dunwiddie , C. T. ; Persinger , C. C. ; Munos , B. H. ; Lindborg , S. R. ; Schacht , A. L. How to improve R&D productivity: the pharmaceutical industry’s grand challenge . Nature Reviews Drug Discovery 2010 , 9 , 203 – 214 . OpenUrl CrossRef PubMed Web of Science (2). ↵ DiMasi , J. A. ; Grabowski , H. G. ; Hansen , R. W. Innovation in the pharmaceutical industry: new estimates of R&D costs . Journal of Health Economics 2016 , 47 , 20 – 33 . OpenUrl CrossRef PubMed (3). ↵ Cournia , Z. ; Chipot , C. ; Roux , B. ; York , D. M. ; Sherman , W. Free energy methods in drug discovery: Current state and future directions ; ACS Publications , 2021 ; pp 1 – 38 . (4). ↵ Li , Z. ; Li , X. ; Huang , Y.-Y. ; Wu , Y. ; Liu , R. ; Zhou , L. ; Lin , Y. ; Wu , D. ; Zhang , L. ; Liu , H. ; others Identify potent SARS-CoV-2 main protease inhibitors via accelerated free energy perturbation-based virtual screening of existing drugs . Proceedings of the National Academy of Sciences 2020 , 117 , 27381 – 27387 . OpenUrl Abstract / FREE Full Text (5). ↵ Alibay , I. ; Magarkar , A. ; Seeliger , D. ; Biggin , P. C. Evaluating the use of absolute binding free energy in the fragment optimisation process . Communications Chemistry 2022 , 5 , 1 – 13 . OpenUrl PubMed (6). ↵ Qian , R. ; Xue , J. ; Xu , Y. ; Huang , J. Alchemical Transformations and Beyond: Recent Advances and Real-World Applications of Free Energy Calculations in Drug Discovery . Journal of Chemical Information and Modeling 2024 , 64 , 7214 – 7237 . OpenUrl CrossRef PubMed (7). ↵ Zwanzig , R. W. High-temperature equation of state by a perturbation method. I. Nonpolar gases . The Journal of Chemical Physics 1954 , 22 , 1420 – 1426 . OpenUrl CrossRef PubMed Web of Science (8). ↵ Zwanzig , R. W. High-Temperature Equation of State by a Perturbation Method. I. Nonpolar Gases . Journal of Chemical Physics 1954 , 22 , 1420 – 1426 . OpenUrl CrossRef PubMed Web of Science (9). ↵ Kollman , P. A. ; Massova , I. ; Reyes , C. ; Kuhn , B. ; Huo , S. ; Chong , L. ; Lee , M. ; Lee , T. ; Duan , Y. ; Wang , W. ; Donini , O. ; Cieplak , P. ; Srinivasan , J. ; Case , D. A. ; Cheatham , T. E. Calculating structures and free energies of complex molecules: Combining molecular mechanics and continuum models . Accounts of Chemical Research 2000 , 33 , 889 – 897 . OpenUrl CrossRef PubMed Web of Science (10). ↵ Genheden , S. ; Ryde , U. The MM/PBSA and MM/GBSA methods to estimate ligandbinding affinities . Expert Opinion on Drug Discovery 2015 , 10 , 449 – 461 . OpenUrl CrossRef PubMed (11). ↵ Muegge , I. ; Hu , Y. Recent Advances in Alchemical Binding Free Energy Calculations for Drug Discovery . ACS Medicinal Chemistry Letters 2023 , 14 , 244 – 250 . OpenUrl PubMed (12). ↵ Mey , A. S. ; Allen , B. K. ; Macdonald , H. E. B. ; Chodera , J. D. ; Hahn , D. F. ; Kuhn , M. ; Michel , J. ; Mobley , D. L. ; Naden , L. N. ; Prasad , S. ; Rizzi , A. ; Scheen , J. ; Shirts , M. R. ; Tresadern , G. ; Xu , H. Best Practices for Alchemical Free Energy Calculations [Article v1.0] . Living Journal of Computational Molecular Science 2020 , 2 , 1 – 52 . OpenUrl (13). ↵ Deflorian , F. ; Perez-Benito , L. ; Lenselink , E. B. ; Congreve , M. ; van Vlijmen , H. W. ; Mason , J. S. ; Graaf , C. d. ; Tresadern , G. Accurate prediction of GPCR ligand binding affinity with free energy perturbation . Journal of Chemical Information and Modeling 2020 , 60 , 5563 – 5579 . OpenUrl PubMed (14). ↵ Aldeghi , M. ; Gapsys , V. ; de Groot , B. L. Accurate estimation of ligand binding affinity changes upon protein mutation . ACS Central Science 2018 , 4 , 1708 – 1718 . OpenUrl PubMed (15). ↵ Wang , L. et al. Accurate and reliable prediction of relative ligand binding potency in prospective drug discovery by way of a modern free-energy calculation protocol and force field . Journal of the American Chemical Society 2015 , 137 , 2695 – 2703 . OpenUrl CrossRef PubMed (16). Kuhn , M. ; Firth-Clark , S. ; Tosco , P. ; Mey , A. S. ; Mackey , M. ; Michel , J. Assessment of binding affinity via alchemical free-energy calculations . Journal of Chemical Information and Modeling 2020 , 60 , 3120 – 3130 . OpenUrl CrossRef PubMed (17). ↵ Lin , Z. ; Zou , J. ; Liu , S. ; Peng , C. ; Li , Z. ; Wan , X. ; Fang , D. ; Yin , J. ; Gobbo , G. ; Chen , Y. ; others A cloud computing platform for scalable relative and absolute binding free energy predictions: New opportunities and challenges for drug discovery . Journal of Chemical Information and Modeling 2021 , 61 , 2720 – 2732 . OpenUrl CrossRef PubMed (18). ↵ Moore , J. H. ; Margreitter , C. ; Janet , J. P. ; Engkvist , O. ; de Groot , B. L. ; Gapsys , V. Automated relative binding free energy calculations from SMILES to ΔΔG . Communications Chemistry 2023 , 6 , 1 – 12 . OpenUrl PubMed (19). Martins , L. C. ; Cino , E. A. ; Ferreira , R. S. PyAutoFEP: An Automated Free Energy Perturbation Workflow for GROMACS Integrating Enhanced Sampling Methods . Journal of Chemical Theory and Computation 2021 , 17 , 4262 – 4273 . OpenUrl (20). Zavitsanou , S. ; Tsengenes , A. ; Papadourakis , M. ; Amendola , G. ; Chatzigoulas , A. ; Dellis , D. ; Cosconati , S. ; Cournia , Z. FEPrepare: A Web-Based Tool for Automating the Setup of Relative Binding Free Energy Calculations . Journal of Chemical Information and Modeling 2021 , 61 , 4131 – 4138 . OpenUrl CrossRef PubMed (21). Gowers , R. J. ; Alibay , I. ; Swenson , D. W. ; Henry , M. M. ; Ries , B. ; Baumann , H. M. ; Eastwood , J. R. B. ; Mitchell , J. A. ; Dotson , D. ; Horton , J. T. ; Thompson , M. The Open Free Energy library . 2024 ; https://github.com/openfreeEnergy/openfe . (22). Lee , T.-S. ; Allen , B. K. ; Giese , T. J. ; Guo , Z. ; Li , P. ; Lin , C. ; McGee Jr , T. D. ; Pearlman , D. A. ; Radak , B. K. ; Tao , Y. ; others Alchemical binding free energy calculations in AMBER20: Advances and best practices for drug discovery . Journal of Chemical Information and Modeling 2020 , 60 , 5595 – 5623 . OpenUrl CrossRef PubMed (23). ↵ Heinzelmann , G. ; Gilson , M. K. Automation of absolute protein-ligand binding free energy calculations for docking refinement and compound evaluation . Scientific Reports 2021 , 11 , 1 – 18 . OpenUrl PubMed (24). Homeyer , N. ; Gohlke , H. FEW: A workflow tool for free energy calculations of ligand binding . Journal of Computational Chemistry 2013 , 34 , 965 – 973 . OpenUrl CrossRef PubMed (25). Jespers , W. ; Esguerra , M. ; Åqvist , J. ; Gutiérrez-De-Terán , H. Qligfep: An automated workflow for small molecule free energy calculations in Q . Journal of Cheminformatics 2019 , 11 , 1 – 16 . OpenUrl PubMed (26). Matricon , P. ; Ranganathan , A. ; Warnick , E. ; Gao , Z. G. ; Rudling , A. ; Lambertucci , C. ; Marucci , G. ; Ezzati , A. ; Jaiteh , M. ; Ben , D. D. ; Jacobson , K. A. ; Carlsson , J. Fragment optimization for GPCRs by molecular dynamics free energy calculations: Probing druggable subpockets of the A2A adenosine receptor binding site . Scientific Reports 2017 , 7 , 1 – 12 . OpenUrl CrossRef PubMed (27). Lenselink , E. B. et al. Predicting Binding Affinities for GPCR Ligands Using Free-Energy Perturbation . ACS Omega 2016 , 1 , 293 – 304 . OpenUrl PubMed (28). Hedges , L. ; Mey , A. ; Laughton , C. ; Gervasio , F. ; Mulholland , A. ; Woods , C. ; Michel , J. BioSimSpace: An interoperable Python framework for biomolecular simulation . Journal of Open Source Software 2019 , 4 , 1831 . OpenUrl (29). ↵ Kim , S. ; Oshima , H. ; Zhang , H. ; Kern , N. R. ; Re , S. ; Lee , J. ; Roux , B. ; Sugita , Y. ; Jiang , W. ; Im , W. CHARMM-GUI Free Energy Calculator for Absolute and Relative Ligand Solvation and Binding Free Energy Simulations . Journal of Chemical Theory and Computation 2020 , 16 , 7207 – 7218 . OpenUrl (30). ↵ Hahn , D. F. ; Bayly , C. I. ; Boby , M. L. ; Macdonald , H. E. B. ; Chodera , J. D. ; Gapsys , V. ; Mey , A. S. J. S. ; Mobley , D. L. ; Benito , L. P. ; Schindler , C. E. M. ; Tresadern , G. ; Warren , G. L. Best Practices for Constructing, Preparing, and Evaluating Protein-Ligand Binding Affinity Benchmarks [Article v1.0] . Living Journal of Computational Molecular Science 2022 , 4 , 1 – 35 . OpenUrl (31). ↵ Zhang , C. H. et al. Potent Noncovalent Inhibitors of the Main Protease of SARSCoV-2 from Molecular Sculpting of the Drug Perampanel Guided by Free Energy Perturbation Calculations . ACS Central Science 2021 , 7 , 467 – 475 . OpenUrl PubMed (32). Tresadern , G. ; Velter , I. ; Trabanco , A. A. ; Keybus , F. V. D. ; MacDonald , G. J. ; Somers , M. V. ; Vanhoof , G. ; Leonard , P. M. ; Lamers , M. B. ; Roosbroeck , Y. E. V. ; Buijnsters , P. J. [1,2,4]Triazolo[1,5-a]pyrimidine Phosphodiesterase 2A Inhibitors: Structure and Free-Energy Perturbation-Guided Exploration . Journal of Medicinal Chemistry 2020 , 63 , 12887 – 12910 . OpenUrl PubMed (33). O’Donovan , D. H. ; Gregson , C. ; Packer , M. J. ; Greenwood , R. ; Pike , K. G. ; Kawatkar , S. ; Bloecher , A. ; Robinson , J. ; Read , J. ; Code , E. ; others Free energy perturbation in the design of EED ligands as inhibitors of polycomb repressive complex 2 (PRC2) methyltransferase . Bioorganic & Medicinal Chemistry Letters 2021 , 39 , 127904 . OpenUrl PubMed (34). ↵ Papadourakis , M. ; Sinenka , H. ; Matricon , P. ; Hénin , J. ; Brannigan , G. ; Pérez-Benito , L. ; Pande , V. ; van Vlijmen , H. ; de Graaf , C. ; Deflorian , F. ; others Alchemical Free Energy Calculations on Membrane-Associated Proteins . Journal of Chemical Theory and Computation 2023 , 19 , 7437 – 7458 . OpenUrl (35). Majellaro , M. ; Jespers , W. ; Crespo , A. ; Núñez , M. J. ; Novio , S. ; Azuaje , J. ; Prieto-Diaz , R. ; Gioé , C. ; Alispahic , B. ; Brea , J. ; others 3, 4-Dihydropyrimidin-2 (1 H)-ones as antagonists of the human A2B adenosine receptor: optimization, structure–activity relationship studies, and enantiospecific recognition . Journal of Medicinal Chemistry 2020 , 64 , 458 – 480 . OpenUrl PubMed (36). Keränen , H. ; de Terán , H. G. ; Åqvist , J. Structural and energetic effects of A 2A adenosine receptor mutations on agonist and antagonist binding . PLoS ONE 2014 , 9 , 1 – 12 . OpenUrl CrossRef PubMed (37). ↵ Xue , B. et al. Development and Comprehensive Benchmark of a High-Quality AMBER-Consistent Small Molecule Force Field with Broad Chemical Space Coverage for Molecular Modeling and Free Energy Calculation . Journal of Chemical Theory and Computation 2024 , 20 , 799 – 818 . OpenUrl (38). ↵ Boresch , S. ; Tettinger , F. ; Leitgeb , M. ; Karplus , M. Absolute binding free energies: A quantitative approach for their calculation . Journal of Physical Chemistry B 2003 , 107 , 9535 – 9551 . OpenUrl (39). ↵ Feng , M. ; Heinzelmann , G. ; Gilson , M. K. Absolute binding free energy calculations improve enrichment of actives in virtual compound screening . Scientific Reports 2022 , 12 , 13640 . OpenUrl PubMed (40). ↵ Michino , M. ; Beautrait , A. ; Boyles , N. A. ; Nadupalli , A. ; Dementiev , A. ; Sun , S. ; Ginn , J. ; Baxt , L. ; Suto , R. ; Bryk , R. ; Jerome , S. V. ; Huggins , D. J. ; Vendome , J. Shape-Based Virtual Screening of a Billion-Compound Library Identifies Mycobacterial Lipoamide Dehydrogenase Inhibitors . ACS Bio and Med Chem Au 2023 , 3 , 507 – 515 . OpenUrl (41). ↵ Andersen , L. ; Rausch-Dupont , M. ; Martínez León , A. ; Volkamer , A. ; Hub , J. S. ; Klakow , D. Accelerating ligand discovery by combining Bayesian optimization with MMGBSA-based binding affinity calculations . bioRxiv 2025 , 1 – 40 . (42). ↵ Aldeghi , M. ; Heifetz , A. ; Bodkin , M. J. ; Knapp , S. ; Biggin , P. C. Accurate calculation of the absolute free energy of binding for drug molecules . Chemical Science 2016 , 7 , 207 – 218 . OpenUrl PubMed (43). ↵ Vernuccio , R. et al. Structural insights into tecovirimat antiviral activity and poxvirus resistance . Nature Microbiology 2025 , 10 , 734 – 748 . OpenUrl PubMed (44). ↵ Aldeghi , M. ; Heifetz , A. ; Bodkin , M. J. ; Knapp , S. ; Biggin , P. C. Predictions of ligand selectivity from absolute binding free energy calculations . Journal of the American Chemical Society 2017 , 139 , 946 – 957 . OpenUrl CrossRef PubMed (45). ↵ Lindahl , E. ; Friedman , R. Estimation of Absolute Binding Free Energies for Drugs That Bind Multiple Proteins . Journal of Chemical Information and Modeling 2025 , 65 , 3431 – 3438 . OpenUrl PubMed (46). ↵ Chen , W. ; Cui , D. ; Jerome , S. V. ; Michino , M. ; Lenselink , E. B. ; Huggins , D. J. ; Beautrait , A. ; Vendome , J. ; Abel , R. ; Friesner , R. A. ; Wang , L. Enhancing Hit Discovery in Virtual Screening through Absolute Protein-Ligand Binding Free-Energy Calculations . Journal of Chemical Information and Modeling 2023 , 63 , 3171 – 3185 . OpenUrl CrossRef PubMed (47). Araujo , S. C. ; Maltarollo , V. G. ; Almeida , M. O. ; Ferreira , L. L. ; Andricopulo , A. D. ; Honorio , K. M. Structure-based virtual screening, molecular dynamics and binding free energy calculations of HIT candidates as ALK-5 inhibitors . Molecules 2020 , 25 , 1 – 14 . OpenUrl CrossRef PubMed (48). ↵ Li , Z. ; Huang , Y. ; Wu , Y. ; Chen , J. ; Wu , D. ; Zhan , C. G. ; Luo , H. B. Absolute Binding Free Energy Calculation and Design of a Subnanomolar Inhibitor of Phosphodiesterase-10 . Journal of Medicinal Chemistry 2019 , 62 , 2099 – 2111 . OpenUrl PubMed (49). ↵ Ries , B. ; Alibay , I. ; Anand , N. M. ; Biggin , P. C. ; Magarkar , A. Automated Absolute Binding Free Energy Calculation Workflow for Drug Discovery . Journal of Chemical Information and Modeling 2024 , 64 , 5357 – 5364 . OpenUrl PubMed (50). ↵ Clark , F. ; Robb , G. R. ; Cole , D. J. ; Michel , J. Automated Adaptive Absolute Binding Free Energy Calculations . Journal of Chemical Theory and Computation 2024 , 20 , 7806 – 7828 . OpenUrl (51). ↵ Heinzelmann , G. ; Huggins , D. J. ; Gilson , M. K. BAT2: an Open-Source Tool for Flexible, Automated, and Low Cost Absolute Binding Free Energy Calculations . Journal of Chemical Theory and Computation 2024 , 20 , 6518 – 6530 . OpenUrl (52). ↵ Santiago-McRae , E. ; Ebrahimi , M. ; Sandberg , J. W. ; Brannigan , G. ; Hénin , J. Computing absolute binding affinities by Streamlined Alchemical Free Energy Perturbation [Article v1.0] . Living Journal of Computational Molecular Science 2023 , 5 , 1 – 25 . OpenUrl (53). ↵ Fu , H. ; Chen , H. ; Blazhynska , M. ; de Lacam , E. G. C. ; Szczepaniak , F. ; Pavlova , A. ; Shao , X. ; Gumbart , J. C. ; Dehez , F. ; Roux , B. ; Cai , W. ; Chipot , C. Accurate determination of protein:ligand standard binding free energies from molecular dynamics simulations . Nature Protocols 2022 , 17 , 1114 – 1141 . OpenUrl PubMed (54). ↵ Fu , H. ; Chen , H. ; Cai , W. ; Shao , X. ; Chipot , C. BFEE2: Automated, Streamlined, and Accurate Absolute Binding Free-Energy Calculations . Journal of Chemical Information and Modeling 2021 , 61 , 2116 – 2123 . OpenUrl PubMed (55). ↵ Fu , H. ; Gumbart , J. C. ; Chen , H. ; Shao , X. ; Cai , W. ; Chipot , C. BFEE: A User-Friendly Graphical Interface Facilitating Absolute Binding Free-Energy Calculations . Journal of Chemical Information and Modeling 2018 , 58 , 556 – 560 . OpenUrl CrossRef PubMed (56). ↵ Dakka , J. ; Farkas-Pall , K. ; Turilli , M. ; Wright , D. W. ; Coveney , P. V. ; Jha , S. Concurrent and adaptive extreme scale binding free energy calculations . Proceedings -IEEE 14th International Conference on eScience, e-Science 2018 . 2018 ; pp 189 – 200 . (57). ↵ Loeffler , H. H. ; Michel , J. ; Woods , C. FESetup: Automating Setup for Alchemical Free Energy Simulations . Journal of Chemical Information and Modeling 2015 , 55 , 2485 – 2490 . OpenUrl CrossRef PubMed (58). ↵ Fu , H. ; Zhou , Y. ; Jing , X. ; Shao , X. ; Cai , W. Meta-Analysis Reveals That Absolute Binding Free-Energy Calculations Approach Chemical Accuracy . Journal of Medicinal Chemistry 2022 , 65 , 12970 – 12978 . OpenUrl CrossRef PubMed (59). ↵ Kirkwood , J. G. Statistical mechanics of fluid mixtures . The Journal of Chemical Physics 1935 , 3 , 300 – 313 . OpenUrl CrossRef (60). ↵ Shirts , M. R. ; Chodera , J. D. Statistically optimal analysis of samples from multiple equilibrium states . Journal of Chemical Physics 2008 , 129 , 1 – 10 . OpenUrl (61). ↵ Wang , J. ; Wolf , R. M. ; Caldwell , J. W. ; Kollman , P. A. ; Case , D. A. Development and testing of a general Amber force field . Journal of Computational Chemistry 2004 , 25 , 1157 – 1174 . OpenUrl CrossRef PubMed Web of Science (62). ↵ Behara , P. K. et al. Development and Benchmarking of Open Force Field 2.0.0: The Sage Small Molecule Force Field . Journal of Chemical Theory and Computation 2023 , 19 , 3251 – 3275 . OpenUrl (63). ↵ Wang , Y. ; Fass , J. ; Kaminow , B. ; Herr , J. E. ; Rufa , D. ; Zhang , I. ; Pulido , I. ; Henry , M. ; Macdonald , H. E. B. ; Takaba , K. ; Chodera , J.D. End-to-end differentiable construction of molecular mechanics force fields . Chemical Science 2022 , 13 , 12016 – 12033 . OpenUrl PubMed (64). ↵ Chipot , C. ; Pohorille , A. Free energy calculations ; Springer , 2007 ; Vol. 86 . (65). ↵ Valdés-Tresanco , M. S. ; Valdés-Tresanco , M. E. ; Valiente , P. A. ; Moreno , E. gmx_MMPBSA: a new tool to perform end-state free energy calculations with GROMACS . Journal of chemical theory and computation 2021 , 17 , 6281 – 6291 . OpenUrl (66). ↵ Wang , E. ; Sun , H. ; Wang , J. ; Wang , Z. ; Liu , H. ; Zhang , J. Z. ; Hou , T. End-point binding free energy calculation with MM/PBSA and MM/GBSA: strategies and applications in drug design . Chemical Reviews 2019 , 119 , 9478 – 9508 . OpenUrl CrossRef PubMed (67). ↵ Lee , B. ; Richards , F. M. The interpretation of protein structures: estimation of static accessibility . Journal of molecular biology 1971 , 55 , 379 – IN4 . OpenUrl CrossRef PubMed Web of Science (68). ↵ Duan , L. ; Liu , X. ; Zhang , J. Z. Interaction entropy: A new paradigm for highly efficient and reliable computation of protein–ligand binding free energy . Journal of the American Chemical Society 2016 , 138 , 5722 – 5728 . OpenUrl CrossRef PubMed (69). ↵ Ekberg , V. ; Ryde , U. On the Use of Interaction Entropy and Related Methods to Estimate Binding Entropies . Journal of Chemical Theory and Computation 2021 , 17 , 5379 – 5391 . OpenUrl (70). ↵ Mölder , F. ; Jablonski , K. P. ; Letcher , B. ; Hall , M. B. ; Tomkins-Tinch , C. H. ; Sochat , V. ; Forster , J. ; Lee , S. ; Twardziok , S. O. ; Kanitz , A. ; Wilm , A. ; Holtgrewe , M. ; Rahmann , S. ; Nahnsen , S. ; Köster , J. Sustainable data analysis with Snakemake . F1000Research 2021 , 10 , 33 . OpenUrl (71). ↵ Abraham , M. J. ; Murtola , T. ; Schulz , R. ; Páll , S. ; Smith , J. C. ; Hess , B. ; Lindahl , E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers . SoftwareX 2015 , 1 , 19 – 25 . OpenUrl CrossRef (72). ↵ Yoo , A. B. ; Jette , M. A. ; Grondona , M. Slurm: Simple linux utility for resource management . Workshop on job scheduling strategies for parallel processing . 2003 ; pp 44 – 60 . (73). ↵ Lindorff-Larsen , K. ; Piana , S. ; Palmo , K. ; Maragakis , P. ; Klepeis , J. L. ; Dror , R. O. ; Shaw , D. E. Improved side-chain torsion potentials for the Amber ff99SB protein force field . Proteins: Structure, Function, and Bioinformatics 2010 , 78 , 1950 – 1958 . OpenUrl (74). ↵ Grote , F. ; Lyubartsev , A. P. Optimization of Slipids Force Field Parameters Describing Headgroups of Phospholipids . Journal of Physical Chemistry B 2020 , 124 , 8784 – 8793 . OpenUrl (75). ↵ Jorgensen , W. L. ; Chandrasekhar , J. ; Madura , J. D. ; Impey , R. W. ; Klein , M. L. Comparison of simple potential functions for simulating liquid water . The Journal of Chemical Physics 1983 , 79 , 926 – 935 . OpenUrl CrossRef PubMed Web of Science (76). ↵ Martínez León , A. TOFF: Topology from Open Force Fields . 2023 ; doi: 10.5281/zenodo.8189649 . OpenUrl CrossRef (77). ↵ Rizzi , A. ; Murkli , S. ; McNeill , J. N. ; Yao , W. ; Sullivan , M. ; Gilson , M. K. ; Chiu , M. W. ; Isaacs , L. ; Gibb , B. C. ; Mobley , D. L. ; Chodera , J. D. Overview of the SAMPL6 host-guest binding affinity prediction challenge . Journal of Computer-Aided Molecular Design 2018 , 32 , 937 – 963 . OpenUrl CrossRef PubMed (78). ↵ Isik , M. ; Rizzi , A. ; Mobley , D. L. ; Bergazin , T. D. ; Shirts , M. SAMPL6 . 2024 . (79). ↵ Khalak , Y. ; Tresadern , G. ; Aldeghi , M. ; Baumann , H. M. ; Mobley , D. L. ; de Groot , B. L. ; Gapsys , V. Alchemical absolute protein-ligand binding free energies for drug design . Chemical Science 2021 , 12 , 13958 – 13971 . OpenUrl PubMed (80). ↵ Croux , C. ; Dehon , C. Influence functions of the Spearman and Kendall correlation measures . Statistical Methods and Applications 2010 , 19 , 497 – 515 . OpenUrl CrossRef Web of Science (81). ↵ Ross , G. A. ; Lu , C. ; Scarabelli , G. ; Albanese , S. K. ; Houang , E. ; Abel , R. ; Harder , E. D. ; Wang , L. The maximal and current accuracy of rigorous protein-ligand binding free energy calculations . Communications Chemistry 2023 , 6 , 1 – 12 . OpenUrl PubMed (82). ↵ Lu , C. ; Wu , C. ; Ghoreishi , D. ; Chen , W. ; Wang , L. ; Damm , W. ; Ross , G. A. ; Dahlgren , M. K. ; Russell , E. ; Von Bargen , C. D. ; others OPLS4: improving force field accuracy on challenging regimes of chemical space . Journal of chemical theory and computation 2021 , 17 , 4291 – 4300 . OpenUrl (83). ↵ Su , P. C. ; Tsai , C. C. ; Mehboob , S. ; Hevener , K. E. ; Johnson , M. E. Comparison of radii sets, entropy, QM methods, and sampling on MM-PBSA, MM-GBSA, and QM/MM-GBSA ligand binding energies of F. tularensis enoyl-ACP reductase (FabI) . Journal of Computational Chemistry 2015 , 36 , 1859 – 1873 . OpenUrl CrossRef PubMed (84). ↵ Wang , K. ; Chodera , J. D. ; Yang , Y. ; Shirts , M. R. Identifying ligand binding sites and poses using GPU-accelerated Hamiltonian replica exchange molecular dynamics . Journal of Computer-Aided Molecular Design 2013 , 27 , 989 – 1007 . OpenUrl CrossRef PubMed (85). Chodera , J. D. ; Shirts , M. R. Replica exchange and expanded ensemble simulations as Gibbs sampling: Simple improvements for enhanced mixing . The Journal of Chemical Physics 2011 , 135 , 1 – 15 . OpenUrl CrossRef (86). ↵ Sugita , Y. ; Kitao , A. ; Okamoto , Y. Multidimensional replica-exchange method for free-energy calculations . The Journal of Chemical Physics 2000 , 113 , 6042 – 6051 . OpenUrl CrossRef Web of Science (87). ↵ Darden , T. ; York , D. ; Pedersen , L. Particle mesh Ewald: An N·log(N) method for Ewald sums in large systems . The Journal of Chemical Physics 1998 , 98 , 10089 – 10092 . OpenUrl CrossRef (88). ↵ Essmann , U. ; Perera , L. ; Berkowitz , M. L. ; Darden , T. ; Lee , H. ; Pedersen , L. G. A smooth particle mesh Ewald method . The Journal of Chemical Physics 1998 , 103 , 8577 – 8593 . OpenUrl CrossRef (89). ↵ Wu , Z. ; Biggin , P. C. Correction Schemes for Absolute Binding Free Energies Involving Lipid Bilayers . Journal of Chemical Theory and Computation 2022 , 18 , 2657 – 2672 . OpenUrl (90). ↵ Chen , W. ; Deng , Y. ; Russell , E. ; Wu , Y. ; Abel , R. ; Wang , L. Accurate calculation of relative binding free energies between ligands with different net charges . Journal of chemical theory and computation 2018 , 14 , 6346 – 6358 . OpenUrl (91). ↵ Rocklin , G. J. ; Mobley , D. L. ; Dill , K. A. ; Hünenberger , P. H. Calculating the binding free energies of charged species based on explicit-solvent simulations employing latticesum methods: An accurate correction scheme for electrostatic finite-size effects . The Journal of Chemical Physics 2013 , 139 , 1 – 32 . OpenUrl (92). ↵ Öhlknecht , C. ; Lier , B. ; Petrov , D. ; Fuchs , J. ; Oostenbrink , C. Correcting electrostatic artifacts due to net-charge changes in the calculation of ligand binding free energies . Journal of Computational Chemistry 2020 , 41 , 986 – 999 . OpenUrl PubMed (93). ↵ Ross , G. A. ; Bodnarchuk , M. S. ; Essex , J. W. Water sites, networks, and free energies with grand canonical Monte Carlo . Journal of the American Chemical Society 2015 , 137 , 14930 – 14943 . OpenUrl CrossRef PubMed (94). ↵ Bruce Macdonald , H. E. ; Cave-Ayland , C. ; Ross , G. A. ; Essex , J. W. Ligand binding free energies with adaptive water networks: two-dimensional grand canonical alchemical perturbations . Journal of Chemical Theory and Computation 2018 , 14 , 6586 – 6597 . OpenUrl (95). ↵ Haider , K. ; Cruz , A. ; Ramsey , S. ; Gilson , M. K. ; Kurtzman , T. Solvation structure and thermodynamic mapping (SSTMap): an open-source, flexible package for the analysis of water in molecular dynamics trajectories . Journal of Chemical Theory and Computation 2018 , 14 , 418 – 425 . OpenUrl (96). ↵ Grubmüller , H. ; Groll , V. SOLVATE . https://www.mpinat.mpg.de/grubmueller/solvate , 2025 ; Max Planck Institute for Multidisciplinary Sciences . (97). ↵ Behera , S. ; Hahn , D. F. ; Wilson , C. J. ; Marsili , S. ; Tresadern , G. ; Gapsys , V. ; de Groot , B. L. Quantification of the Impact of Structure Quality on Predicted Binding Free Energy Accuracy . Journal of Chemical Information and Modeling 2025 , 65 , 6927 – 6938 . OpenUrl PubMed (98). ↵ Jo , S. ; Kim , T. ; Iyer , V. G. ; Im , W. CHARMM-GUI: A web-based graphical user interface for CHARMM . Journal of Computational Chemistry 2008 , 29 , 1859 – 1865 . OpenUrl CrossRef PubMed (99). ↵ Maier , J. A. ; Martinez , C. ; Kasavajhala , K. ; Wickstrom , L. ; Hauser , K. E. ; Simmerling , C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB . Journal of Chemical Theory and Computation 2015 , 11 , 3696 – 3713 . OpenUrl (100). ↵ Eastman , P. ; Swails , J. ; Chodera , J. D. ; McGibbon , R. T. ; Zhao , Y. ; Beauchamp , K. A. ; Wang , L. P. ; Simmonett , A. C. ; Harrigan , M. P. ; Stern , C. D. ; Wiewiora , R. P. ; Brooks , B. R. ; Pande , V. S. OpenMM 7: Rapid development of high performance algorithms for molecular dynamics . PLoS Computational Biology 2017 , 13 , 1 – 17 . OpenUrl (101). ↵ Shirts , M. R. ; Klein , C. ; Swails , J. M. ; Yin , J. ; Gilson , M. K. ; Mobley , D. L. ; Case , D. A. ; Zhong , E. D. Lessons learned from comparing molecular dynamics engines on the SAMPL5 dataset . Journal of Computer-Aided Molecular Design 2017 , 31 , 147 – 161 . OpenUrl CrossRef PubMed (102). ↵ Hess , B. ; Bekker , H. ; Berendsen , H. J. ; Fraaije , J. G. LINCS: A Linear Constraint Solver for molecular simulations . Journal of Computational Chemistry 1997 , 18 , 1463 – 1472 . OpenUrl CrossRef PubMed Web of Science (103). ↵ Miyamoto , S. ; Kollman , P. A. Settle: An analytical version of the SHAKE and RATTLE algorithm for rigid water models . Journal of Computational Chemistry 1992 , 13 , 952 – 962 . OpenUrl CrossRef PubMed Web of Science (104). ↵ Bernetti , M. ; Bussi , G. Pressure control using stochastic cell rescaling . Journal of Chemical Physics 2020 , 153 , 1 – 12 . OpenUrl CrossRef (105). ↵ Berendsen , H. J. C. ; Postma , J. P. M. ; van Gunsteren , W. F. ; DiNola , A. ; Haak , J. R. Molecular dynamics with coupling to an external bath . The Journal of Chemical Physics 1984 , 81 , 3684 – 3690 . OpenUrl CrossRef PubMed Web of Science (106). ↵ Parrinello , M. ; Rahman , A. Crystal structure and pair potentials: A moleculardynamics study . Physical Review Letters 1980 , 45 , 1196 – 1199 . OpenUrl CrossRef Web of Science (107). ↵ Alibay , I. IAlibay/MDRestraintsGenerator: MDRestraintsGenerator v0.2.0 . 2022 ; doi: 10.5281/zenodo.6972482 . OpenUrl CrossRef (108). ↵ Wu , Z. et al. alchemlyb: the simple alchemistry library . Journal of Open Source Software 2024 , 9 , 6934 . OpenUrl PubMed View the discussion thread. Back to top Previous Next Posted September 27, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following BindFlow: a free, user-friendly pipeline for absolute binding free energy calculations using free energy perturbation or MM(PB/GB)SA Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share BindFlow: a free, user-friendly pipeline for absolute binding free energy calculations using free energy perturbation or MM(PB/GB)SA Alejandro Martínez León , Lucas Andersen , Jochen S. Hub bioRxiv 2025.09.25.678545; doi: https://doi.org/10.1101/2025.09.25.678545 Share This Article: Copy Citation Tools BindFlow: a free, user-friendly pipeline for absolute binding free energy calculations using free energy perturbation or MM(PB/GB)SA Alejandro Martínez León , Lucas Andersen , Jochen S. Hub bioRxiv 2025.09.25.678545; doi: https://doi.org/10.1101/2025.09.25.678545 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Biophysics Subject Areas All Articles Animal Behavior and Cognition (7629) Biochemistry (17660) Bioengineering (13881) Bioinformatics (41911) Biophysics (21436) Cancer Biology (18578) Cell Biology (25482) Clinical Trials (138) Developmental Biology (13371) Ecology (19887) Epidemiology (2067) Evolutionary Biology (24302) Genetics (15599) Genomics (22483) Immunology (17728) Microbiology (40364) Molecular Biology (17163) Neuroscience (88537) Paleontology (666) Pathology (2830) Pharmacology and Toxicology (4821) Physiology (7637) Plant Biology (15129) Scientific Communication and Education (2045) Synthetic Biology (4290) Systems Biology (9817) Zoology (2269)