Polarizable Model Research Articles

Polarizable force fields for accurate molecular simulations of aqueous solutions of electrolytes, crystalline salts, and solubility: Li+, Na+, K+, Rb+, F−, Cl−, Br−, I−

We develop and study polarizable microscopic models, force fields, for molecular simulations of alkali-halide electrolyte aqueous solutions, their crystals, and phase equilibria. We start from the AH/BK3 force fields of Kiss and Baranyai (P. T. Kiss and A. A. Baranyai, J. Chem. Phys. 2014, 141, 114501), which we refine using an approach for determining ion-ion interaction parameters, directly targeting experimental values of the lattice energy, pressure at a given density, and bulk and shear moduli of anhydrous electrolyte crystals. We apply the approach to 16 alkali-halide salts crystallizing in the face-centered cubic rock salt structure. As a result, we obtain force fields which predict rather accurately properties of crystals including their chemical potentials, properties of aqueous solutions including their static permittivity, and aqueous solubility. The force fields can thus find applications in molecular simulations of alkali-halide aqueous electrolytes, their interfaces, phase equilibria, metastable states, and other cases where less accurate models may exhibit various undesirable features such as excessive ion pairing, spurious precipitation, or incorrect low mobility.

A ΔSCF model for excited states within a polarisable embedding

Hybrid TDDFT/MM approaches are very popular methods for describing electronic transitions of molecules in solution or embedded in more complex (bio)matrices. However, when combined with a polarisable force field some problems can appear depending on the type of environment response scheme that is used. In particular, specific a posteriori corrections are generally needed to accurately describe charge-transfer states implying a large reorganisation of the electron density in the excited state. Here, we present a possible strategy to solve this issue by introducing a ΔSCF formulation. As the ΔSCF strategy has the advantage of being intrinsically state specific, its coupling to a polarisable model is expected to be particularly suited to describe all cases where the standard, linear response, formulation of a polarisable TDDFT/MM approach is not sufficient.

A Unified Approach to Derive Atomic Partial Charges and Polarizabilities of Ionic Liquids.

We propose a unified approach to fit simultaneously a set of atomic partial charges and polarizabilities of the polarizable model against the ab initio electrostatic potential (ESP) and polarizability. The polarizable model is represented with interactive atomic dipoles with distance-dependent attenuation. For the polarizable model employed in this study, the internal electric field on the polarization sites is fully turned on, and thus allows self-induced dipoles, which persist even for an isolated molecule/ion. By such treatment, the contribution of ESP stems not only from the partial charges but also from the self-induced dipoles, and the atomic partial charges and polarizabilities can be fitted simultaneously against ESP in a unified manner. The fitting with 1-ethyl-3-methylimidazolium (EMIM+) and nitrate (NO3-), a prototypical organic cation and inorganic anion, respectively, that can form ionic liquid, demonstrates that allowance of the self-induced dipoles gives much better fitness. Moreover, test on the total dipole of an EMIM+/NO3- ion pair shows that the agreement with the ab initio dipole is also much improved for the polarizable model, which highlights the importance of the polarization effects of ionic liquids.

A Polarizable Cationic Dummy Metal Ion Model.

A novel locally polarizable multisite model based on the original cation dummy atom (CDA) model is described for molecular dynamics simulations of ions in condensed phases. Polarization effects are introduced by the electronegativity equalization model (EEM) method where charges on the metal ion and its dummy atoms can fluctuate to respond to the environment. This model includes explicit polarization and ion-induced interactions and can be coupled with nonpolarizable or polarizable water models, making it more transferable to simpler force fields. This approach allows us to enhance the original fixed charge CDA model where the charge distribution cannot adapt to the local solvent structure. To illustrate the new CDApol model, we examined properties of the Zn2+, Al3+, and Zr4+ ions in aqueous solution. The polarizable model and Lennard-Jones parameters were refined for octahedrally coordinated Zn2+, Al3+, and Zr4+ CDAs to reproduce thermodynamic and geometrical properties. Using this locally polarizable model, we were able to obtain the experimental hydration free energy, ion-oxygen distance, and coordination number coupled with the standard 12-6 Lennard-Jones model. This model can be used in myriad additional applications where local polarization and charge transfer is important.

A fast method for electronic couplings in embedded multichromophoric systems

Electronic couplings are key to understanding exciton delocalization and transport in natural and artificial light harvesting processes. We develop a method to compute couplings in multichromophoric aggregates embedded in complex environments without running expensive quantum chemical calculations. We use a transition charge approximation to represent the quantum mechanical transition densities of the chromophores and an atomistic and polarizable classical model to describe the environment atoms. We extend our framework to estimate transition charges directly from the chromophore geometry, i.e., bypassing completely the quantum mechanical calculations using a regression approach. The method allows to rapidly compute accurate couplings for a large number of geometries along molecular dynamics trajectories.

Thermal conductance of the water-gold interface: The impact of the treatment of surface polarization in non-equilibrium molecular simulations.

Interfacial thermal conductance (ITC) quantifies heat transport across material-fluid interfaces. It is a property of crucial importance to study heat transfer processes at both macro- and nanoscales. Therefore, it is essential to accurately model the specific interactions between solids and liquids. Here, we investigate the thermal conductance of gold-water interfaces using polarizable and non-polarizable models. Both models have been fitted to reproduce the interfacial tension of the gold-water interface, but they predict significantly different ITCs. We demonstrate that the treatment of polarization using Drude-like models, widely employed in molecular simulations, leads to a coupling of the solid and liquid vibrational modes that give rise to a significant overestimation of the ITCs. We analyze the dependence of the vibrational coupling with the mass of the Drude particle and propose a solution to the artificial enhancement of the ITC, preserving at the same time the polarization response of the solid. Based on our calculations, we estimate ITCs of 200 MW/(m2 K) for the water-gold interface. This magnitude is comparable to that reported recently for gold-water interfaces [279 ± 16 MW/(m2 K)] using atomic fluctuating charges to account for the polarization contribution.

Conical Shape Fluctuations Determine the Rate of Ion Evaporation and the Emitted Cluster Size Distribution from Multicharged Droplets.

The ion evaporation mechanism (IEM) is perceived to be a major pathway for disintegration of multi-ion charged droplets found in atmospheric and sprayed aerosols. However, the precise mechanism of IEM and the effect of the nature of the ions in the emitted cluster size distribution have not yet been established despite its broad use in mass spectrometry and atmospheric chemistry over the past half century. Here, we present a systematic study of the emitted ion cluster distribution in relation to their spatial distribution in the parent droplet using atomistic modeling. It is found that in the parent droplet, multiple kosmotropic and weakly polarizable chaotropic ions (Cs+) are buried deeper within the droplet than polarizable chaotropic ions (Cl-, I-). This differentiation in the ion location is only captured by a polarizable model. It is demonstrated that the emitted cluster size distribution is determined by dynamic conical deformations and not by the equilibrium ion depth within the parent droplet as the IEM models assume. Critical factors that determine the cluster size distribution such as the charge sign asymmetry that have not been considered in models and in experiments are presented. We argue that the existing IEM analytical models do not establish a clear difference between IEM and Rayleigh fission. We propose a shift in the existing view for IEM from the equilibrium properties of the parent droplet to the chemistry in the conical shape fluctuations that serve as the centers for ion emission. Consequently, chemistry in the conical fluctuations may also be a key element to explain charge states of macromolecules in mass spectrometry and may have potential applications in catalysis due to the electric field in the conical region.

Isotropic and anisotropic collision-induced light scattering spectra of krypton gas simultaneously fitted by ground state interatomic potential

Isotropic and anisotropic binary collision-induced light scattering spectra of gaseous krypton obtained at room temperature are analyzed in terms of different interatomic potentials, both reference and recent, and trace and anisotropy models of interaction-induced polarizability. The spectral intensities obtained numerically at low frequencies are determined by bound and free transitions. At intermediate and high frequencies, the spectra are sensitive to both the attractive part of the potential and to short-range values of the trace or anisotropy. An empirical interatomic potential for the krypton gas interaction is developed by simultaneously fitting the Barker et al. (BFW) and modified Tang-Toennies (MTT) potentials to thermophysical and transport properties over a wide temperature range. The quality of the present potentials was checked by comparison between the calculated and experimental vibrational energy levels. The results show that these are the most accurate potentials reported to date for this system, both for the reproduction of spectral line shapes and spectral moments.

Ion Correlation in Choline Chloride-Urea Deep Eutectic Solvent (Reline) from Polarizable Molecular Dynamics Simulations.

In recent years, deep eutectic solvents (DESs) emerged as highly tunable and environmentally friendly alternatives to common ionic liquids and organic solvents. In this work, a polarizable model based on the CHARMM Drude polarizable force field is developed for a 1:2 ratio mixture of choline chloride/urea (reline) DES. To successfully reproduce the structure of the liquid as compared to first-principles molecular dynamics simulations, a damping factor was introduced to correct the observed over-binding between the chloride and the hydrogen bonding site of choline. Investigated radial distributions reveal the formation of hydrogen bonds between all the constituents of reline and similar interactions for chloride and urea's oxygen atoms, which could contribute to the melting point depression of the mixture. Predicted dynamic properties from our polarizable force field are in good agreement with experiments, showing significant improvements over nonpolarizable models. Similar to some ionic liquids, an oscillatory behavior in the velocity autocorrelation function of the anion is visible, which can be interpreted as a rattling motion of the lighter anion surrounded by the heavier cations. The obtained results for ionic conductivity of reline show some degree of correlated ion motion in this DES. However, a joint diffusion of ion pairs cannot be observed during the simulations.

Modified Poisson–Boltzmann equations and macroscopic forces in inhomogeneous ionic fluids

We propose a field-theoretical approach based on the thermodynamic perturbation theory and within it derive a grand thermodynamic potential of the inhomogeneous ionic fluid as a functional of electrostatic potential for an arbitrary reference fluid system. We obtain a modified Poisson–Boltzmann (PB) equation as the Euler–Lagrange equation for the obtained functional. Applying Noether’s theorem to this functional, we derive a general mean-field expression for the stress tensor consistent with the respective modified PB equation. We derive a general expression for the macroscopic force acting on the dielectric or conductive body immersed in an ionic fluid. In particular, we derive a general mean-field expression for the disjoining pressure of an ionic fluid in a slit pore. We apply the developed formalism to describe three ionic fluid models of practical importance: nonpolarizable models (including the well-known PB and Poisson–Fermi equations), polarizable models (ions carry nonzero permanent dipole or static polarizability), and models of ion-dipole mixtures (including the well-known PB–Langevin equation). For these models, we obtain modified PB equations and respective stress tensors, which could be valuable for different applications, where it is necessary to estimate the macroscopic forces acting on the dielectric or conductive bodies (electrodes, colloids, membranes, etc) together with the local electrostatic potential (field) and ionic concentrations.

Influence of Polarizability on the Structure, Dynamic Characteristics, and Ion-Transport Mechanisms in Polymeric Ionic Liquids.

We used atomistic simulations and compared the prediction of three different implementations of force fields, namely, the original full partial charge system, the scaled partial charge system, and the Drude oscillator polarizable force field and its effect on the structural and dynamic properties of a polymeric ionic liquid, poly(1-butyl-3-methyl-imidazolium hexafluorophosphate). We found that both the scaled and the polarizable force field models yield comparable predictions of structural and dynamic properties, although the scaled charge model artificially lowers the first-neighbor peak of the radial distribution function and therefore leads to a slight reduction in density. The full charge model was not accurate in its prediction of the dynamic properties but could reproduce the structural properties. With a refined analysis method for the ion-hopping mechanisms, we found that all three methods produce very similar conclusions, namely, that the mobile anion is associated with three cations from two distinct polymer chains and that the fractions of inter- and intramolecular hopping events are comparable. Our results demonstrate that the scaled charge force fields provide a computationally efficient means to capture polarizability effects on both the structural and dynamic properties of polymeric ionic liquid systems.

Stress tensor and constant pressure simulation for polarizable Gaussian multipole model.

Our previous article has established the theory of molecular dynamics (MD) simulations for systems modeled with the polarizable Gaussian multipole (pGM) electrostatics [Wei et al., J. Chem. Phys. 153(11), 114116 (2020)]. Specifically, we proposed the covalent basis vector framework to define the permanent multipoles and derived closed-form energy and force expressions to facilitate an efficient implementation of pGM electrostatics. In this study, we move forward to derive the pGM internal stress tensor for constant pressure MD simulations with the pGM electrostatics. Three different formulations are presented for the flexible, rigid, and short-range screened systems, respectively. The analytical formulations were implemented in the SANDER program in the Amber package and were first validated with the finite-difference method for two different boxes of pGM water molecules. This is followed by a constant temperature and constant pressure MD simulation for a box of 512 pGM water molecules. Our results show that the simulation system stabilized at a physically reasonable state and maintained the balance with the externally applied pressure. In addition, several fundamental differences were observed between the pGM and classic point charge models in terms of the simulation behaviors, indicating more extensive parameterization is necessary to utilize the pGM electrostatics.

Molecular dynamics calculations of partial molar volumes of amino acids in aqueous solutions

Partial molar volumes of amino acids in their zwitterionic and molecular forms were calculated using molecular dynamics simulations of these systems in aqueous solutions. Calculations performed with the TIP4P, SPC (rigid and flexible), SPC/E, and polarizable water models showed that the choice of water model can have a significant impact on the calculated volumes. The effect of treatment of long-range electrostatic interactions on the calculated results was also investigated. The volumes obtained in simulations with a properly chosen water model fit well the experimental data for both molecular and zwitterionic forms of amino acids.

High-Dimensional Parameter Search Method to Determine Force Field Mixing Terms in Molecular Simulations.

Molecular dynamics (MD) force fields for lipids and ions are typically developed independently of one another. In simulations consisting of both lipids and ions, lipid-ion interaction energies are estimated using a predefined set of mixing rules for Lennard-Jones (LJ) interactions. This, however, does not guarantee their reliability. In fact, compared to the quantum mechanical reference data, Lorentz-Berthelot mixing rules substantially underestimate the binding energies of Na+ ions with small-molecule analogues of lipid headgroups, yielding errors on the order of 80 and 130 kJ/mol, respectively, for methyl acetate and diethyl phosphate. Previously, errors associated with mixing force fields have been reduced using approaches such as "NB-fix" in which LJ interactions are computed using explicit cross terms rather than those from mixing rules. Building on this idea, we derive explicit lipid-ion cross terms that also may implicitly include many-body cooperativity effects. Additionally, to account for the interdependency between cross terms, we optimize all cross terms simultaneously by performing high-dimensional searches using our ParOpt software. The cross terms we obtain reduce the errors due to mixing rules to below 10 kJ/mol. MD simulation of the lipid bilayer conducted using these optimized cross terms resolves the structural discrepancies between our previous simulations and small-angle X-ray and neutron scattering experiments. These results demonstrate that simulations of lipid bilayers with ions that are accurate up to structural data from scattering experiments can be performed without explicit polarization terms. However, it is worth noting that such NB-fix cross terms are not based on any physical principle; a polarizable lipid model would be more realistic and is still desired. Our approach is generic and can be applied to improve the accuracies of simulations employing mixed force fields.

A transferable explicit-solvent polarizable coarse-grained model for proteins

The application of classical molecular dynamics (MD) simulations at atomic resolution (fine-grained (FG) level), to the majority of biomolecular processes, remains limited because of the associated computational complexity of representing all the atoms. This problem is magnified in presence of protein-based biomolecular systems that have a very large conformational space and MD simulations with FG resolution have slow dynamics to explore this space. Current transferrable coarse-grained (CG) force fields in literature are either limited to only peptides with the environment encoded in an implicit form (cannot study environmental heterogeneity) or cannot capture transitions into secondary/tertiary peptide structures from a primary sequence of amino acids.

PyRESP: A flexible program for polarizable force field parameterizations

Molecular modeling on atomic level have been applied in a wide range of biological systems. The widely adopted additive force fields typically use fixed atom-centered partial charges to model electrostatic interactions. However, the additive force fields cannot accurately model polarization effects, leading to unrealistic simulations in polarization sensitive processes. Numerous attempts have been directed to developing induced dipole based polarizable models, including the Applequist point dipole model, the Thole damped dipole model, and the recently proposed polarizable Gaussian multipole (pGM) model.

Molecular understanding of aqueous electrolyte properties and dielectric effect in a CDI system

Capacitive deionization (CDI) provides performance of ion removal via electrodes with high porosity and confined nanopores. Hitherto, the main emphasis has been on the enhancement of removal efficiency based on a lab-scale setup, neglecting anomalous dynamic and dielectric properties in confinement. This work was completed to reveal structural and dynamic properties of hydrated ions in a 3-nm nanopore of a CDI electrode using molecular dynamics simulation based on the nonpolarizable (SPC/E) and polarizable (SWM4-DP/-NDP) force fields. It was shown that all the force fields caused marginal differences in hydration structure but greater differences in ionic dynamics. Furthermore, dielectric constants obtained from all the force fields were lower in confinement than in bulk, although they were slightly higher in the polarizable models. Angular distributions of water bonds and dipole vectors revealed this decrease was caused by the restriction of water rotational dynamics and hence a depression of dipolar fluctuations. Hydrogen bond lifetimes indicated that higher availability of hydrogen bonding sites also facilitated increasing the dielectric constant as it allowed a more frequent “switch” of the fast hydrogen bond. Both the angular distribution and the hydrogen bonding contributed to the slightly higher dielectric constant in the polarizable model than in the nonpolarizable model. Our local sensitivity analysis of a ion transport model revealed that the dielectric constant greatly determined ion transport in the electrode, although the diffusion coefficients of typical ions had limited influences. This work provides fundamental theories and great implications for the theoretical analysis of ion electrosorption in a CDI electrode and enhances its application in water treatment.

Advances in Molecular Modeling of Ion-Protein Interaction Systems Towards Accurate Electrostatics: Methods and Applications

Metal ions are ubiquitous in complex with proteins and play key roles in protein structure and function. The ion-protein interactions are electrostatic delicate in nature, however, the description of electrostatic interactions could be problematic in conventional additive fixed-charge force fields. With empowered computational sources going beyond the common approximations, many efforts have been done to take account in more elaborate electrostatic description in molecular modeling using more sophisticated physical models and dilated algorithms/implementations. Here we review rencent progress in advanced polarizable models and new impletment approaches towards accurate electrostatics and highlight some successful application cases in ion-protein interaction systems in recent years.

Computational and Experimental Characterization of rDNA and rRNA G-Quadruplexes.

DNA G-quadruplexes in human telomeres and gene promoters are being extensively studied for their role in controlling the growth of cancer cells. G-quadruplexes have been unambiguously shown to exist both in vitro and in vivo, including in the guanine (G)-rich DNA genes encoding pre-ribosomal RNA (pre-rRNA), which is transcribed in the cell's nucleolus. Recent studies strongly suggest that these DNA sequences ("rDNA"), and the transcribed rRNA, are a potential anticancer target through the inhibition of RNA polymerase I (Pol I) in ribosome biogenesis, but the structures of ribosomal G-quadruplexes at atomic resolution are unknown and very little biophysical characterization has been performed on them to date. In the present study, circular dichroism (CD) spectroscopy is used to show that two putative rDNA G-quadruplex sequences, NUC 19P and NUC 23P and their counterpart rRNAs, predominantly adopt parallel topologies, reminiscent of the analogous telomeric quadruplex structures. Based on this information, we modeled parallel topology atomistic structures of the putative ribosomal G-quadruplexes. We then validated and refined the modeled ribosomal G-quadruplex structures using all-atom molecular dynamics (MD) simulations with the CHARMM36 force field in the presence and absence of stabilizing K+. Motivated by preliminary MD simulations of the telomeric parallel G-quadruplex (TEL 24P) in which the K+ ion is expelled, we used updated CHARMM36 force field K+ parameters that were optimized, targeting the data from quantum mechanical calculations and the polarizable Drude model force field. In subsequent MD simulations with optimized CHARMM36 parameters, the K+ ions are predominantly in the G-quadruplex channel and the rDNA G-quadruplexes have more well-defined, predominantly parallel-topology structures as compared to rRNA. In addition, NUC 19P is more structured than NUC 23P, which contains extended loops. Results from this study set the structural foundation for understanding G-quadruplex functions and the design of novel chemotherapeutics against these nucleolar targets and can be readily extended to other DNA and RNA G-quadruplexes.

OS100: A Benchmark Set of 100 Digitized UV-Visible Spectra and Derived Experimental Oscillator Strengths.

Excited-state quantum chemical calculations usually report excitation energies and oscillator strengths, f, for each electronic transition. On the other hand, UV-visible spectrophotometric experiments measure energy-dependent molar extinction/attenuation coefficients, ε(v), that give absorption band line shapes when plotted. ε(v) and f are related, but this relation is complicated by broadening and solvation effects. We fitted and integrated 100 experimental UV-visible spectra to obtain 164 fexp values for absorption bands appearing in these spectra. The 100 UV-visible spectra belong to solvated organic molecules ranging in size from 6-34 atoms. We estimated uncertainties in the fitting to indicate confidence level in the reported fexp values. The corresponding computed oscillator strengths (fcomp) were obtained with time-dependent density functional theory and a polarizable continuum solvent model. By expressing experimental and computed absorption strengths using a common quantity, we directly compared fcomp and fexp. Although fcomp and fexp are well correlated (linear regression R2 = 0.921), fcomp in most cases overestimated fexp (regression slope = 1.34). The agreement between absolute fcomp and fexp values was substantially improved by accounting for a solvent refractive index factor, as suggested in some derivations in the literature. The 100 digitized UV-visible spectra are included as plain text files in the Supporting Information to aid in benchmarking computational or machine learning methods that aim to simulate realistic UV-visible absorption spectra.