Logo: to the web site of Uppsala University

Most general force fields only implement a harmonic potential to model covalent bonds. In addition, in some force fields, all or a selection of the covalent bonds are constrained in molecular dynamics simulations. Nevertheless, it is possible to implement accurate bond potentials for a relatively small computational cost. Such potentials may be important for spectroscopic applications, free energy perturbation calculations or for studying reactions using empirical valence bond theory. Here, we evaluate different bond potentials for diatomic molecules. Based on quantum-chemical scans around the equilibrium distance of 71 molecules using the MP2/aug-cc-pVTZ level of theory as well as CCSD(T) with the same basis-set, we determine the quality of fit to the data of 28 model potentials. As expected, a large spread in accuracies of the potentials is found and more complex potentials generally provide a better fit. As a second and more challenging test, five spectroscopic parameters (ω_e, ω_ex_e, α_e, B_e and D_e) predicted based on quantum chemistry as well as the fitted potentials are compared to experimental data. A handful of the 28 potentials tested are found to be accurate. Of these, we suggest that the potential due to Hua (Phys. Rev. A, 42 (1990), 2524) could be a suitable choice for implementation in molecular simulations codes, since it is considerably more accurate than the well-known Morse potential (Phys. Rev., 34 (1929), 57) at a very similar computational cost.

This article presents an evaluation of the environmental impact of forestry based on landscape theory. It has been argued that this type of forestry offers a positive impact on the climate because there is a balance between the amount of greenhouse gas emissions and the absorption of these gases within an entire forested area. However, this analysis will demonstrate that the arrangement and composition of managed forests are linked to a significant carbon debt. This debt represents the disparity between the carbon that would typically be stored in a natural forest and the actual amount of carbon stored in the managed forest. While this excess carbon remains in the atmosphere rather than being sequestered, it contributes to the greenhouse effect. Using Swedish forestry as an example, the carbon debt is estimated to be comparable in scale to the total accumulated fossil fuel emissions of the country.

The systems of noble gases are particularly instructive for molecular modeling due to the elemental nature of their interactions. They do not normally form bonds nor possess a (permanent) dipole moment, and the only forces determining their bonding/clustering stems from van der Waals forces─dispersion and Pauli repulsion, which can be modeled by empirical potential functions. Combination rules, that is, formulas to derive parameters for pair potentials of heterodimers from parameters of corresponding homodimers, have been studied at length for the Lennard-Jones 12-6 potentials but not in great detail for other, more accurate, potentials. In this work, we examine the usefulness of nine empirical potentials in their ability to reproduce quantum mechanical (QM) benchmark dissociation curves of noble gas dimers (He, Ne, Ar, Kr, and Xe homo- and heterodimers), and we systematically study the efficacy of different permutations of combination relations for each parameter of the potentials. Our QM benchmark comprises dissociation curves computed by several different coupled cluster implementations as well as symmetry-adapted perturbation theory. The two-parameter Lennard-Jones potentials were decisively outperformed by more elaborate potentials that sport a 25–30 times lower root-mean-square error (RMSE) when fitted to QM dissociation curves. Very good fits to the QM dissociation curves can be achieved with relatively inexpensive four- or even three-parameter potentials, for instance, the damped 14-7 potential (Halgren, J. Am. Chem. Soc. 1992, 114, 7827–7843), a four-parameter Buckingham potential (Werhahn et al., Chem. Phys. Lett. 2015, 619, 133–138), or the three-parameter Morse potential (Morse, Phys. Rev. 1929, 34, 57–64). Potentials for heterodimers that are generated from combination rules have an RMSE that is up to 20 times higher than potentials that are directly fitted to the QM dissociation curves. This means that the RMSE, in particular, for light atoms, is comparable in magnitude to the well-depth of the potential. Based on a systematic permutation of combination rules, we present one or more combination rules for each potential tested that yield a relatively low RMSE. Two new combination rules are introduced that perform well, one for the van der Waals radius σ_ij as (1/2(𝜎³_𝑖+𝜎³_𝑗))^1/3 and one for the well-depth ϵ_ij as (1/2(𝜖⁻²_𝑖+𝜖⁻²_𝑗))^−1/2. The QM data and the fitted potentials were evaluated in the gas phase against experimental second virial coefficients for homo- and heterodimers, the latter of which allowed evaluation of the combination rules. The fitted models were used to perform condensed phase molecular dynamics simulations to verify the melting points, liquid densities at the melting point, and the enthalpies of vaporization produced by the models for pure substances. Subtle differences in the benchmark potentials, in particular, the well-depth, due to the level of theory used were found here to have a profound effect on the macroscopic properties of noble gases: second virial coefficients or the bulk properties in simulations. By explicitly including three-body dispersion in molecular simulations employing the best pair potential, we were able to obtain accurate melting points as well as satisfactory densities and enthalpies of vaporization.

Computational chemistry is an important tool in numerous scientific disciplines, including drug discovery and structural biology. Coarse-grained models offer simple representations of molecular systems that enable simulations of large-scale systems. Because there has been an increase in the adoption of such models for simulations of biomolecular systems, critical evaluation is warranted. Here, the stability of the amyloid peptide and organic crystals is evaluated using the Martini 3 coarse-grained force field. The crystals change shape drastically during the simulations. Radial distribution functions show that the distance between backbone beads in β-sheets increases by ∼1 Å, breaking the crystals. The melting points of organic compounds are much too low in the Martini force field. This suggests that Martini 3 lacks the specific interactions needed to accurately simulate peptides or organic crystals without imposing artificial restraints. The problems may be exacerbated by the use of the 12-6 potential, suggesting that a softer potential could improve this model for crystal simulations.

In some compounds, exchange repulsion is orientation dependent. However, in contrast to quantum chemical methods that treat exchange explicitly, empirical models assume exchange to be spherically symmetric, yielding an average description only. Here we quantify the anisotropy of exchange and dispersion energy for hydrogen halides and water by probing these compounds with a helium atom using the symmetry-adapted perturbation theory (SAPT). The exchange interaction is reduced by up to 33% due to the sigma-hole in hydrogen iodide, depending on the location of the probe. We demonstrate how this anisotropy can be modeled in empirical force fields either using an angle-dependent potential or by introducing virtual sites, reducing the error in the empirical model by a factor of 5 compared to isotropic atoms. Lone-pairs on water, positioned close to perpendicular to the plane of the molecule, on a line with the oxygen atom, and, surprisingly, sigma-holes on water both modulate the exchange interaction strongly. Both lone-pairs and sigma-holes can be modeled by virtual sites, leading to an 80% reduced error.

There are many problems in biochemistry that are difficult to study experimentally. Simulation methods are appealing due to direct availability of atomic coordinates as a function of time. However, direct molecular simulations are challenged by the size of systems and the time scales needed to describe relevant motions. In theory, enhanced sampling algorithms can help to overcome some of the limitations of molecular simulations. Here, we discuss a problem in biochemistry that offers a significant challenge for enhanced sampling methods and that could, therefore, serve as a benchmark for comparing approaches that use machine learning to find suitable collective variables. In particular, we study the transitions LacI undergoes upon moving between being non-specifically and specifically bound to DNA. Many degrees of freedom change during this transition and that the transition does not occur reversibly in simulations if only a subset of these degrees of freedom are biased. We also explain why this problem is so important to biologists and the transformative impact that a simulation of it would have on the understanding of DNA regulation.

A profound understanding of the physicochemical properties of organic crystals is crucial for topics from material science to drug discovery. Using molecular dynamics (MD) simulations with a sufficiently accurate force field, microscopic insight into structure and dynamics can be obtained of materials, including liquids and biomolecules. They are a valuable complement to experimental investigations that are used routinely in drug design, but not very often for studies of organic crystals. Indeed, the often delicate interactions in organic crystals act as a sensitive probe to investigate the accuracy of force fields. Here, we study the structural, dynamic, and thermodynamic properties of 30 organic crystals using the popular general AMBER force field (GAFF). In particular, we investigate both solid-solid and solid-liquid phase transitions. Melting points were determined using extensive solid-liquid coexistence simulations. For many compounds, we detect a phase transition from an ordered to a plastic crystal in the simulations. Based on the translational and rotational dynamics of the compounds, we can rationalize the properties of the plastic crystal phase. MD simulations can therefore help to answer the important question of whether or not organic crystals have a plastic crystal phase, and if so, what are the underlying factors in the molecular structure determining that.

There still is little treatment available for amyloid diseases, despite their significant impact on individuals and the social and economic implications for society. One reason for this is that the physical nature of amyloid formation is not understood sufficiently well. Therefore, fundamental research at the molecular level remains necessary to support the development of therapeutics. A few structures of short peptides from amyloid-forming proteins have been determined. These can in principle be used as scaffolds for designing aggregation inhibitors. Attempts to this end have often used the tools of computational chemistry, in particular molecular simulation. However, few simulation studies of these peptides in the crystal state have been presented so far. Hence, to validate the capability of common force fields (AMBER19SB, CHARMM36m, and OPLS-AA/M) to yield insight into the dynamics and structural stability of amyloid peptide aggregates, we have performed molecular dynamics simulations of twelve different peptide crystals at two different temperatures. From the simulations, we evaluate the hydrogen bonding patterns, the isotropic B-factors, the change in energy, the Ramachandran plots, and the unit cell parameters and compare the results with the crystal structures. Most crystals are stable in the simulations but for all force fields there is at least one that deviates from the experimental crystal, suggesting more work is needed on these models.

The human genome codes only a few thousand druggable proteins, mainly receptors and enzymes. While this pool of available drug targets is limited, there is an untapped potential for discovering new drug-binding mechanisms and modes. For example, enzymes with long binding cavities offer numerous prerequisite binding sites that may be visited by an inhibitor during migration from a bulk solution to the destination site. Drug design can use these prerequisite sites as new structural targets. However, identifying these ephemeral sites is challenging. Here, we introduce a new method called NetBinder for the systematic identification and classification of prerequisite binding sites at atomic resolution. NetBinder is based on atomistic simulations of the full inhibitor binding process and provides a networking framework on which to select the most important binding modes and uncover the entire binding mechanism, including previously undiscovered events. NetBinder was validated by a study of the binding mechanism of blebbistatin (a potent inhibitor) to myosin 2 (a promising target for cancer chemotherapy). Myosin 2 is a good test enzyme because, like other potential targets, it has a long internal binding cavity that provides blebbistatin with numerous potential prerequisite binding sites. The mechanism proposed by NetBinder of myosin 2 structural changes during blebbistatin binding shows excellent agreement with experimentally determined binding sites and structural changes. While NetBinder was tested on myosin 2, it may easily be adopted to other proteins with long internal cavities, such as G-protein-coupled receptors or ion channels, the most popular current drug targets. NetBinder provides a new paradigm for drug design by a network-based elucidation of binding mechanisms at an atomic resolution.

Equilibrium simulations of molecular systems allow to extract many physicochemical properties. Given an "accurate enough" model, a "large enough" simulation system and "long enough" simulations, such calculations should yield accurate predictions of properties that can be tested by experimental measurements. Non-equilibrium simulations can be used as a tool to obtain specific properties like viscosity or conductivity, but they have the drawback that in general only one property per simulation is produced. In addition, a range of methods is available for computing free energy differences. We here review the state of the art of using classical simulation models for generating quantitative predictions. Popular force fields have significant predictive power already but there is room for improvement. Bonded force potentials may need to be replaced by more accurate ones to better reproduce vibrational frequencies. Simplification of non-bonded force terms, such as cut-offs for electrostatic or dispersion interactions, should be avoided. Routine usage of force field methods will therefore require some tuning of parameters. Despite the extensive toolbox that is available for producing quantitative results, the computational cost of explicit solvent simulation is significant and therefore, approximate methods like implicit solvent models remain popular and are still being developed. Based on fundamental arguments as well as on examples of solvation free energies, host-guest complexation and non-covalent association of molecules in solution, we conclude that implicit solvents as well as algorithmic simplifications are most useful when validation using experimental data or rigorous theoretical treatments is possible.