Molecular Mechanics Force Field Research Articles

Multiscale modeling and simulation of surface‐enhanced spectroscopy and plasmonic photocatalysis

AbstractPlasmonic metal nanoparticles (PMNPs) are capable of localized surface plasmon resonance (LSPR) and have become an important component in many experimental settings, such as the surface‐enhanced spectroscopy and plasmonic photocatalysts, in which PMNPs are used to regulate the nearby molecular photophysical and photochemical behaviors by means of the complex interplay between the plasmon and molecular quantum transitions. Building computational models of these coupled plasmon‐molecule systems can help us better understand the bound molecular properties and reactivity, and make better decisions to design and control such systems. Ab initio modeling the nanosystem remains highly challenging. Many hybrid quantum‐classical (or ‐quantum) computing models have thus been developed to model the coupled systems, in which the molecular system of interest is designated as the quantum mechanical (QM) sub‐region and treated by the excited‐state electronic structure approaches such as the time‐dependent density functional theory (TDDFT), while the electromagnetic response of PMNPs is usually described using either a computational/classical electrodynamic (CED) model, polarizable continuum model(PCM), a polarizable molecular mechanics (MM) force field, or a collective of optical oscillators in QED model, leading to many hybrid approaches, such as QM/CED, QM/PCM, QM/MM or ab initio QED. In this review, we summarize recent advances in the development of these hybrid models as well as their advantages and limitations, with a specific emphasis on the TDDFT‐based approaches. Some numerical simulations on the plasmon‐enhanced absorption and Raman spectroscopy, plasmon‐driven water splitting reaction and interfacial electronic injection dynamics in dye‐sensitized solar cell are demonstrated.This article is categorized under: Electronic Structure Theory > Ab Initio Electronic Structure Methods Theoretical and Physical Chemistry > Spectroscopy Software > Quantum Chemistry Electronic Structure Theory > Combined QM/MM Methods

GAFF-AIC: reoptimisation of the GAFF force field for realistic densities and viscosities in aromatic isocyanates

ABSTRACT Molecular mechanical force field methods, which help to investigate material properties beyond molecular scales, describe the intra- and intermolecular interactions of molecules assuming a classical force field comprised of several terms. Force fields should accurately describe the energy interactions of all moleculer elements to, eventually, accurately reproduce macromolecular properties of that materials. In the case of aromatic isocyanates, general purpose force fields failed to properly account for molecular interactions and, consequently, at the prediction of macromolecular properties of isocyanate-based materials. The present study presents a new force field based on the General AMBER Force Field (GAFF), namely GAFF-AIC, which shows improved results for the prediction of macromolecular properties, in particular density and viscosities, of materials based on aromatic isocyanates. For that, a new set of parameters for intramolecular and intermolecular interactions has been optimised by means of a least square fitting. It is shown that GAFF-AIC is able to accurately reproduce the experimental density and viscosity of liquid 1-naphtyl isocyanate. Moreover, these newly proposed force field parameters were evaluated to see how efficient they were at predicting the liquid densities and viscosities for several considered aromatic isocyanates, resulting in an accurate estimation that is in excellent agreement with experimental results.

On hidden anisotropy of formally charged fragments.

The proper and precise reproduction of the molecular electrostatic potential (MEP) is crucial to describe correctly electrostatic interactions in molecular modeling. Most of the classical molecular mechanics force fields for biomolecules and drug-like molecules use the atom-centered point charges to describe MEP. However, it has been systematically pointed out in literature that such an approximation is not always enough, and some groups, like amino group or heavy halogens, require the use of anisotropic model for better description of their MEP. At the same time, the formally charged groups have not been as extensively and systematically studied as their neutral counterparts. In this report, we demonstrate that the anisotropic models for formally charged groups do bring improvements in the reference MEP reproduction, that are comparable in magnitude to those for neutral groups.

Integrated Multiscale Multilevel Approach to Open Shell Molecular Systems.

We present a novel multiscale approach to study the electronic structure of open shell molecular systems embedded in an external environment. The method is based on the coupling of multilevel Hartree-Fock (MLHF) and Density Functional Theory (MLDFT), suitably extended to the unrestricted formalism, to Molecular Mechanics (MM) force fields (FF). Within the ML region, the system is divided into active and inactive parts, thus describing the most relevant interactions (electrostatic, polarization, and Pauli repulsion) at the quantum level. The surrounding MM part, which is formulated in terms of nonpolarizable or polarizable FFs, permits a physically consistent treatment of long-range electrostatics and polarization effects. The approach is extended to the calculation of hyperfine coupling constants and applied to selected nitroxyl radicals in an aqueous solution.

CHARMM-GUI QM/MM interfacer for the hybrid QM/MM molecular dynamics simulations.

Molecular mechanical (MM) force fields are constructed based on the empirical potentials, describing intramolecular and intermolecular potentials. Molecular dynamics (MD) simulation utilizing MM force fields is a robust tool to investigate large biomolecular systems, but the changes in electronic structures during a chemical process are overlooked in classical mechanics. Quantum mechanical (QM) treatments are needed for such changes including bond forming/breaking and charge transfers, but a high computational cost limits its usages with up to hundreds of atoms or less. A hybrid QM/MM approach, where the reactive part of the system is treated with QM and the remainder by cheap MM potentials, can effectively capture the reaction without significantly increasing the computational cost. For the different nature of the QM versus MM theories, however, the setup of biomolecular QM/MM simulation is known to be time consuming even to QM/MM experts. We present CHARMM-GUI QM/MM Interfacer that prepares QM/MM simulations with a user-selected QM region and various QM theories in multiple simulation packages. A non-trivial and error-prone QM/MM simulation setup involving link-atom insertion and charge distribution can be seamlessly prepared in QM/MM Interfacer. A few representative systems are simulated to delineate the robustness of QM/MM Interfacer. We hope both QM/MM beginners and experts to find the module's quick and accurate QM/MM simulation setup useful for their important applications.

Improved polarizable force field for modeling Mg2+ binding to proteins and nucleotides.

Magnesium (Mg2+) is an indispensable regulatory cation for numerous cellular enzymes and membrane proteins that participate in a variety of biological processes, including cell signaling, glycolysis, cell differentiation, and genome stability. Molecular mechanics (MM) simulations have the potential to provide detailed insights into the mechanisms of these Mg-dependent proteins. However, MM force fields continue to suffer from large inaccuracies when dealing with divalent cations. The inclusion of explicit polarization in MM models, such as in CHARMM Drude and AMOEBA models, has been promising, but despite substantial improvements over non-polarizable models, errors in interactions of proteins with divalent cations remain in excess of 10 kcal/mol. Such errors are present even with ‘nonbonded-fix (NB-fix)’ corrections. Modeling becomes even more problematic when dealing with enzymes like kinases and phosphatases, where Mg2+ ions interact directly with both proteins and nucleotides. Toward addressing this long-standing challenge, we recently presented a recalibration of the AMOEBA model (Delgado et al. J Chem. Info. Model. 2022). Notably, polarization terms of cation-coordinating groups were improved to respond better to the high electric fields present near cations. This systematically reduced errors in interactions of monovalent cations with proteins. Here we show that this revised model, along with our many body NB-fix (MB-NB-fix) corrections, reduces errors in Mg2+-protein interaction energies to 5 kcal/mol. MB-NB-fix corrections also improve predictions of Mg2+-ATP interactions substantially. Errors are computed with respect to estimates from vdW-inclusive density functional theory that we benchmark against "gold standard" CCSD(T) calculations and experiments. We will also present the ramifications of these improvements in local interactions on Mg2+-ATP binding free energies and binding modes in bulk water, and on complexation of Mg2+ and ATP with proteins.

MELD-DNA: Accurate structural prediction of protein-DNA interactions.

Protein-DNA interaction is one of the most important aspects of biology. It is a main step of gene regulation in all organisms and is a potential target for therapeutic and biotech applications. Protein-DNA interactions are distinct from protein-protein interactions due to several factors: 1) DNA sequence and specific interactions (base readout); 2) DNA shape and its deformability (shape readout); 3) Exploration and searching on the DNA for alternative binding states. Such factors dictate affinity and specificity of transcription factors. Traditionally, these interactions were only predicted at the sequence level (e.g. using position weight matrices) and more recently using docking and machine learning models. Here we present MELD-DNA, a too that combines cheap information from experiment (which can be noisy, ambiguous and sparse) with physics-based molecular mechanics force fields using bayesian inference. MELD will accelerate sampling of the binding process by temperature and hamiltonian replica exchanging that is guided by experimental data. We have tested the tool on a diverse set of 15 systems from various TF families and we observed native-like bound structures most of them with room for more improvements. We also compare and benchmark with the state-of-the-art in docking (HADDOCK) and ML (RosettaFold2NA). In this study we discuss several issues in the available methods and techniques that need to be addressed for more accurate predictions. We believe MELD-DNA can be an essential tool in the community from drug design, regulation and fundamental research.

Transferability of the Electrostatic Parameters of the Polarizable Gaussian Multipole Model.

Accuracy and transferability are the two highly desirable properties of molecular mechanical force fields. Compared with the extensively used point-charge additive force fields that apply fixed atom-centered point partial charges to model electrostatic interactions, polarizable force fields are thought to have the advantage of modeling the atomic polarization effects. Previous works have demonstrated the accuracy of the recently developed polarizable Gaussian multipole (pGM) models. In this work, we assessed the transferability of the electrostatic parameters of the pGM models with (pGM-perm) and without (pGM-ind) atomic permanent dipoles in terms of reproducing the electrostatic potentials surrounding molecules/oligomers absent from electrostatic parameterizations. Encouragingly, both the pGM-perm and pGM-ind models show significantly improved transferability than the additive model in the tests (1) from water monomer to water oligomer clusters; (2) across different conformations of amino acid dipeptides and tetrapeptides; (3) from amino acid tetrapeptides to longer polypeptides; and (4) from nucleobase monomers to Watson-Crick base pair dimers and tetramers. Furthermore, we demonstrated that the double-conformation fittings using amino acid tetrapeptides in the αR and β conformations can result in good transferability not only across different tetrapeptide conformations but also from tetrapeptides to polypeptides with lengths ranging from 1 to 20 repetitive residues for both the pGM-ind and pGM-perm models. In addition, the observation that the pGM-ind model has significantly better accuracy and transferability than the point-charge additive model, even though they have an identical number of parameters, strongly suggest the importance of intramolecular polarization effects. In summary, this and previous works together show that the pGM models possess both accuracy and transferability, which are expected to serve as foundations for the development of next-generation polarizable force fields for modeling various polarization-sensitive biological systems and processes.

QMrebind: incorporating quantum mechanical force field reparameterization at the ligand binding site for improved drug-target kinetics through milestoning simulations.

Understanding the interaction of ligands with biomolecules is an integral component of drug discovery and development. Challenges for computing thermodynamic and kinetic quantities for pharmaceutically relevant receptor-ligand complexes include the size and flexibility of the ligands, large-scale conformational rearrangements of the receptor, accurate force field parameters, simulation efficiency, and sufficient sampling associated with rare events. Our recently developed multiscale milestoning simulation approach, SEEKR2 (Simulation Enabled Estimation of Kinetic Rates v.2), has demonstrated success in predicting unbinding (koff) kinetics by employing molecular dynamics (MD) simulations in regions closer to the binding site. The MD region is further subdivided into smaller Voronoi tessellations to improve the simulation efficiency and parallelization. To date, all MD simulations are run using general molecular mechanics (MM) force fields. The accuracy of calculations can be further improved by incorporating quantum mechanical (QM) methods into generating system-specific force fields through reparameterizing ligand partial charges in the bound state. The force field reparameterization process modifies the potential energy landscape of the bimolecular complex, enabling a more accurate representation of the intermolecular interactions and polarization effects at the bound state. We present QMrebind (Quantum Mechanical force field reparameterization at the receptor-ligand binding site), an ORCA-based software that facilitates reparameterizing the potential energy function within the phase space representing the bound state in a receptor-ligand complex. With SEEKR2 koff estimates and experimentally determined kinetic rates, we compare and interpret the receptor-ligand unbinding kinetics obtained using the newly reparameterized force fields for model host-guest systems and HSP90-inhibitor complexes. This method provides an opportunity to achieve higher accuracy in predicting receptor-ligand koff rate constants.

Simulation Reveals the Chameleonic Behavior of Macrocycles.

Conformational analysis is central to the design of bioactive molecules. It is particularly challenging for macrocycles due to noncovalent transannular interactions, steric interactions, and ring strain that are often coupled. Herein, we simulated the conformations of five macrocycles designed to express a progression of increasing complexity in environment-dependent intramolecular interactions and verified the results against NMR measurements in chloroform and dimethyl sulfoxide. Molecular dynamics using an explicit solvent model, but not the Monte Carlo method with implicit solvation, handled both solvents correctly. Refinement of conformations at the ab initio level was fundamental to reproducing the experimental observations─standard state-of-the-art molecular mechanics force fields were insufficient. Our simulations correctly predicted the intramolecular interactions between side chains and the macrocycle and revealed an unprecedented solvent-induced conformational switch of the macrocyclic ring. Our results provide a platform for the rational, prospective design of molecular chameleons that adapt to the properties of the environment.

Molecular mechanics simulations of lattice dynamical properties of the spin crossover complex [Fe(pyrazine)][Ni(CN)4

An empirical molecular mechanics force field has been built for the [Fe(pyrazine)][Ni(CN)4] spin crossover complex using Raman spectroscopic and nuclear inelastic scattering data. Based on this force field, molecular dynamics simulations were conducted to investigate the lattice dynamical properties of the complex in the high-spin and low-spin states. The extracted thermodynamic and mechanical parameters match well with the corresponding experimental data. Most importantly, the force field can successfully reproduce the low-frequency acoustic part of the partial densities of vibrational states, which is vital for the understanding the thermodynamic analysis of the spin crossover phenomenon in crystalline solids and which would be difficult to tackle by quantum mechanical methods.

Enhancing Biomolecular Simulations with Hybrid Potentials Incorporating NMR Data.

Some recent advances in biomolecular simulation and global optimization have used hybrid restraint potentials, where harmonic restraints that penalize conformations inconsistent with experimental data are combined with molecular mechanics force fields. These hybrid potentials can be used to improve the performance of molecular dynamics, structure prediction, energy landscape sampling, and other computational methods that rely on the accuracy of the underlying force field. Here, we develop a hybrid restraint potential based on NapShift, an artificial neural network trained to predict protein nuclear magnetic resonance (NMR) chemical shifts from sequence and structure. In addition to providing accurate predictions of experimental chemical shifts, NapShift is fully differentiable with respect to atomic coordinates, which allows us to use it for structural refinement. By employing NapShift to predict chemical shifts from the protein conformation at each simulation step, we can compute an energy penalty and the corresponding hybrid restraint forces based on the difference between the predicted values and the experimental chemical shifts. The performance of the hybrid restraint potential was benchmarked using both basin-hopping global optimization and molecular dynamics simulations. In each case, the NapShift hybrid potential improved the accuracy, leading to better structure prediction via basin-hopping and increased local stability in molecular dynamics simulations. Our results suggest that neural network hybrid potentials based on NMR observables can enhance a broad range of molecular simulation methods, and the prediction accuracy will improve as more experimental training data become available.

Atomic-Resolution Experimental Structural Biology and Molecular Dynamics Simulations of Hyaluronan and Its Complexes.

This review summarizes the atomic-resolution structural biology of hyaluronan and its complexes available in the Protein Data Bank, as well as published studies of atomic-resolution explicit-solvent molecular dynamics simulations on these and other hyaluronan and hyaluronan-containing systems. Advances in accurate molecular mechanics force fields, simulation methods and software, and computer hardware have supported a recent flourish in such simulations, such that the simulation publications now outnumber the structural biology publications by an order of magnitude. In addition to supplementing the experimental structural biology with computed dynamic and thermodynamic information, the molecular dynamics studies provide a wealth of atomic-resolution information on hyaluronan-containing systems for which there is no atomic-resolution structural biology either available or possible. Examples of these summarized in this review include hyaluronan pairing with other hyaluronan molecules and glycosaminoglycans, with ions, with proteins and peptides, with lipids, and with drugs and drug-like molecules. Despite limitations imposed by present-day computing resources on system size and simulation timescale, atomic-resolution explicit-solvent molecular dynamics simulations have been able to contribute significant insight into hyaluronan's flexibility and capacity for intra- and intermolecular non-covalent interactions.

In Silico Simulations and Functional Cell Studies Evidence Similar Potency and Distinct Binding of Pacific and Caribbean Ciguatoxins

Ciguatoxins (CTX) cause ciguatera poisoning, which is the most common reported human food poisoning related to natural marine toxins. Pacific ciguatoxins are the most abundant and studied CTX analogues; however, the growing distribution of Caribbean analogues and the limited data available on their biological effects make necessary to re-evaluate their relative potency. For decades, the guidelines established by regulatory agencies have assumed that the potency of the Caribbean CTXs were tenfold lower than the Pacific CTXs. We present here an integrated study involving Neuro-2a cells (the method used worldwide to test ciguatoxins), electrophysiological assays, and in silico simulations that evidence the similar cytotoxicity of Caribbean and Pacific ciguatoxins and their asymmetry binding within sodium channels. The binding mode of the toxins was first explored by molecular docking using the GOLD program and the resulting binary complexes were further studied by Molecular Dynamics simulation studies using the molecular mechanics force field AMBER. The simulation studies explain their distinct impact on the activation potential of the channel as experimentally observed and provide a detailed picture of the effects caused by these toxins on an atomic scale.Graphical

Parametrization of Force Field Bonded Terms under Structural Inconsistency.

Parametrization of the bonded part of molecular mechanics (MM) force fields (FFs) is typically done by fitting reference quantum mechanical (QM) energies or forces of representative structures. FFs for small molecules are constructed in incremental parametrization procedures, where parameters developed previously are retained for novel molecules, followed by optimization of missing, not previously optimized parameters. Equilibrium QM and MM geometries of molecules can deviate due to parameters transferred from existing molecules in the FF. In this work, we demonstrate that conventional parametrization methods based on fitting QM energies and/or forces to derive parameters for bond and angle terms produce largely suboptimal force constants when MM and QM equilibrium structures deviate. We further developed and tested a new method to derive CHARMM FF parameters based on the potential energy surface scans where a structural deviation between QM and MM optimized geometries is explicitly allowed during parametrization. The test of the new method was performed on a diverse set of 32 molecules. The results show that without any need for additional restraints, the new method produces robust and largely transferable parameters for bond and angle terms. The new method also improves the agreement for the normal modes for all molecules in the test set, reducing the average error in the reproduction of QM normal mode frequencies from 9.5% computed with CGenFF parameters to 6.8% computed with the new parameters. The new method will allow parametrization of molecules under structural deviations, common for force fields for small molecules, producing robust and transferable parameters.

Best Practices in Constant pH MD Simulations: Accuracy and Sampling.

Various approaches have been proposed to include the effect of pH in molecular dynamics (MD) simulations. Among these, the λ-dynamics approach proposed by Brooks and co-workers [Kong, X.; Brooks III, C. L. J. Chem. Phys.1996, 105, 2414−2423] can be performed with little computational overhead and hfor each typeence be used to routinely perform MD simulations at microsecond time scales, as shown in the accompanying paper [Aho, N. et al. J. Chem. Theory Comput.2022, DOI: 10.1021/acs.jctc.2c00516]. At such time scales, however, the accuracy of the molecular mechanics force field and the parametrization becomes critical. Here, we address these issues and provide the community with guidelines on how to set up and perform long time scale constant pH MD simulations. We found that barriers associated with the torsions of side chains in the CHARMM36m force field are too high for reaching convergence in constant pH MD simulations on microsecond time scales. To avoid the high computational cost of extending the sampling, we propose small modifications to the force field to selectively reduce the torsional barriers. We demonstrate that with such modifications we obtain converged distributions of both protonation and torsional degrees of freedom and hence consistent pKa estimates, while the sampling of the overall configurational space accessible to proteins is unaffected as compared to normal MD simulations. We also show that the results of constant pH MD depend on the accuracy of the correction potentials. While these potentials are typically obtained by fitting a low-order polynomial to calculated free energy profiles, we find that higher order fits are essential to provide accurate and consistent results. By resolving problems in accuracy and sampling, the work described in this and the accompanying paper paves the way to the widespread application of constant pH MD beyond pKa prediction.

A Quantum-Guided Molecular Mechanics Force Field for the Ferrocene Scaffold.

Ferrocene derivatives have a wide range of applications, including as ligands in asymmetric catalysis, due to their chemical stability, rigid backbone, steric bulk, and ability to encode stereochemical information via planar chirality. Unfortunately, few of the available molecular mechanics force fields incorporate parameters for the accurate study of this important building block. Here, we present a MM3* force field for ferrocenyl ligands, which was generated using the quantum-guided molecular mechanics (Q2MM) method. Detailed validation by comparison to DFT calculations and crystal structures demonstrates the accuracy of the parameters and uncovers the physical origin of deviations through excess energy analysis. Combining the ferrocene force field with a force field for Pd-allyl complexes and comparing the crystal structures shows the compatibility with previously developed MM3* force fields. Finally, the ferrocene force field was combined with a previously published transition-state force field to predict the stereochemical outcomes of the aminations of Pd-allyl complexes with different amines and different chiral ferrocenyl ligands, with an R2 of ∼0.91 over 10 examples.

Extension of the CHARMM Classical Drude Polarizable Force Field to N- and O-Linked Glycopeptides and Glycoproteins.

Molecular dynamic simulations are an effective tool to study complex molecular systems and are contingent upon the availability of an accurate and reliable molecular mechanics force field. The Drude polarizable force field, which allows for the explicit treatment of electronic polarization in a computationally efficient fashion, has been shown to reproduce experimental properties that were difficult or impossible to reproduce with the CHARMM additive force field, including peptide folding cooperativity, RNA hairpin structures, and DNA base flipping. Glycoproteins are essential components of glycoconjugate vaccines, antibodies, and many pharmaceutically important molecules, and an accurate polarizable force field that includes compatibility between the protein and carbohydrate aspect of the force field is essential to study these types of systems. In this work, we present an extension of the Drude polarizable force field to glycoproteins, including both N- and O-linked species. Parameter optimization focused on the dihedral terms using a reweighting protocol targeting NMR solution J-coupling data for model glycopeptides. Validation of the model include eight model glycopeptides and four glycoproteins with multiple N- and O-linked glycosylations. The new glycoprotein carbohydrate force field can be used in conjunction with the remainder of Drude polarizable force field through a variety of MD simulation programs including GROMACS, OPENMM, NAMD, and CHARMM and may be accessed through the Drude Prepper module in the CHARMM-GUI.

(Invited) Theoretical Study on the Excited and Charge-Separated States of Donor-Acceptor-Linked Molecules

Photoinduced charge separation process plays important roles in various fields such as natural and artificial photosynthesis, solar cell, and light-emitting diode. Toward an understanding of the molecular mechanism of photoinduced charge separation, we have been studying the excited and charge-separated states of donor-acceptor-linked molecules with theoretical methods such as QM/MM method and MD simulation. In this presentation, we will mainly introduce two topics of our research. The first one is porphyrin-fullerene linked molecules. It is experimentally known that the yields of charge-separated states of these molecules is not linear with the bridge length, but reaches a maximum at a specific length. It is also known that the energy of the singlet excited state of this compound is lower in benzonitrile solvent in the following order: local excited state of porphyrin, local excited state of fullerene, and charge separation state. However, by simply using the stable structure of the ground state in the theoretical calculations, the energies of each excited state are in the exact opposite order, and thus the experimental results could not be reproduced. In order to reproduce the experimental results, the structural relaxation of the molecule in the excited state and the response of the solvent are important. However, a simple polarization continuum model of the solvent is insufficient to stabilize the charge separation state. Therefore, we used the QM/MM RWFE-SCF method, which efficiently combines the QM/MM method with MD simulations to optimize the geometry of the solute molecule (QM region) with quantum chemical calculations including the effects of the thermal fluctuations of the surrounding environment such as the solvent molecules (MM region) treated by the molecular mechanics force field. The energy differences calculated with the QM/MM RWFE-SCF method are in good agreement with the experimental results. The charge-separated state is greatly stabilized by solvent molecules due to the large dipole moment. The next topic is acceptor-donor-acceptor type non-fullerene acceptor called TACIC. The excited state of TACIC is known to have a longer lifetime in the film state than in solution. This is in contrast to a similar non-fullerene acceptor molecule, ITIC, which has a shorter lifetime in the film state. However, the molecular mechanism is unknown. We analyzed the structures and excited states of the TACIC and ITIC by using quantum chemical calculations and MD simulation. It is found that the packing structure of TACIC in the film state is much different from that of ITIC, which could be related to the difference in lifetime of the excited state. Further analysis is currently in progress.

Evaluating Geometric Definitions of Stacking for RNA Dinucleoside Monophosphates Using Molecular Mechanics Calculations.

RNA modulation via small molecules is a novel approach in pharmacotherapies, where the determination of the structural properties of RNA motifs is considered a promising way to develop drugs capable of targeting RNA structures to control diseases. However, due to the complexity and dynamic nature of RNA molecules, the determination of RNA structures using experimental approaches is not always feasible, and computational models employing force fields can provide important insight. The quality of the force field will determine how well the predictions are compared to experimental observables. Stacking in nucleic acids is one such structural property, originating mainly from London dispersion forces, which are quantum mechanical and are included in molecular mechanics force fields through nonbonded interactions. Geometric descriptions are utilized to decide if two residues are stacked and hence to calculate the stacking free energies for RNA dinucleoside monophosphates (DNMPs) through statistical mechanics for comparison with experimental thermodynamics data. Here, we benchmark four different stacking definitions using molecular dynamics (MD) trajectories for 16 RNA DNMPs produced by two different force fields (RNA-IL and ff99OL3) and show that our stacking definition better correlates with the experimental thermodynamics data. While predictions within an accuracy of 0.2 kcal/mol at 300 K were observed in RNA CC, CU, UC, AG, GA, and GG, stacked states of purine-pyrimidine and pyrimidine-purine DNMPs, respectively, were typically underpredicted and overpredicted. Additionally, population distributions of RNA UU DNMPs were poorly predicted by both force fields, implying a requirement for further force field revisions. We further discuss the differences predicted by each RNA force field. Finally, we show that discrete path sampling (DPS) calculations can provide valuable information and complement the MD simulations. We propose the use of experimental thermodynamics data for RNA DNMPs as benchmarks for testing RNA force fields.