A generalized reactive force field for nonlinear hydrogen bonds: hydrogen dynamics and transfer in malonaldehyde

A Bond-Based Machine Learning Model for Molecular Polarizabilities and A Priori Raman Spectra

The use of machine learning (ML) algorithms in molecular simulations has become commonplace in recent years. There now exists, for instance, a multitude of ML force field algorithms that have enabled simulations approaching ab initio level accuracy at time scales and system sizes that significantly exceed what is otherwise possible with traditional methods. Far fewer algorithms exist for predicting rotationally equivariant, tensorial properties such as the electric polarizability. Here, we introduce a kernel ridge regression algorithm for machine learning of the polarizability tensor. This algorithm is based on the bond polarizability model and allows prediction of the tensor components at the cost similar to that of scalar quantities. We subsequently show the utility of this algorithm by simulating gas phase Raman spectra of biphenyl and malonaldehyde using classical molecular dynamics simulations of these systems performed with the recently developed MACE-OFF23 potential. The calculated spectra are shown to agree very well with the experiments and thus confirm the expediency of our algorithm as well as the accuracy of the used force field. More generally, this work demonstrates the potential of physics-informed approaches to yield simple yet effective machine learning algorithms for molecular properties.

Investigating proton tunneling dynamics in the time-dependent Schrödinger equation

Understanding the temporal evolution of the wave function in the proton tunneling reactions allows us to make theoretical predictions on the possible femtosecond spectroscopy patterns. However, the analytical solution of the time-dependent Schrödinger equation of complex molecules is challenging and requires a high computational cost. In this study, we solve the time-dependent Schrödinger equation using the Fourier grid Hamiltonian method, highlighting its simplicity of calculation, even for multidimensional tunneling reactions. The obtained model is applied to studying malonaldehyde's two-dimensional intramolecular proton tunneling, comparing the results with those obtained using other computational methods.

Intramolecular Proton Transfer in the Hydrogen Oxalate Anion and the Cooperativity Effects of the Low-Frequency Vibrations: A Driven Molecular Dynamics Study

We report first-principles molecular dynamics (MD) and dipole-driven molecular dynamics (μ-DMD) simulations of the hydrogen oxalate anion at the MP2/aug-cc-pVDZ level of theory. We examine the role of vibrational coupling between the OH stretching bands, that is, the fundamental and a few combination bands spanning the 2900-3100 cm^-1 range, and several of the low-frequency bending and stretching fundamental modes. The low-frequency modes between 300 and 825 cm^-1 play a crucial role in the proton-transfer motion. Strong involvement of CO₂ and CCO bending and the CC stretching vibrations indicate that these large amplitude motions cause the shortening of the O···O distance and thus promote H⁺ transfer to the other oxygen by bringing it over the 3.4 kcal/mol barrier. Analysis of resonant μ-DMD trajectories shows that the complex spectral feature near 825 cm^-1, closely corresponding to both an overtone of two quanta of 425 cm^-1 and a combination band of low-frequency CO₂ rocking (300 cm^-1) and CCO bending (575 cm^-1) modes, is involved in the proton transfer. μ-DMD shows that exciting the system at these mode combinations leads to faster barrier activation than exciting at the OH fundamental mode.

A simplistic computational procedure for tunneling splittings caused by proton transfer

In this manuscript, we present an approach for computing tunneling splittings for large amplitude motions. The core of the approach is a solution of an effective one-dimensional Schrödinger equation with an effective mass and an effective potential energy surface composed of electronic and harmonic zero-point vibrational energies of small amplitude motions in the molecule. The method has been shown to work in cases of three model motions: nitrogen inversion in ammonia, single proton transfer in malonaldehyde, and double proton transfer in the formic acid dimer. In the current work, we also investigate the performance of different DFT and post-Hartree–Fock methods for prediction of the proton transfer tunneling splittings, quality of the effective Schrödinger equation parameters upon the isotopic substitution, and possibility of a complete basis set (CBS) extrapolation for the resulting tunneling splittings.

Impact of the Characteristics of Quantum Chemical Databases on Machine Learning Prediction of Tautomerization Energies

An essential aspect for adequate predictions of chemical properties by machine learning models is the database used for training them. However, studies that analyze how the content and structure of the databases used for training impact the prediction quality are scarce. In this work, we analyze and quantify the relationships learned by a machine learning model (Neural Network) trained on five different reference databases (QM9, PC9, ANI-1E, ANI-1, and ANI-1x) to predict tautomerization energies from molecules in Tautobase. For this, characteristics such as the number of heavy atoms in a molecule, number of atoms of a given element, bond composition, or initial geometry on the quality of the predictions are considered. The results indicate that training on a chemically diverse database is crucial for obtaining good results and also that conformational sampling can partly compensate for limited coverage of chemical diversity. The overall best-performing reference database (ANI-1x) performs on average by 1 kcal/mol better than PC9, which, however, contains about 2 orders of magnitude fewer reference structures. On the other hand, PC9 is chemically more diverse by a factor of ∼5 as quantified by the number of atom-in-molecule-based fragments (amons) it contains compared with the ANI family of databases. A quantitative measure for deficiencies is the Kullback-Leibler divergence between reference and target distributions. It is explicitly demonstrated that when certain types of bonds need to be covered in the target database (Tautobase) but are undersampled in the reference databases, the resulting predictions are poor. Examples of this include the poor performance of all databases analyzed to predict C(sp²)-C(sp²) double bonds close to heteroatoms and azoles containing N-N and N-O bonds. Analysis of the results with a Tree MAP algorithm provides deeper understanding of specific deficiencies in predicting tautomerization energies by the reference datasets due to inadequate coverage of chemical space. Capitalizing on this information can be used to either improve existing databases or generate new databases of sufficient diversity for a range of machine learning (ML) applications in chemistry.

Vibrational band-structures caused by internal rotations of the boron Wankel rotor B11. RSC Adv 2021;11:3613-3621. [PMID: 35424268 PMCID: PMC8694174 DOI: 10.1039/d0ra08821h]

Nuclear quantum effects are often neglected for systems without hydrogen atoms. However some planar boron rotors turn out to exhibit remarkable nuclear quantum effects. Recent experiment on infrared spectroscopy of B₁₃⁺ shows unexpected spectral broadening which still awaits physical explanation. Here we present quantitative investigations of the vibrational energy levels of B₁₁⁻ up to full dimension. A harmonic-bath averaged Hamiltonian suitable for planar boron rotors is constructed and used to predict typical types of vibrational states of B₁₁⁻. Band structures caused by internal rotations are found for all the investigated vibrational states. The experimental phenomenon of spectral broadening is thus due to the band structures of the corresponding vibrational levels. The detailed information of the relevant vibrational states reported in the present work may provide valuable references for future investigations of high resolution spectroscopy of B₁₁⁻.

The band structures of the vibrational energy levels of B₁₁⁻ lead to corresponding spectral broadening. The vibrational band-structures of planar boron rotors are caused by internal rotations.

Reactive Molecular Dynamics at Constant Pressure via Nonreactive Force Fields: Extending the Empirical Valence Bond Method to the Isothermal-Isobaric Ensemble

The Empirical Valence Bond (EVB) method offers a suitable framework to obtain reactive potentials through the coupling of nonreactive force fields. In this formalism, most of the implemented coupling terms are built using functional forms that depend on spatial coordinates, while parameters are fitted against reference data to model the change of chemistry between the participating nonreactive states. In this work, we demonstrate that the use of such coupling terms precludes the computation of the stress tensor for condensed phase systems and prevents the possibility to carry out EVB molecular dynamics in the isothermal-isobaric (NPT) ensemble. Alternatively, we make use of coupling terms that depend on the energy gaps, defined as the energy differences between the participating nonreactive force fields, and derive a general expression for the EVB stress tensor suitable for computation. Implementation of this new methodology is tested for a model of a single reactive malonaldehyde solvated in nonreactive water. Mass densities and probability distributions for the values of the energy gaps computed in the NPT ensemble reveal a negligible role of the reactive potential in the limit of low concentrated solutions, thus corroborating for the first time the validity of approximations based on the canonical NVT ensemble, customarily adopted for EVB simulations. The presented formalism also aims to contribute to future implementations and extensions of the EVB method to research the limit of highly concentrated solutions.

Observation of photoassociation spectroscopy of ultralong 37D5/2 + 6S1/2Cs2 Rydberg molecules

We present an experimental observation of 37D_5/2 + 6S_1/2Cs₂ Rydberg-ground molecules by employing a two-photon photoassociation method. Two distinct Rydberg-ground molecular signals, deep and shallow bound molecules, are observed at the red detuning of atomic line. In theory, the model of scattering interaction between the Rydberg electron and ground-state atom is used to simulate the experiments. Two potential energy curves with energy minimum, deep pure triplet ³Σ and shallow hyperfine-mixed singlet-triplet ^1,3Σ potentials, refer to the attained Rydberg-ground molecular signals, respectively. Calculations of the binding energy of triplet ³Σ and mixed ^1,3Σv = 0 states are compared with the measurements. The agreement between the calculated and measured values of the binding energy yields zero-energy scattering lengths a_s ^T(0) = -19.2a₀ and a_s ^S(0) = -1.3a₀, respectively.

Multistate Reactive Molecular Dynamics Simulations of Proton Diffusion in Water Clusters and in the Bulk

The molecular mechanics with proton transfer (MMPT) force field is combined with multistate adiabatic reactive molecular dynamics (MS-ARMD) to describe proton transport in the condensed phase. Parametrization for small protonated water clusters based on electronic structure calculations at the MP2/6-311+G(2d,2p) level of theory and refinement by comparing with infrared spectra for a protonated water tetramer yields a force field which faithfully describes the minimum energy structures of small protonated water clusters. In protonated water clusters up to (H₂O)₁₀₀H⁺, the proton hopping rate is around 100 hops/ns. This rate converges for 21 ≤ n ≤ 31, and no further speedup in bulk water is found. This indicates that bulklike behavior requires the solvation of a Zundel motif by ∼25 water molecules, which corresponds to the second solvation sphere. For smaller cluster sizes, the number of available states (i.e., the number of proton acceptors) is too small and slows down proton-transfer rates. The cluster simulations confirm that the excess proton is typically located on the surface. The free-energy surface as a function of the weights of the two lowest states and a configurational parameter suggests that the "special pair" plays a central role in rapid proton transport. The barriers between this minimum-energy structure and the Zundel and Eigen minima are sufficiently low (∼1 kcal/mol, consistent with recent experiments and commensurate with a hopping rate of ∼100/ns or 1 every 10 ps), leading to a highly dynamic environment. These findings are also consistent with recent experiments which find that Zundel-type hydration geometries are prevalent in bulk water.

Infrared Spectra of Deprotonated Dicarboxylic Acids: IRMPD Spectroscopy and Empirical Valence‐Bond Modeling

Experimental infrared multiple-photon dissociation (IRMPD) spectra recorded for a series of deprotonated dicarboxylic acids, HO₂ (CH₂ )_n CO 2 - (n=2-4), are interpreted using a variety of computational methods. The broad bands centered near 1600 cm^-1 can be reproduced neither by static vibrational calculations based on quantum chemistry nor by a dynamical description of individual structures using the many-body polarizable AMOEBA force field, strongly suggesting that these molecules experience dynamical proton sharing between the two carboxylic ends. To confirm this assumption, AMOEBA was combined with a two-state empirical valence-bond (EVB) model to allow for proton transfer in classical molecular dynamics simulations. Upon suitable parametrization based on ab initio reference data, the EVB-AMOEBA model satisfactorily reproduces the experimental infrared spectra, and the finite temperature dynamics reveals a significant amount of proton sharing in such systems.

A reactive, scalable, and transferable model for molecular energies from a neural network approach based on local information

Despite the ever-increasing computer power, accurate ab initio calculations for large systems (thousands to millions of atoms) remain infeasible. Instead, approximate empirical energy functions are used. Most current approaches are either transferable between different chemical systems, but not particularly accurate, or they are fine-tuned to a specific application. In this work, a data-driven method to construct a potential energy surface based on neural networks is presented. Since the total energy is decomposed into local atomic contributions, the evaluation is easily parallelizable and scales linearly with system size. With prediction errors below 0.5 kcal mol^-1 for both unknown molecules and configurations, the method is accurate across chemical and configurational space, which is demonstrated by applying it to datasets from nonreactive and reactive molecular dynamics simulations and a diverse database of equilibrium structures. The possibility to use small molecules as reference data to predict larger structures is also explored. Since the descriptor only uses local information, high-level ab initio methods, which are computationally too expensive for large molecules, become feasible for generating the necessary reference data used to train the neural network.

Spectroscopy of Proton Coordination with Ethylenediamine

Protonated ethylenediamine monomer, dimer, and trimer were produced in the gas phase by an electrical discharge/supersonic expansion of argon seeded with ethylenediamine (C₂H₈N₂, en) vapor. Infrared spectra of these ions were measured in the region from 1000 to 4000 cm^-1 using laser photodissociation and argon tagging. Computations at the CBS-QB3 level were performed to explore possible isomers and understand the infrared spectra. The protonated monomer exhibits a gauche conformation and an intramolecular hydrogen bond. Its parallel shared proton vibration occurs as a broad band around 2785 cm^-1, despite the formally equivalent proton affinities of the two amino groups involved, which usually leads to low frequency bands. The barrier to intramolecular proton transfer is 2.2 kcal mol^-1 and does not vanish upon addition of the zero-point energy, unlike the related protonated ammonia dimer. The structure of the dimer is formed by chelation of the monomer's NH₃⁺ group, thereby localizing the excess proton and increasing the frequency of the intramolecular shared proton vibration to 3157 cm^-1. Other highly fluxional dimer structures with facile intermolecular proton transfer and concomitant structural reorganization were computed to lie within 2 kcal mol^-1 of the experimentally observed structure. The spectrum of the trimer is rather diffuse, and a clear assignment is not possible. However, an isomer with an intramolecular proton transfer like that of the monomer is most consistent with the experimental spectrum.

2-Chloromalonaldehyde, a model system of resonance-assisted hydrogen bonding: vibrational investigation

The chelated enol isomer of 2-chloromalonaldehyde (2-ClMA) is experimentally characterized for the first time by IR and Raman spectroscopies. The spectra are obtained by trapping the molecule in cryogenic matrices and analyzed with the assistance of theoretical calculations. Experiments were performed in argon, neon and para-hydrogen matrices. The results highlight puzzling matrix effects, beyond site effects, which are interpreted as due to a tunneling splitting of the vibrational levels related to the proton transfer along the internal hydrogen bond (IHB). 2-ClMA is thus one of the very few molecules in which the H tunneling has been observed in cryogenic matrices. The comparison with its parent molecule (malonaldehyde) shows experimentally and theoretically the weakening of the IHB upon chlorination, with a reduced cooperative effect in the resonance assisted hydrogen bond. In addition, the Cl substitution induces an important stabilization of two open enol conformers. These two open forms appear in the spectra of as-deposited samples, meaning that, in contrast with other well-studied molecules of the same family (β-dialdehydes and β-diketones), they are present in the gas phase at room temperature.

Kinetic isotope effects and how to describe them

We review several methods for computing kinetic isotope effects in chemical reactions including semiclassical and quantum instanton theory. These methods describe both the quantization of vibrational modes as well as tunneling and are applied to the ⋅H + H2 and ⋅H + CH4 reactions. The absolute rate constants computed with the semiclassical instanton method both using on-the-fly electronic structure calculations and fitted potential-energy surfaces are also compared directly with exact quantum dynamics results. The error inherent in the instanton approximation is found to be relatively small and similar in magnitude to that introduced by using fitted surfaces. The kinetic isotope effect computed by the quantum instanton is even more accurate, and although it is computationally more expensive, the efficiency can be improved by path-integral acceleration techniques. We also test a simple approach for designing potential-energy surfaces for the example of proton transfer in malonaldehyde. The tunneling splittings are computed, and although they are found to deviate from experimental results, the ratio of the splitting to that of an isotopically substituted form is in much better agreement. We discuss the strengths and limitations of the potential-energy surface and based on our findings suggest ways in which it can be improved.

Ultrafast dynamics induced by the interaction of molecules with electromagnetic fields: Several quantum, semiclassical, and classical approaches

Several strategies for simulating the ultrafast dynamics of molecules induced by interactions with electromagnetic fields are presented. After a brief overview of the theory of molecule-field interaction, we present several representative examples of quantum, semiclassical, and classical approaches to describe the ultrafast molecular dynamics, including the multiconfiguration time-dependent Hartree method, Bohmian dynamics, local control theory, semiclassical thawed Gaussian approximation, phase averaging, dephasing representation, molecular mechanics with proton transfer, and multipolar force fields. In addition to the general overview, some focus is given to the description of nuclear quantum effects and to the direct dynamics, in which the ab initio energies and forces acting on the nuclei are evaluated on the fly. Several practical applications, performed within the framework of the Swiss National Center of Competence in Research "Molecular Ultrafast Science and Technology," are presented: These include Bohmian dynamics description of the collision of H with H2, local control theory applied to the photoinduced ultrafast intramolecular proton transfer, semiclassical evaluation of vibrationally resolved electronic absorption, emission, photoelectron, and time-resolved stimulated emission spectra, infrared spectroscopy of H-bonding systems, and multipolar force fields applications in the condensed phase.

Vibrational Spectroscopy and Proton Transfer Dynamics in Protonated Oxalate

The dynamics and infrared spectroscopic signatures of proton transfer in protonated oxalate (p-Oxa) are studied using classical and quantum dynamics. The intermolecular interactions are described by a force field suitable to follow proton transfer. This allows to carry out multiple extended classical molecular dynamics (MD) and ring polymer MD simulations from which the infrared spectrum is determined. Simulations at 600 K sample the quantum mechanical ground state probability distribution and best reproduce the experimentally observed maximum absorption wavelength and part of the line shape. Comparison with the experimentally measured spectrum provides an estimate for the barrier height for proton transfer which can not be determined directly from experiment. A barrier of 4.2 kcal/mol is found to best reproduce the position and width of the infrared absorption of the transferring proton in p-Oxa and also leads to an infrared (IR) spectrum in good agreement with experiment for the deuterated species d-Oxa. A novel means to capture the two resonance forms of oxalate depending on the localization of the excess proton on either CO moiety is found to yield improved results for the spectroscopy in the framework region between 1000 and 2000 cm^-1.

Spectroscopy and dynamics of double proton transfer in formic acid dimer

We present the isolated gas phase infrared spectra of formic acid dimer, (HCOOH)2, and its deuterated counterpart formic-d acid, (DCOOH)2, at room temperature. The formic acid dimer spectrum was obtained by spectral subtraction of a spectrum of formic acid vapor recorded at low pressure from that recorded at a higher pressure. The spectra of formic acid vapor contain features from both formic acid monomer and formic acid dimer, but at low and high pressures of formic acid, the equilibrium is pushed towards the monomer and dimer, respectively. A similar approach was used for the formic-d acid dimer. Building on the previous development of the Molecular Mechanics with Proton Transfer (MMPT) force field for simulating proton transfer reactions, molecular dynamics (MD) simulations were carried out to interpret the experimental spectra in the OH-stretching region. Within the framework of MMPT, a combination of symmetric single and double minimum potential energy surfaces (PESs) provides a good description of the double proton transfer PES. In a next step, potential morphing together with electronic structure calculations at the B3LYP and MP2 level of theory was used to align the computed and experimentally observed spectral features in the OH-stretching region. From this analysis, a barrier for double proton transfer between 5 and 7 kcal mol(-1) was derived, which compares with a CCSD(T)/aug-cc-pVTZ calculated barrier of 7.9 kcal mol(-1). Such a combination of experimental and computational techniques for estimating barriers for proton transfer in gas phase systems is generic and holds promise for further improved PESs and energetics of these important systems. Additional MD simulations at the semi-empirical DFTB level of theory agree quite well for the center band position but underestimate the width of the OH-stretching band.

Theoretical study of the C-H/O-H stretching vibrations in malonaldehyde

IR and Raman spectra of the malonaldehyde molecule and its deuterated analogues were calculated in the B3LYP/cc-pVQZ approximation. Anharmonicity effects were taken into account both in the context of a standard model of the second order perturbation theory and by constructing the potential energy surfaces (PES) with a limited number of dimensions using the Cartesian coordinates of the hydroxyl hydrogen atom and the stretching coordinates of С-Н, C-D, O-H, and O-D bonds. It was shown that in each of the two equivalent forms of the molecule, besides the global minimum, an additional local minimum at the PES is formed with the energy more than 3,000 cm(-1) higher than the energy in the global minimum. Calculations carried out by constructing the 2D and 3D PESs indicate a high anharmonicity level and multiple manifestations of the stretching О-Н vibrations, despite the fact that the model used does not take into account the splitting of the ground-state and excited vibrational energy levels. In particular, the vibration with the frequency 3,258 cm(-1) may be associated with proton transfer to the region of a local minimum of energy. Comparing the results obtained with the experimental data presented in the literature allowed us to propose a new variant of bands assignments in IR and Raman spectra of the molecule in the spectral region 2,500-3,500 cm(-1).

Symmetry breaking in the axial symmetrical configurations of enolic propanedial, propanedithial, and propanediselenal: pseudo Jahn–Teller effect versus the resonance-assisted hydrogen bond theory

The origin of the symmetry breaking in the axial symmetrical configurations of enolic propanedial (1), propanedithial (2), and propanediselenal (3) have been investigated by means of time-dependence density functional theory and natural bond orbital interpretations. The results obtained at the quantum chemistry composite (G2MP2, CBS-QB3), ab initio molecular orbital (MP2/6-311++G**), and hybrid density functional theory (B3LYP/6-311++G**) levels of theory showed that the hydrogen-centered synchronous axial symmetrical (C2v) configurations of compounds 1–3 possessing the maximum π-electron delocalization within the M1=C2–C3=C4–M5–H6 keto-enol groups are less stable than their corresponding plane symmetrical (Cs) forms. Importantly, the symmetry breaking in the C2v configurations of the enol forms of compounds 1–3 to their corresponding plane symmetrical Cs configurations is due to the pseudo Jahn–Teller effect (PJTE) by mixing the ground A1 and excited B2 electronic states resulting in a PJT (A1 + B2) ⊗ b2 problem. We may expect that by the decrease of the energy gaps between reference states in the C2v forms that are involved in the PJTE decrease from compound 1 to compound 3, the PJT stabilization energy (PJTSE) may increase but the results obtained showed that the corresponding PJTSEs decrease. This fact can be justified by the increase of the electron delocalizations from the nonbonding orbitals of the C=M moieties to the antibonding orbitals of the H–M bonds, which leads to an increase of the π-electron delocalization within the M1=C2–C3=C4–M5–H6 keto-enol groups. In confrontation between the impacts of the resonance-assisted hydrogen bond and PJTE in the structural and configurational properties of compounds 1–3, PJTE has an overwhelming contribution and causes the symmetry breaking of the C2v configurations to their corresponding Cs forms. The correlations between the structural parameters, synchronicity indices, natural charges, PJTSEs, electron delocalizations, and the hardness of compounds 1–3 have been investigated.

Harmonic bath averaged Hamiltonian: an efficient tool to capture quantum effects of large systems

Starting from a reaction path Hamiltonian, a suitably reduced harmonic bath averaged Hamiltonian is derived by averaging over all the normal mode coordinates. Generalization of the harmonic bath averaged Hamiltonian to any dimensions are performed and the feasibility to use a linear reaction path/surface are investigated and discussed. By use of a harmonic bath averaged Hamiltonian, the tunneling splitting and proton transfer dynamics of malonaldehyde is briefly discussed and shows that the harmonic bath averaged Hamiltonian is an efficient tool to capture quantum effects in larger systems.