The interactions between biomolecules (BMs) and nanomaterials (NMs) are receiving growing research interest due to their diverse applications in fields such as nanomedicine, biosensing, and biocatalysis. BMs, which include nucleic acids, proteins, and carbohydrates [1], interact with various types of NMs, such as 0D quantum dots, 1D nanotubes, 2D nanosheets, and 3D nanocomposites [2]. For example, C₃N, an ultra-small nanodot, could effectively disassemble mature Aβ fibrils and hence reduce aggregation-related neuron cytotoxicity both in vitro and in vivo [3]. An ultra-thin MoS₂-graphene heterostructure nanopore was shown to prolong the translocation time of λDNA and bovine serum albumin (BSA), while acquiring detailed information about these single molecules [4]. Candida rugosa lipase (CRL) adsorbed onto multi-walled carbon nanotubes has enhanced the production of geranyl propionate, doubling enzyme activity compared to free CRL [5]. Despite these successes, a deeper understanding of molecular interactions often requires the aid of an important in silico method–molecular simulation. This area has rapidly expanded over the years, driven by significant advances in computer power and revolutionized by machine learning techniques. Molecular simulations are especially adept at elucidating underlying molecular mechanisms and predicting unknown molecular properties. They have been extensively applied to study BM-NM interactions, in which free energy serves as an essential quantity to estimate the interaction strength between molecules and evaluate the spontaneity of a physical process. Unfortunately, although free energy calculations are very popular in computer-aided drug design (CADD) to accurately predict the drug-target interactions (DTIs) [6], they generally receive much less attention in BM-NM interactions, due to limited availability of experimental data, lack of standardization in interaction types, and complexity and diversity of NMs [7]. This mini-review focuses on the most widely used computational techniques for free energy calculations specifically for the BM-NM systems, providing readers with a succinct guidance on computing various kinds of free energy.

There are many classifications concerning free energy concepts. First, free energy changes help determine the direction of spontaneous processes. Whilst Helmholtz free energy difference (ΔA) is defined for systems at constant volume, Gibbs free energy difference (ΔG) applies to systems at constant pressure. In molecular simulations of condensed matter systems, the NPT ensemble is generally adopted for production runs, i.e., the number of molecules (N), the pressure (P), and the temperature (T) of the system remain constant during simulations. Consequently, this paper only takes Gibbs free energy into consideration instead of Helmholtz free energy. However, it should be noted that, in such simulations, the volume change ΔV is usually small, making the two free energy measures generally comparable. Second, the binding affinity between BMs and NMs is described using binding free energy, which comes in two forms: absolute binding free energy (ABFE) and relative binding free energy (RBFE). ABFE offers an exact value of the binding free energy between the bound state (e.g., protein-graphene complex) and the unbound states (e.g., free protein and free graphene), to be compared directly with an experimental value. By contrast, RBFE measures the difference in binding free energy between two or more ligands (e.g., proteins) binding to the same receptor (e.g., graphene), to make a comparison between different ligands. The former is typically achieved by molecular dynamics (MD) simulations with enhanced sampling techniques, while the latter is typically achieved by alchemical methods. Thirdly, to calculate free energies, some methods only consider the final bound state, some require the initial and final states, whereas others compute the free energy changes along a reaction coordinate. Accordingly, we categorize free energy calculation methods into three main types: endpoint methods, alchemical methods, and pathway methods.

This paper is organized as follows. In Section 2, we introduce two similar endpoint methods: molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) and molecular mechanics/generalized Born surface area (MM/GBSA). In Section 3, we discuss in detail three alchemical methods: Free Energy Perturbation (FEP), Thermodynamic Integration (TI), and Bennett Acceptance Ratio (BAR). In Section 4, pathway methods are elaborated including Umbrella Sampling (US), Jarzynski Equality (JE), and Metadynamics (MtD). Their relative advantages and disadvantages are compared, and some examples are illustrated especially for BM-NM systems. Finally, we outline future directions for improvement in free energy calculations.

MM/PBSA or MM/GBSA In these two similar methods, the binding free energy can be calculated as ΔG = ΔG_BN-ΔG_B-ΔG_N and divided into different contributions using Equation 1: ∆ G = ∆ E M M + ∆ G s o l − T ∆ S = ∆ E int + ∆ E v d W + ∆ G s o l _ p o l a r + ∆ G s o l _ n o n p o l a r − T ∆ S (1) where ΔE _MM is the gas-phase molecular mechanics (MM) energy, which can be decomposed into contributions of the internal energy ΔE _int (the sum of ΔE _bond, ΔE _angle, and ΔE _diheral), the electrostatic energy ΔE _ele, and the van der Waals energy ΔE _vdW; ΔG _sol is the solvation free energy, which can be decomposed into the polar solvation energy ΔG _{sol_polar} calculated from either the Poisson-Boltzmann (PB) or generalized Born (GB) model, and the nonpolar contribution ΔG _{sol_nonpolar} estimated by the solvent accessible surface area (SASA); and −TΔS denotes the change in conformational entropy of BMs, normally derived from normal mode analysis. For more detailed information about these two methods, readers are referred to the original works [8–11] and recent review articles [12–14]. Figure 1 provides a typical example of these endpoint methods, illustrating the binding of a peptide drug (ID: 7mll) to a (16,16) carbon nanotube (CNT). Several computational programs specifically designed for these methods are available, such as “g_mmpbsa” [15] and “gmx mmpbsa” [16] for Gromacs, which are popular software tools for modelling BM-NM systems [17]. These programs are user-friendly and facilitate quick access to free energy calculations. Due to their simplicity, straightforward application, and low computational cost, MM/PBSA and MM/GBSA methods have been extensively used in BM-NM systems [18–25]. These methods strike a good balance between accuracy and computational efficiency, making them among the most popular choices for calculating binding free energies. Moreover, the breakdown of energy contributions also provides useful information for further energy analysis.

It is important to note that the entropic term −TΔS is not included in binding free energy calculations using “g_mmpbsa” and “gmx mmpbsa”, due to its high computational cost, trivial contributions to the total binding energy between proteins and drugs, and the substantial standard error compared to enthalpic and solvation terms [15, 16]. Consequently, strictly speaking, the binding free energies calculated by these programs are not true ABFE, and their accuracy may vary depending on the system. For instance, when a BM binds to a rigid NM, it may undergo varying degrees of conformational changes, and in some cases, its structure could be significantly altered or even disrupted [26–28]. Nevertheless, these programs are well-suited for calculating RBFE, such as comparing the binding affinities of different amino acid residues to the same material [19]. Recent developments in entropy calculation methods, such as quasi-harmonic approximation [29], interaction entropy [30], and multiscale cell correlation [31], might help address this limitation [32], though further investigation is still required [33]. Another concern over these end-point methods is the use of implicit solvent models for evaluating solvation energy, whereas MD simulations employ explicit water models that significantly influence BM structure and dynamics [34, 35].

Free Energy Perturbation In 1954, Zwanzig developed the equation for the FEP theory in which the binding free energy was associated with the potential energy difference of two ensemble states 1 and 0 using Equation 2 [36]: ∆ G = G 1 − G 0 = − k B T I n 〈 exp − U 1 − U 0 k B T 〉 0 (2) where U ₁ and U ₀ are the potential energy of states 1 and 0, respectively; k _B is the Boltzmann constant, and T is the temperature. If the initial and final states are only influenced by a small perturbation, state 1 can be considered as the final state and state 0 as the initial state. However, in most situations, state 1 and state 0 have little overlap, and therefore a series of alchemical (nonphysical) intermediate states are created via a coupling parameter λ, which varies from 0 to 1. The potential function for these intermediate states is then achieved using U_λ = λU₀ + (1-λ)U₁. The total free energy is finally calculated through adding all the free energy differences between successive states through constructed thermodynamic cycle. For more information about FEP, the readers can refer to more recent advances [37–41], review papers [6, 42], as well as real applications [43–47]. Of all the tools to implement FEP, perhaps the most famous is FEP+, which was developed by Schrödinger Inc. [48].

Thermodynamics Integration In 1935, Kirkwood proposed a theory of fluid mixtures, which later became the foundation for the TI method [49]. In this method, the free energy change between two states of the investigated system is connected by a coupling variable λ, being 0 for the initial state and 1 for the final state. By continuously changing λ, a series of intermediate states between the initial and final states are obtained. By sampling these intermediate states, the free energy change is finally calculated using Equation 3: ∆ G = ∫ 0 1 〈 ∂ U λ ∂ λ 〉 λ d λ (3) where U _λ is the potential energy function of the system as a function of λ, and the angle brackets denote an ensemble average at each λ value. The major difference between FEP and TI is that the adjacent states sampled by FEP need to be overlapped in phase space whereas those sampled by TI may not be overlapped [50]. For more information about this method, readers can further refer to many of the classical works [50–55]. To handle the complicated procedure for FEP and TI, some tools were developed based on alchemical transformations specifically for Gromacs [56] and Amber [57].

Bennett Acceptance Ratio In 1976, Bennett introduced an implicit equation obtaining free energy difference by an iterative numerical approach using Equation 4 [58]: Δ G = G 1 − G 0 = − k B T I n 〈 1 1 + ⁡ exp U 1 − U 0 + ∆ G k B T 〉 0 〈 1 1 + ⁡ exp U 1 − U 0 − ∆ G k B T 〉 1 (4) where U ₁ and U ₀ are the potential energy of states 1 and 0, respectively. Unlike FEP, BAR requires gathering samples of potential energies for configurations at both states 0 and 1. This method was later extended to the multistate Bennett acceptance ratio (MBAR) method in 2008, which calculates the free energy from multiple states rather than two states [59]. For more information about this method, readers can refer to method comparative studies [60, 61], relevant review papers [62–64], application-focused works [65], and recent developments [66–68].

Figure 1 also illustrates the application of alchemical methods for calculating the binding energies between an α-helical protein (ID: 1fmh) and a (26, 26) CNT. Although alchemical methods are appealing for calculating protein-ligand binding free energies and solvation free energies, they are typically limited to small ligands [35]. Many BMs, however, are macromolecules, and their interactions with NMs can involve multiple binding modes and orientations. This complexity makes it computationally expensive to adequately sample the conformational space needed for converged free energy estimates. Additionally, NMs often have complex and heterogeneous surfaces, and their properties such as hydrophobicity and charge distribution can also complicate the accurate modeling of interactions with BMs. These factors, alongside the unavailability of relevant experimental data, hinder the application of alchemical methods to BM-NM systems. Consequently, we can find a limited quantity of literature reporting the application of FEP, TI, and BAR to BM-NM systems. Shen et al. used TI to investigate the binding free energies of single nucleotides adenine (A) and thymine (T) on CNT inner walls [69]. In this scenario, A and T are relatively small ligands, which may potentially mitigate some of the major limitations associated with using alchemical methods in BM-NM systems. Additionally, alchemical methods could provide insights into how mutations in BMs or changes in NM surface chemistry influence binding and interactions. These methods could also be applied to explore the interactions between small nanoparticles (e.g., fullerene) and BMs during the formation of hybrid functional assemblies [70].

Umbrella Sampling In statistical thermodynamics, the free energy of a system is related to the probability distribution of the system. The free energy difference from state 0 to state 1 can be obtained from the probability distribution of the system in the two states using Equation 5: ∆ G = G 1 − G 0 = − k B T I n ρ 0 ρ 1 (5) where ρ ₀ and ρ ₁ are the probability distributions of the system in states 0 and 1. However, for transitions between states that have high energy barriers, the regions near these barriers are poorly sampled. To overcome these challenges, umbrella sampling employs a series of biased simulations, where a harmonic umbrella potential is added to sample configurations within each window along the reaction coordinate ξ [71]. The unbiased free energy profile in the form of potential of mean force (PMF) [72, 73] can be constructed from data in all the windows using the weighted histogram analysis method (WHAM) [74–76]. This method was proposed by Torrie and Valleau in 1977. A typical example using US for protein-ligand interactions can be found here [77], and recent advances are also provided here [72, 78, 79]. In BM-NM systems, the US method can be widely used to compute PME versus ξ with good accuracy despite its high computational cost [26, 80–83]. It can be applied to obtain free energy differences as well as free energy profiles. For example, it was reported that a protein can spontaneously enter carbon nanotubes (CNTs) in aqueous solutions, while releasing it from the tubes is a rare event [84]. Using US, one can easily understand the spontaneity of this encapsulation process by calculating the PMF along the reaction coordinate, where a wide potential well can be clearly observed [26, 80].

Jarzynski Equality JE seeks to relate the free energy differences between two states to the work done during a non-equilibrium transformation between those states. To visit the high free energy regions between states, a possible method is steered molecular dynamics (SMD), in which the harmonic potential guides the system from low free energy regions to high free energy regions. It should be noted that US also uses SMD to build starting configurations along ξ. However, unlike in US sampling where additional harmonic potential is static, the harmonic potential in JE sampling moves at a constant speed over time. By a sufficient number of samples of irreversible work, the free energy difference of the reaction process is directly related to the irreversible work done to the system using Equation 6 [85]: e − ∆ G k B T = 〈 e − W k B T 〉 (6) where W is the irreversible work. This Jarzynski equality, which was proposed in 1997, provides a way to calculate equilibrium free energy differences using non-equilibrium measurements. Further improvements on JE can be found in the literature [86–89]. Figure 1 depicts an example where SMD simulations are used to pull a protein (ID: 1jmq) out of a CNT. In this process, US extracts a series of configurations along the reaction pathway as sampling windows for further biased simulations, while the JE method directly calculates the work done on the protein by integrating the force over the pulling distance. Numerous studies have shown that JE can yield results comparable in accuracy to US when determining free energy for simple systems [54, 90–93]. However, there are some limitations to the JE method. First, the SMD pulling might not always align with the most physically favorable separation pathway required to accurately reproduce the binding affinity. Second, for slightly more complex systems, energy dissipation and the impact of the applied harmonic potential on the conformations of BMs can significantly affect the accuracy of free energy calculations [94–98]. Zhang et al. employed JE to calculate the adsorption free energy of biomolecules on the surface of hydroxyapatite, indicating that the structured water layer at the solid-liquid interface causes SMD simulations to produce a large amount of dissipative work [98]. They further proposed a hybrid scheme combining JE with US to address this issue.

Metadynamics MtD, similar to US and JE, uses an additional potential to sample high free energy regions of the system. However, unlike the other two methods, MtD operates by applying a time-dependent biasing potential to selected collective variables (CVs). The biasing potential is periodically introduced through the periodic addition of Gaussian-shaped repulsive potentials to the overall potential energy of the system. This prevents the system from revisiting areas of the reaction coordinate, and continually drives the system towards higher free energy regions. Ultimately, free energy can be calculated based on the series of repulsive potentials using Equation 7: ∆ G = − ∑ i = 1 N G i ⁡ exp − s − s i 2 σ 2 (7) where G _i is the height of the ith Gaussian, s _i is the position of the ith Gaussian in CV space, and σ is the width of the Gaussian functions. This method, proposed by Laio and Parrinello in 2002 [99], offers several advantages for calculating free energy. Unlike US and JE, MtD does not require the prior determination of the reaction pathway. Instead, it accelerates the sampling of high free energy events by gradually pushing the system from low to high free energy regions through a series of added repulsive potentials. MtD is especially useful for high-dimensional free energy calculations, such as two-dimensional free energy surfaces. Consequently, it has become widely used to map the free energy landscape of protein-ligand interactions, and can be easily extended to BM-NM systems, even though the number of such studies has still been small so far. MtD can be implemented using the portable plugin PLUMED [100]. For more details, the readers are asked to see more review papers [101–104]. NM-BM systems often involve rare events, such as the dissociation of a protein from its binding NM. MtD can enhance the sampling of these events, which might be difficult to capture using conventional MD simulations. An example of using MtD to calculate the binding free energy between P-glycoprotein and two drugs, paclitaxel and doxorubicin, is presented here [105], which can also be extended to BM-NM interactions. However, the success of MtB depends on the selection of CVs, biasing potentials, and the complex interpretation of data, making it more suitable for experienced users. Also, the application of MtD to large BMs demands substantial computational resources. These drawbacks underscore the necessity for continued improvements in this method.

Free energy methods can accurately predict the interaction strength between BMs and NMs. Table 1 summarizes highlights and limitations of all the methods discussed within this paper. Future research should address several key areas to enhance these methods. First, improving and optimizing force fields are crucial for achieving accurate free energy calculations. Second, the development of more user-friendly software and programs is needed to make these methods more accessible and commercially viable. Third, while many end-to-end DTI models can precisely predict protein-ligand interactions [106–108], there is a need for extensive investigation into the application of machine learning to BM-NM interactions.