Quantum Signature of Anisotropic Singularities in Hydrogen Bond Breaking of Water Dimer


 Molecular simulations from small molecules to large bio-macromolecules and polymer systems are routinely used to simulate thermodynamics properties of interests by molecular mechanics-based potentials. In a recent paper, via three different semi-empirical methods, we reported quantum singularities in molecular mechanics torsion potentials as signature of chemical bond break-up process revealed under experimental X-ray as broken chemical moieties. In this present work, applying first principle methods of Hartree-Fock, Density Functional as well as Moller-Plesset techniques we have reconfirmed the previous general predictions of singularities in the torsion potential for the case of water dimer that connects two water monomers by weak hydrogen bond. Due to quantum nature of chemical bond breaking process leading to break-point conditions in otherwise connected molecular topology aided by molecular mechanics based potentials, the singularities are also suggestive of large forces as onset in the bond-breaking process. We have presented the details of these novel interesting findings in this paper. These results of quantum singularities can have significant impacts to improve current force fields and can open up new areas we define as “Fracture Molecular Mechanics” or “Fracture Force Field” in overlap regions of molecular and quantum mechanics based approaches to explore and account for chemical bond-breaking mechanisms in molecular simulation techniques.


Introduction
In our recent paper we have reported for the rst time quantum singularities in molecular mechanics torsion potential as signature of chemical bond break-up process revealed in the X-ray structure 1 , 2 . In this present work, via several different ab initio techniques like Hartree-Fock 3 , Density Functional 4 as well as Moller-Plesset 5 methods and with higher basis set we have computed the torsion potential of dissociative water dimer and have consistently observed our previous general predictions of singularities in breaking chemical bonds. Due to the quantum nature of chemical bond break-up process leading to break-point conditions in connected molecular topology or in other expression the singularities in molecular mechanics based potential energy are simultaneously also indicative of large classical force at onset of bond-breaking by torsion of a molecule. Therefore, the effect can be used to enhance current molecular simulation techniques. Although purely quantum mechanical phenomenon, due to weak spectroscopic feature, tracking individual chemical bond break-up experimentally in dissociative or other chemical reaction mechanisms can be challenging both in gas or solution phase by traditional spectroscopies 6 . Spectroscopic detection of such torsion singularities will require advanced state-of-art single molecule based force or torque spectroscopy experiments 7 , 8 , 9 , 10 .
Development of physics-based force eld potentials for molecular modeling of bio-macromolecules and large polymer systems have their origin in the rst successful molecular dynamics (MD) simulation of liquid argon carried out by implementing the Lennard Jones potential 11 , 12 . Based on the Born-Oppenheimer theory via assigning an atomic charge at nuclear co-ordinates with van der Waals radii as soft excluded volume, the MD scheme uses Boltzmann equilibrium assembly properties in correlating ensemble averages with macroscopic properties 13 . In MD simulation force eld categories of both BONDED and NON-BONDED terms do connect the atoms by a topology of network and in practice NON-BONDED terms practically link the atoms at longer geometrical distance. In molecular simulations the terms are usually truncated to 10-12 Ȧ in order to save processor time. However, all the force eld terms including traditional non-bonded terms are basically a representation of connectivity of atoms in a bonded network with smooth potential as MD inherently does not simulate any phenomena that is based on electronic orbital structure both for nearest or long distance atomic neighbors. However, as routine methods via using high-performance computing, MD simulations have signi cant successes in prediction of X-ray, Cryo-EM structure within reasonable atomic accuracy 14 , 15 , 16 . The connectivity of atoms by molecular mechanics potential reduces accessible phase space particularly to capture critical phenomena when molecular potential might be rede ned in critical transition process; in recent years efforts have been applied in enhanced sampling methods with adaptive polarization and energetics in the assembly 17 , 18 . Another simulation eld used the Monte Carlo (MC) technique mainly based on the pioneering work of Metropolis et al. 19 so far has reported and reproduced some feature following MD trait 20 , 21 . However, all reported classical MD/MC methods in bio-molecular simulation elds so far have been based on connected molecular topology with distinctive phase space sampling technique with no categorical difference in molecular mechanics atomic connectivity based on bond breaking or making information by wave mechanics techniques. As also with the case with ab initio MD where re ned forces are computed quantum mechanically rather than using molecular mechanics for enhanced phase space sampling but does fall short to capture actual quantized singularity effects in bond-breaking or bondmaking 22 . Due to above inherent limitations, it is not surprising that normal-mode analysis in MD trajectory analysis does not have the expected success in correlating experimental spectroscopy compared to the ab initio methods 23 , 24 .
Since its discovery water dimer 25 has been an extensively studied chemical system by various experimental and theoretical techniques over the last few decades 26 , 27 , 28 , 29 . It is the simplest form in water cluster formation as has been con rmed by subsequent numerous theoretical and experimental investigations 30 , 31 , 32 . Since its rst reporting more than 60 years ago the ab initio based studies have been carried out to study and correlate water dimer geometry, dynamics and spectroscopic feature by several methods 33 , 34 . In this current report we have repeated and reproduced some of the standard ab initio water dimer results previously reported related to optimized geometry and vibrational spectra correlated with experimental methods 35 , 36 . Finally results of H-bond breaking as quantum signature in torsion potential have been presented as further conclusive evidence previously reported with semiempirical techniques. In the following sections we have discussed the methods for computational techniques, discussion, results analysis and conclusions as well as some quantum computational results with scopes for new experiment design based on reported singularity feature as logical sequence of the theoretical observations reported in this work.

Computational Methodologies
Ground State Geometry of Water Dimer Ground state geometry of water dimer has been computed by several rst principle techniques as implemented in Gaussian 16 37 and SPARTAN 18 platform (URL: https://www.wavefun.com; Wavefunction, INC., Irvine, CA). Geometry optimization, bond length, bond angle, HOMO-LUMO orbital, energy gap difference, electrostatic potential map of water dimer have been computed by Moller-Plesset (MP2) with aug-ccPVDz basis, Hartree-Fock (HF) with 6-31G* basis, Density Functional Theory (DFT) at theory level wB97X and B3LYP with 6-31G* basis in gas phase, water and polar medium. DFT with wB97X is a relatively new functional form compared to B3LYP 38 that uses a long-range corrected hybrid density functional with atom-atom damped dispersion correction 39 . Besides DFT, post Hartree-Fock ab initio technique with two types of basis set HF/6-31G* 40 , 41 and HF/3-21G have been applied for comparison of optimized geometry in gas phase, water and polar medium 42 , 43 .

Water Dimer Binding Energy and Vibrational Spectra
Hydrogen bonds are represented by the well-known 12-10 potential 44 for bond length range between 1.65 to 3.00 Å and acceptor-donor angle 90° < θ< 180°.
Theoretical binding energy of water dimer formation was estimated from the hydrogen bond well depth between O4 and H2 ( Fig.1) atoms linking water monomers. Likewise for geometry optimization several rst principle methods have been applied to estimate the binding energy that has reported a value of about 5.2 kcal/mol without zero point energy correction 45 , 46 , 47 . As shown in Table 4, the binding energy value obtained by the MP2 technique is in good consistency with reported values after basis-set superposition corrections 48 . Each energy scan for a speci c method has been carried out with the structure initially obtained from geometry optimization. The spectral calculation procedure as mentioned in above section and discussed below used the same rst principle technique. The scan was performed with O4-H2 distance as variable at 0.1 Å resolution ranging from 1.5 Å to 10 Å as shown in Figs 3 and 4.
The dissociation energy has been estimated from the well depth from global minima and the asymptotic unbounded energy threshold computed from the potential plot for each technique. The above method was repeated with O1-O4 49 , 50 interaction co-ordinates as well to compare dimer binding energy results from H-bonded O4-H2 separation distance method. The results have been shown in Table. 4.
For vibrational spectral studies several methods with different basis sets have been used for benchmarking and best reproducibility of reported theoretical and experimental spectra 51 , 52 , 53 . Infrared vibration peaks were assigned to corresponding different vibrational modes of water dimer atomic vibration groups 54 . Computational vibrational spectra computed by Gaussian 16 have also been reported in the supplementary section. For Spartan 18 based calculations, the systematic error primarily due to the harmonic approximation is corrected by uniformly scaling the amplitude and spectra lines are broadened to account for nite temperature due primarily to rotational structure.

Torsion-Dependent H-bond Energy by MM and QM
Torsion potential energy for dissociative water dimer molecule around O4-H2 bond (Figs 9) has been calculated by three separate ab initio methods. A set of three representative ab initio techniques of Hartree-Fock (basis set HF/6-31G*), Moller Plesset (basis MP2/aug-ccPDVZ) and Density Functional Theory (basis set B3YLP/6-31G*) have been applied to map out the torsional energy in the entire conformation span for O4-H2 hydrogen bonded atoms as well as O4-O1 interaction atoms forming reaction co-ordinates also shown in Figs 1 and 8. Torsion potential energy is given by following equation (symbols have their standard meanings).
Prior to the energy scan each structure was geometry optimized via the corresponding technique for ground state parameters and vibrational spectra as reported below. Dihedral energy scans with two different torsion planes of H5O4H2O1 (O4-H2 hydrogen bonded co-ordinates) and H5O4O1H3 (O4-O1 interaction co-ordinates) have been carried out for entire conformational span of 0o to 360o with 1o resolution. For comparison of the torsion potential with molecular mechanics potential, the process for both cases were repeated with the Merck Force Field popularly known as MMFF 55 .

Anisotropic Singularites in Torsion Dependent H-bond Energy Surface
Anisotropic torsion energy of water dimer has been explored both in gas phase and in water medium by semi-empirical AM1 and PM3 techniques which previously we reported can record quantum chemical signature of bond break-up mechanisms. As seen in Fig 9, with xed O4-H2 or O4-O1 lengths and varied torsion angle for the entire torsional space 0o to 360o, nature of energetics of water dimer can be explored in the critical range where H-bond interaction should be weak to hold two water monomers as integral chemical entity. Since both O1 and O4 atoms are heavier compared to the four hydrogen atoms (H2, H3, H5 and H6), although not bonded by Lewis de nition the interaction distance parameters O1-O4 can serve in good approximation as center-of-mass distance between two water monomers as has been pointed out in respect to assessing water dimer binding energy. Any singularities in energetics of such reaction co-ordinates in (r, f) does imply break-up of molecular mechanics based connected topology of atoms in otherwise chemicaly stable molecules. For a xed r, that is either O4-H2 or O1-O4 atomic distance; torsion angle f has been scanned for entire conformational space from 0o to 360o at 1o

Optimized Water Dimer Geometry in Gas Phase, Water and Polar Medium
As presented in Tables 1,2, 3 and Supplementary Table 1, prior to vibrational spectra calculations, water dimer geometry was optimized by several ab initio methods in gas phase, water and polar medium. The results for geometric parameters of bond length, bond angle in gas phase are shown in Table 1. For each method subsequent vibrational spectra computations have been performed with the optimized geometry as the starting structure as discussed in next section. As observed below for optimized geometry parameters, all three applied techniques of DFT, MP2 and HF with basis set of 6-31G*, 6-311+G** and 6-311+G(2df, 2p) are in general in good consistency as carried out under Spartan 18 software tools with earlier reported works 59 , 60 . Also in the Supplementary Table 1, some of the calculations presented have been repeated and reproduced by Gaussian 16 with similar theory level and basis set and all are in good agreement. The optimized geometry has been also assessed in water and polar medium shown in Table   2    The above Table 5. shows results for vibrational spectra mode calculations of water dimer. The OHvibration mode is split into symmetric and asymmetric stretch which has distinctively also split to donor and acceptor mode as observed in The hydrogen bond linking two water molecules as reported and computed above has weak binding energy of around 5.2 kcal/mol without zero point energy corrections. The range of energy makes water dimer structure very unstable under any bond torsion motion. That dissociative bond-breaking tendency makes the water dimer an ideal system to explore for the quantum signature in torsion potential that we reported earlier with semi-empirical techniques. For the dihedral H5O4H2O1 plane in Fig 1 and Fig 9, dihedral rotation does strain the weak O4H2 hydrogen bond and reaches critical conformation that breaks the water dimer structure eventually. As observed in each case of HF/6-31G*, B3LYP/ B97X-D/6-31G*, MP2/aug-cc/PVDZ based calculations in Figs 11, 12 and 13, except for narrow isolated angle ranges of discontinuities, water dimer H bond has mostly isotropic at energy feature as revealed under with varied torsion angle scans. This common feature in rst principle calculations is quite in contrast to molecular mechanics based torsion potentials computed by MMFF method shown in Fig. 10. In contrast to QM potential as reconstructed from the orbital basis set, MM potential curve does not show any singularities in torsion energy curve as also continuous rst derivative has been observed in entire the conformation space in contrast to each of the ab initio case for exploring weak H bond in water dimer. The amount of energy transition of around ~ 5.5 eV is well above H-bond binding energy of 0.22 eV or 5.2 Kcal/mol. The values are clear indicative of bond-breaking phenomena under torsion for weak H bonds in water dimer as the jumps have large energetic feature for its dissociation. Based on that observations, ab initio based torsion potential can be applied as good exploratory probe to theoretically check molecular stability and estimate of bond incision at different locations of large molecules provided potential pro le is reconstructed from quantum functional basis. Reproducibility results by Gaussian 16 also by rst principle techniques as seen in Fig. 17 and 18 also have shown instability of water dimer for torsional motion around O4-O1 reaction co-ordinates. The estimates of torsion energy slopes in bond-breaking region by B3YLP/6-31G+ and wBxd4/6-31G+ methods are respectively 2x10-3 and 8.5x10- 3 Hartree/Bohr which are equivalent to 0. 16

Conclusions
The scheme of connected topology of atoms by molecular mechanics potentials, although has inherent limitations to underpin bond-breaking feature via quantum mechanics, still has advantages to de ne a collection of atoms to de ne molecular entity from small organic molecules to large macro-molecules, like protein, DNA/RNA, carbohydrates, lipids, polymers etc. Although classical force eld theory allows explicit bonding between neighboring atoms and by a networks of those linkages; quantum mechanics does not eliminate any long distance instantaneous interaction between atoms not bonded via Lewis de nition of proximity 76 , 77 . Given the quantum nature of non-localization of electrons it will not be surprising if distant electrons in a large macro-molecule are effectively coupled as Born-Oppenheimer centers affecting observable macro-molecular properties in rather detectable ways. In addressing such problems, the bottlenecks are how big a system in number of atoms current state-of-art quantum tools can handle so that some new types of reaction co-ordinates can be de ned over traditional methods.
Towards that goal we have de ned a new area of molecular simulation studies by "FRACTURE MOLECULAR MECHANICS" or "FRACTURE FORCE FIELD" that can be an interesting overlap of molecular mechanics and quantum mechanics incorporating chemical bond making and breaking processes in time scale based on electronic orbital properties rather than the current practices of so called hybrid applications of MM and QM in partitioned regions in molecules. Such quantum assembly of atoms and molecules in different physics conditions can capture many spontaneous reaction processes compared to current available state-of-art molecular simulation tools mostly evolved from the paradigm of classical molecular dynamics. Obviously, such simulation results will needed to be rigorously bench-marked against reaction-based fast dynamic laser spectroscopy. While the above should be the enhancement of current state-of-art force elds that can potentially mimic the quantum dynamic nature of bond breaking and making events with time scale with collection of large number of molecules, bio-macromolecules or polymers in the scheme of Fracture Molecular Mechanics/Fracture Force Field, in this paper we have underlined the expected outcomes for single molecules torsion or torque spectroscopy experiments with the water dimer. In above section, we have reported the expected values for bond rupture force in H-bond breaking of water dimer that in fact exceeds the current reported experimental values to actuate bond breaking process. In future publications we expect to report further generalized experimental aspects of the results.