All backbone torsion angles, except the recently corrected α and γ4 (link), were parameterized using representative model compounds (Supplementary Fig. 23), for which torsional profiles were obtained at the MP2/aug-cc-pVDZ level using B3LYP/6-31++G(d,p)-optimized geometries21 -23 . Single-point calculations at crucial points of the conformational space were performed at the CCSD(T)/complete basis set (CBS) level24 -26 . Solvent effects were introduced using our MST28 method as implemented in Gaussian (http://www.gaussian.com). See Supplementary Notes for additional details on QM calculations.
Parameters were fitted using a flexible Monte-Carlo procedure4 (link), which minimizes the error between QM reference profiles (in solution) and classical potentials of mean force calculations in aqueous solution obtained from umbrella sampling calculations29 . By default we used gas phase-fitted values as first guess, and always limited the torsional representation to a three Fourier expansion terms, while reinforcing in the fitting the weight of the points described at the highest level of theory and those geometrical regions that are specially populated in experimental structures. Around 5–10 acceptable solutions of the Monte Carlo refinement were tested on short MD simulations (around 50–100 ns) for one small duplex d(CGATCG)2 rejecting those leading to distorted structures. The optimum parameter set, without additional refinement was extensively tested against experimental results. Additional details (and references) on the parameterization procedure are given in the Supplementary Notes. Note that the way in which the parameters were derived does not guarantee their validity for RNA simulations. The use of others, already validated, RNA force-fields are recommended.
As shown in the Online Methods Table 1, refined parmbsc1 parameters fit very well high-level QM data. The syn-anti equilibrium, which was non-optimal in parmbsc0, is now well reproduced (Supplementary Fig. 24). The fitting to sugar puckering profile was improved by increasing the East barrier, and by displacing the North and South minima to more realistic regions (Online Methods Table 1 and Supplementary Fig. 25). Additionally, parmbsc1 provides ε and ζ conformational map almost indistinguishable from the CCSD(T)/CBS results in solution (Supplementary Fig. 26), with errors in the estimates of relative BI/BII stability and transition barrier equal to 0.2 and 0.0 kcal mol−1 respectively.