Analysis of Nuclear Quantum Effects on Hydrogen Bonding

Nuclear-Electronic Orbital QM/MM Approach: Geometry Optimizations and Molecular Dynamics

Hybrid quantum mechanical/molecular mechanical (QM/MM) methods allow simulations of chemical reactions in atomistic solvent and heterogeneous environments such as proteins. Herein, the nuclear-electronic orbital (NEO) QM/MM approach is introduced to enable the quantization of specified nuclei, typically protons, in the QM region using a method such as NEO-density functional theory (NEO-DFT). This approach includes proton delocalization, polarization, anharmonicity, and zero-point energy in geometry optimizations and dynamics. Expressions for the energies and analytical gradients associated with the NEO-QM/MM method, as well as the previously developed polarizable continuum model (NEO-PCM), are provided. Geometry optimizations of small organic molecules hydrogen bonded to water in either dielectric continuum solvent or explicit atomistic solvent illustrate that aqueous solvation can strengthen hydrogen-bonding interactions for the systems studied, as indicated by shorter intermolecular distances at the hydrogen-bond interface. We then performed a real-time direct dynamics simulation of a phenol molecule in explicit water using the NEO-QM/MM method. These developments and initial examples provide the foundation for future studies of nuclear-electronic quantum dynamics in complex chemical and biological environments.

A Tug of War between the Self- and Cross-associating Hydrogen Bonds in Neutral Ammonia-Water Clusters: Energetic Insights by Molecular Tailoring Approach

In the present work, the energies of various types of individual HBs observed in neutral (NH₃ )_m (H₂ O)_n , (m+n=2 to 7) clusters were estimated using the molecular tailoring approach (MTA)-based method. The calculated individual HB energies suggest that the O-H…N HBs are the strongest (1.21 to 12.49 kcal mol^-1 ). The next ones are the O-H…O (3.97 to 9.30 kcal mol^-1 ) HBs. The strengths of N-H…N (1.09 to 5.29 kcal mol^-1 ) and N-H…O (2.85 to 5.56 kcal mol^-1 ) HBs are the weakest. The HB energies in dimers also follow this rank ordering. However, the HB energies in dimers are much smaller than those obtained by the MTA-based method due to the loss in cooperativity contribution in the dimers. Thus, the calculated cooperativity contributions, for different types of HBs, fall in the range 0.64 to 5.73 kcal mol^-1 . We wish to emphasize based on the energetic rank ordering obtained by the MTA-based method that the O-H of water is a better HB donor than the N-H of ammonia. The reasons for the observed energetic rank ordering are two folds: (i) intrinsically stronger O-H…N HBs than the O-H…O ones as revealed by dimer energies and (ii) the higher cooperativity contribution in the former than the later ones. Indeed, the MTA-based method is useful in providing the missing energetic rank ordering of various type of HBs in neutral (NH₃ )_m (H₂ O)_n clusters, in the literature.

Energetic Ordering of Hydrogen Bond Strengths in Methanol-Water Clusters: Insights via Molecular Tailoring Approach

In this work, we examine the strength of various types of individual hydrogen bond (HB) in mixed methanol-water M_n W_m , (n+m=2 to 7) clusters, with an aim to understand the relative order of their strength, using our recently proposed molecular tailoring-based approach (MTA). Among all the types of HB, it is observed that the O_M -H…O_W HBs are the strongest (6.9 to 12.4 kcal mol^-1 ). The next ones are O_M -H…O_M HBs (6.5 to 11.6 kcal mol^-1 ). The O_W -H…O_W (0.2 to 10.9 kcal mol^-1 ) and O_W -H…O_M HBs (0.3 to 10.3 kcal mol^-1 ) are the weakest ones. This energetic ordering of HBs is seen to be different from the respective HB energies in the dimer i. e., O_M -H…O_M (5.0 to 6.0 kcal mol^-1 )>O_W -H…O_M (1.5 to 6.0 kcal mol^-1 )>O_M -H…O_W (3.8 to 5.6 kcal mol^-1 )>O_W -H…O_W (1.2 to 5.0 kcal mol^-1 ). The plausible reason for the difference in the HB energy ordering may be attributed to the increase or decrease in HB strengths due to the formation of cooperative or anti-cooperative HB networks. For instance, the cooperativity contribution towards the different types of HB follows: O_M -H…O_W (2.4 to 8.6 kcal mol^-1 )>O_M -H…O_M (1.3 to 6.3 kcal mol^-1 )>O_W -H…O_W (-1.0 to 6.5 kcal mol^-1 )>O_W -H…O_M (-1.2 to 5.3 kcal mol^-1 ). This ordering of cooperativity contribution is similar to the HB energy ordering obtained by the MTA-based method. It is emphasized here that, the interplay between the cooperative and anti-cooperative contributions are indispensable for the correct energetic ordering of these HBs.

Assessment of hydrogen bond strengths and cooperativity in self- and cross-associating cyclic (HF)m(H2O)n (m + n = 2 to 8) clusters

In this work, we investigated the strengths of various self- and cross-associating hydrogen bonds (HBs) in mixed hydrogen fluoride–water cyclic (HF)m(H2O)n (m + n = 2 to 8) clusters, employing a molecular tailoring approach (MTA)-based method.

Unusually Large Hydrogen-Bond Cooperativity in Hydrogen Fluoride Clusters, (HF)n, n = 3 to 8, Revealed by the Molecular Tailoring Approach

In this work, our recently proposed molecular tailoring approach (MTA)-based method is employed for the evaluation of individual hydrogen-bond (HB) energies in linear (L) and cyclic (C) hydrogen fluoride clusters, (HF)_n (n = 3 to 8). The estimated individual HB energies calculated at the MP2(full)/aug-cc-pVTZ level for the L-(HF)_n are between 6.2 to 9.5 kcal/mol and those in the C-(HF)_n lie between 7.9 to 11.4 kcal/mol. The zero-point energy corrections and basis set superposition corrections are found to be very small (less than 0.6 and 1.2 kcal/mol, respectively). The cooperativity contribution toward individual HBs is seen to fall between 1.0 to 4.8 kcal/mol and 3.2 to 6.9 kcal/mol for linear and cyclic clusters, respectively. Interestingly, the HB energies in dimers, cleaved from these clusters, lie in a narrow range (4.4 to 5.2 kcal/mol) suggesting that the large HB strength in (HF)_n clusters is mainly due to the large cooperativity contribution, especially for n ≥ 5 (50 to 62% of the HBs energy). Furthermore, the HB energies in these clusters show a good qualitative correlation with geometrical parameters (H···F distance and F-H···F angles), stretching frequencies of F-H bonds, and electron density values at the (3, -1) bond critical points.

Quantum Stabilization of the Frustrated Hydrogen Bonding Structure in the Hydrogen Fluoride Trimer

We performed ab initio path integral molecular dynamics (PIMD) and molecular dynamics (MD) simulations to discuss the thermal and nuclear quantum effects on the stabilities of hydrogen bonding network in a hydrogen fluoride trimer (HF)₃ cluster. By the conventional molecular orbital calculation, the (HF)₃ cluster has an equilateral triangle shape, which has a frustration in the chemical structure of the hydrogen bonds, whereas the hydrogen bonding structure of a hydrogen fluoride dimer (HF)₂ cluster is nearly perpendicular to the acceptor molecule. The ratio of the triangular structures with the three hydrogen bondings in the PIMD simulation is larger than that in the MD one, whereas nonhydrogen bonding conformations such as a dimerlike structure are often found in MD simulation. The nuclear quantum effect stabilizes the frustrated hydrogen bonding network of the triangular (HF)₃ cluster.

Analysis of the basis set superposition error in molecular dynamics of hydrogen-bonded liquids: application to methanol

An efficient protocol is presented to compensate for the basis set superposition error (BSSE) in DFT molecular dynamics (MD) simulations using localized Gaussian basis sets. We propose a classical correction term that can be added a posteriori to account for BSSE. It is tested to what extension this term will improve radial distribution functions (RDFs). The proposed term is pairwise between certain atoms in different molecules and was calibrated by fitting reference BSSE data points computed with the counterpoise method. It is verified that the proposed exponential decaying functional form of the model is valid. This work focuses on hydrogen-bonded liquids, i.e., methanol, and more specific on the intermolecular hydrogen bond, but in principle the method is generally applicable on any type of interaction where BSSE is significant. We evaluated the relative importance of the Grimme-dispersion versus BSSE and found that they are of the same order of magnitude, but with an opposite sign. Upon introduction of the correction, the relevant RDFs, obtained from MD, have amplitudes equal to experiment.

Quantum nature of the hydrogen bond

Hydrogen bonds are weak, generally intermolecular bonds, which hold much of soft matter together as well as the condensed phases of water, network liquids, and many ferroelectric crystals. The small mass of hydrogen means that they are inherently quantum mechanical in nature, and effects such as zero-point motion and tunneling must be considered, though all too often these effects are not considered. As a prominent example, a clear picture for the impact of quantum nuclear effects on the strength of hydrogen bonds and consequently the structure of hydrogen bonded systems is still absent. Here, we report ab initio path integral molecular dynamics studies on the quantum nature of the hydrogen bond. Through a systematic examination of a wide range of hydrogen bonded systems we show that quantum nuclear effects weaken weak hydrogen bonds but strengthen relatively strong ones. This simple correlation arises from a competition between anharmonic intermolecular bond bending and intramolecular bond stretching. A simple rule of thumb is provided that enables predictions to be made for hydrogen bonded materials in general with merely classical knowledge (such as hydrogen bond strength or hydrogen bond length). Our work rationalizes the influence of quantum nuclear effects, which can result in either weakening or strengthening of the hydrogen bonds, and the corresponding structures, across a broad range of hydrogen bonded materials. Furthermore, it highlights the need to allow flexible molecules when anharmonic potentials are used in force field-based studies of quantum nuclear effects.

Properties of liquid water from a systematic refinement of a high-rank multipolar electrostatic potential

We build on previous work [S. Y. Liem and P. L. A. Popelier, J. Chem. Theory Comput. 4, 353 (2008)], where for the first time, a high-rank multipolar electrostatic potential was used in molecular dynamics simulations of liquid water at a wide range of pressures and temperatures, and using a multipolar Ewald summation. Water is represented as a rigid body, with atomic multipole moments defined by quantum chemical topology partitioning its gas phase electron density. The effect of the level of theory on the local structure of liquid water is systematically addressed. Values for Lennard-Jones (LJ) parameters are optimized, for both oxygen and hydrogen atoms, against bulk properties. The best LJ parameters were then used in a set of simulations at 30 different temperatures (1 atm) and another set at 11 different pressures (at 298 K). Inclusion of the hydrogen LJ parameters significantly increases the self-diffusion coefficient. The behavior of bulk properties was studied and the local water structure analyzed by both radial and spatial distribution functions. Comparisons with familiar point-charge potentials, such as TIP3P, TIP4P, TIP5P, and simple point charge, show the benefits of multipole moments.