Ab-initio force prediction for single molecule force spectroscopy made simple

1. Problem Statement

Single Molecule Force Spectroscopy (SMFS) is a primary experimental technique for studying mechanochemistry, where mechanical forces induce bond rupture or chemical reactions. While experimental force-extension curves can be reproduced using ab-initio calculations, a significant gap exists in predicting the specific rupture forces measured in experiments directly from first principles.

• The Discrepancy: Previous computational strategies, specifically the COGEF (Constrained Geometry Simulates Forces) method, calculate the maximal force a bond can withstand ($F_{max}$). However, these calculations typically overestimate experimental rupture forces by an order of magnitude.
• The Root Cause: The overestimation arises because COGEF calculations are often performed at 0 K, neglecting thermal fluctuations. In reality, bond rupture is a stochastic process driven by temperature fluctuations that allow the system to overcome energy barriers before the maximal theoretical force is reached.
• The Challenge: There is a need for a closed-form analytical expression that connects ab-initio calculated bond properties (dissociation energy and maximal force) with experimental parameters (temperature and loading rate) to accurately predict the most probable rupture force.

2. Methodology

The authors developed a theoretical framework combining statistical mechanics with ab-initio quantum chemical calculations.

A. Theoretical Framework

1. Probabilistic Model: The bond breaking process is treated as a probabilistic event (following Evans and Ritchie). The probability of the bond remaining intact ($P$) decreases over time as the external force increases.
2. Force-Dependent Barrier: The authors propose that the energy barrier to rupture, $\Delta H^\ddagger(F)$, under an external force $F$ can be approximated by a quadratic form:
   $$\Delta H^\ddagger(f) = \Delta U^\ddagger (1 - f)^2$$
   Where:
   • $\Delta U^\ddagger$ is the dissociation energy (barrier height at zero force).
   • $f = F/F_{max}$ is the normalized force.
   • $F_{max}$ is the maximal force the bond can withstand (the point where the barrier vanishes).
3. Derivation of Most Probable Force ($F^*$):
   • Assuming a constant loading rate ($\alpha = dF/dt$) and neglecting bond reformation (valid for strong bonds), the authors derived the distribution of rupture forces.
   • They solved for the peak of this distribution (the most probable force) using the Lambert W function, resulting in a closed-form expression (Eq. 11):
     $$f^* = 1 - \sqrt{\frac{1}{2\beta\Delta U^\ddagger} W\left( \frac{1}{2\beta\Delta U^\ddagger} \frac{t_{max}^2}{h^2 \beta^2} \right)}$$
   • This expression accounts for the interplay between the barrier height, thermal energy ($k_B T$), and the time scale of force application ($t_{max}$).
   • Note: In the low-force limit, this expression reduces to the famous Bell model, but the quadratic form provides higher accuracy across a wider range of loading rates.

B. Computational Strategy (Ab-initio)

To obtain the two key parameters ($\Delta U^\ddagger$ and $F_{max}$) required for the prediction:

1. Dissociation Energy ($\Delta U^\ddagger$): Calculated using Nudged Elastic Band (NEB) methods to find the transition state in the absence of force. This ensures the correct thermal activation barrier is used, avoiding artifacts from force-induced pathway changes (e.g., violating Woodward-Hoffmann rules).
2. Maximal Force ($F_{max}$): Calculated using the COGEF method. This involves constraining the geometry of the molecule to simulate stretching until the bond ruptures, identifying the peak force.
3. Software: Calculations were performed using Density Functional Theory (DFT) with the GPAW package (PBE functional, PAW method).

3. Key Contributions

• Unified Predictive Model: The paper establishes a simple, closed-form analytical equation (Eq. 11) that predicts experimental rupture forces using only two ab-initio parameters ($\Delta U^\ddagger$ and $F_{max}$) and experimental conditions ($T$ and $\alpha$).
• Resolution of the COGEF Discrepancy: The authors demonstrate that the large overestimation in previous COGEF studies is due to the neglect of thermal activation. By incorporating the loading rate and temperature into the statistical model, the predictions align with experimental data.
• Validation of the Quadratic Barrier: The study validates that the force-dependent barrier for many mechanochemical systems follows a quadratic dependence on force, justifying the use of the derived analytical expression over more complex numerical integrations for many cases.
• General Applicability: The method is shown to work for diverse systems, including covalent bond breaking in ring-opening reactions and non-covalent transformations.

4. Results

The authors tested their model against experimental data from the literature for three distinct classes of molecules:

1. Cyclopropanes (CP): Ring-opening reactions involving the rearrangement of multiple bonds.
2. Benzocyclobutenes (BCB): Similar ring-opening mechanisms.
3. Diarylethenes (DAE): Non-covalent transformations with lever-arm effects.

Findings:

• High Accuracy: The predicted most probable forces ($F^*$) showed excellent agreement with experimental measurements across a wide range of forces (0 to 2.5 nN).
• Comparison with Bell Model: While the Bell model (linear barrier approximation) works for low loading rates, the authors' quadratic model (Eq. 11) is significantly more accurate, particularly at higher loading rates where the Bell model fails to capture the distribution width and peak position correctly.
• Error Analysis: The predicted values fell within the experimental error bars and the theoretical uncertainty derived from the probability distribution width.
• Parameter Sensitivity: The model correctly predicted that increasing the loading rate increases the rupture force, while increasing temperature decreases it, matching experimental trends.

5. Significance

• Simplification of Mechanochemical Prediction: The paper provides a "simple" yet rigorous workflow. Researchers no longer need to perform complex, computationally expensive molecular dynamics simulations with explicit solvent or time-dependent pulling to predict rupture forces.
• Design of Mechanophores: This method enables the rational design of new mechanophores (molecules that respond to force) by allowing chemists to predict rupture forces directly from the molecular structure using standard DFT tools.
• Bridging Theory and Experiment: It successfully bridges the gap between static ab-initio calculations and dynamic experimental observations, clarifying that the "missing link" was the statistical treatment of thermal fluctuations combined with the specific loading rate.
• Limitations: The authors note that the quadratic approximation may break down for extremely complex barriers (e.g., retro-Diels-Alder reactions) that cannot be described by a simple cusp-like potential, but it is highly effective for the majority of standard mechanochemical bond ruptures.