Stability and Vibrations of Proteins in Vacuum and Water: Bridging Quantum Accuracy and Force-Field Efficiency

TL;DR

The work tackles the challenge of predicting biomolecular vibrations and metastable energies with quantum accuracy at force-field scale. It introduces SO3LR, a general-purpose MLFF trained on $PBE0+MBD$ data, and demonstrates its ability to reproduce DFT-level PES, vibrational densities of states, and mode eigenvectors for small molecules and complex biomolecular assemblies in vacuum and water. Key findings include near-DFT accuracy for frequencies and IR spectra, faithful capture of anharmonicity and environment-driven effects, and strong quantitative agreement for solvent- and protein-protein interaction-induced shifts, surpassing traditional MMFFs. This approach offers a scalable, transferable, and environment-aware framework for predictive biomolecular dynamics, with potential to enable quantum-accurate simulations of IDPs and non-natural amino acids at reduced computational cost.

Abstract

Predicting biomolecular thermodynamics and spectroscopy requires accurate relative energies of metastable states and local curvatures on the potential-energy surface. We show that the general-purpose SO3LR machine-learned force field (MLFF) reproduces PBE0+MBD density-functional theory with unprecedented fidelity across bio-relevant molecules spanning sizes and complexities far beyond its training dataset. For 23 small molecules, SO3LR captures harmonic and anharmonic vibrational features, including frequencies, displacement patterns, and IR spectra. We perform detailed dynamical studies of the amino acid oF-Phe+, folding of the alanine-15 peptide, and assembly of monomeric p53 transactivation domains into tetramers, in vacuum and water. SO3LR consistently reproduces DFT-level potential-energy surfaces, vibrational densities of states, and mode eigenvectors, capturing anharmonicity, polarization, and medium-range environment-driven interactions crucial for proteins. Our results establish that MLFF-driven dynamics provide a quantum-accurate picture of metastable minima and vibrational properties, achieving DFT-level accuracy at force-field computational cost and opening new possibilities for the computational study of biomolecules.

Figures (5)

• Figure 1: Vibrational benchmark of SO3LR and GAFF2 using DFT at PBE0+MBD as reference calculation.a Aspartame. b NMA frequencies of SO3LR (left) and GAFF2 (right) against PBE0+MBD, and MAE based on 1 to 1 fitting. c–d Absolute values of correlation matrices of vibrational eigenvectors from SO3LR c and GAFF2 d versus PBE0+MBD (eq. \ref{['eq.corr']}). Spread is reported (eq. \ref{['eq.spread']}). For visualisation purposes, a 4x4 max pool filter was applied to the matrix. e Whole set comparison of SO3LR (blue) and GAFF2 (orange) versus PBE0+MBD. Error bars represent the standard deviation per molecule. Molecules: 2-Chloroethanol ($\alpha$), 2-Fluoroethanol ($\beta$), 4-Chlorotoluene ($\gamma$), Acetic acid ($\delta$), Acetone ($\epsilon$), Acetonitrile ($\zeta$), Acetophenone ($\eta$), Aspartame ($\theta$), CO$_2$ ($\iota$), DMSO ($\kappa$), Ethane ($\lambda$), Ethanol ($\mu$), Ethoxyethane ($\nu$), Formaldehyde ($\xi$), Haloperidol ($o$), Hexane ($\pi$), Ibuprofen ($\rho$), Methanol ($\sigma$), Methylcyclohexane ($\tau$), Phenol ($\upsilon$), Tetrahydrofuran ($\phi$), Toluene ($\chi$) and tert-Butylmercaptan ($\psi$). f-i Comparison of IR spectra versus experiment (red) using GAFF2 (orange) and SO3LR (blue) from dipole autocorrelation and from harmonic approximation SO3LR (purple) and PBE0+MBD (green). Experimental infrared spectra of Acetic acid f and Ethoxyethane g are extracted from HITRAN2016_XSCHITRAN2020. Ethane h and Aspartame i are extracted from NISTWebBook2025Aspartame. Source data are provided as a Source Data File.
• Figure 2: Vibrational analysis of L-$o$F-phenylalanine+H+.a Four stable conformational isomers, ordered from most to least stable, used as initial geometries for harmonic approximation and molecular dynamics simulations. b Experimental IR spectrum Safferthal2023 (red) compared with GAFF2 (orange) and SO3LR (blue) spectra obtained from dipole autocorrelation, alongside harmonic approximation results from PBE0+MBD (black)*. c Potential energy surface mapped as a function of $\phi_1$ (N–C$_\alpha$–C$_\beta$–C$_\gamma$) and $\phi_2$ (H$_{\beta^1}$–C$_\beta$–C$_\gamma$–C$_{\sigma^2}$) dihedral angles. Letters indicate conformers defined in panel A. d Four independent dihedral sampling during 100 ps simulations initiated from conformer A at 300K. *SO3LR and PBE0+MBD spectrum used a correction factor of 0.966 as in the reference. Source data are provided as a Source Data File.
• Figure 3: AceAla15NMe analysis across folding structures.a Three-dimensional structures of the four conformers analyzed: canonical $\alpha$-helix, canonical 310-helix, intermediate state, and coil (fully extended). b Eigenvector and frequency analysis comparing SO3LR and MMFFs to PBE0+MBD. Mean standard error of frequencies and the average spread (see Eq. \ref{['eq.spread']}) for the four conformers. Error bars represent standard deviations. c Energy differences (kcal mol-1) of the 310-helix, intermediate state, and coil relative to the $\alpha$-helix. Potential energy and vibrational entropy contributions ($-TS_{Truhlar}$) are computed using the harmonic approximation at 300K and $\upsilon_0$ is 50 cm-1. Source data are provided as a Source Data File.
• Figure 4: High frequency stretching modes in p53 monomer. Normalized spectral profiles of C-H, N-H and O-H modes from 2700 to 4200 cm-1 in p53 monomer (PDBid: 1SAE) Clore1995. Different computational methods are shown: PBE+MBD (dotted black), SO3LR (blue), AMOEBA (orange), FF14SB (green) and CHARMM36m (red). Source data are provided as a Source Data File.
• Figure 5: Environmental effects on stretching modes in p53 multimeric protein. C–H (green), N–H (blue), and O–H (red) stretching modes in the 2700-4200 cm-1 region are studied under different conditions. a Vibrational shifts between dimer and monomer configurations (x‑axis) for the capped fragment of residues 327–333 (green), computed using PBE0+MBD, SO3LR, AMOEBA, and FF14SB. b Vibrational shifts plotted against bond‑length differences between solvated tetramer and monomer configurations, including the linear correlation $R^2$c Potential energy profiles for bond‑length perturbations around equilibrium for the N–H stretching mode of the Ile332 backbone, shown for the monomer (blue) and tetramer (orange). d Energy comparison for the same atoms as a function of scaled bond distance (0.9–1.3× equilibrium) in the monomer. The shaded region represents the area predicted by a fully harmonic potential from 1.0 to 1.3× the equilibrium distance (dotted). e Normalized area distributions (scaled to 1 at maximum) as an anharmonicity measure comparing SO3LR and MMFFs, grouped by bond type. In scatter plots, each point corresponds to the most similar mode matched by the atom index. Source data are provided as a Source Data File.