Molecular Simulation Strategies for Understanding the Degradation Mechanisms of Acrylic Polymers.

Acrylic polymers, commonly used in paints, can degrade over time by several different chemical and physical mechanisms, depending on structure and exposure conditions. While exposure to UV light and temperature results in irreversible chemical damage, acrylic paint surfaces in museums can also accumulate pollutants, such as volatile organic compounds (VOCs) and moisture, that affect their material properties and stability. In this work, we studied the effects of different degradation mechanisms and agents on properties of acrylic polymers found in artists' acrylic paints for the first time using atomistic molecular dynamics simulations. Through the use of enhanced sampling methods, we investigated how pollutants are absorbed into thin acrylic polymer films from the environment around the glass transition temperature. Our simulations suggest that the absorption of VOCs is favorable (-4 to -7 kJ/mol depending on VOCs), and the pollutants can easily diffuse and be emitted back into the environment slightly above glass transition temperature when the polymer is soft. However, typical environmental fluctuations in temperature (<16 °C) can lead for these acrylic polymers to transition to glassy state, in which case the trapped pollutants act as plasticizers and cause a loss of mechanical stability in the material. This type of degradation results in disruption of polymer morphology, which we investigate through calculation of structural and mechanical properties. In addition, we also investigate the effects of chemical damage, such as backbone bond scission and side-chain cross-linking reactions on polymer properties.

Sampling rare event energy landscapes via birth-death augmented dynamics.

A common problem that affects simulations of complex systems within the computational physics and chemistry communities is the so-called sampling problem or rare event problem where proper sampling of energy landscapes is impeded by the presences of high kinetic barriers that hinder transitions between metastable states on typical simulation time scales. Many enhanced sampling methods have been developed to address this sampling problem and more efficiently sample rare event systems. An interesting idea, coming from the field of statistics, was introduced in a recent work [Lu, Lu, and Nolen, Accelerating Langevin sampling with birth-death, arXiv:1905.09863] in the form of a novel sampling algorithm that augments overdamped Langevin dynamics with a birth-death process. In this work, we expand on this idea and show that this birth-death sampling scheme can efficiently sample prototypical rare event energy landscapes, and that the speed of equilibration is independent of the barrier height. We amend a crucial shortcoming of the original algorithm that leads to incorrect sampling of barrier regions by introducing an alternative approximation of the birth-death term. We establish important theoretical properties of the modified algorithm and prove mathematically that the relevant convergence results still hold. We investigate via numerical simulations the effect of various parameters, and we investigate ways to reduce the computational effort of the sampling scheme. We show that the birth-death mechanism can be used to accelerate sampling in the more general case of underdamped Langevin dynamics that is more commonly used in simulating physical systems. Our results show that this birth-death scheme is a promising method for sampling rare event energy landscapes.

The Crucial Role of Solvation Forces in the Steric Stabilization of Nanoplatelets.

The steric stability of inorganic colloidal particles in an apolar solvent is usually described in terms of the balance between three contributions: the van der Waals attraction, the free energy of mixing, and the ligand compression. However, in the case of nanoparticles, the discrete nature of the ligand shell and the solvent has to be taken into account. Cadmium selenide nanoplatelets are a special case. They combine a weak van der Waals attraction and a large facet to particle size ratio. We use coarse grained molecular dynamics simulations of nanoplatelets in octane to demonstrate that solvation forces are strong enough to induce the formation of nanoplatelet stacks and by that have a crucial impact on the steric stability. In particular, we demonstrate that for sufficiently large nanoplatelets, solvation forces are proportional to the interacting facet area, and their strength is intrinsically tied to the softness of the ligand shell.

Reweighted Manifold Learning of Collective Variables from Enhanced Sampling Simulations.

Enhanced sampling methods are indispensable in computational chemistry and physics, where atomistic simulations cannot exhaustively sample the high-dimensional configuration space of dynamical systems due to the sampling problem. A class of such enhanced sampling methods works by identifying a few slow degrees of freedom, termed collective variables (CVs), and enhancing the sampling along these CVs. Selecting CVs to analyze and drive the sampling is not trivial and often relies on chemical intuition. Despite routinely circumventing this issue using manifold learning to estimate CVs directly from standard simulations, such methods cannot provide mappings to a low-dimensional manifold from enhanced sampling simulations, as the geometry and density of the learned manifold are biased. Here, we address this crucial issue and provide a general reweighting framework based on anisotropic diffusion maps for manifold learning that takes into account that the learning data set is sampled from a biased probability distribution. We consider manifold learning methods based on constructing a Markov chain describing transition probabilities between high-dimensional samples. We show that our framework reverts the biasing effect, yielding CVs that correctly describe the equilibrium density. This advancement enables the construction of low-dimensional CVs using manifold learning directly from the data generated by enhanced sampling simulations. We call our framework reweighted manifold learning. We show that it can be used in many manifold learning techniques on data from both standard and enhanced sampling simulations.

Improving the Efficiency of Variationally Enhanced Sampling with Wavelet-Based Bias Potentials.

Collective variable-based enhanced sampling methods are routinely used on systems with metastable states, where high free energy barriers impede the proper sampling of the free energy landscapes when using conventional molecular dynamics simulations. One such method is variationally enhanced sampling (VES), which is based on a variational principle where a bias potential in the space of some chosen slow degrees of freedom, or collective variables, is constructed by minimizing a convex functional. In practice, the bias potential is taken as a linear expansion in some basis function set. So far, primarily basis functions delocalized in the collective variable space, like plane waves, Chebyshev, or Legendre polynomials, have been used. However, there has not been an extensive study of how the convergence behavior is affected by the choice of the basis functions. In particular, it remains an open question if localized basis functions might perform better. In this work, we implement, tune, and validate Daubechies wavelets as basis functions for VES. The wavelets construct orthogonal and localized bases that exhibit an attractive multiresolution property. We evaluate the performance of wavelet and other basis functions on various systems, going from model potentials to the calcium carbonate association process in water. We observe that wavelets exhibit excellent performance and much more robust convergence behavior than all other basis functions, as well as better performance than metadynamics. In particular, using wavelet bases yields far smaller fluctuations of the bias potential within individual runs and smaller differences between independent runs. Based on our overall results, we can recommend wavelets as basis functions for VES.

Acrylic Paints: An Atomistic View of Polymer Structure and Effects of Environmental Pollutants.

Most of the artwork and cultural heritage objects are stored in museums under conditions that are difficult to monitor. While advanced technologies aim to control and prevent the degradation of cultural heritage objects in line with preventive conservation measures, there is much to be learned in terms of the physical processes that lead to the degradation of the synthetic polymers that form the basis of acrylic paints largely used in contemporary art. In museums, stored objects are often exposed to temperature and relative humidity fluctuations as well as airborne pollutants such as volatile organic compounds (VOCs). The glass transition of acrylic paints is below room temperature; while low temperatures may cause cracking, at high temperatures the sticky surface of the paint becomes vulnerable to pollutants. Here we develop fully atomistic models to understand the structure of two types of acrylic copolymers and their interactions with VOCs and water. The structure and properties of acrylic copolymers are slighlty modified by incorporation of a monomer with a longer side chain. With favorable solvation free energies, once absorbed, VOCs and water interact with the polymer side chains to form hydrogen bonds. The cagelike structure of the polymers prevents the VOCs and water to diffuse freely below the glass transition temperature. In addition, our model forms the foundation for developing mesoscopic and continuum models that will allow us to access longer time and length scales to further our understanding of the degradation of artwork.

Multiscale Reweighted Stochastic Embedding: Deep Learning of Collective Variables for Enhanced Sampling.

Machine learning methods provide a general framework for automatically finding and representing the essential characteristics of simulation data. This task is particularly crucial in enhanced sampling simulations. There we seek a few generalized degrees of freedom, referred to as collective variables (CVs), to represent and drive the sampling of the free energy landscape. In theory, these CVs should separate different metastable states and correspond to the slow degrees of freedom of the studied physical process. To this aim, we propose a new method that we call multiscale reweighted stochastic embedding (MRSE). Our work builds upon a parametric version of stochastic neighbor embedding. The technique automatically learns CVs that map a high-dimensional feature space to a low-dimensional latent space via a deep neural network. We introduce several new advancements to stochastic neighbor embedding methods that make MRSE especially suitable for enhanced sampling simulations: (1) weight-tempered random sampling as a landmark selection scheme to obtain training data sets that strike a balance between equilibrium representation and capturing important metastable states lying higher in free energy; (2) a multiscale representation of the high-dimensional feature space via a Gaussian mixture probability model; and (3) a reweighting procedure to account for training data from a biased probability distribution. We show that MRSE constructs low-dimensional CVs that can correctly characterize the different metastable states in three model systems: the Müller-Brown potential, alanine dipeptide, and alanine tetrapeptide.

Free-energy landscape of polymer-crystal polymorphism.

Polymorphism rationalizes how processing can control the final structure of a material. The rugged free-energy landscape and exceedingly slow kinetics in the solid state have so far hampered computational investigations. We report for the first time the free-energy landscape of a polymorphic crystalline polymer, syndiotactic polystyrene. Coarse-grained metadynamics simulations allow us to efficiently sample the landscape at large. The free-energy difference between the two main polymorphs, α and ß, is further investigated by quantum-chemical calculations. The results of the two methods are in line with experimental observations: they predict ß as the more stable polymorph under standard conditions. Critically, the free-energy landscape suggests how the α polymorph may lead to experimentally observed kinetic traps. The combination of multiscale modeling, enhanced sampling, and quantum-chemical calculations offers an appealing strategy to uncover complex free-energy landscapes with polymorphic behavior.

Erratum: Enhancing Entropy and Enthalpy Fluctuations to Drive Crystallization in Atomistic Simulations [Phys. Rev. Lett. 119, 015701 (2017)].

This corrects the article DOI: 10.1103/PhysRevLett.119.015701.

Mg2+ Sensing by an RNA Fragment: Role of Mg2+-Coordinated Water Molecules.

RNA molecules selectively bind to specific metal ions to populate their functional active states, making it important to understand their source of ion selectivity. In large RNA systems, metal ions interact with the RNA at multiple locations, making it difficult to decipher the precise role of ions in folding. To overcome this complexity, we studied the role of different metal ions (Mg2+, Ca2+, and K+) in the folding of a small RNA hairpin motif (5'-ucCAAAga-3') using unbiased all-atom molecular dynamics simulations. The advantage of studying this system is that it requires specific binding of a single metal ion to fold to its native state. We find that even for this small RNA, the folding free energy surface (FES) is multidimensional as different metal ions present in the solution can simultaneously facilitate folding. The FES shows that specific binding of a metal ion is indispensable for its folding. We further show that in addition to the negatively charged phosphate groups, the spatial organization of electronegative nucleobase atoms drives the site-specific binding of the metal ions. Even though the binding site cannot discriminate between different metal ions, RNA folds efficiently only in a Mg2+ solution. We show that the rigid network of Mg2+-coordinated water molecules facilitates the formation of important interactions in the transition state. The other metal ions such as K+ and Ca2+ cannot facilitate the formation of such interactions. These results allow us to hypothesize possible metal-sensing mechanisms in large metalloriboswitches and also provide useful insights into the design of appropriate collective variables for studying large RNA molecules using enhanced sampling methods.

Finding multiple reaction pathways of ligand unbinding.

Searching for reaction pathways describing rare events in large systems presents a long-standing challenge in chemistry and physics. Incorrectly computed reaction pathways result in the degeneracy of microscopic configurations and inability to sample hidden energy barriers. To this aim, we present a general enhanced sampling method to find multiple diverse reaction pathways of ligand unbinding through nonconvex optimization of a loss function describing ligand-protein interactions. The method successfully overcomes large energy barriers using an adaptive bias potential and constructs possible reaction pathways along transient tunnels without the initial guesses of intermediate or final states, requiring crystallographic information only. We examine the method on the T4 lysozyme L99A mutant which is often used as a model system to study ligand binding to proteins, provide a previously unknown reaction pathway, and show that by using the bias potential and the tunnel widths, it is possible to capture heterogeneity of the unbinding mechanisms between the found transient protein tunnels.

Chemical potential calculations in non-homogeneous liquids.

The numerical computation of chemical potential in dense non-homogeneous fluids is a key problem in the study of confined fluid thermodynamics. To this day, several methods have been proposed; however, there is still need for a robust technique, capable of obtaining accurate estimates at large average densities. A widely established technique is the Widom insertion method, which computes the chemical potential by sampling the energy of insertion of a test particle. Non-homogeneity is accounted for by assigning a density dependent weight to the insertion points. However, in dense systems, the poor sampling of the insertion energy is a source of inefficiency, hampering a reliable convergence. We have recently presented a new technique for the chemical potential calculation in homogeneous fluids. This novel method enhances the sampling of the insertion energy via well-tempered metadynamics, reaching accurate estimates at very large densities. In this paper, we extend the technique to the case of non-homogeneous fluids. The method is successfully tested on a confined Lennard-Jones fluid. In particular, we show that, thanks to the improved sampling, our technique does not suffer from a systematic error that affects the classic Widom method for non-homogeneous fluids, providing a precise and accurate result.

Frequency adaptive metadynamics for the calculation of rare-event kinetics.

The ability to predict accurate thermodynamic and kinetic properties in biomolecular systems is of both scientific and practical utility. While both remain very difficult, predictions of kinetics are particularly difficult because rates, in contrast to free energies, depend on the route taken. For this reason, specific enhanced sampling methods are needed to calculate long-time scale kinetics. It has recently been demonstrated that it is possible to recover kinetics through the so-called "infrequent metadynamics" simulations, where the simulations are biased in a way that minimally corrupts the dynamics of moving between metastable states. This method, however, requires the bias to be added slowly, thus hampering applications to processes with only modest separations of time scales. Here we present a frequency-adaptive strategy which bridges normal and infrequent metadynamics. We show that this strategy can improve the precision and accuracy of rate calculations at fixed computational cost and should be able to extend rate calculations for much slower kinetic processes.

Conformational Entropy as Collective Variable for Proteins.

Many enhanced sampling methods rely on the identification of appropriate collective variables. For proteins, even small ones, finding appropriate descriptors has proven challenging. Here we suggest that the NMR S2 order parameter can be used to this effect. We trace the validity of this statement to the suggested relation between S2 and conformational entropy. Using the S2 order parameter and a surrogate for the protein enthalpy in conjunction with metadynamics or variationally enhanced sampling, we are able to reversibly fold and unfold a small protein and draw its free energy at a fraction of the time that is needed in unbiased simulations. We also use S2 in combination with the free energy flooding method to compute the unfolding rate of this peptide. We repeat this calculation at different temperatures to obtain the unfolding activation energy.

Enhancing Entropy and Enthalpy Fluctuations to Drive Crystallization in Atomistic Simulations.

Crystallization is a process of great practical relevance in which rare but crucial fluctuations lead to the formation of a solid phase starting from the liquid. As in all first order first transitions, there is an interplay between enthalpy and entropy. Based on this idea, in order to drive crystallization in molecular simulations, we introduce two collective variables, one enthalpic and the other entropic. Defined in this way, these collective variables do not prejudge the structure into which the system is going to crystallize. We show the usefulness of this approach by studying the cases of sodium and aluminum that crystallize in the bcc and fcc crystalline structures, respectively. Using these two generic collective variables, we perform variationally enhanced sampling and well tempered metadynamics simulations and find that the systems transform spontaneously and reversibly between the liquid and the solid phases.

Coarse graining from variationally enhanced sampling applied to the Ginzburg-Landau model.

A powerful way to deal with a complex system is to build a coarse-grained model capable of catching its main physical features, while being computationally affordable. Inevitably, such coarse-grained models introduce a set of phenomenological parameters, which are often not easily deducible from the underlying atomistic system. We present a unique approach to the calculation of these parameters, based on the recently introduced variationally enhanced sampling method. It allows us to obtain the parameters from atomistic simulations, providing thus a direct connection between the microscopic and the mesoscopic scale. The coarse-grained model we consider is that of Ginzburg-Landau, valid around a second-order critical point. In particular, we use it to describe a Lennard-Jones fluid in the region close to the liquid-vapor critical point. The procedure is general and can be adapted to other coarse-grained models.

Variational Flooding Study of a SN2 Reaction.

We have studied the reaction dynamics of a prototypical organic reaction using a variationally optimized truncated bias to accelerate transitions between educt and product reactant states. The asymmetric SN2 nucleophilic substitution reaction of fluoromethane and chloromethane CH3F + Cl- ⇌ CH3Cl + F- is considered, and many independent biased molecular dynamics simulations have been performed at 600, 900, and 1200 K, collecting several hundred transitions at each temperature. The transition times and relative rate constants have been obtained for both reaction directions. The activation energies extracted from an Arrhenius plot compare well with standard static calculations.

Electrostatic versus Resonance Interactions in Photoreceptor Proteins: The Case of Rhodopsin.

Light sensing in photoreceptor proteins is subtly modulated by the multiple interactions between the chromophoric unit and its binding pocket. Many theoretical and experimental studies have tried to uncover the fundamental origin of these interactions but reached contradictory conclusions as to whether electrostatics, polarization, or intrinsically quantum effects prevail. Here, we select rhodopsin as a prototypical photoreceptor system to reveal the molecular mechanism underlying these interactions and regulating the spectral tuning. Combining a multireference perturbation method and density functional theory with a classical but atomistic and polarizable embedding scheme, we show that accounting for electrostatics only leads to a qualitatively wrong picture, while a responsive environment can successfully capture both the classical and quantum dominant effects. Several residues are found to tune the excitation by both differentially stabilizing ground and excited states and through nonclassical "inductive resonance" interactions. The results obtained with such a quantum-in-classical model are validated against both experimental data and fully quantum calculations.

Hierarchical Protein Free Energy Landscapes from Variationally Enhanced Sampling.

In recent work, we demonstrated that it is possible to obtain approximate representations of high-dimensional free energy surfaces with variationally enhanced sampling ( Shaffer, P.; Valsson, O.; Parrinello, M. Proc. Natl. Acad. Sci. , 2016 , 113 , 17 ). The high-dimensional spaces considered in that work were the set of backbone dihedral angles of a small peptide, Chignolin, and the high-dimensional free energy surface was approximated as the sum of many two-dimensional terms plus an additional term which represents an initial estimate. In this paper, we build on that work and demonstrate that we can calculate high-dimensional free energy surfaces of very high accuracy by incorporating additional terms. The additional terms apply to a set of collective variables which are more coarse than the base set of collective variables. In this way, it is possible to build hierarchical free energy surfaces, which are composed of terms that act on different length scales. We test the accuracy of these free energy landscapes for the proteins Chignolin and Trp-cage by constructing simple coarse-grained models and comparing results from the coarse-grained model to results from atomistic simulations. The approach described in this paper is ideally suited for problems in which the free energy surface has important features on different length scales or in which there is some natural hierarchy.