
Communication: The Al + CO_{2} → AlO + CO reaction: Experiment vs. theory J. Chem. Phys. (IF 2.965) Pub Date : 20171102
Zhi Sun, Kevin B. MooreIII, Henry F. SchaeferIIIBased on their highly sophisticated crossedbeam experimental studies of the Al + CO2 → AlO + CO reaction, Honma and Hirata have directly challenged the results of earlier theoretical studies of this system. We report high level theoretical studies of this system. It is shown that, consistent with HonmaHirata experimental conclusions, the previous theoretical prediction of a substantial barrier height for this reaction was incorrect. However, for the structures of the possible intermediates, in agreement with the 1992 theoretical study of Sakai, we find striking disagreement with the experimental conclusion that the O–C–O moiety is nearly linear. The energies of the three entrance channel intermediates lie 14.4, 15.2, and 16.4 kcal mol−1 below separated Al + CO2.

Performing SELEX experiments in silico J. Chem. Phys. (IF 2.965) Pub Date : 20171101
J. A. J. Wondergem, H. Schiessel, M. TompitakDue to the sequencedependent nature of the elasticity of DNA, many proteinDNA complexes and other systems in which DNA molecules must be deformed have preferences for the type of DNA sequence they interact with. SELEX (Systematic Evolution of Ligands by EXponential enrichment) experiments and similar sequence selection experiments have been used extensively to examine the (indirect readout) sequence preferences of, e.g., nucleosomes (protein spools around which DNA is wound for compactification) and DNA rings. We show how recently developed computational and theoretical tools can be used to emulate such experiments in silico. Opening up this possibility comes with several benefits. First, it allows us a better understanding of our models and systems, specifically about the roles played by the simulation temperature and the selection pressure on the sequences. Second, it allows us to compare the predictions made by the model of choice with experimental results. We find agreement on important features between predictions of the rigid basepair model and experimental results for DNA rings and interesting differences that point out open questions in the field. Finally, our simulations allow application of the SELEX methodology to systems that are experimentally difficult to realize because they come with high energetic costs and are therefore unlikely to form spontaneously, such as very short or overwound DNA rings.

Theory of molecular nonadiabatic electron dynamics in condensed phases J. Chem. Phys. (IF 2.965) Pub Date : 20171101
Kazuo TakatsukaIn light of the rapid progress of ultrafast chemical dynamics driven by the pulse lasers having width as short as several tens of attoseconds, we herein develop a theory of nonadiabatic electron wavepacket dynamics in condensed phases, with which to directly track the dynamics of electronicstate mixing such as electron transfer in liquid solvents. Toward this goal, we combine a theory of pathbranching representation for nonadiabatic electron wavepacket dynamics in vacuum {a mixed quantumclassical representation, Yonehara and Takatsuka [J. Chem. Phys. 129, 134109 (2008)]} and a theory of entropy functional to treat chemical dynamics in condensed phases {a mixed dynamicalstatistical representation, Takatsuka and Matsumoto [Phys. Chem. Chem. Phys. 18, 1771 (2016)]}. Difficulty and complexity in the present theoretical procedure arise in embedding the Schrödinger equation into classically treated statistical environment. Nevertheless, the resultant equations of motion for electronicstate mixing due to the intrinsic nonadiabatic interactions and solutesolvent interactions, along with the force matrix that drives nuclear branching paths, both turn out to be clear enough to make it possible to comprehend the physical meanings behind. We also discuss briefly the nonvalidness of naive application of the notion of nonadiabatic transition dynamics among free energy surfaces.

Benchmark CCSDSAPT study of rare gas dimers with comparison to MPSAPT and DFTSAPT J. Chem. Phys. (IF 2.965) Pub Date : 20171102
Leonid Shirkov, Vladimir SladekSymmetryadapted perturbation theory (SAPT) based on coupled cluster approach with single and double excitations (CCSD) treatment of intramonomer electron correlation effects was applied to study rare gas homodimers from He2 to Kr2. The obtained benchmark CCSDSAPT energies, including cumulant contributions to first order exchange and secondorder exchangeinduction terms, were then compared to their counterparts found using other methods—MøllerPlessetSAPT based on manybody MøllerPlesset perturbation theory and DFTSAPT based on density functional theory. The SAPT terms up to the secondorder were calculated with the basis sets close to the complete basis set at the large range of interatomic distances R. It was shown that overestimation of the binding energies De found with DFTSAPT reported in the work of Shirkov and Makarewicz [J. Chem. Phys. 142, 064102 (2015)] for Ar2 and Kr2 is mostly due to underestimation of the exchange energy Eexch(1) when comparing to the CCSDSAPT benchmark. The CCSDSAPT potentials were found to give the following values of the dissociation energies D0: 0.0006 cm−1 for He2, 16.71 cm−1 for Ne2, 85.03 cm−1 for Ar2, and 129.81 cm−1 for Kr2, which agree well with the values found from previously reported highly accurate ab initio supermolecular potentials and experimental data. The longrange dispersion coefficients C2n up to n = 6 that give the dispersion energy asymptotically equivalent to its SAPT counterpart were calculated from dynamic multipole polarizabilities at different levels of theory.

Similarity transformed equation of motion coupledcluster theory based on an unrestricted HartreeFock reference for applications to highspin openshell systems J. Chem. Phys. (IF 2.965) Pub Date : 20171102
Lee M. J. Huntington, Martin Krupička, Frank Neese, Róbert IzsákThe similarity transformed equation of motion coupledcluster approach is extended for applications to highspin openshell systems, within the unrestricted HartreeFock (UHF) formalism. An automatic active space selection scheme has also been implemented such that calculations can be performed in a blackbox fashion. It is observed that both the canonical and automatic active space selecting similarity transformed equation of motion (STEOM) approaches perform about as well as the more expensive equation of motion coupledcluster singles doubles (EOMCCSD) method for the calculation of the excitation energies of doublet radicals. The automatic active space selecting UHF STEOM approach can therefore be employed as a viable, lower scaling alternative to UHF EOMCCSD for the calculation of excited states in highspin openshell systems.

Optimal updating magnitude in adaptive flatdistribution sampling J. Chem. Phys. (IF 2.965) Pub Date : 20171103
Cheng Zhang, Justin A. Drake, Jianpeng Ma, B. Montgomery PettittWe present a study on the optimization of the updating magnitude for a class of free energy methods based on flatdistribution sampling, including the WangLandau (WL) algorithm and metadynamics. These methods rely on adaptive construction of a bias potential that offsets the potential of mean force by histogrambased updates. The convergence of the bias potential can be improved by decreasing the updating magnitude with an optimal schedule. We show that while the asymptotically optimal schedule for the singlebin updating scheme (commonly used in the WL algorithm) is given by the known inversetime formula, that for the Gaussian updating scheme (commonly used in metadynamics) is often more complex. We further show that the singlebin updating scheme is optimal for very long simulations, and it can be generalized to a class of bandpass updating schemes that are similarly optimal. These bandpass updating schemes target only a few longrange distribution modes and their optimal schedule is also given by the inversetime formula. Constructed from orthogonal polynomials, the bandpass updating schemes generalize the WL and LangfeldLuciniRago algorithms as an automatic parameter tuning scheme for umbrella sampling.

Nonorthogonal internally contracted multiconfigurational perturbation theory (NICPT): Dynamic electron correlation for large, compact active spaces J. Chem. Phys. (IF 2.965) Pub Date : 20171106
Sven Kähler, Jeppe OlsenA computational method is presented for systems that require highlevel treatments of static and dynamic electron correlation but cannot be treated using conventional complete active space selfconsistent fieldbased methods due to the required size of the active space. Our method introduces an efficient algorithm for perturbative dynamic correlation corrections for compact nonorthogonal MCSCF calculations. In the algorithm, biorthonormal expansions of orbitals and CIwave functions are used to reduce the scaling of the performance determining step from quadratic to linear in the number of configurations. We describe a hierarchy of configuration spaces that can be chosen for the active space. Potential curves for the nitrogen molecule and the chromium dimer are compared for different configuration spaces. Already the most compact spaces yield qualitatively correct potentials that with increasing size of configuration spaces systematically approach complete active space results.

Delayed Slater determinant update algorithms for high efficiency quantum Monte Carlo J. Chem. Phys. (IF 2.965) Pub Date : 20171107
T. McDaniel, E. F. D’Azevedo, Y. W. Li, K. Wong, P. R. C. KentWithin ab initio Quantum Monte Carlo simulations, the leading numerical cost for large systems is the computation of the values of the Slater determinants in the trial wavefunction. Each Monte Carlo step requires finding the determinant of a dense matrix. This is most commonly iteratively evaluated using a rank1 ShermanMorrison updating scheme to avoid repeated explicit calculation of the inverse. The overall computational cost is, therefore, formally cubic in the number of electrons or matrix size. To improve the numerical efficiency of this procedure, we propose a novel multiple rank delayed update scheme. This strategy enables probability evaluation with an application of accepted moves to the matrices delayed until after a predetermined number of moves, K. The accepted events are then applied to the matrices en bloc with enhanced arithmetic intensity and computational efficiency via matrixmatrix operations instead of matrixvector operations. This procedure does not change the underlying Monte Carlo sampling or its statistical efficiency. For calculations on large systems and algorithms such as diffusion Monte Carlo, where the acceptance ratio is high, order of magnitude improvements in the update time can be obtained on both multicore central processing units and graphical processing units.

Relativistic effects on the NMR parameters of Si, Ge, Sn, and Pb alkynyl compounds: Scalar versus spinorbit effects J. Chem. Phys. (IF 2.965) Pub Date : 20171101
Taye B. DemissieThe NMR chemical shifts and indirect spinspin coupling constants of 12 molecules containing 29Si, 73Ge, 119Sn, and 207Pb [X(CCMe)4, Me2X(CCMe)2, and Me3XCCH] are presented. The results are obtained from nonrelativistic as well as two and fourcomponent relativistic density functional theory (DFT) calculations. The scalar and spin–orbit relativistic contributions as well as the total relativistic corrections are determined. The main relativistic effect in these molecules is not due to spin–orbit coupling but rather to the scalar relativistic contraction of the sshells. The correlation between the calculated and experimental indirect spin–spin coupling constants showed that the fourcomponent relativistic density functional theory (DFT) approach using the Perdew’s hybrid scheme exchangecorrelation functional (PBE0; using the PerdewBurkeErnzerhof exchange and correlation functionals) gives results in good agreement with experimental values. The indirect spinspin coupling constants calculated using the spinorbit zeroth order regular approximation together with the hybrid PBE0 functional and the specially designed Jcoupling (JCPL) basis sets are in good agreement with the results obtained from the fourcomponent relativistic calculations. For the coupling constants involving the heavy atoms, the relativistic corrections are of the same order of magnitude compared to the nonrelativistically calculated results. Based on the comparisons of the calculated results with available experimental values, the best results for all the chemical shifts and nonexisting indirect spin–spin coupling constants for all the molecules are reported, hoping that these accurate results will be used to benchmark future DFT calculations. The present study also demonstrates that the fourcomponent relativistic DFT method has reached a level of maturity that makes it a convenient and accurate tool to calculate indirect spin–spin coupling constants of “large” molecular systems involving heavy atoms.

Correlated electron dynamics in strongfield nonsequential double ionization of Mg J. Chem. Phys. (IF 2.965) Pub Date : 20171106
Ning Li, Yueming Zhou, Xiaomeng Ma, Min Li, Cheng Huang, Peixiang LuUsing the classical ensemble model, we systematically investigate strongfield nonsequential double ionization (NSDI) of Mg by intense elliptically polarized laser pulses with different wavelengths. Different from the noble atoms, NSDI occurs for Mg driven by elliptically and circularly polarized laser fields. Our results show that in elliptically and circularly polarized laser fields, the NSDI yield is sharply suppressed as the wavelength increases. Interestingly, the correlated behavior in the electron momentum spectra depends sensitively on the wavelengths. The corresponding electron dynamics is revealed by back tracing the classical trajectory.

A study on the plasticity of sodalime silica glass via molecular dynamics simulations J. Chem. Phys. (IF 2.965) Pub Date : 20171102
Shingo Urata, Yosuke SatoMolecular dynamics (MD) simulations were applied to construct a plasticity model, which enables one to simulate deformations of sodalime silica glass (SLSG) by using continuum methods. To model the plasticity, stress induced by uniaxial and a variety of biaxial deformations was measured by MD simulations. We found that the surfaces of yield and maximum stresses, which are evaluated from the equivalent stressstrain curves, are reasonably represented by the MohrCoulomb ellipsoid. Comparing a finite element model using the constructed plasticity model to a large scale atomistic model on a nanoindentation simulation of SLSG reveals that the empirical method is accurate enough to evaluate the SLSG mechanical responses. Furthermore, the effect of ionexchange on the SLSG plasticity was examined by using MD simulations. As a result, it was demonstrated that the effects of the initial compressive stress on the yield and maximum stresses are anisotropic contrary to our expectations.

The accretion of the new ice layer on the surface of hexagonal ice crystal and the influence of the local electric field on this process J. Chem. Phys. (IF 2.965) Pub Date : 20171107
Joanna Grabowska, Anna Kuffel, Jan ZielkiewiczThe process of creation of a new layer of ice on the basal plane and on the prism plane of a hexagonal ice crystal is analyzed. It is demonstrated that the ordering of water molecules in the already existing crystal affects the freezing. On the basal plane, when the orientations of water molecules in the ice block are random, the arrangement of the new layer in a cubic manner is observed more frequently—approximately 1.7 times more often than in a hexagonal manner. When the water molecules in the ice block are more ordered, it results in the predominance of the oxygen atoms or the hydrogen atoms on the most outer part of the surface of the ice block. In this case, the hexagonal structure is formed more frequently when the supercooling of water exceeds 10 K. This phenomenon is explained by the influence of the oriented electric field, present as a consequence of the ordering of the dipoles of water molecules in the ice block. This field modifies the structure of solvation water (i.e., the layer of water in the immediate vicinity of the ice surface). We showed that the structure of solvation water predetermines the kind of the newly created layer of ice. This effect is temperaturedependent: when the temperature draws nearer to the melting point, the cubic structure becomes the prevailing form. The temperature at which the cubic and the hexagonal structures are formed with the same probabilities is equal to about 260 K. In the case of the prism plane, the new layer that is formed is always the hexagonal one, which is independent of the arrangement of water molecules in the ice block and is in agreement with previous literature data. For the basal plane, as well as for the prism plane, no evident dependence on the ordering of water molecules that constitute the ice block on the rate of crystallization can be observed.

Going beyond the standard line tension: Sizedependent contact angles of water nanodroplets J. Chem. Phys. (IF 2.965) Pub Date : 20171101
Matej KandučThe dependence of the contact angle on the size of a nanoscopic droplet residing on a flat substrate is traditionally ascribed solely to line tension. Other contributions, stemming from the droplet geometry dependence of the surface tension and line tension, are typically ignored. Here, we perform molecular dynamics simulations of water droplets of cylindrical morphology on surfaces of a wide range of polarities. In the cylindrical geometry, where the line tension is not operative directly, we find significant contact angle dependence on the droplet size. The effect is most pronounced on hydrophilic surfaces, with the contact angle increase of up to 1 0 ° with a decreasing droplet size. On hydrophobic surfaces, the trend is reversed and considerably weaker. Our analysis suggests that these effects can be attributed to the Tolman correction due to the curved water–vapor interface and to a generalized line tension that possesses a contact angle dependence. The latter is operative also in the cylindrical geometry and yields a comparable contribution to the contact angle as the line tension itself in case of spherical droplets.

Demixing of active particles in the presence of external fields J. Chem. Phys. (IF 2.965) Pub Date : 20171102
Sunita Kumari, André S. Nunes, Nuno A. M. Araújo, Margarida M. Telo da GamaSelfpropelled active particles are inherently out of equilibrium as they collect energy from their surroundings and transform it into directed motion. A recent theoretical study suggests that binary mixtures of active particles with distinct effective diffusion coefficients exhibit dynamical demixing when their diffusion coefficients differ by more than one order of magnitude. Here, we show that this difference may be reduced drastically in the presence of external fields even when the response to the field is the same for both species. We investigate this demixing as a function of the ratio of the diffusion coefficients and discuss the implications of the results for active systems.

MnNiO_{3} revisited with modern theoretical and experimental methods J. Chem. Phys. (IF 2.965) Pub Date : 20171103
Allison L. Dzubak, Chandrima Mitra, Michael Chance, Stephen Kuhn, Gerald E. JellisonJr., Athena S. Sefat, Jaron T. Krogel, Fernando A. ReboredoMnNiO3 is a strongly correlated transition metal oxide that has recently been investigated theoretically for its potential application as an oxygenevolution photocatalyst. However, there is no experimental report on critical quantities such as the band gap or bulk modulus. Recent theoretical predictions with standard functionals such as LDA+U and HSE show large discrepancies in the band gaps (about 1.23 eV), depending on the nature of the functional used. Hence there is clearly a need for an accurate quantitative prediction of the band gap to gauge its utility as a photocatalyst. In this work, we present a diffusion quantum Monte Carlo study of the bulk properties of MnNiO3 and revisit the synthesis and experimental properties of the compound. We predict quasiparticle band gaps of 2.0(5) eV and 3.8(6) eV for the majority and minority spin channels, respectively, and an equilibrium volume of 92.8 Å3, which compares well to the experimental value of 94.4 Å3. A bulk modulus of 217 GPa is predicted for MnNiO3. We rationalize the difficulty for the formation of ordered ilmenitetype structure with specific sites for Ni and Mn to be potentially due to the formation of antisite defects that form during synthesis, which ultimately affects the physical properties of MnNiO3.

Retention of contaminants Cd and Hg adsorbed and intercalated in aluminosilicate clays: A first principles study J. Chem. Phys. (IF 2.965) Pub Date : 20171103
F. D. Crasto de Lima, R. H. Miwa, Caetano R. MirandaLayered clay materials have been used to incorporate transition metal (TM) contaminants. Based on firstprinciples calculations, we have examined the energetic stability and the electronic properties due to the incorporation of Cd and Hg in layered clay materials, kaolinite (KAO) and pyrophyllite (PYR). The TM can be (i) adsorbed on the clay surface as well as (ii) intercalated between the clay layers. For the intercalated case, the contaminant incorporation rate can be optimized by controlling the interlayer spacing of the clay, namely, pillared clays. Our total energy results reveal that the incorporation of the TMs can be maximized through a suitable tuning of vertical distance between the clay layers. Based on the calculated TM/clay binding energies and the Langmuir absorption model, we estimate the concentrations of the TMs. Further kinetic properties have been examined by calculating the activation energies, where we found energy barriers of ∼20 and ∼130 meV for adsorbed and intercalated cases, respectively. The adsorption and intercalation of ionized TM adatoms were also considered within the deprotonated KAO surface. This also leads to an optimal interlayer distance which maximizes the TM incorporation rate. By mapping the total charge transfers at the TM/clay interface, we identify a net electronic charge transfer from the TM adatoms to the topmost clay surface layer. The effect of such a charge transfer on the electronic structure of the clay (host) has been examined through a set of Xray absorption near edge structure (XANES) simulations, characterizing the changes of the XANES spectra upon the presence of the contaminants. Finally, for the pillared clays, we quantify the Cd and Hg Kedge energy shifts of the TMs as a function of the interlayer distance between the clay layers and the Al Kedge spectra for the pristine and pillared clays.

Chargeregularized swelling kinetics of polyelectrolyte gels: Elasticity and diffusion J. Chem. Phys. (IF 2.965) Pub Date : 20171101
Swati Sen, Arindam KundagramiWe apply a recently developed method [S. Sen and A. Kundagrami, J. Chem. Phys. 143, 224904 (2015)], using a phenomenological expression of osmotic stress, as a function of polymer and charge densities, hydrophobicity, and network elasticity for the swelling of spherical polyelectrolyte (PE) gels with fixed and variable charges in a saltfree solvent. This expression of stress is used in the equation of motion of swelling kinetics of spherical PE gels to numerically calculate the spatial profiles for the polymer and free ion densities at different time steps and the time evolution of the size of the gel. We compare the profiles of the same variables obtained from the classical linear theory of elasticity and quantitatively estimate the bulk modulus of the PE gel. Further, we obtain an analytical expression of the elastic modulus from the linearized expression of stress (in the small deformation limit). We find that the estimated bulk modulus of the PE gel decreases with the increase of its effective charge for a fixed degree of deformation during swelling. Finally, we match the gelfront locations with the experimental data, taken from the measurements of charged reversible additionfragmentation chain transfer gels to show an increase in gelsize with charge and also match the same for PNIPAM (uncharged) and imidazoliumbased (charged) minigels, which specifically confirms the decrease of the gel modulus value with the increase of the charge. The agreement between experimental and theoretical results confirms general diffusive behaviour for swelling of PE gels with a decreasing bulk modulus with increasing degree of ionization (charge). The new formalism captures large deformations as well with a significant variation of charge content of the gel. It is found that PE gels with large deformation but same initial size swell faster with a higher charge.

A lattice model for the impact of volume fraction fluctuations upon percolation by cylinders J. Chem. Phys. (IF 2.965) Pub Date : 20171103
Avik P. Chatterjee, Claudio GrimaldiA lattice model for continuum percolation by cylindrical rods is generalized to account for inhomogeneities in the volume fraction that are indicative of particle clustering or aggregation. The percolation threshold is evaluated from a formalism that uses two different categories of occupied sites (denoting particles) with different occupation probabilities that represent large and small local volume fractions. Our modeling framework enables independent variations in (i) the strength of the correlation that adjacent particles experience high (or low) effective volume fractions, (ii) the disparity between the macroscopically averaged volume fraction and (say) the volume fraction characterizing the regions with high effective particle concentrations, and (iii) the overall proportion of particles that are located in regions with either high or low volume fraction. Calculations performed for monodisperse cylinders show that enhancement in each of the above factors leads to reduction in the macroscopically averaged volume fraction at the percolation threshold.

Nonconstant link tension coefficient in the tumblingsnake model subjected to simple shear J. Chem. Phys. (IF 2.965) Pub Date : 20171106
Pavlos S. Stephanou, Martin KrögerThe authors of the present study have recently presented evidence that the tumblingsnake model for polymeric systems has the necessary capacity to predict the appearance of pronounced undershoots in the timedependent shear viscosity as well as an absence of equally pronounced undershoots in the transient two normal stress coefficients. The undershoots were found to appear due to the tumbling behavior of the director u when a rotational Brownian diffusion term is considered within the equation of motion of polymer segments, and a theoretical basis concerning the use of a link tension coefficient given through the nematic order parameter had been provided. The current work elaborates on the quantitative predictions of the tumblingsnake model to demonstrate its capacity to predict undershoots in the timedependent shear viscosity. These predictions are shown to compare favorably with experimental rheological data for both polymer melts and solutions, help us to clarify the microscopic origin of the observed phenomena, and demonstrate in detail why a constant link tension coefficient has to be abandoned.

Note: MSM lag time cannot be used for variational model selection J. Chem. Phys. (IF 2.965) Pub Date : 20171103
Brooke E. Husic, Vijay S. PandeThe variational principle for conformational dynamics has enabled the systematic construction of Markov state models through the optimization of hyperparameters by approximating the transfer operator. In this note, we discuss why the lag time of the operator being approximated must be held constant in the variational approach.

Preface: Special Topic: From Quantum Mechanics to Force Fields J. Chem. Phys. (IF 2.965) Pub Date : 20171024
JeanPhilip Piquemal, Kenneth D. JordanThis Special Topic issue entitled “From Quantum Mechanics to Force Fields” is dedicated to the ongoing efforts of the theoretical chemistry community to develop a new generation of accurate force fields based on data from highlevel electronic structure calculations and to develop faster electronic structure methods for testing and designing force fields as well as for carrying out simulations. This issue includes a collection of 35 original research articles that illustrate recent theoretical advances in the field. It provides a timely snapshot of recent developments in the generation of approaches to enable more accurate molecular simulations of processes important in chemistry, physics, biophysics, and materials science.

Cheap but accurate calculation of chemical reaction rate constants from ab initio data, via systemspecific, blackbox force fields J. Chem. Phys. (IF 2.965) Pub Date : 20170414
Julien Steffen, Bernd HartkeBuilding on the recently published quantummechanically derived force field (QMDFF) and its empirical valence bond extension, EVBQMDFF, it is now possible to generate a reliable potential energy surface for any given elementary reaction step in an essentially black box manner. This requires a limited and predefined set of reference data near the reaction path and generates an accurate approximation of the reference potential energy surface, on and off the reaction path. This intermediate representation can be used to generate reaction rate data, with far better accuracy and reliability than with traditional approaches based on transition state theory (TST) or variational extensions thereof (VTST), even if those include sophisticated tunneling corrections. However, the additional expense at the reference level remains very modest. We demonstrate all this for three arbitrarily chosen example reactions.

Mapping the Drude polarizable force field onto a multipole and induced dipole model J. Chem. Phys. (IF 2.965) Pub Date : 20170605
Jing Huang, Andrew C. Simmonett, Frank C. PickardIV, Alexander D. MacKerellJr., Bernard R. BrooksThe induced dipole and the classical Drude oscillator represent two major approaches for the explicit inclusion of electronic polarizability into force fieldbased molecular modeling and simulations. In this work, we explore the equivalency of these two models by comparing condensed phase properties computed using the Drude force field and a multipole and induced dipole (MPID) model. Presented is an approach to map the electrostatic model optimized in the context of the Drude force field onto the MPID model. Condensed phase simulations on water and 15 small model compounds show that without any reparametrization, the MPID model yields properties similar to the Drude force field with both models yielding satisfactory reproduction of a range of experimental values and quantum mechanical data. Our results illustrate that the Drude oscillator model and the point induced dipole model are different representations of essentially the same physical model. However, results indicate the presence of small differences between the use of atomic multipoles and offcenter charge sites. Additionally, results on the use of dispersion particle mesh Ewald further support its utility for treating longrange Lennard Jones dispersion contributions in the context of polarizable force fields. The main motivation in demonstrating the transferability of parameters between the Drude and MPID models is that the more than 15 years of development of the Drude polarizable force field can now be used with MPID formalism without the need for dualthermostat integrators nor selfconsistent iterations. This opens up a wide range of new methodological opportunities for polarizable models.

A comparison of sodium and hydrogen halides at the airwater interface J. Chem. Phys. (IF 2.965) Pub Date : 20170605
Collin D. WickNew molecular models, parameterized to ab initio calculations, were developed to describe HBr and HI at the airwater interface. These were used to compare how the airwater interface influenced dissociation of NaX and HX, with X being Cl, Br, or I, and also their propensity for the interface. The polarizable multistate empirical valence bond method, which explicitly describes proton sharing, was used to model HX. Results showed that the airwater interface suppressed HX dissociation from a contact ion pair to a solvent separated to a greater degree than NaX dissociation. Furthermore, HX had a greater propensity for the interface than NaX, which was a consequence of the hydronium ion having a greatest interfacial activity of all species studied. As a consequence of this, the average configuration of dissociated HX, while in both contact ion and solvent separated ion pairs near the airwater interface, is with the dissociated hydrogen oriented more towards the air than the X atom.

Intermolecular interactions in the condensed phase: Evaluation of semiempirical quantum mechanical methods J. Chem. Phys. (IF 2.965) Pub Date : 20170613
Anders S. Christensen, Jimmy C. Kromann, Jan H. Jensen, Qiang CuiTo facilitate further development of approximate quantum mechanical methods for condensed phase applications, we present a new benchmark dataset of intermolecular interaction energies in the solution phase for a set of 15 dimers, each containing one charged monomer. The reference interaction energy in solution is computed via a thermodynamic cycle that integrates dimer binding energy in the gas phase at the coupled cluster level and solutesolvent interaction with density functional theory; the estimated uncertainty of such calculated interaction energy is ±1.5 kcal/mol. The dataset is used to benchmark the performance of a set of semiempirical quantum mechanical (SQM) methods that include DFTB3D3, DFTB3/CPED3, OM2D3, PM6D3, PM6D3H+, and PM7 as well as the HF3c method. We find that while all tested SQM methods tend to underestimate binding energies in the gas phase with a rootmeansquared error (RMSE) of 25 kcal/mol, they overestimate binding energies in the solution phase with an RMSE of 34 kcal/mol, with the exception of DFTB3/CPED3 and OM2D3, for which the systematic deviation is less pronounced. In addition, we find that HF3c systematically overestimates binding energies in both gas and solution phases. As most approximate QM methods are parametrized and evaluated using data measured or calculated in the gas phase, the dataset represents an important first step toward calibrating QM based methods for application in the condensed phase where polarization and exchange repulsion need to be treated in a balanced fashion.

Explicit treatment of hydrogen bonds in the universal force field: Validation and application for metalorganic frameworks, hydrates, and hostguest complexes J. Chem. Phys. (IF 2.965) Pub Date : 20170615
Damien E. Coupry, Matthew A. Addicoat, Thomas HeineA straightforward means to include explicit hydrogen bonds within the Universal Force Field (UFF) is presented. Instead of treating hydrogen bonds as nonbonded interaction subjected to electrostatic and LennardJones potentials, we introduce an explicit bond with a negligible bond order, thus maintaining the structural integrity of the Hbonded complexes and avoiding the necessity to assign arbitrary charges to the system. The explicit hydrogen bond changes the coordination number of the acceptor site and the approach is thus most suitable for systems with undercoordinated atoms, such as many metalorganic frameworks; however, it also shows an excellent performance for other systems involving a hydrogenbonded framework. In particular, it is an excellent means for creating starting structures for molecular dynamics and for investigations employing more sophisticated methods. The approach is validated for the hydrogen bonded complexes in the S22 dataset and then employed for a set of metalorganic frameworks from the ComputationReady Experimental database and several hydrogen bonded crystals including water ice and clathrates. We show that the direct inclusion of hydrogen bonds reduces the maximum error in predicted cell parameters from 66% to only 14%, and the mean unsigned error is similarly reduced from 14% to only 4%. We posit that with the inclusion of hydrogen bonding, the solventmediated breathing of frameworks such as MIL53 is now accessible to rapid UFF calculations, which will further the aim of rapid computational scanning of metalorganic frameworks while providing better starting points for electronic structure calculations.

Interpolation of intermolecular potentials using Gaussian processes J. Chem. Phys. (IF 2.965) Pub Date : 20170626
Elena Uteva, Richard S. Graham, Richard D. Wilkinson, Richard J. WheatleyA procedure is proposed to produce intermolecular potential energy surfaces from limited data. The procedure involves generation of geometrical configurations using a Latin hypercube design, with a maximin criterion, based on inverse internuclear distances. Gaussian processes are used to interpolate the data, using overspecified inverse molecular distances as covariates, greatly improving the interpolation. Symmetric covariance functions are specified so that the interpolation surface obeys all relevant symmetries, reducing prediction errors. The interpolation scheme can be applied to many important molecular interactions with trivial modifications. Results are presented for three systems involving CO2, a system with a deep energy minimum (HF−HF), and a system with 48 symmetries (CH4−N2). In each case, the procedure accurately predicts an independent test set. Training this method with highprecision ab initio evaluations of the CO2−CO interaction enables a parameterfree, firstprinciples prediction of the CO2−CO cross virial coefficient that agrees very well with experiments.

On the development of polarizable and LennardJones force fields to study hydration structure and dynamics of actinide(III) ions based on effective ionic radii J. Chem. Phys. (IF 2.965) Pub Date : 20170705
Riccardo Spezia, Valentina Migliorati, Paola D’AngeloIn this contribution, we show how it is possible to develop polarizable and nonpolarizable force fields to study hydration properties of a whole chemical series based on atomic properties such as ionic radii. In particular, we have addressed the actinide(III) ion series, from U3+ to Cf3+, for which Xray absorption data and effective ionic radii are available. A polarizable force field has been reparameterized improving the original one [M. Duvail et al., J. Chem. Phys. 135, 044503 (2011)] which was based on solid state ionic radii. The new force field does not depend on solid state properties but directly on the liquid phase ones, and it can be used to study these ions in liquid water without any ambiguity. Furthermore, we have shown that it is possible to parameterize also a nonpolarizable potential using standard LennardJones and Coulombic forces, which can be transferred to other systems in condensed phase. The structural and dynamical properties of these two force fields are compared to each other and with data available in the literature, providing a good agreement. Moreover, we show the comparison with experimental Xray absorption data that are very well reproduced by both force fields.

A general intermolecular force field based on tightbinding quantum chemical calculations J. Chem. Phys. (IF 2.965) Pub Date : 20170711
Stefan Grimme, Christoph Bannwarth, Eike Caldeweyher, Jana Pisarek, Andreas HansenA blackbox type procedure is presented for the generation of a moleculespecific, intermolecular potential energy function. The method uses quantum chemical (QC) information from our recently published extended tightbinding semiempirical scheme (GFNxTB) and can treat noncovalently bound complexes and aggregates with almost arbitrary chemical structure. The necessary QC information consists of the equilibrium structure, Mulliken atomic charges, charge centers of localized molecular orbitals, and also of frontier orbitals and orbital energies. The molecular pair potential includes model density dependent Pauli repulsion, penetration, as well as point charge electrostatics, the newly developed D4 dispersion energy model, Drude oscillators for polarization, and a chargetransfer term. Only one elementspecific and about 20 global empirical parameters are needed to cover systems with nuclear charges up to radon (Z = 86). The method is tested for standard small molecule interaction energy benchmark sets where it provides accurate intermolecular energies and equilibrium distances. Examples for structures with a few hundred atoms including charged systems demonstrate the versatility of the approach. The method is implemented in a standalone computer code which enables rigidbody, global minimum energy searches for molecular aggregation or alignment.

Li^{+} solvation and kinetics of Li^{+}–BF_{4}^{−}/PF_{6}^{−} ion pairs in ethylene carbonate. A molecular dynamics study with classical rate theories J. Chem. Phys. (IF 2.965) Pub Date : 20170719
TsunMei Chang, Liem X. DangUsing our polarizable forcefield models and employing classical rate theories of chemical reactions, we examine the ethylene carbonate (EC) exchange process between the first and second solvation shells around Li+ and the dissociation kinetics of ion pairs Li+–[BF4] and Li+–[PF6] in this solvent. We calculate the exchange rates using transition state theory and correct them with transmission coefficients computed by the reactive flux, Impey, Madden, and McDonald approaches, and GroteHynes theory. We found that the residence times of EC around Li+ ions varied from 60 to 450 ps, depending on the correction method used. We found that the relaxation times changed significantly from Li+–[BF4] to Li+–[PF6] ion pairs in EC. Our results also show that, in addition to affecting the free energy of dissociation in EC, the anion type also significantly influences the dissociation kinetics of ion pairing.

Structure and polarization near the Li^{+} ion in ethylene and propylene carbonates J. Chem. Phys. (IF 2.965) Pub Date : 20170719
Travis P. Pollard, Thomas L. BeckResearch on fundamental interactions in Liion batteries is accelerating due to the importance of developing batteries with enhanced energy and power densities while maintaining safety. Improving electrode materials and controlling the formation of the solid electrolyte interphase during the first battery charge have been the main focus areas for research. Ionsolvent interactions in the electrolyte are also of great importance in tuning solvation and transport properties, however. Here we present ab initio density functional theory simulations of a Li+ ion in ethylene and propylene carbonates. The aim is to obtain a detailed analysis of local solvation structure and solvent polarization near the ion and in the bulk. The results indicate the significance of molecular polarization for developing accurate solvation models. The simulations illustrate the substantial differences between ion solvation in water and in organic materials.

Structural study of Na_{2}O–B_{2}O_{3}–SiO_{2} glasses from molecular simulations using a polarizable force field J. Chem. Phys. (IF 2.965) Pub Date : 20170719
Fabien Pacaud, JeanMarc Delaye, Thibault Charpentier, Laurent Cormier, Mathieu SalanneSodium borosilicate glasses Na2O–B2O3–SiO2 (NBS) are complex systems from a structural point of view. Three main building units are present: tetrahedral SiO4 and BO4 (BIV) and triangular BO3 (BIII). One of the salient features of these compounds is the change of the BIII/BIV ratio with the alkali concentration, which is very difficult to capture in force fieldsbased molecular dynamics simulations. In this work, we develop a polarizable force field that is able to reproduce the boron coordination and more generally the structure of several NBS systems in the glass and in the melt. The parameters of the potential are fitted from density functional theory calculations only, in contrast with the existing empirical potentials for NBS systems. This ensures a strong improvement on the transferability of the parameters from one composition to another. Using this new force field, the structure of NBS systems is validated against neutron diffraction and nuclear magnetic resonance experiments. A special focus is given to the distribution of BIII/BIV with respect to the composition and the temperature.

Minimal distributed charges: Multipolar quality at the cost of point charge electrostatics J. Chem. Phys. (IF 2.965) Pub Date : 20170721
Oliver T. Unke, Mike Devereux, Markus MeuwlyMost empirical force fields use atomcentered point charges (PCs) to represent the electrostatic potential (ESP) around molecules. While such PC models are computationally efficient, they are unable to capture anisotropic electronic features, such as σ holes or lone pairs. These features are better described using atomic multipole (MTP) moments, which significantly improve the quality of the resulting ESP. However, the improvement comes at the expense of a considerably increased computational complexity and cost for calculating the interaction energies and forces. In the present work, a novel minimal distributed charge model (MDCM) based on offcentered point charges is presented and the quality of the resulting ESP is compared to the performance of MTPs and atomcentered PC models for several test molecules. All three models are fitted using the same algorithm based on differential evolution, which is available as a Fortran90 program from the authors upon request. We show that the MDCM is capable of approximating the reference ab initio ESP with an accuracy as good as, or better than, MTPs without the need for computationally expensive higher order multipoles. Further it is demonstrated that the MDCM is numerically stable in molecular dynamics simulations and is able to reproduce electrostatic interaction energies and thermodynamic quantities with the same accuracy as MTPs at reduced computational cost.

Combining configurational energies and forces for molecular force field optimization J. Chem. Phys. (IF 2.965) Pub Date : 20170721
Lukas Vlcek, Weiwei Sun, Paul R. C. KentWhile quantum chemical simulations have been increasingly used as an invaluable source of information for atomistic model development, the high computational expenses typically associated with these techniques often limit thorough sampling of the systems of interest. It is therefore of great practical importance to use all available information as efficiently as possible, and in a way that allows for consistent addition of constraints that may be provided by macroscopic experiments. Here we propose a simple approach that combines information from configurational energies and forces generated in a molecular dynamics simulation to increase the effective number of samples. Subsequently, this information is used to optimize a molecular force field by minimizing the statistical distance similarity metric. We illustrate the methodology on an example of a trajectory of configurations generated in equilibrium molecular dynamics simulations of argon and water and compare the results with those based on the force matching method.

Twocomponent, ab initio potential energy surface for CO_{2}—H_{2}O, extension to the hydrate clathrate, CO_{2}@(H_{2}O)_{20}, and VSCF/VCI vibrational analyses of both J. Chem. Phys. (IF 2.965) Pub Date : 20170721
Qingfeng (Kee) Wang, Joel M. BowmanWe report an ab initio, fulldimensional, potential energy surface (PES) for CO2—H2O, in which twobody interaction energies are fit using a basis of permutationally invariant polynomials and combined with accurate potentials for the noninteracting monomers. This approach which we have termed “plug and play” is extended here to improve the precision of the 2body fit in the long range. This is done by combining two separate fits. One is a fit to 47 593 2body energies in the region of strong interaction and approaching the long range, and the second one is a fit to 6244 2body energies in the long range. The two fits have a region of overlap which permits a smooth switch from one to the other. All energies are obtained at the CCSD(T)F12b/augccpVTZ level of theory. Properties of the full PES, i.e., stationary points, harmonic frequencies of the global minimum, etc., are shown to be in excellent agreement with direct CCSD(T)F12b/augccpVTZ results. Diffusion Monte Carlo calculations of the dimer zeropoint energy (ZPE) are performed, and a dissociation energy, D0, of 787 cm−1 is obtained using that ZPE, De, and the rigorous ZPEs of the monomers. Using a benchmark De, D0 is 758 cm−1. Vibrational selfconsistent field (VSCF)/virtual state configuration interaction (VCI) MULTIMODE calculations of intramolecular fundamentals are reported and are in good agreement with available experimental results. Finally, the full dimer PES is combined with an existing ab initio water potential to develop a potential for the CO2 hydrate clathrate CO2(H2O)20(512 water cage). A full normalmode analysis of this hydrate clathrate is reported as are localmonomer VSCF/VCI calculations of the fundamentals of CO2.

Toward chemical accuracy in the description of ion–water interactions through manybody representations. Alkaliwater dimer potential energy surfaces J. Chem. Phys. (IF 2.965) Pub Date : 20170724
Marc Riera, Narbe Mardirossian, Pushp Bajaj, Andreas W. Götz, Francesco PaesaniThis study presents the extension of the MBnrg (ManyBody energy) theoretical/computational framework of transferable potential energy functions (PEFs) for molecular simulations of alkali metal ionwater systems. The MBnrg PEFs are built upon the manybody expansion of the total energy and include the explicit treatment of onebody, twobody, and threebody interactions, with all higherorder contributions described by classical induction. This study focuses on the MBnrg twobody terms describing the fulldimensional potential energy surfaces of the M+(H2O) dimers, where M+ = Li+, Na+, K+, Rb+, and Cs+. The MBnrg PEFs are derived entirely from “first principles” calculations carried out at the explicitly correlated coupledcluster level including single, double, and perturbative triple excitations [CCSD(T)F12b] for Li+ and Na+ and at the CCSD(T) level for K+, Rb+, and Cs+. The accuracy of the MBnrg PEFs is systematically assessed through an extensive analysis of interaction energies, structures, and harmonic frequencies for all five M+(H2O) dimers. In all cases, the MBnrg PEFs are shown to be superior to both polarizable force fields and ab initio models based on density functional theory. As previously demonstrated for halidewater dimers, the MBnrg PEFs achieve higher accuracy by correctly describing shortrange quantummechanical effects associated with electron density overlap as well as longrange electrostatic manybody interactions.

Electrostatic solvation free energies of charged hard spheres using molecular dynamics with density functional theory interactions J. Chem. Phys. (IF 2.965) Pub Date : 20170726
Timothy T. Duignan, Marcel D. Baer, Gregory K. Schenter, Chistopher J. MundyDetermining the solvation free energies of single ions in water is one of the most fundamental problems in physical chemistry and yet many unresolved questions remain. In particular, the ability to decompose the solvation free energy into simple and intuitive contributions will have important implications for models of electrolyte solution. Here, we provide definitions of the various types of single ion solvation free energies based on different simulation protocols. We calculate solvation free energies of charged hard spheres using density functional theory interaction potentials with molecular dynamics simulation and isolate the effects of charge and cavitation, comparing to the Born (linear response) model. We show that using uncorrected Ewald summation leads to unphysical values for the single ion solvation free energy and that charging free energies for cations are approximately linear as a function of charge but that there is a small nonlinearity for small anions. The charge hydration asymmetry for hard spheres, determined with quantum mechanics, is much larger than for the analogous real ions. This suggests that real ions, particularly anions, are significantly more complex than simple charged hard spheres, a commonly employed representation.

Implementation of analytical gradients and of a mixed real and momentum space DVR method for excess electron systems described by a selfconsistent polarization model J. Chem. Phys. (IF 2.965) Pub Date : 20170731
Tae Hoon Choi, Tijo Vazhappilly, Kenneth D. JordanThis work presents two extensions of our selfconsistent polarization model for treating nonvalence excess electron systems. The first extension is the implementation of analytical gradients, and the second extension is the implementation of a mixed real space plus momentum space approach combined with fast Fourier transforms to reduce the computational time compared to a purely real space discrete variable representation approach. The performance of the new algorithms is assessed in calculations of the excess electron states of various size water clusters and of the nonvalence correlationbound anion of the C240 fullerene.

Geometrydependent atomic multipole models for the water molecule J. Chem. Phys. (IF 2.965) Pub Date : 20170731
O. Loboda, C. MillotModels of atomic electric multipoles for the water molecule have been optimized in order to reproduce the electric potential around the molecule computed by ab initio calculations at the coupled cluster level of theory with up to noniterative triple excitations in an augmented triplezeta quality basis set. Different models of increasing complexity, from atomic charges up to models containing atomic charges, dipoles, and quadrupoles, have been obtained. The geometry dependence of these atomic multipole models has been investigated by changing bond lengths and HOH angle to generate 125 molecular structures (reduced to 75 symmetryunique ones). For several models, the atomic multipole components have been fitted as a function of the geometry by a Taylor series of fourth order in monomer coordinate displacements.

Development of reactive force fields using ab initio molecular dynamics simulation minimally biased to experimental data J. Chem. Phys. (IF 2.965) Pub Date : 20170801
Chen Chen, Christopher Arntsen, Gregory A. VothIncorporation of quantum mechanical electronic structure data is necessary to properly capture the physics of many chemical processes. Proton hopping in water, which involves rearrangement of chemical and hydrogen bonds, is one such example of an inherently quantum mechanical process. Standard ab initio molecular dynamics (AIMD) methods, however, do not yet accurately predict the structure of water and are therefore less than optimal for developing force fields. We have instead utilized a recently developed method which minimally biases AIMD simulations to match limited experimental data to develop novel multiscale reactive molecular dynamics (MSRMD) force fields by using relative entropy minimization. In this paper, we present two new MSRMD models using such a parameterization: one which employs water with harmonic internal vibrations and another which uses anharmonic water. We show that the newly developed MSRMD models very closely reproduce the solvation structure of the hydrated excess proton in the target AIMD data. We also find that the use of anharmonic water increases proton hopping, thereby increasing the proton diffusion constant.

Organic ion association in aqueous phase and ab initiobased force fields: The case of carboxylate/ammonium salts J. Chem. Phys. (IF 2.965) Pub Date : 20170818
Céline Houriez, Valérie Vallet, Florent Réal, Michael MeotNer (Mautner), Michel MasellaWe performed molecular dynamics simulations of carboxylate/methylated ammonium ion pairs solvated in bulk water and of carboxylate/methylated ammonium salt solutions at ambient conditions using an ab initiobased polarizable force field whose parameters are assigned to reproduce only high end quantum computations, at the MøllerPlesset secondorder perturbation theory/complete basis set limit level, regarding single ions and ion pairs as isolated and microhydrated in gas phase. Our results agree with the available experimental results regarding carboxylate/ammonium salt solutions. For instance, our force field approach predicts the percentage of acetate associated with ammonium ions in CH3COO−/CH3NH3+ solutions at the 0.2–0.8M concentration scale to range from 14% to 35%, in line with the estimates computed from the experimental ion association constant in liquid water. Moreover our simulations predict the number of water molecules released from the ion first hydration shell to the bulk upon ion association to be about 2.0 ± 0.6 molecules for acetate/protonated amine ion pairs, 3.1 ± 1.5 molecules for the HCOO−/NH4+ pair and 3.3 ± 1.2 molecules for the CH3COO−/(CH3)4N+ pair. For protonated aminebased ion pairs, these values are in line with experiment for alkali/halide pairs solvated in bulk water. All these results demonstrate the promising feature of ab initiobased force fields, i.e., their capacity in accurately modeling chemical systems that cannot be readily investigated using available experimental techniques.

Assessing manybody contributions to intermolecular interactions of the AMOEBA force field using energy decomposition analysis of electronic structure calculations J. Chem. Phys. (IF 2.965) Pub Date : 20170824
Omar Demerdash, Yuezhi Mao, Tianyi Liu, Martin HeadGordon, Teresa HeadGordonIn this work, we evaluate the accuracy of the classical AMOEBA model for representing manybody interactions, such as polarization, charge transfer, and Pauli repulsion and dispersion, through comparison against an energy decomposition method based on absolutely localized molecular orbitals (ALMOEDA) for the water trimer and a variety of ionwater systems. When the 2 and 3body contributions according to the manybody expansion are analyzed for the ionwater trimer systems examined here, the 3body contributions to Pauli repulsion and dispersion are found to be negligible under ALMOEDA, thereby supporting the validity of the pairwiseadditive approximation in AMOEBA’s 147 van der Waals term. However AMOEBA shows imperfect cancellation of errors for the missing effects of charge transfer and incorrectness in the distance dependence for polarization when compared with the corresponding ALMOEDA terms. We trace the larger 2body followed by 3body polarization errors to the Thole damping scheme used in AMOEBA, and although the width parameter in Thole damping can be changed to improve agreement with the ALMOEDA polarization for points about equilibrium, the correct profile of polarization as a function of intermolecular distance cannot be reproduced. The results suggest that there is a need for reexamining the damping and polarization model used in the AMOEBA force field and provide further insights into the formulations of polarizable force fields in general.

From dimers to the solidstate: Distributed intermolecular forcefields for pyridine J. Chem. Phys. (IF 2.965) Pub Date : 20170828
Alexander A. Aina, Alston J. Misquitta, Sarah L. PriceAn anisotropic atomatom forcefield for pyridine, using distributed atomic multipoles, polarizabilities, and dispersion coefficients and an anisotropic atomatom repulsion model derived from symmetryadapted perturbation theory (density functional theory) dimer calculations, is used to model pyridine crystal structures. Here we show that this distributed intermolecular forcefield (DIFF) models the experimental crystal structures as accurately as modelling all but the electrostatic term with an isotropic repulsiondispersion potential that has been fitted to experimental crystal structures. In both cases, the differences are comparable to the changes in the crystal structure with temperature, pressure, or neglect of zeropoint vibrational effects. A crystal structure prediction study has been carried out, and the observed polymorphs contrasted with hypothetical thermodynamically competitive crystal structures. The DIFF model was able to identify the structure of an unreported high pressure phase of pyridine, unlike the empirically fitted potential. The DIFF model approach therefore provides a model of the underlying pair potential energy surface that we have transferred to the crystalline phase with a considerable degree of success, though the treatment of the manybody terms needs improvement and the pair potential is slightly overbinding. Furthermore, this study of a system that exhibits isotopic polymorphism highlights that the use of an empirical potential has partially absorbed temperature and zeropoint motion effects as well as the intermolecular forces not explicitly represented in the functional form. This study therefore highlights the complexity in modelling crystallization phenomena from a realistic pair potential energy surface.

Analytical gradients for tensor hypercontracted MP2 and SOSMP2 on graphical processing units J. Chem. Phys. (IF 2.965) Pub Date : 20170829
Chenchen Song, Todd J. MartínezAnalytic energy gradients for tensor hypercontraction (THC) are derived and implemented for secondorder MøllerPlesset perturbation theory (MP2), with and without the scaledoppositespin (SOS)MP2 approximation. By exploiting the THC factorization, the formal scaling of MP2 and SOSMP2 gradient calculations with respect to system size is reduced to quartic and cubic, respectively. An efficient implementation has been developed that utilizes both graphics processing units and sparse tensor techniques exploiting spatial sparsity of the atomic orbitals. THCMP2 has been applied to both geometry optimization and ab initio molecular dynamics (AIMD) simulations. The resulting energy conservation in microcanonical AIMD demonstrates that the implementation provides accurate nuclear gradients with respect to the THCMP2 potential energy surfaces.

The truncated conjugate gradient (TCG), a noniterative/fixedcost strategy for computing polarization in molecular dynamics: Fast evaluation of analytical forces J. Chem. Phys. (IF 2.965) Pub Date : 20170829
Félix Aviat, Louis Lagardère, JeanPhilip PiquemalIn a recent paper [F. Aviat et al., J. Chem. Theory Comput. 13, 180–190 (2017)], we proposed the Truncated Conjugate Gradient (TCG) approach to compute the polarization energy and forces in polarizable molecular simulations. The method consists in truncating the conjugate gradient algorithm at a fixed predetermined order leading to a fixed computational cost and can thus be considered “noniterative.” This gives the possibility to derive analytical forces avoiding the usual energy conservation (i.e., drifts) issues occurring with iterative approaches. A key point concerns the evaluation of the analytical gradients, which is more complex than that with a usual solver. In this paper, after reviewing the present state of the art of polarization solvers, we detail a viable strategy for the efficient implementation of the TCG calculation. The complete cost of the approach is then measured as it is tested using a multitime step scheme and compared to timings using usual iterative approaches. We show that the TCG methods are more efficient than traditional techniques, making it a method of choice for future long molecular dynamics simulations using polarizable force fields where energy conservation matters. We detail the various steps required for the implementation of the complete method by software developers.

Improving the accuracy of MøllerPlesset perturbation theory with neural networks J. Chem. Phys. (IF 2.965) Pub Date : 20170906
Robert T. McGibbon, Andrew G. Taube, Alexander G. Donchev, Karthik Siva, Felipe Hernández, Cory Hargus, KaHei Law, John L. Klepeis, David E. ShawNoncovalent interactions are of fundamental importance across the disciplines of chemistry, materials science, and biology. Quantum chemical calculations on noncovalently bound complexes, which allow for the quantification of properties such as binding energies and geometries, play an essential role in advancing our understanding of, and building models for, a vast array of complex processes involving molecular association or selfassembly. Because of its relatively modest computational cost, secondorder MøllerPlesset perturbation (MP2) theory is one of the most widely used methods in quantum chemistry for studying noncovalent interactions. MP2 is, however, plagued by serious errors due to its incomplete treatment of electron correlation, especially when modeling van der Waals interactions and πstacked complexes. Here we present spinnetworkscaled MP2 (SNSMP2), a new semiempirical MP2based method for dimer interactionenergy calculations. To correct for errors in MP2, SNSMP2 uses quantum chemical features of the complex under study in conjunction with a neural network to reweight terms appearing in the total MP2 interaction energy. The method has been trained on a new data set consisting of over 200 000 complete basis set (CBS)extrapolated coupledcluster interaction energies, which are considered the gold standard for chemical accuracy. SNSMP2 predicts goldstandard binding energies of unseen test compounds with a mean absolute error of 0.04 kcal mol−1 (rootmeansquare error 0.09 kcal mol−1), a 6 to 7fold improvement over MP2. To the best of our knowledge, its accuracy exceeds that of all extant density functional theory and wavefunctionbased methods of similar computational cost, and is very close to the intrinsic accuracy of our benchmark coupledcluster methodology itself. Furthermore, SNSMP2 provides reliable perconformation confidence intervals on the predicted interaction energies, a feature not available from any alternative method.

Performing the Millikan experiment at the molecular scale: Determination of atomic MillikanThomson charges by computationally measuring atomic forces J. Chem. Phys. (IF 2.965) Pub Date : 20170907
T. Ryan Rogers, Feng WangAn atomic version of the Millikan oil drop experiment is performed computationally. It is shown that for planar molecules, the atomic version of the Millikan experiment can be used to define an atomic partial charge that is free from charge flow contributions. We refer to this charge as the MillikanThomson (MT) charge. Since the MT charge is directly proportional to the atomic forces under a uniform electric field, it is the most relevant charge for force field developments. The MT charge shows good stability with respect to different choices of the basis set. In addition, the MT charge can be easily calculated even at postHartreeFock levels of theory. With the MT charge, it is shown that for a planar water dimer, the charge transfer from the proton acceptor to the proton donor is about −0.052 e. While both planar hydrated cations and anions show signs of charge transfer, anions show a much more significant charge transfer to the hydration water than the corresponding cations. It might be important to explicitly model the ion charge transfer to water in a force field at least for the anions.

The BioFragment Database (BFDb): An opendata platform for computational chemistry analysis of noncovalent interactions J. Chem. Phys. (IF 2.965) Pub Date : 20170911
Lori A. Burns, John C. Faver, Zheng Zheng, Michael S. Marshall, Daniel G. A. Smith, Kenno Vanommeslaeghe, Alexander D. MacKerellJr., Kenneth M. MerzJr., C. David SherrillAccurate potential energy models are necessary for reliable atomistic simulations of chemical phenomena. In the realm of biomolecular modeling, large systems like proteins comprise very many noncovalent interactions (NCIs) that can contribute to the protein’s stability and structure. This work presents two highquality chemical databases of common fragment interactions in biomolecular systems as extracted from highresolution Protein DataBank crystal structures: 3380 sidechainsidechain interactions and 100 backbonebackbone interactions that inaugurate the BioFragment Database (BFDb). Absolute interaction energies are generated with a computationally tractable explicitly correlated coupled cluster with perturbative triples [CCSD(T)F12] “silver standard” (0.05 kcal/mol average error) for NCI that demands only a fraction of the cost of the conventional “gold standard,” CCSD(T) at the complete basis set limit. By sampling extensively from biological environments, BFDb spans the natural diversity of protein NCI motifs and orientations. In addition to supplying a thorough assessment for lower scaling forcefield (2), semiempirical (3), density functional (244), and wavefunction (45) methods (comprising >1M interaction energies), BFDb provides interactive tools for running and manipulating the resulting large datasets and offers a valuable resource for potential energy model development and validation.

Quasichemical theory of F^{−}(aq): The “no split occupancies rule” revisited J. Chem. Phys. (IF 2.965) Pub Date : 20170912
Mangesh I. Chaudhari, Susan B. Rempe, Lawrence R. PrattWe use ab initio molecular dynamics (AIMD) calculations and quasichemical theory (QCT) to study the innershell structure of F−(aq) and to evaluate that singleion free energy under standard conditions. Following the “no split occupancies” rule, QCT calculations yield a free energy value of −101 kcal/mol under these conditions, in encouraging agreement with tabulated values (−111 kcal/mol). The AIMD calculations served only to guide the definition of an effective innershell constraint. QCT naturally includes quantum mechanical effects that can be concerning in more primitive calculations, including electronic polarizability and induction, electron density transfer, electron correlation, molecular/atomic cooperative interactions generally, molecular flexibility, and zeropoint motion. No direct assessment of the contribution of dispersion contributions to the internal energies has been attempted here, however. We anticipate that other aqueous halide ions might be treated successfully with QCT, provided that the structure of the underlying statistical mechanical theory is absorbed, i.e., that the “no split occupancies” rule is recognized.

Understanding the manybody expansion for large systems. III. Critical role of fourbody terms, counterpoise corrections, and cutoffs J. Chem. Phys. (IF 2.965) Pub Date : 20170913
KuanYu Liu, John M. HerbertPapers I and II in this series [R. M. Richard et al., J. Chem. Phys. 141, 014108 (2014); K. U. Lao et al., ibid. 144, 164105 (2016)] have attempted to shed light on precision and accuracy issues affecting the manybody expansion (MBE), which only manifest in larger systems and thus have received scant attention in the literature. Manybody counterpoise (CP) corrections are shown to accelerate convergence of the MBE, which otherwise suffers from a mismatch between how basisset superposition error affects subsystem versus supersystem calculations. In water clusters ranging in size up to (H2O)37, fourbody terms prove necessary to achieve accurate results for both total interaction energies and relative isomer energies, but the sheer number of tetramers makes the use of cutoff schemes essential. To predict relative energies of (H2O)20 isomers, two approximations based on a lower level of theory are introduced and an ONIOMtype procedure is found to be very well converged with respect to the appropriate MBE benchmark, namely, a CPcorrected supersystem calculation at the same level of theory. Results using an energybased cutoff scheme suggest that if reasonable approximations to the subsystem energies are available (based on classical multipoles, say), then the number of requisite subsystem calculations can be reduced even more dramatically than when distancebased thresholds are employed. The end result is several accurate fourbody methods that do not require charge embedding, and which are stable in large basis sets such as augccpVTZ that have sometimes proven problematic for fragmentbased quantum chemistry methods. Even with aggressive thresholding, however, the fourbody approach at the selfconsistent field level still requires roughly ten times more processors to outmatch the performance of the corresponding supersystem calculation, in test cases involving 1500–1800 basis functions.

Links between the charge model and bonded parameter force constants in biomolecular force fields J. Chem. Phys. (IF 2.965) Pub Date : 20171004
David S. Cerutti, Karl T. Debiec, David A. Case, Lillian T. ChongThe ff15ipq protein force field is a fixed charge model built by automated tools based on the two charge sets of the implicitly polarized charge method: one set (appropriate for vacuum) for deriving bonded parameters and the other (appropriate for aqueous solution) for running simulations. The duality is intended to treat waterinduced electronic polarization with an understanding that fitting data for bonded parameters will come from quantum mechanical calculations in the gas phase. In this study, we compare ff15ipq to two alternatives produced with the same fitting software and a further expanded data set but following more conventional methods for tailoring bonded parameters (harmonic angle terms and torsion potentials) to the charge model. First, ff15ipqQsolv derives bonded parameters in the context of the ff15ipq solution phase charge set. Second, ff15ipqVac takes ff15ipq’s bonded parameters and runs simulations with the vacuum phase charge set used to derive those parameters. The IPolQ charge model and associated protocol for deriving bonded parameters are shown to be an incremental improvement over protocols that do not account for the material phases of each source of their fitting data. Both force fields incorporating the polarized charge set depict stable globular proteins and have varying degrees of success modeling the metastability of short (5–19 residues) peptides. In this particular case, ff15ipqQsolv increases stability in a number of αhelices, correctly obtaining 70% helical character in the K19 system at 275 K and showing appropriately diminishing content up to 325 K, but overestimating the helical fraction of AAQAA3 by 50% or more, forming longlived αhelices in simulations of a βhairpin, and increasing the likelihood that the disordered p53 Nterminal peptide will also form a helix. This may indicate a systematic bias imparted by the ff15ipqQsolv parameter development strategy, which has the hallmarks of strategies used to develop other popular force fields, and may explain some of the need for manual corrections in this force fields’ evolution. In contrast, ff15ipqVac incorrectly depicts globular protein unfolding in numerous systems tested, including Trp cage, villin, lysozyme, and GB3, and does not perform any better than ff15ipq or ff15ipqQsolv in tests on short peptides. We analyze the free energy surfaces of individual amino acid dipeptides and the electrostatic potential energy surfaces of each charge model to explain the differences.

Computational and experimental characterization of a pyrrolidiniumbased ionic liquid for electrolyte applications J. Chem. Phys. (IF 2.965) Pub Date : 20171004
Hedieh Torabifard, Luke Reed, Matthew T. Berry, Jason E. Hein, Erik Menke, G. Andrés CisnerosThe development of Liion batteries for energy storage has received significant attention. The synthesis and characterization of electrolytes in these batteries are an important component of this development. Ionic liquids (ILs) have been proposed as possible electrolytes in these devices. Thus, the accurate determination of thermophysical properties for these solvents becomes important for determining their applicability as electrolytes. In this contribution, we present the synthesis and experimental/computational characterization of thermodynamic and transport properties of a pyrrolidinium based ionic liquid as a first step to investigate the possible applicability of this class of ILs for Liion batteries. A quantum mechanicalbased force field with manybody polarizable interactions has been developed for the simulation of spirocyclic pyrrolidinium, [sPyr+], with BF4− and Li+. Molecular dynamics calculations employing intramolecular polarization predicted larger heat of vaporization and selfdiffusion coefficients and smaller densities in comparison with the model without intramolecular polarization, indicating that the inclusion of this term can significantly effect the interionic interactions. The calculated properties are in good agreement with available experimental data for similar IL pairs and isothermal titration calorimetry data for [sPyr+][BF4−].

Internal force corrections with machine learning for quantum mechanics/molecular mechanics simulations J. Chem. Phys. (IF 2.965) Pub Date : 20171012
Jingheng Wu, Lin Shen, Weitao YangAb initio quantum mechanics/molecular mechanics (QM/MM) molecular dynamics simulation is a useful tool to calculate thermodynamic properties such as potential of mean force for chemical reactions but intensely time consuming. In this paper, we developed a new method using the internal force correction for lowlevel semiempirical QM/MM molecular dynamics samplings with a predefined reaction coordinate. As a correction term, the internal force was predicted with a machine learning scheme, which provides a sophisticated force field, and added to the atomic forces on the reaction coordinate related atoms at each integration step. We applied this method to two reactions in aqueous solution and reproduced potentials of mean force at the ab initio QM/MM level. The saving in computational cost is about 2 orders of magnitude. The present work reveals great potentials for machine learning in QM/MM simulations to study complex chemical processes.

Study of interactions between metal ions and protein model compounds by energy decomposition analyses and the AMOEBA force field J. Chem. Phys. (IF 2.965) Pub Date : 20171017
Zhifeng Jing, Rui Qi, Chengwen Liu, Pengyu RenThe interactions between metal ions and proteins are ubiquitous in biology. The selective binding of metal ions has a variety of regulatory functions. Therefore, there is a need to understand the mechanism of proteinion binding. The interactions involving metal ions are complicated in nature, where shortrange chargepenetration, charge transfer, polarization, and manybody effects all contribute significantly, and a quantitative description of all these interactions is lacking. In addition, it is unclear how well current polarizable force fields can capture these energy terms and whether these polarization models are good enough to describe the manybody effects. In this work, two energy decomposition methods, absolutely localized molecular orbitals and symmetryadapted perturbation theory, were utilized to study the interactions between Mg2+/Ca2+ and model compounds for amino acids. Comparison of individual interaction components revealed that while there are significant chargepenetration and chargetransfer effects in Ca complexes, these effects can be captured by the van der Waals (vdW) term in the AMOEBA force field. The electrostatic interaction in Mg complexes is well described by AMOEBA since the charge penetration is small, but the distancedependent polarization energy is problematic. Manybody effects were shown to be important for proteinion binding. In the absence of manybody effects, highly charged binding pockets will be overstabilized, and the pockets will always favor Mg and thus lose selectivity. Therefore, manybody effects must be incorporated in the force field in order to predict the structure and energetics of metalloproteins. Also, the manybody effects of charge transfer in Ca complexes were found to be nonnegligible. The absorption of chargetransfer energy into the additive vdW term was a main source of error for the AMOEBA manybody interaction energies.

Fast and reliable ab initio calculation of crystal field splittings in lanthanide complexes J. Chem. Phys. (IF 2.965) Pub Date : 20171023
P. P. Hallmen, C. Köppl, G. Rauhut, H. Stoll, J. van SlagerenAb initio calculations of crystal field splittings and magnetic properties of lanthanide complexes are usually performed using stateaveraged complete active space selfconsistent field (CASSCF) calculations and a subsequent spinorbit calculation mixing the CASSCF wave functions (CASSCF/state interaction with spinorbit coupling). Because this approach becomes very timeconsuming for large molecules, simplifications have been proposed in the literature to determine the stateaveraged orbitals by configurationaveraged HartreeFock (CAHF) instead of CASSCF. We present an approach which is an extension of the CAHF method. We combine the techniques of local density fitting with CAHF and achieve a significant speedup compared to CASSCF without loss in accuracy. To assess the performance of our method, we apply it to three wellknown molecules, namely, Er[N(SiMe3)2]3, Er(trensal), and the doubledecker (NBu4)+ [Er(Pc)2]−.

FokkerPlanck equation for the nonMarkovian Brownian motion in the presence of a magnetic field J. Chem. Phys. (IF 2.965) Pub Date : 20171023
Joydip Das, Shrabani Mondal, Bidhan Chandra BagIn the present study, we have proposed the FokkerPlanck equation in a simple way for a Langevin equation of motion having ordinary derivative (OD), the Gaussian random force and a generalized frictional memory kernel. The equation may be associated with or without conservative force field from harmonic potential. We extend this method for a charged Brownian particle in the presence of a magnetic field. Thus, the present method is applicable for a Langevin equation of motion with OD, the Gaussian colored thermal noise and any kind of linear force field that may be conservative or not. It is also simple to apply this method for the colored Gaussian noise that is not related to the damping strength.

Acceleration and sensitivity analysis of lattice kinetic Monte Carlo simulations using parallel processing and rate constant rescaling J. Chem. Phys. (IF 2.965) Pub Date : 20171023
M. Núñez, T. Robie, D. G. VlachosKinetic Monte Carlo (KMC) simulation provides insights into catalytic reactions unobtainable with either experiments or meanfield microkinetic models. Sensitivity analysis of KMC models assesses the robustness of the predictions to parametric perturbations and identifies rate determining steps in a chemical reaction network. Stiffness in the chemical reaction network, a ubiquitous feature, demands lengthy run times for KMC models and renders efficient sensitivity analysis based on the likelihood ratio method unusable. We address the challenge of efficiently conducting KMC simulations and performing accurate sensitivity analysis in systems with unknown time scales by employing two acceleration techniques: rate constant rescaling and parallel processing. We develop statistical criteria that ensure sufficient sampling of nonequilibrium steady state conditions. Our approach provides the twofold benefit of accelerating the simulation itself and enabling likelihood ratio sensitivity analysis, which provides further speedup relative to finite difference sensitivity analysis. As a result, the likelihood ratio method can be applied to real chemistry. We apply our methodology to the watergas shift reaction on Pt(111).

Meshfree hierarchical clustering methods for fast evaluation of electrostatic interactions of point multipoles J. Chem. Phys. (IF 2.965) Pub Date : 20171023
H. A. BoatengElectrostatic interactions involving point multipoles are being increasingly implemented to achieve higher accuracy in molecular simulations. A major drawback of multipolar electrostatics is the increased computational cost. Here we develop and compare two Cartesian tree algorithms which employ Taylor approximations and hierarchical clustering to speed up the evaluation of point multipole interactions. We present results from applying the algorithms to compute the free space Coulomb potential and forces of different sets of interacting point multipoles with different densities. The methods achieve high accuracy and speedup of more than an order of magnitude over direct sum calculations and scale well in parallel.

Crossing conditions in coupled cluster theory J. Chem. Phys. (IF 2.965) Pub Date : 20171024
Eirik F. Kjønstad, Rolf H. Myhre, Todd J. Martínez, Henrik KochWe derive the crossing conditions at conical intersections between electronic states in coupled cluster theory and show that if the coupled cluster Jacobian matrix is nondefective, two (three) independent conditions are correctly placed on the nuclear degrees of freedom for an inherently real (complex) Hamiltonian. Calculations using coupled cluster theory on a 21A′/31A′ conical intersection in hypofluorous acid illustrate the nonphysical artifacts associated with defects at accidental samesymmetry intersections. In particular, the observed intersection seam is folded about a space of the correct dimensionality, indicating that minor modifications to the theory are required for it to provide a correct description of conical intersections in general. We find that an accidental symmetry allowed 11A″/21A″ intersection in hydrogen sulfide is properly described, showing no artifacts as well as linearity of the energy gap to first order in the branching plane.

Localmetrics errorbased Shepard interpolation as surrogate for highly nonlinear material models in high dimensions J. Chem. Phys. (IF 2.965) Pub Date : 20171024
Juan M. Lorenzi, Thomas Stecher, Karsten Reuter, Sebastian MateraMany problems in computational materials science and chemistry require the evaluation of expensive functions with locally rapid changes, such as the turnover frequency of first principles kinetic Monte Carlo models for heterogeneous catalysis. Because of the high computational cost, it is often desirable to replace the original with a surrogate model, e.g., for use in coupled multiscale simulations. The construction of surrogates becomes particularly challenging in highdimensions. Here, we present a novel version of the modified Shepard interpolation method which can overcome the curse of dimensionality for such functions to give faithful reconstructions even from very modest numbers of function evaluations. The introduction of local metrics allows us to take advantage of the fact that, on a local scale, rapid variation often occurs only across a small number of directions. Furthermore, we use local error estimates to weigh different local approximations, which helps avoid artificial oscillations. Finally, we test our approach on a number of challenging analytic functions as well as a realistic kinetic Monte Carlo model. Our method not only outperforms existing isotropic metric Shepard methods but also stateoftheart Gaussian process regression.
Some contents have been Reproduced by permission of The Royal Society of Chemistry.