Pairwise Alternatives to Ewald Summation for Calculating Long-Range Electrostatics in Ionic Liquids

Monte Carlo QM/MM simulation studies of the Cannizzaro reaction in ionic liquids for improved biofuel production

The conversion of biomass to 5-hydroxymethylfurfural (HMF) holds substantial promise as a renewable energy source. Notably, HMF can be transformed into 2,5-bis(hydroxymethyl)furan (BHMF), a crucial reactant in biofuel production, but requires harsh operating conditions, H2, and precious metal catalysts. A recently reported Cannizzaro reaction of HMF to BHMF, characterized by its efficiency, mild conditions, and eco-friendliness, instead employed ionic liquids (ILs) to achieve high yields. In this study, combined quantum mechanical and molecular mechanical (QM/MM) simulations in conjunction with Metropolis Monte Carlo statistical mechanics and free-energy perturbation theory utilized M06-2X/6-31+G(d), PDDG/PM3, and the OPLS-VSIL force field to uncover important solute-solvent interactions present in the HMF to BHMF reaction pathway. The Cannizzaro reaction was examined in water and in five ILs composed of the 1-butyl-3-methylimidazolium [BMIM] cation coupled to hexafluorophosphate, tetrafluoroborate, thiocyanate, chloride, and bromide. Energetic and structural analysis of the rate-determining hydride transfer between HMF and the hydride-donor anion HMFOH- attributed the enhanced reactivity to highly organized solvent interactions featuring (1) hydrogen bonding between the ring protons of [BMIM] and the negatively charged carbonyl oxygen atoms on the transition structure, (2) favorable electrostatic interactions between the IL anions and solute hydroxyl groups, and (3) beneficial π-π stacking interactions between [BMIM] and the two aromatic rings present in HMF and HMFOH-.

Unification of Ewald and shifted force methods to calculate Coulomb interactions in molecular simulations

Three new Ewald series are derived using a new strategy that does not start with a proposed charge spreading function. Of these, the Ewald series produced using shifted potential interactions for the Ewald real space series converges relatively slowly, while the corresponding expression using a shifted force (SF) interaction does not converge. A comparison is made between several approximations of the Ewald method and the SF route to include Coulomb interactions in molecular dynamics (MD) computer simulations. MD simulations of a model bulk molten salt and water were carried out. The recently derived α' variant of Ewald, by K. D. Hammonds and D. M. Heyes [J. Chem. Phys. 157, 074108 (2022)], has been developed analytically and found to be more accurate and computationally efficient than SF in part due to the smaller real space truncation distance that can be used. In addition, with α', the number of reciprocal lattice vectors required is reduced considerably compared with the standard Ewald implementations to give the same accuracy. The invention of the α' method shifts the computational balance back toward using an Ewald construction. The SF method shows greater errors in the Coulomb pressure and time dependent fluctuation properties compared to α'. It does not conserve the shadow Hamiltonian in a microcanonical MD simulation, whereas the α' method does, which facilitates long time stability and insignificant drift of properties over time. The speed of the Ewald computer code is improved by using a new lookup table method.

Analytic Gradients for the Electrostatic Embedding QM/MM Model in Periodic Boundary Conditions Using Particle-Mesh Ewald Sums and Electrostatic Potential Fitted Charge Operators

Long-range electrostatic effects are fundamental for describing chemical reactivity in the condensed phase. Here, we present the methodology of an efficient quantum mechanical/molecular mechanical (QM/MM) model in periodic boundary conditions (PBC) compatible with QM/MM boundaries at chemical bonds. The method combines electrostatic potential fitted charge operators and electrostatic potentials derived from the smooth particle-mesh Ewald (PME) sum approach. The total energy and its analytic first derivatives with respect to QM, MM, and lattice vectors allow QM/MM molecular dynamics (MD) in the most common thermodynamic ensembles. We demonstrate the robustness of the method by performing a QM/MM MD equilibration of methanol in water. We simulate the cis/trans isomerization free-energy profiles in water of proline amino acid and a proline-containing oligopeptide, showing a correct description of the reaction barrier. Our PBC-compatible QM/MM model can efficiently be used to study the chemical reactivity in the condensed phase and enzymatic catalysis.

Non-Ewald methods for evaluating the electrostatic interactions of charge systems: similarity and difference

In molecular simulations, it is essential to properly calculate the electrostatic interactions of particles in the physical system of interest. Here we consider a method called the non-Ewald method, which does not rely on the standard Ewald method with periodic boundary conditions, but instead relies on the cutoff-based techniques. We focus on the physicochemical and mathematical conceptual aspects of the method in order to gain a deeper understanding of the simulation methodology. In particular, we take into account the reaction field (RF) method, the isotropic periodic sum (IPS) method, and the zero-multipole summation method (ZMM). These cutoff-based methods are based on different physical ideas and are completely distinguishable in their underlying concepts. The RF and IPS methods are "additive" methods that incorporate information outside the cutoff region, via dielectric medium and isotropic boundary condition, respectively. In contrast, the ZMM is a "subtraction" method that tries to remove the artificial effects, generated near the boundary, from the cutoff sphere. Nonetheless, we find physical and/or mathematical similarities between these methods. In particular, the modified RF method can be derived by the principle of neutralization utilized in the ZMM, and we also found a direct relationship between IPS and ZMM.

DFT-based QM/MM with Particle-Mesh Ewald for Direct, Long-Range Electrostatic Embedding

We present a DFT-based, QM/MM implementation with long-range electrostatic embedding achieved by direct real-space integration of the particle mesh Ewald (PME) computed electrostatic potential. The key transformation is the interpolation of the electrostatic potential from the PME grid to the DFT quadrature grid, from which integrals are easily evaluated utilizing standard DFT machinery. We provide benchmarks of the numerical accuracy with choice of grid size and real-space corrections, and demonstrate that good convergence is achieved while introducing nominal computational overhead. Furthermore, the approach requires only small modification to existing software packages, as is demonstrated with our implementation in the OpenMM and Psi4 software. After presenting convergence benchmarks, we evaluate the importance of long-range electrostatic embedding in three solute/solvent systems modeled with QM/MM. Water and BMIM/BF4 ionic liquid were considered as ``simple' and ``complex' solvents respectively, with water and p-phenylenediamine (PPD) solute molecules treated at QM level of theory. While electrostatic embedding with standard real-space truncation may introduce negligible error for simple systems such as water solute in water solvent, errors become more significant when QM/MM is applied to complex solvents such as ionic liquids. An extreme example is the electrostatic embedding energy for oxidized PPD in BMIM/BF4 for which real-space truncation produces severe error even at 2-3 nm cutoff distances. This latter example illustrates that utilization of QM/MM to compute redox potentials within concentrated electrolytes/ionic media requires carefully chosen long-range electrostatic embedding algorithms, with our presented algorithm providing a general and robust approach.

Development and application of quantum mechanics/molecular mechanics methods with advanced polarizable potentials

Quantum mechanics/molecular mechanics (QM/MM) simulations are a popular approach to study various features of large systems. A common application of QM/MM calculations is in the investigation of reaction mechanisms in condensed-phase and biological systems. The combination of QM and MM methods to represent a system gives rise to several challenges that need to be addressed. The increase in computational speed has allowed the expanded use of more complicated and accurate methods for both QM and MM simulations. Here, we review some approaches that address several common challenges encountered in QM/MM simulations with advanced polarizable potentials, from methods to account for boundary across covalent bonds and long-range effects, to polarization and advanced embedding potentials.

Assessment of the Wolf method using the Stillinger-Lovett sum rules: From strong electrolytes to weakly charged colloidal dispersions

The Ewald method has been the cornerstone in molecular simulations for modeling electrostatic interactions of charge-stabilized many-body systems. In the late 1990s, Wolf and collaborators developed an alternative route to describe the long-range nature of electrostatic interactions; from a computational perspective, this method provides a more efficient and straightforward way to implement long-range electrostatic interactions than the Ewald method. Despite these advantages, the validity of the Wolf potential to account for the electrostatic contribution in charged fluids remains controversial. To alleviate this situation, in this contribution, we implement the Wolf summation method to both electrolyte solutions and charged colloids with moderate size and charge asymmetries in order to assess the accuracy and validity of the method. To this end, we verify that the proper selection of parameters within the Wolf method leads to results that are in good agreement with those obtained through the standard Ewald method and the theory of integral equations of simple liquids within the so-called hypernetted chain approximation. Furthermore, we show that the results obtained with the original Wolf method do satisfy the moment conditions described by the Stillinger-Lovett sum rules, which are directly related to the local electroneutrality condition and the electrostatic screening in the Debye-Hückel regime. Hence, the fact that the solution provided by the Wolf method satisfies the first and second moments of Stillinger-Lovett proves, for the first time, the reliability of the method to correctly incorporate the electrostatic contribution in charge-stabilized fluids. This makes the Wolf method a powerful alternative compared to more demanding computational approaches.

Current Status of AMOEBA-IL: A Multipolar/Polarizable Force Field for Ionic Liquids

Computational simulations of ionic liquid solutions have become a useful tool to investigate various physical, chemical and catalytic properties of systems involving these solvents. Classical molecular dynamics and hybrid quantum mechanical/molecular mechanical (QM/MM) calculations of IL systems have provided significant insights at the atomic level. Here, we present a review of the development and application of the multipolar and polarizable force field AMOEBA for ionic liquid systems, termed AMOEBA–IL. The parametrization approach for AMOEBA–IL relies on the reproduction of total quantum mechanical (QM) intermolecular interaction energies and QM energy decomposition analysis. This approach has been used to develop parameters for imidazolium– and pyrrolidinium–based ILs coupled with various inorganic anions. AMOEBA–IL has been used to investigate and predict the properties of a variety of systems including neat ILs and IL mixtures, water exchange reactions on lanthanide ions in IL mixtures, IL–based liquid–liquid extraction, and effects of ILs on an aniline protection reaction.

Speeding up the study of diluted dipolar systems

We study the regimes of a diluted dipolar system through Monte Carlo numerical simulations in the NVT ensemble. To accelerate the dynamics, several approximations and speed-up algorithms are proposed and tested. In particular, it turns out that "cluster move Monte Carlo" algorithm speeds-up to two decades faster than the traditional Monte Carlo, depending on temperature and density. We find simple-fluid, chain-fluid, ring-fluid, gel, and antiparallel columnar regimes, which are studied and characterized through positional, orientational, and thermodynamical observables.

LICHEM 1.1: Recent Improvements and New Capabilities

The QM/MM method has become a useful tool to investigate various properties of complex systems. We previously introduced the layered interacting chemical models (LICHEM) package to enable QM/MM simulations with advanced potentials by combining various (unmodified) QM and MM codes ( J. Comp. Chem., 2016, 37, 1019). LICHEM provides several capabilities such as the ability to use polarizable force fields, such as AMOEBA, for the MM environment. Here, we describe an updated version of LICHEM (v1.1), which includes several new functionalities including a new method to account for long-range electrostatic effects in QM/MM (QM/MM-LREC), a new implementation for QM/MM with the Gaussian electrostatic model (GEM), and new capabilities for path optimizations using the quadratic string model (QSM) coupled with restrained MM environment optimization.

On short-ranged pair-potentials for long-range electrostatics

Fast and accurate summation of long-range electrostatics by using a short-ranged pair-potential that ensures moment cancellation in the cutoff sphere.

A Modified Shifted Force Approach to the Wolf Summation

The Wolf method for calculation of electrostatic interactions in molecular simulations is known to describe the energy well, whereas the forces have discontinuities. For a more reliable description of the forces this method can be extended with a shifted force approach. This leads to a good description of the forces and precise molecular dynamics simulation, but the description of the energy becomes poorer. In this study we propose a modification of a shifted force extension to describe the energy as well as the forces in better agreement to reference data as determined from the Ewald summation. We show that vapor-liquid phase equilibria (VLE) calculated with Monte Carlo simulations in the grand canonical ensemble and dynamic properties calculated with molecular dynamics simulations can be calculated reliably using this modification to describe the electrostatic interactions.

Molecular Simulation of Vapor-Liquid Equilibria Using the Wolf Method for Electrostatic Interactions

The applicability of the Wolf method for calculating electrostatic interactions is verified for simulating vapor-liquid equilibria of hydrogen sulfide, methanol, and carbon dioxide. Densities, chemical potentials, and critical properties are obtained with Monte Carlo simulations using the Continuous Fractional Component version of the Gibbs Ensemble. Saturated vapor pressures are obtained from NPT simulations. Excellent agreement is found between simulation results and data from literature (simulations using the Ewald summation). It is also shown how to choose the optimal parameters for the Wolf method. Even though the Wolf method requires a large simulation box in the gas phase, due to the lack of screening of electrostatics, one can consider the Wolf method as a suitable alternative to the Ewald summation in VLE calculations.

Single-pass incremental force updates for adaptively restrained molecular dynamics

Adaptively restrained molecular dynamics (ARMD) allows users to perform more integration steps in wall-clock time by switching on and off positional degrees of freedoms. This article presents new, single-pass incremental force updates algorithms to efficiently simulate a system using ARMD. We assessed different algorithms for speedup measurements and implemented them in the LAMMPS MD package. We validated the single-pass incremental force update algorithm on four different benchmarks using diverse pair potentials. The proposed algorithm allows us to perform simulation of a system faster than traditional MD in both NVE and NVT ensembles. Moreover, ARMD using the new single-pass algorithm speeds up the convergence of observables in wall-clock time. © 2017 Wiley Periodicals, Inc.

Treating electrostatics with Wolf summation in combined quantum mechanical and molecular mechanical simulations

The Wolf summation approach [D. Wolf et al., J. Chem. Phys. 110, 8254 (1999)], in the damped shifted force (DSF) formalism [C. J. Fennell and J. D. Gezelter, J. Chem. Phys. 124, 234104 (2006)], is extended for treating electrostatics in combined quantum mechanical and molecular mechanical (QM/MM) molecular dynamics simulations. In this development, we split the QM/MM electrostatic potential energy function into the conventional Coulomb r(-1) term and a term that contains the DSF contribution. The former is handled by the standard machinery of cutoff-based QM/MM simulations whereas the latter is incorporated into the QM/MM interaction Hamiltonian as a Fock matrix correction. We tested the resulting QM/MM-DSF method for two solution-phase reactions, i.e., the association of ammonium and chloride ions and a symmetric SN2 reaction in which a methyl group is exchanged between two chloride ions. The performance of the QM/MM-DSF method was assessed by comparing the potential of mean force (PMF) profiles with those from the QM/MM-Ewald and QM/MM-isotropic periodic sum (IPS) methods, both of which include long-range electrostatics explicitly. For ion association, the QM/MM-DSF method successfully eliminates the artificial free energy drift observed in the QM/MM-Cutoff simulations, in a remarkable agreement with the two long-range-containing methods. For the SN2 reaction, the free energy of activation obtained by the QM/MM-DSF method agrees well with both the QM/MM-Ewald and QM/MM-IPS results. The latter, however, requires a greater cutoff distance than QM/MM-DSF for a proper convergence of the PMF. Avoiding time-consuming lattice summation, the QM/MM-DSF method yields a 55% reduction in computational cost compared with the QM/MM-Ewald method. These results suggest that, in addition to QM/MM-IPS, the QM/MM-DSF method may serve as another efficient and accurate alternative to QM/MM-Ewald for treating electrostatics in condensed-phase simulations of chemical reactions.

A new smoothing function to introduce long-range electrostatic effects in QM/MM calculations

A new method to account for long range electrostatic contributions is proposed and implemented for quantum mechanics/molecular mechanics long range electrostatic correction (QM/MM-LREC) calculations. This method involves the use of the minimum image convention under periodic boundary conditions and a new smoothing function for energies and forces at the cutoff boundary for the Coulomb interactions. Compared to conventional QM/MM calculations without long-range electrostatic corrections, the new method effectively includes effects on the MM environment in the primary image from its replicas in the neighborhood. QM/MM-LREC offers three useful features including the avoidance of calculations in reciprocal space (k-space), with the concomitant avoidance of having to reproduce (analytically or approximately) the QM charge density in k-space, and the straightforward availability of analytical Hessians. The new method is tested and compared with results from smooth particle mesh Ewald (PME) for three systems including a box of neat water, a double proton transfer reaction, and the geometry optimization of the critical point structures for the rate limiting step of the DNA dealkylase AlkB. As with other smoothing or shifting functions, relatively large cutoffs are necessary to achieve comparable accuracy with PME. For the double-proton transfer reaction, the use of a 22 Å cutoff shows a close reaction energy profile and geometries of stationary structures with QM/MM-LREC compared to conventional QM/MM with no truncation. Geometry optimization of stationary structures for the hydrogen abstraction step by AlkB shows some differences between QM/MM-LREC and the conventional QM/MM. These differences underscore the necessity of the inclusion of the long-range electrostatic contribution.

Assessing the accuracy of the isotropic periodic sum method through Madelung energy computation

We tested the isotropic periodic sum (IPS) method for computing Madelung energies of ionic crystals. The performance of the method, both in its nonpolar (IPSn) and polar (IPSp) forms, was compared with that of the zero-charge and Wolf potentials [D. Wolf, P. Keblinski, S. R. Phillpot, and J. Eggebrecht, J. Chem. Phys. 110, 8254 (1999)]. The results show that the IPSn and IPSp methods converge the Madelung energy to its reference value with an average deviation of ∼10(-4) and ∼10(-7) energy units, respectively, for a cutoff range of 18-24a (a/2 being the nearest-neighbor ion separation). However, minor oscillations were detected for the IPS methods when deviations of the computed Madelung energies were plotted on a logarithmic scale as a function of the cutoff distance. To remove such oscillations, we introduced a modified IPSn potential in which both the local-region and long-range electrostatic terms are damped, in analogy to the Wolf potential. With the damped-IPSn potential, a smoother convergence was achieved. In addition, we observed a better agreement between the damped-IPSn and IPSp methods, which suggests that damping the IPSn potential is in effect similar to adding a screening potential in IPSp.

Simulating chemical reactions in ionic liquids using QM/MM methodology

The use of ionic liquids as a reaction medium for chemical reactions has dramatically increased in recent years due in large part to the numerous reported advances in catalysis and organic synthesis. In some extreme cases, ionic liquids have been shown to induce mechanistic changes relative to conventional solvents. Despite the large interest in the solvents, a clear understanding of the molecular factors behind their chemical impact is largely unknown. This feature article reviews our efforts developing and applying mixed quantum and molecular mechanical (QM/MM) methodology to elucidate the microscopic details of how these solvents operate to enhance rates and alter mechanisms for industrially and academically important reactions, e.g., Diels-Alder, Kemp eliminations, nucleophilic aromatic substitutions, and β-eliminations. Explicit solvent representation provided the medium dependence of the activation barriers and atomic-level characterization of the solute-solvent interactions responsible for the experimentally observed "ionic liquid effects". Technical advances are also discussed, including a linear-scaling pairwise electrostatic interaction alternative to Ewald sums, an efficient polynomial fitting method for modeling proton transfers, and the development of a custom ionic liquid OPLS-AA force field.

Quantum and Molecular Mechanical (QM/MM) Monte Carlo Techniques for Modeling Condensed-Phase Reactions

A recent review (Acc. Chem. Res. 2010, 43:142-151) examined our use and development of a combined quantum and molecular mechanical (QM/MM) technique for modelling organic and enzymatic reactions. Advances included the PDDG/PM3 semiempirical QM (SQM) method, computation of multi-dimensional potentials of mean force (PMF), incorporation of on-the-fly QM in Monte Carlo simulations, and a polynomial quadrature method for rapidly treating proton-transfer reactions. The current article serves as a follow up on our progress. Highlights include new reactions, alternative SQM methods, a polarizable OPLS force field, and novel solvent environments, e.g., "on water" and room temperature ionic liquids. The methodology is strikingly accurate across a wide range of condensed-phase and antibody-catalyzed reactions including substitution, decarboxylation, elimination, isomerization, and pericyclic classes. Comparisons are made to systems treated with continuum-based solvents and ab initio or density functional theory (DFT) methods. Overall, the QM/MM methodology provides detailed characterization of reaction paths, proper configurational sampling, several advantages over implicit solvent models, and a reasonable computational cost.