Polarizable Model Research Articles

Coupled cluster theory with the polarizable continuum model of solvation

AbstractThe environment may significantly affect molecular properties. Thus, it is desirable to account explicitly for these effects on the wave function and its derivatives, especially when the latter are evaluated with accurate methods, such as those belonging to coupled cluster (CC) theory. In this tutorial review, we discuss how to combine CC methods with the polarizable continuum model of solvation (PCM). We describe useful approximations that include the solvent response to the correlation and excited state equations while maintaining the computational cost comparable to in vacuo calculations. Although applied to PCM, the theoretical framework presented in this review is general and can be used with any polarizable embedding model. Representative applications of the CC‐PCM method to ground and excited state properties of solvated molecules are presented, and comparisons with experiment, and between the full and approximate schemes are discussed.

Calculation of the surface tension of water: 40 years of molecular simulations

ABSTRACTWe focus here on the calculation of the surface tension of water by molecular simulations of planar interfaces. We provide a chronological list of the key dates in the calculation of the surface tension of water by using atomistic empirical models, polarisable, many-body and coarse-grained (CG) models. We investigate how the SPC/E, TIP4P and TIP4P/2005 atomistic force fields perform in the reproduction of the surface tension along the coexistence curve. We also study the temperature dependence of the surface tension of water for CG models. We also discuss the ability of ab initio molecular dynamics simulations to reproduce the surface tension of water at room temperature.

Polarized Raman spectroscopy and lattice dynamics of potassic-magnesio-arfvedsonite

We report polarized Raman spectra from potassic-magnesio-arfvedsonite in all informative scattering configurations. On the basis of the polarization selection rules, several Ag vibrational modes have been identified. The Bg modes, however, are below the detection limits of the Raman spectrometer. The OH stretching band is situated between 3630 and 3750 cm−1, and its spectral shape is typical of amphiboles with high occupancy of the A site. It is composed of seven overlapping but resolvable subbands, which stem from occupied A-site configurations M(1)M(1)M(3)–OH–A(K/Na)–WOH and M(1)M(1)M(3)–OH–A(K/Na)–WF, as well as from vacant A-site configurations M(1)M(1)M(3)–OH–A□–WOH, with different Mg and Fe occupancy of the M(1) and M(3) sites. The experimental Raman spectra are compared with the results of theoretical calculations based on a shell-model force-field and a bond polarizability model. The simulated partial Raman spectra allowed us to assign many low-frequency Raman bands to stretching vibrations involving specific cation-oxygen bonds, as well as the higher-frequency modes of the Si–O skeleton. On the basis of our calculations we hypothesize that the Raman bands at 467, 540 and 589 cm−1 are related to a superposition of M(2)Fe3+–O bond stretching and Si–O–Si bending vibrations.

Structural study of Na2O-B2O3-SiO2-La2O3 glasses from molecular simulations using a polarizable force field

The properties of Na2O-B2O3-SiO2+La2O3 glasses (NBS_La) are not well understood. In particular the influence of the lanthanides on the polymerized network structure and organization are not well known. In this work, we use a Polarizable Ion Model (PIM) that was fitted on electronic structure calculations to simulate a series of NBS_La glasses using molecular dynamics. The simulations account well for the main structural characteristics of NBS_La glasses such as the structure factors and the ionic local environments despite minor discrepancies with experiments about the Si-O distance and the B coordination. In particular, we examine the impact of La2O3 addition to Na2O-B2O3-SiO2 glasses in term of structural changes and in the Na migration paths which become more homogeneous.

The hydration structure of methylthiolate from QM/MM molecular dynamics.

Thiols are widely present in biological systems, most notably as the side chain of cysteine amino acids in proteins. Thiols can be deprotonated to form a thiolate which affords a diverse range of enzymatic activity and modes for chemical modification of proteins. Parameters for modeling thiolates using molecular mechanical force fields have not yet been validated, in part due to the lack of structural data on thiolate solvation. Here, the CHARMM36 and Amber models for thiolates in aqueous solutions are assessed using free energy perturbation and hybrid quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations. The hydration structure of methylthiolate was calculated from 1 ns of QM/MM MD (PBE0-D3/def2-TZVP//TIP3P), which shows that the water-S- distances are approximately 2 Å with a coordination number near 6. The CHARMM thiolate parameters predict a thiolate S radius close to the QM/MM value and predict a hydration Gibbs energy of -329.2 kJ/mol, close to the experimental value of -318 kJ/mol. The cysteine thiolate model in the Amber force field underestimates the thiolate radius by 0.2 Å and overestimates the thiolate hydration energy by 119 kJ/mol because it uses the same Lennard-Jones parameters for thiolates as for thiols. A recent Drude polarizable model for methylthiolate with optimized thiolate parameters also performs well. SAPT2+ [Symmetry Adapted Perturbation Theory (SAPT)] analysis indicates that exchange repulsion is larger for the methylthiolate, consistent with it having a more diffuse electron density distribution in comparison with the parent thiol. These data demonstrate that it is important to define distinct non-bonded parameters for the protonated/deprotonated states of amino acid side chains in molecular mechanical force fields.

An efficient water force field calibrated against intermolecular THz and Raman spectra.

A polarizable water model is presented which has been calibrated against experimental THz and Raman spectra of bulk water. These low-frequency spectra directly probe the dynamics, and thereby intermolecular interactions, on time scales relevant to molecular motions. The model is based on the TL4P force field developed recently by Tavan and co-workers [J. Phys. Chem. B 117, 9486 (2013)], which has been designed to be transferable between different environments; in particular, to correctly describe the electrostatic properties of both the isolated water molecule in the gas-phase and the liquid water at ambient conditions. Following this design philosophy, TL4P was amended with charge transfer across hydrogen-bonded dimers as well as an anisotropic polarizability in order to correctly reproduce the THz and Raman spectra. The thermodynamic and structural properties of the new model are of equal quality as those of TL4P, and at the same time, an almost quantitative agreement with the spectroscopic data could be achieved. Since TL4P is a rigid model with a single polarizable site, it is computationally very efficient, while the numerical overhead for the addition of charge transfer and the anisotropic polarizability is minor. Overall, the model is expected to be well suited for, e.g., large scale simulations of 2D-Raman-THz spectra or biomolecular simulations.

Polarizable Molecular Simulations Reveal How Silicon-Containing Functional Groups Govern the Desalination Mechanism in Nanoporous Graphene.

We report a molecular dynamics (MD) simulation study of reverse osmosis desalination using nanoporous monolayer graphene passivated by SiH2 and Si(OH)2 functional groups. A highly accurate and detailed polarizable molecular mechanics force field model was developed for simulating graphene nanopores of various sizes and geometries. The simulated water fluxes and ion rejection percentages are explained using detailed atomistic mechanisms derived from analysis of the simulation trajectories. Our main findings are as follows: (1) The Si(OH)2 pores possess superior ion rejection rates due to selective electrostatic repulsion of Cl- ions, but Na+ ions are attracted to the pore and block water transfer. (2) By contrast, the SiH2 pores operate via a steric mechanism that excludes ions based on the size and flexibility of their hydration layers. (3) In the absence of ions, water flux is directly proportional to the solvent accessible area within the pore; however, simulated fluxes are lower than those inferred from recent experimental work. We also provide some hypotheses that could resolve the differences between simulation and experiment.

Quantifying the Short-Range Order in Amorphous Silicon by Raman Scattering.

Quantification of the short-range order in amorphous silicon has been formulized using Raman scattering by taking into account established frameworks for studying the spectral line-shape and size dependent Raman peak shift. A theoretical line-shape function has been proposed for representing the observed Raman scattering spectrum from amorphous-Si-based on modified phonon confinement model framework. While analyzing modified phonon confinement model, the term "confinement size" used in the context of nanocrystalline Si was found analogous to the short-range order distance in a-Si thus enabling one to quantify the same using Raman scattering. Additionally, an empirical formula has been proposed using bond polarizability model for estimating the short-range order making one capable to quantify the distance of short-range order by looking at the Raman peak position alone. Both the proposals have been validated using three different data sets reported by three different research groups from a-Si samples prepared by three different methods making the analysis universal.

Molecular Treatment of Nano-Kaolinite Generations.

A procedure is developed for defining a compositionally and structurally realistic, atomic-scale description of exfoliated clay nanoparticles from the kaolinite family of phylloaluminosilicates. By use of coordination chemical principles, chemical environments within a nanoparticle can be separated into inner, outer, and peripheral spheres. The edges of the molecular models of nanoparticles were protonated in a validated manner to achieve charge neutrality. Structural optimizations using semiempirical methods (NDDO Hamiltonians and DFTB formalism) and ab initio density functionals with a saturated basis set revealed previously overlooked molecular origins of morphological changes as a result of exfoliation. While the use of semiempirical methods is desirable for the treatment of nanoparticles composed of tens of thousands of atoms, the structural accuracy is rather modest in comparison to DFT methods. We report a comparative survey of our infrared data for untreated crystalline and various exfoliated states of kaolinite and halloysite. Given the limited availability of experimental techniques for providing direct structural information about nano-kaolinite, the vibrational spectra can be considered as an essential tool for validating structural models. The comparison of experimental and calculated stretching and bending frequencies further justified the use of the preferred level of theory. Overall, an optimal molecular model of the defect-free, ideal nano-kaolinite can be composed with respect to stationary structure and curvature of the potential energy surface using the PW91/SVP level of theory with empirical dispersion correction (PW91+D) and polarizable continuum solvation model (PCM) without the need for a scaled quantum chemical force field. This validated theoretical approach is essential in order to follow the formation of exfoliated clays and their surface reactivity that is experimentally unattainable.

Pseudo-3D constrained inversion of transient electromagnetic data for a polarizable SMS hydrothermal system in the Deep Sea

Seafloor massive sulfide (SMS) deposits are generated by high-temperature hydrothermal systems. Their precious resources have attracted global interest. A number of investigations with controlled-source electromagnetic (CSEM) methods have been implemented in recent years. There are three major problems with SMS surveying using EM methods. First, SMS imaging techniques for hydrothermal systems have a limited range. Simulations and applications have validated only simple layered models. Second, their inversion efficiencies must be improved further. Laterally constrained inversions and spatially constrained inversions are usually used to map geological structures. However, choosing their suitable weighting parameters is inefficient. Third, the effects of induced polarization (IP) on ore deposits are not considered in such inversions. A non-polarizable model is unable to accurately depict a polarizable model. To resolve these problems, an advanced strategy is used to improve the efficiency of the pseudo-3D inversion process. The proposed imaging method has the ability to map complex 3D geoelectrical structures, and therefore, it can both obtain information regarding surface ore deposits and distinguish between active and inactive hydrothermal systems. However, this method can also be used to depict the distributions of alteration zones and buried deposits. Furthermore, the influences of IP on the inversion are discussed with respect to the Cole- Cole model, and it is shown that the effects of IP on polarizable deposits cannot be ignored during the inversion.

Insights into Stabilizing Forces in Amyloid Fibrils of Differing Sizes from Polarizable Molecular Dynamics Simulations

Pathological aggregation of amyloid-forming proteins is a hallmark of a number of human diseases, including Alzheimer's, type 2 diabetes, Parkinson's, and more. Despite having very different primary amino acid sequences, these amyloid proteins form similar supramolecular, fibril structures that are highly resilient to physical and chemical denaturation. To better understand the structural stability of disease-related amyloids and to gain a greater understanding of factors that stabilize functional amyloid assemblies, insights into tertiary and quaternary interactions are needed. We performed molecular dynamics simulations on human tau, amyloid-β, and islet amyloid polypeptide fibrils to determine key physicochemical properties that give rise to their unique characteristics and fibril structures. These simulations are the first of their kind in employing a polarizable force field to explore properties of local electric fields on dipole properties and other electrostatic forces that contribute to amyloid stability. Across these different amyloid fibrils, we focused on how the underlying forces stabilize fibrils to elucidate the driving forces behind the protein aggregation. The polarizable model allows for an investigation of how side-chain dipole moments, properties of structured water molecules in the fibril core, and the local environment around salt bridges contribute to the formation of interfaces essential for fibril stability. By systematically studying three amyloidogenic proteins of various fibril sizes for key structural properties and stabilizing forces, we shed light on properties of amyloid structures related to both diseased and functional states at the atomistic level.

Classical Polarizable Force Field to Study Hydrated Hectorite: Optimization on DFT Calculations and Validation against XRD Data

Following our previous works on dioctahedral clays, we extend the classical Polarizable Ion Model (PIM) to trioctahedral clays, by considering dry Na-, Cs-, Ca- and Sr-hectorites as well as hydrated Na-hectorite. The parameters of the force field are determined by optimizing the atomic forces and dipoles on density functional theory calculations. The simulation results are validated by comparison with experimental X-ray diffraction (XRD) data. The XRD patterns calculated from classical molecular dynamics simulations performed with the PIM force field are in very good agreement with experimental results. In the bihydrated state, the less structured electronic density profile obtained with PIM compared to the one from the state-of-the-art non-polarizable force field clayFF explains the slightly better agreement between the PIM results and experiments.

CHARMM Drude Polarizable Force Field for Glycosidic Linkages Involving Pyranoses and Furanoses.

We present an extension of the CHARMM Drude polarizable force field to enable modeling of polysaccharides containing pyranose and furanose monosaccharides. The new force field parameters encompass 1↔2, 1→3, 1→4, and 1→6 pyranose-furanose linkages, 2→1 and 2→6 furanose-furanose linkages, 2→2, 2→3, and 2→4 furanose-pyranose, and 1↔1, 1→2, 1→3, 1→4, and 1→6 pyranose-pyranose linkages. For the glycosidic linkages, both simple model compounds and the full disaccharides with methylation at the reducing end were used for parameter optimization. The model compounds were chosen to be monomers or glycosidic-linked dimers of tetrahydropyran (THP) and tetrahydrofuran (THF). Target data for optimization included one- and two-dimensional potential energy scans of ω and the Φ/Ψ glycosidic dihedral angles in the model compounds and full disaccharides computed by quantum mechanical (QM) RIMP2/cc-pVQZ single point energies on MP2/6-31G(d) optimized structures. Also included in the target data are extensive sets of QM gas phase monohydrate water-saccharide interactions, dipole moments, and molecular polarizabilities for both model compounds and full disaccharides. The resulting polarizable model is shown to be in good agreement with a range of QM data, offering a significant improvement over the additive CHARMM36 carbohydrate force field, as well as experimental data including crystal structures and conformational properties of disaccharides and a trisaccharide in aqueous solution.

Quantum Chemical Calculation of p Kas of Environmentally Relevant Functional Groups: Carboxylic Acids, Amines, and Thiols in Aqueous Solution.

Developing accurate quantum chemical approaches for calculating p Kas is of broad interest. Useful accuracy can be obtained by using density functional theory (DFT) in combination with a polarizable continuum solvent model. However, some classes of molecules present problems for this approach, yielding errors greater than 5 p K units. Various methods have been developed to improve the accuracy of the combined strategy. These methods perform well but either do not generalize or introduce additional degrees of freedom, increasing the computational cost. The Solvation Model based on Density (SMD) has emerged as one of the most commonly used continuum solvent models. Nevertheless, for some classes of organic compounds, e.g., thiols, the p Kas calculated with the original SMD model show errors of 6-10 p K units, and we traced these errors to inaccuracies in the solvation free energies of the anions. To improve the accuracy of p Kas calculated with DFT and the SMD model, we developed a scaled solvent-accessible surface approach for constructing the solute-solvent boundary. By using a "direct" approach, in which all quantities are computed in the presence of the continuum solvent, the use of thermodynamic cycles is avoided. Furthermore, no explicit water molecules are required. Three benchmark data sets of experimentally measured p Ka values, including 28 carboxylic acids, 10 aliphatic amines, and 45 thiols, were used to assess the optimized SMD model, which we call SMD with a scaled solvent-accessible surface (SMDsSAS). Of the methods tested, the M06-2X density functional approximation, 6-31+G(d,p) basis set, and SMDsSAS solvent model provided the most accurate p Kas for each set, yielding mean unsigned errors of 0.9, 0.4, and 0.5 p K units, respectively, for carboxylic acids, aliphatic amines, and thiols. This approach is therefore useful for efficiently calculating the p Kas of environmentally relevant functional groups.

Polarizable Force Field for Molecular Ions Based on the Classical Drude Oscillator.

Development of accurate force field parameters for molecular ions in the context of a polarizable energy function based on the classical Drude oscillator is a crucial step toward an accurate polarizable model for modeling and simulations of biological macromolecules. Toward this goal we have undertaken a hierarchical approach in which force field parameter optimization is initially performed for small molecules for which experimental data exists that serve as building blocks of macromolecular systems. Small molecules representative of the ionic moieties of biological macromolecules include the cationic ammonium and methyl substituted ammonium derivatives, imidazolium, guanidinium and methylguanidinium, and the anionic acetate, phenolate, and alkanethiolates. In the present work, parameters for molecular ions in the context of the Drude polarizable force field are optimized and compared to results from the nonpolarizable additive CHARMM general force field (CGenFF). Electrostatic and Lennard-Jones parameters for the model compounds are developed in the context of the polarizable SWM4-NDP water model, with emphasis on assuring that the hydration free energies are consistent with previously reported parameters for atomic ions. The final parameters are shown to be in good agreement with the selected quantum mechanical (QM) and experimental target data. Analysis of the structure of water around the ions reveals substantial differences between the Drude and additive force fields indicating the important role of polarization in dictating the molecular details of aqueous solvation. The presented parameters represent the foundation for the charged functionalities in future generations of the Drude polarizable force field for biological macromolecules as well as for drug-like molecules.

Structure of supercooled water: Polarizable BK3 model versus non-polarizable models

The polarizable BK3 model of water was simulated, using both Monte Carlo and molecular dynamics methods, along isobars p = 1 and 1000 bar for temperatures ranging from supercritical down to 228.15 K with the focus on details of its structure and changes with temperature. Since overwhelming majority of results for supercooled model water have been obtained for non-polarizable models, in most cases the ST2 one, we also simulated at corresponding conditions this model along with another, and qualitatively different non-polarizable model, TIP4P/Ice for comparison. In addition to commonly used various order parameters to characterize the structure, and a detailed analysis of the spatial distribution function over the first three coordination shells, we also employ a geometric construction based on the Delaunay tessellation, the method independent of any a priori defined quantity. All these analyses aim at an examination to what extent qualitative differences in the models (force fields) affect the observed structure. It is found that in most cases the polarizable BK3 model and TIP4P/Ice water yield very similar structural properties. On the contrary, ST2 water exhibits, in majority cases, even qualitatively different properties and these findings are discussed in detail. The obtained results suggest that any conclusion on the behavior of supercooled water based on the ST2 model should be taken, at least, with caution.

Complexation reactions in pyridine and 2,6-dimethylpyridine-water system: The quantum-chemical description and the path to liquid phase separation.

Phase separation in substituted pyridines in water is usually described as an interplay between temperature-driven breakage of hydrogen bonds and the associating interaction of the van der Waals force. In previous quantum-chemical studies, the strength of hydrogen bonding between one water and one pyridine molecules (the 1:1 complex) was assigned a pivotal role. It was accepted that the disassembly of the 1:1 complex at a critical temperature leads to phase separation and formation of the miscibility gap. Yet, for over two decades, notable empirical data and theoretical arguments were presented against that view, thus revealing the need in a revised quantum-mechanical description. In the present study, pyridine-water and 2,6-dimethylpyridine-water systems at different complexation stages are calculated using high level Kohn-Sham theory. The hydrophobic-hydrophilic properties are accounted for by the polarizable continuum solvation model. Inclusion of solvation in free energy of formation calculations reveals that 1:1 complexes are abundant in the organically rich solvents but higher level oligomers (i.e., 2:1 dimers with two pyridines and one water molecule) are the only feasible stable products in the more polar media. At the critical temperature, the dissolution of the external hydrogen bonds between the 2:1 dimer and the surrounding water molecules induces the demixing process. The 1:1 complex acts as a precursor in the formation of the dimers but is not directly involved in the demixing mechanism. The existence of the miscibility gap in one pyridine-water system and the lack of it in another is explained by the ability of the former to maintain stable dimerization. Free energy of formation of several reaction paths producing the 2:1 dimers is calculated and critically analyzed.

Structure and Raman Spectra of C60 and C70 Fullerenes Encased into Single-Walled Boron Nitride Nanotubes: A Theoretical Study

We report the structures and the nonresonant Raman spectra of hybrid systems composed of carbon fullerenes ( C 60 and C 70 ) encased within single walled boron nitride nanotube. The optimal structure of these systems are derived from total energy minimization using a convenient Lennard-Jones expression of the van der Waals intermolecular potential. The Raman spectra have been calculated as a function of nanotube diameter and fullerene concentration using the bond polarizability model combined with the spectral moment method. These results should be useful for the interpretation of the experimental Raman spectra of boron nitride nanotubes encasing C 60 and C 70 fullerenes.

One Solution of the Forward Problem of DC Resistivity Well Logging by the Method of Volume Integral Equations with Allowance for Induced Polarization

For theoretically studying the intensity of the influence exerted by the polarization of the rocks on the results of direct current (DC) well logging, a solution is suggested for the direct inner problem of the DC electric logging in the polarizable model of plane-layered medium containing a heterogeneity by the example of the three-layer model of the hosting medium. Initially, the solution is presented in the form of a traditional vector volume-integral equation of the second kind (IE2) for the electric current density vector. The vector IE2 is solved by the modified iteration–dissipation method. By the transformations, the initial IE2 is reduced to the equation with the contraction integral operator for an axisymmetric model of electrical well-logging of the three-layer polarizable medium intersected by an infinitely long circular cylinder. The latter simulates the borehole with a zone of penetration where the sought vector consists of the radial Jr and Jz axial (relative to the cylinder’s axis) components. The decomposition of the obtained vector IE2 into scalar components and the discretization in the coordinates r and z lead to a heterogeneous system of linear algebraic equations with a block matrix of the coefficients representing 2x2 matrices whose elements are the triple integrals of the mixed derivatives of the second-order Green’s function with respect to the parameters r, z, r', and z'. With the use of the analytical transformations and standard integrals, the integrals over the areas of the partition cells and azimuthal coordinate are reduced to single integrals (with respect to the variable t = cos ϕ on the interval [−1, 1]) calculated by the Gauss method for numerical integration. For estimating the effective coefficient of polarization of the complex medium, it is suggested to use the Siegel–Komarov formula.

Cation-π Interactions between Quaternary Ammonium Ions and Amino Acid Aromatic Groups in Aqueous Solution.

Cation-π interactions play important roles in the stabilization of protein structures and protein-ligand complexes. They contribute to the binding of quaternary ammonium ligands (mainly RNH3+ and RN(CH3)3+) to various protein receptors and are likely involved in the blockage of potassium channels by tetramethylammonium (TMA+) and tetraethylammonium (TEA+). Polarizable molecular models are calibrated for NH4+, TMA+, and TEA+ interacting with benzene, toluene, 4-methylphenol, and 3-methylindole (representing aromatic amino acid side chains) based on the ab initio MP2(full)/6-311++G(d,p) properties of the complexes. Whereas the gas-phase affinity of the ions with a given aromatic follows the trend NH4+ > TMA+ > TEA+, molecular dynamics simulations using the polarizable models show a reverse trend in water, likely due to a contribution from the hydrophobic effect. This reversed trend follows the solubility of aromatic hydrocarbons in quaternary ammonium salt solutions, which suggests a role for cation-π interactions in the salting-in of aromatic compounds in solution. Simulations in water show that the complexes possess binding free energies ranging from -1.3 to -3.3 kcal/mol (compared to gas-phase binding energies between -8.5 and -25.0 kcal/mol). Interestingly, whereas the most stable complexes involve TEA+ (the largest ion), the most stable solvent-separated complexes involve TMA+ (the intermediate-size ion).