Improved Sampling in Ab Initio Free Energy Calculations of Biomolecules at Solid–Liquid Interfaces: Tight-Binding Assessment of Charged Amino Acids on TiO₂ Anatase (101)

Abstract

Atomistic simulations can complement the scarce experimental data on free energies of molecules at bio-inorganic interfaces. In molecular simulations, adsorption free energy landscapes are efficiently explored with advanced sampling methods, but classical dynamics is unable to capture charge transfer and polarization at the solid–liquid interface. Ab initio simulations do not suffer from this flaw, but only at the expense of an overwhelming computational cost. Here, we introduce a protocol for adsorption free energy calculations that improves sampling on the timescales relevant to ab initio simulations. As a case study, we calculate adsorption free energies of the charged amino acids Lysine and Aspartate on the fully hydrated anatase (101) TiO₂ surface using tight-binding forces. We find that the first-principle description of the system significantly contributes to the adsorption free energies, which is overlooked by calculations with previous methods.

1. Introduction

Engineered inorganic materials are important in many technological applications [1,2,3] such as biomimetics [4], optics [5], biosensors [6,7], and smart surfaces [8]. The key to harness the true potentials of these and other nanobio applications lies in a microscopic understanding of biomolecular adsorption on inorganic materials. Nanotoxicity is another area where molecular interactions determine the effects of nanoparticles on biological organisms, and controls the outcome in terms of nanomaterial safety [9,10]. The bionano region, the nanometer-thick boundary near the surfaces of nanoparticles and engineered nanomaterials, is believed to regulate adsorption behavior [9,11]. This region can be probed with modern atomistic simulations, but is difficult to access in experiments due to the weak signal generated by its comparatively small volume.

Titanium dioxide (TiO₂) is a biocompatible semiconductor used in implants and biomedical applications [12,13,14] with a low bandgap that is ideal for water-splitting applications [15]. TiO₂ has become the standard surface model for water interactions with metal oxides in theoretical chemistry due its perceived simplicity, with reactive sites at undercoordinated titanium atoms (Ti${}_{nc}$) and bridging oxygen atoms (O${}_{br}$) [16,17]. However, on this “simple” surface, the adsorption behavior is modified by surface hydroxyl groups from spontaneous water splitting [18,19,20], and adsorbate interactions are indirectly mediated via strongly adsorbed water layers [21,22,23]. Charge transfer and polarization at the surface are other important factors to consider [24].

The binding free energy of a biomolecule on nano-TiO₂, $\mathsf{\Delta }{F}_{\mathrm{ads}}$, dictates the adsorption behavior in equilibrium. This quantity can be measured experimentally [25] and calculated from simulations, but the latter demands careful evaluation of all energy and entropy contributions at the solid–liquid interface. The Lennard–Jones models with fixed partial Coulomb charges employed in classical molecular dynamics (CMD) simulations may not be sensitive to subtle differences in adsorption behavior between similar biomolecules on, e.g., TiO₂. Certainly, much effort has been devoted with such force fields to compute trends for amino acids on various TiO₂ surfaces [26,27,28,29], amorphous nanoparticles [30], and other inorganic materials [31]. The overall assessment is that amino acids bind strongly if they can penetrate the first water layer at the solid–liquid TiO₂ interface [32], but many questions remain unanswered. For example, the force fields used to model TiO₂–adsorbate interactions are parameterized on small and neutral molecules, where polarization and charge transfer take a small role. Further, there are significant differences in the water density profiles near TiO₂ obtained with density functional theory (DFT) [32] and CMD, respectively.

Overestimated surface/water interactions lead to artificial order at the interface, which prevents direct molecular adsorption [29]. Underestimation, on the other hand, smooths out the distinctive adsorption behavior of individual biomolecules [26,27,28]. These differences in force field parameterizations sometimes result in different adsorption behavior for the same molecules [26,27,28,29]. In the absence of experimental evidence, ab initio simulations that account for reactivity and polarization at the solid–liquid interface can resolve these differences. Ab initio molecular dynamics simulations (AIMD), however, are limited to tens of ps of simulation time, which renders proper sampling impossible.

Here, we validate a simulation protocol that substantially improves sampling in free energy calculations with ab initio dynamics. This method can be used to calculate accurate adsorption free energies of small (even charged) molecules on inorganic (and other) surfaces. This approach is based on semi-empirical Tight-Binding (TB) forces, which captures reactivity, polarization, and charge transfer [33], and can be used for orders-of-magnitude longer simulations compared to DFT with the Generalized-Gradient Approximation (GGA) for the forces [34,35]. We use multi-walker metadynamics [36] with an augmented collective variable (CV) that significantly improves sampling of the target CV compared to standard metadynamics. Furthermore, we reconstruct free energies using a “mean force estimator”, which is superior to the traditional way of cumulatively summing the bias potentials [37].

We show that this combination of advanced simulation techniques enables the calculation of ab initio-based, converged, free energy profiles for small molecules on inorganic surfaces. As a case study, we calculate the adsorption free energies of the charged amino acids Lysine (Lys) and Aspartate (Asp) on fully hydrated TiO₂ anatase (101) surfaces. These molecules have been shown to interact strongly with TiO₂ surfaces in previous single-point DFT calculations, with adsorption energies that failed to be explained with classical models [23,33,38].

2. Methods

2.1. System Preparations

The ab initio molecular dynamics (AIMD) simulations using self-consistent density functional tight-binding (DFTB) approach [39] were done with Cp2K [40] with the same setup as in [32]. Briefly, we used PACKMOL [41] to prepare anatase (101) simulation boxes with sizes $10.35\times 11.4\times 43$ Å. The TiO₂ slabs were built by repeating the unit cell four times along the z-direction (same as in our previous simulations [32]), and the rest of the boxes were filled with one amino acid and water amounts that correspond to 1 atm and 310 K [32]. These simulation boxes are large enough to avoid self-interactions. We used analog molecules for Aspartate and Lysine, i.e., side chains of the amino acids with the backbone replaced with a CH₃-group. The removal of the backbone mimics the state of the amino acid in a protein, where the backbone is buried and prevented to contribute to the adsorption except at the terminal groups. We kept the systems neutral with a counterion (OH⁻ or OH₃⁺) initially placed on the opposite surface slab to the amino acid.

Tight-binding (TB) calculations with amino acids have been reported on dry TiO₂ surfaces [33] but, to the best of the authors’ knowledge, not at full hydration. We therefore first tested the Matsci [34] and Mio-Tiorg [35] TB parameterizations for lysine and aspartate at TiO₂–water interfaces. In both cases, we found that the amino acid $C-H$- and $N-H$-groups deprotonated on the TiO₂ anatase (101) surface. There is no physical basis for such behavior, which implies that the overlap integrals of the underlying interactions are underestimated. We added harmonic constraints on the $C-H$- and $N-H$-groups as a simple remedy in the rest of the simulations. Further, the test simulations revealed geometrical distortions at the solid–liquid interface with Mio-Tiorg parameters, but structures were in good agreement to DFT-derived geometries with the Matsci parameters [32]. Based on these results, we used the Matsci parameterization, without long-range dispersions, for all free energy calculations in this work. Furthermore, the final adsorbing modes of Lys and Asp were fully optimized both with TB and DFT approaches [32] to check the consistency using two different levels of theory. For the DFT calculations, we used the BLYP functional [42,43] augmented with the Grimme DFT-D3 dispersion corrections [44], whereas the GTH normconserving pseudopotentials [45,46] and a double-$\zeta$ Gaussian basis set with polarization functions (DZVP) [47] were used to describe the core and valence electrons respectively. The energy cut-off for the electron density expansion in the GPW method was 400 Ry and the minimization was stopped when the total forces were lower than ${10}^{-3}$ Hartree/Bohr.

2.2. Metadynamics Setup

Metadynamics simulations were carried out with the PLUMED [48] module in CP2K. Specifically, we used well-tempered metadynamics with adaptive Gaussians (AWTMetaD) [49,50,51] and a bias factor of 15. We started with Gaussian heights of $3.5$ $\mathrm{kJ}$mol⁻¹, added new bias potentials every 25 fs, and adjusted the widths of the Gaussian every 75 fs. This parameter combination is a reasonable compromise of accuracy vs. speed for ab initio dynamics with limited time sampling [52].

We used the surface separation distance (SSD) as the target collective variable (CV) in the free energy calculations. The SSD is the distance between the outermost layer of Ti surface atoms and the center of mass of a group of concern in the amino acid, NH₃⁺- and COO⁻-groups respectively for Lys and Asp. In addition, we augmented the calculations with a second collective variable with the aim of boosting the exploration of the adsorption landscape along another important degree of freedom—the adsorbent orientation—thus improving the sampling of the target CV [53]. As augmenting variable, we used the angle between the vector of the surface normal and N−C_fn (in Lys) or C_OO^-−C_fn (in Asp), where C_fn is the first carbon neighbor to the atom in question. The augmented variable is integrated out during the calculation of the free energy (see Equation (2)). We restricted the phase space with a wall potential 1 nm away from the outermost Ti atoms, and launched multi-walker [36] AWTMetaD with eight replicas starting at different CV values. The walkers communicated every 25 fs and we simulated 240–300 ps for each walker.

Traditionally, the potential of mean force (PMFs) is reconstructed from the accumulated bias potentials [29,49,51]

$$W_2(z,\theta )=-\underset{t\to \infty }{lim}\left(V(z,\theta ,t)+k_{B}Tln\left(n(z,\theta ,t)\right)+const\right)$$

(1)

where $V(z,\theta ,t)$ is the accumulated bias potential, $\theta$ is the augmenting CV, and $n(z,\theta ,t)$ is the accumulated histogram of the target collective variable. This correction term from the histogram is needed when the widths of the Gaussians bias are adjusted dynamically [51]. Integrating out the augmenting collective variable $\theta$ yields to the potential of mean force along the main collective variable [53]:

$$PMF\left(z\right)=-k_{B}Tln\int e^{-W_2(z,\theta )/k_{B}T}d\theta +const$$

(2)

Forces can be used directly in the calculations to improve convergence in free energy profiles from metadynamics [37]. With this in mind, we also calculated free energy profiles by thermodynamic integration of $F\left(z\right)$, the mean force on the NH₃⁺- and COO⁻-groups in the metadynamics simulations (more exactly, from the derivative of the energy over the SSD collective variable, not including the bias energy). In this formulation the PMF is calculated as

$${\mathrm{PMF}}_{\mathrm{TI}}\left(z\right)=-{\int }_{r_c}^z\langle F\left(z^{\prime }\right)\rangle \mathrm{d}{z}^{\prime }+\mathrm{const},$$

(3)

where z spans the SSD-values from $r_c$ (the onset of the solid surface) to 1 nm, where the potential wall is set. Canonical averaging of the average force is evaluated from all sampled configurations “i” having the main collective variable z in the range $z-\mathsf{\Delta }z/2<z_i<z+\mathsf{\Delta }z/2$:

$$\langle F\left(z\right)\rangle =\frac{{\sum }_{i}F\left(i\right)exp\left(-V(z;{\theta }_i,t)/k_{B}T\right)}{{\sum }_{i}exp\left(-V(z;{\theta }_i,t)/k_{B}T\right)}$$

(4)

In our computations, bin size for force integration was set to $\mathsf{\Delta }z=0.005$ Å. We will refer to this method as MetaDF (metadynamics with force integration) in the rest of this text. Note that MetaDF is not bound to a specific implementation of the metadynamics. It can use well-tempered metadynamics, or constant Gaussian height metadynamics, as it follows from Equations (3) and (4).

The binding free energy is computed from the PMF as [29]

$$\mathsf{\Delta }{F}_{\mathrm{ads}}=-k_{B}Tln\left(\frac{1}{\delta }{\int }_{r_c}^{r_c+\delta }{e}^{-\mathrm{PMF}\left(z\right)/k_{B}T}\mathrm{d}z\right),$$

(5)

where $k_{B}T$ is the product of the Boltzmann constant and the absolute temperature, $\delta =8\AA$ is the thickness of the adsorption layer and $r_c+\delta$ is the start of the liquid bulk. Equation (5) is the thermodynamically correct route to $\mathsf{\Delta }{F}_{\mathrm{ads}}$, but earlier work has also quantified the adsorbate’s binding strength by the difference of the minimum and bulk values of its PMF [26,27,28]

$$\mathsf{\Delta }{F}_{\mathrm{diff}}=\mathrm{PMF}(r_c+\delta )-\mathrm{PMF}\left(r_c\right)$$

(6)

This difference is always larger than $\mathsf{\Delta }{F}_{\mathrm{ads}}$. Note that although $\mathsf{\Delta }{F}_{\mathrm{diff}}$ depends on the specific choice of the CV (e.g., determined by the molecular center of mass or by a specific atomic group on a sorbent molecule), $\mathsf{\Delta }{F}_{\mathrm{ads}}$ does not depend on such choice.

3. Results and Discussion

3.1. Method Validation

To validate the new protocol described so far and the parameters chosen for ab initio Metadynamics, we run extended classical molecular dynamics (CMD) simulations for the lysine–anatase (101) system. The force field describing the interactions was taken from reference [29]. We run classical simulations with the same system size and the same MetaD parameters as we used in ab initio TB computations, but we extended classical simulations up to 200 ns using only the SSD variable for sampling the free energy and comparing the effect of having a single or eight walkers. We tested the effect of the Gaussian height and insertion rate on the PMF convergence. Eight different Metadynamics simulations were run for the single walker simulations with different combinations of these parameters. Gaussians heights and insertion rate were spanning 1 to 3.5 kJ/mol and 25 to 500 fs respectively. In each case for 200 ns of classical simulations no significant difference in the PMF profiles obtained at different Gaussian heights and insertion rates were noticed.

In Figure 1, the convergence test is shown for MetaDF compared to AWTMetaD for Gaussians height $=3.5\mathrm{kJ}/\mathrm{mol}$ and insertion rate $=50$ fs. Two adsorption modes can be distinguished: one $M_{med}$ where the Lysine-TiO₂ interaction is mediated by water at ($\mathrm{SSD}=5$ Å), and a second one $M_{dir}$ where a direct surface contact occurs on the $O_{\mathrm{br}}$ at $\mathrm{SSD}=3.2$ Å. Both MetaDF and AWTMetaD converge to the same PMF profile after 200 ns as expected. The main difference is that although MetaDF needs barely 1 ns to reach the final profile on the $M_{med}$ mode, AWTMetaD requires 30 ns. This effect is even more pronounced for the second adsorption mode $M_{dir}$ where the AWTMetaD convergence is strongly hindered by the long time sampling of the transition between direct-mediated adsorption modes and the cumulative behavior of the bias potential which changes strongly if the CV falls in a previously poorly sampled region. On the contrary, MetaDF converges once the main regions of the free energy landscape have been scanned due to the fact that the average force acting on the CV does not change significantly with long sampling. When eight walkers are implemented (Figure 1b, bottom) the final PMF profile for the $M_{med}$ adsorption mode is reached after barely 300 ps with MetaDF due to the fact that each walker push the remaining walkers to explore other CV values, thus the whole range of CV is sampled simultaneously. The second mode $M_{dir}$ is also sampled faster than with a single walker and results, after 300 ps, in a profile rather close to that obtained in 200 ns calculations. This set-up shows that the convergence can be obtained at least 100 times faster than using a single walker with AWTMetaD. Finally, we note that the adsorption free energy value of Lysine arising uniquely from the $M_{med}$ mode is $F_{\mathrm{ads}}=-2.7\mathrm{kJ}$mol⁻¹, indicating that in the classical description Lysine adsorbs very weakly on the anatase (101) surface.

3.2. Tight-Binding Results

Charged amino acids are driving adsorption at bio-inorganic interfaces [23,38,54,55], and are the most challenging to model since they can induce strong polarization at interfaces and change the surface’s protonation state. In this work, we used Lysine and Aspartate ($+1$ and $-1$ charge, respectively) to investigate how charged amino acids impact the adsorption free energy in bionano systems.

Figure 2c shows the eight walkers exploring the SSD during the Lysine TB simulation. By analyzing contact configurations, we identified two adsorption modes, corresponding to double ($L_I$) and single ($L_{II}$) adsorption on oxygen bridges ($O_{\mathrm{br}}$) on the TiO₂ surface) (Figure 2a,b). The real strength of the multi-walker method is not in the brute force sampling itself, but how walkers in different regions of the free energy landscape communicate so that, although an individual walker may sample a limited area of the configurational space, all walkers together sample the whole relevant region. This effect is extremely important if combined with the force estimator for the calculation of the PMF, as explained in the validation section for CMD. Multi-walker metadynamics also yields linear scaling to reconstruct the free energy landscape with required precision, whereas single-walker metadynamics is limited by slow diffusion [36]. In the case of adsorption at the TiO₂–water interface, individual walkers penetrate into the elusive adsorbed water layer adjacent to the surface. This region is extremely difficult to sample under normal circumstances [56] (which is also illustrated by our CMD simulations described in the previous section), but multiple walkers solve this problem. Figure 2c shows that several adsorption/desorption events are sampled by different walkers during the calculation, as necessary to obtain proper equilibrium statistics. Figure 2d shows the cumulative SSD distribution from all walkers. The peak in the SSD histogram associated to adsorption modes $L_I$ and $L_{II}$ is significantly pronounced after 300 ps of simulation per walker, and is still accentuated when calculated on the last 50 ps of the simulation trajectory. A flat profile in the CV histogram is the signature of a converged free energy profile when calculated with the standard estimator (Equation (2)). This emphasizes how slow PMF convergence can be in standard metadynamics compared to integrating forces (Equation (3)), which does not depend on the accumulated histogram.

Figure 3 shows potentials of mean force (PMFs) along the SSD for Lysine on anatase (101) calculated with MetaDF and AWTMetaD. $\theta$ is also completely sampled (see Figure S1 in the Supporting Information for the 2D map), and it shows a single global minima at the adsorption site. The PMFs are plotted at different times up to 300 ps per walker and show the convergence behaviors of the two methods. The PMFs calculated with the standard estimator (Equation (2)) fluctuate strongly, which is a manifestation of slow convergence after 300 ps of simulated time (Figure 2d and Figure 3d). The PMFs calculated by integrating forces on the NH₃⁺-group converges to a smooth profile with an error of 3 $\mathrm{kJ}$mol⁻¹ (as estimated from the force variance [57]) after 250 ps per walker. This is visible also from Figure 3d, where $\mathsf{\Delta }{F}_{\mathrm{ads}}$ is plotted as a function of time for MetaDF and AWTMetaD. Via Equation (5), we calculate that the binding free energy for Lysine on anatase (101) is $\mathsf{\Delta }{F}_{\mathrm{ads}}=-58.2\mathrm{kJ}$mol⁻¹ (or $\mathsf{\Delta }{F}_{\mathrm{diff}}=-65.4\mathrm{kJ}$mol⁻¹ from Equation (6)). This is substantially larger (3–10 times) than values reported from different CMD simulations [26,27,28,29], where both direct and indirect adsorption modes were found. This difference suggests that electronic degrees of freedom, in particular polarization and charge transfer effects at the interface, cannot be neglected in the free energy calculations.

To further validate the TB parametrization used in our calculations, we performed optimization computations for the $L_I$ mode bound state using TB and DFT-GGA theory. The optimized bond length between the hydrogen of the amino group and the bridging oxygen ( NH−O_br, Figure 2a) was found to be $d_{\mathrm{NH}-O_{\mathrm{br}}}=1.7\AA$ for DFT-GGA and $d_{\mathrm{NH}-O_{\mathrm{br}}}=1.6\AA$ for TB. No further discrepancies or adjustment in the final adsorption configuration were noticed.

For Aspartate on anatase (101), we also found two adsorption modes: $A_I$, which is not in direct contact with the surface but mediated via the first water layer, and $A_{II}$, which is in direct contact with the TiO₂ surface (Figure 4a,b). Figure 5 shows the difference in PMFs calculated with MetaDF compared to the standard free energy estimator. Mode $A_I$ is well-sampled in 180 ps with MetaDF (within an error of 3 $\mathrm{kJ}$mol⁻¹ as estimated from the force variance [57] and inspecting $F_{\mathrm{ads}}$ as a function of time as shown in Figure 5d), while AWTMetaD (Equation (2)) fails to converge after 300 ps of simulation time (Figure 4d and Figure 5d). The $A_{II}$-mode (direct contact) appears when a single walker penetrates the first surface water layer after 60 ps. The $A_{II}$-mode is separated from the $A_I$-mode by a free energy barrier of $\sim 20\mathrm{kJ}$mol⁻¹, and thus much more challenging to sample. The high barrier and narrow region available to $A_{II}$ implies that the main contribution to the binding free energy is coming from $A_I$. The binding free energy is $\mathsf{\Delta }{F}_{\mathrm{ads}}=12.1\mathrm{kJ}$mol⁻¹ from Equation (5) (or $17.2\mathrm{kJ}$mol⁻¹ via Equation (6)). CMD simulations have reported similar values of the binding free energy when the COO⁻-group is in direct contact to the surface [58] and via surface hydroxyls [59]. A water-mediated adsorption mode has not been found due to the weak nature of this interaction in classical models. The optimized bond length CO−H₂O (Figure 4a) for the $A_I$-mode is $d_{\begin{array}{c}\hfill \mathrm{CO}-{\mathrm{H}}_2\mathrm{O}\end{array}}=1.8\AA$ and it coincides if calculated with TB and DFT-GGA. The present work emphasizes the importance of an atomistic model of the solid–liquid interface that simultaneously reproduces the interfacial water structure/reactivity and the underlying quantum nature of semiconductor materials with electronic correlation effects such as polarization and charge transfer.

The present work has showed that the force estimator, Equation (3), is superior to the traditional Equation (2) for calculating PMFs with metadynamics. Although advantages of the force estimator have been discussed previously [37], these becomes critically important in quantum–chemical simulations with limited sampling time. The slow convergence of the PMFs with the accumulated bias potential as an estimator of the free energy is due to that AWTMetaD is diffusion-limited. Standard thermodynamic integration (including umbrella sampling), on the other hand, is inefficient when the important regions are unknown prior to the simulation, as is usually the case. MetaDF improves sampling by combining these techniques so that the bias potentials provide nearly uniform sampling along the collective variable and create an optimal configuration set for the thermodynamic integration. The multiple walkers boost the sampling of the relevant regions of the CV thus representing the optimal tool for the force estimator. As well-tempered MetaD can suffer from poor sampling if the Gaussian heights became too low before the bias potential reaches the optimal shape, using fixed bias for the metadynamics could eventually speed up even more the convergence. Furthermore, the key point of MetaDF is to collect the forces in the whole relevant region of the free energy landscape, but not to provide the bias potential which exactly compensates the free energy profile. The augmenting collective variable (orientation of the adsorbate molecule) and the multiple walkers further improve sampling of orthogonal (“slow” or “hidden”) degrees of freedom. The mean force as a function of the CV is not affected by the bias, so the forces in a small region (bin) of the CV can be averaged in time. Gaussian bias potentials are inserted every 25 fs, but forces are collected at each time step. Our case study estimates that once the bias potential has explored all values of the CV in the range of interest, force integration converges at least 100 times faster than AWTMetaD with the standard free energy estimator. The MetaDF method can therefore be used in all situations where sampling is limited by diffusion or by strongly bound states.

4. Conclusions

There is a pressing need for simulation protocols that target strongly adsorbing molecules, e.g., amino acids on TiO₂ surfaces, so that adsorption free energies can be accurately determined for such situations. MetaDF is a sampling method that combines metadynamics and thermodynamic integration to accelerate the convergence of PMF calculations. Multiple walkers and augmented collective variables further improves sampling to a point where free energies can be determined with high accuracy, even in cases with strong adsorption. We tested two Tight-Binding parameterizations for adsorption of the charged Lysine and Aspartate amino acids on the anatase (101) surface. For both molecules, we found large adsorption free energies compared to previous CMD studies. We hypothesize that the quantum nature of these systems strongly influence the adsorption behavior due to polarization and charge transfer at the interface. For the present case, the PMFs converge within 180 to 250 ps (per walker) in MetaDF simulations with eight walkers. This time scale is accessible to large-scale ab initio molecular dynamics simulations with GGA-level density functional theory, which opens the possibility to move beyond the simplifications of tight-binding DFT and calculate adsorption PMFs at bio-inorganic interfaces with full electronic treatment. MetaDF is particularly useful in situations with limited sampling, such as ab initio simulations, or when diffusion is hindered, e.g., by bound states or barriers along hidden degrees of freedom. Hopefully, MetaDF will be extensively used incoming systematic investigations of various kinds of bionano interfaces.