Dynamical strengthening of covalent and non-covalent molecular interactions by nuclear quantum effects at finite temperature

Nuclear quantum effects (NQE) tend to generate delocalized molecular dynamics due to the inclusion of the zero point energy and its coupling with the anharmonicities in interatomic interactions. Here, we present evidence that NQE often enhance electronic interactions and, in turn, can result in dynamical molecular stabilization at finite temperature. The underlying physical mechanism promoted by NQE depends on the particular interaction under consideration. First, the effective reduction of interatomic distances between functional groups within a molecule can enhance the n → π* interaction by increasing the overlap between molecular orbitals or by strengthening electrostatic interactions between neighboring charge densities. Second, NQE can localize methyl rotors by temporarily changing molecular bond orders and leading to the emergence of localized transient rotor states. Third, for noncovalent van der Waals interactions the strengthening comes from the increase of the polarizability given the expanded average interatomic distances induced by NQE. The implications of these boosted interactions include counterintuitive hydroxyl–hydroxyl bonding, hindered methyl rotor dynamics, and molecular stiffening which generates smoother free-energy surfaces. Our findings yield new insights into the versatile role of nuclear quantum fluctuations in molecules and materials.

The Dunning's correlation-consistent basis set is a well established basis, then not much detail has to be included in this sections to describe their accuracy and capabilities. On the other hand, numerically tabulated atomcentered orbitals (NAOs), as implemented the FHIaims package, are less known. NAOs basis sets in the FHIaims * sauceda@tu-berlin.de † klaus-robert.mueller@tu-berlin.de ‡ alexandre.tkatchenko@uni.lu package are constructed to be transferable and hierarchical basis sets to systematically reach basis set convergence limit (sub-meV-level accuracy on the total energy). The basis sets are constructed starting from a minimal basis set (minimal free-atom basis) and then systematically adding new energetically-favorable functions from a pool of basis functions (hydrogen-like, cation-like, and atom-like functions with a variable confinement potential). Then, the basis used in this work, named "really tight" in the code, consists in, for example the hydro-   Classical MD PIMD B) sGDML@{CCSD/cc-pVTZ} SUPPLEMENTARY FIGURE 1. Classical (MD) and path integral molecular dynamics (PIMD) simulations at room temperature of aspirin described by the sGDML@CCSD molecular force field using the basis set cc-pVDZ (A) and cc-pVTZ (B). The plots are projections of the dynamics to the two main degrees of freedom of aspirin: carboxyl and ester dihedral angles.

Supplementary Note 3. Calculation of En→π *
Natural Bond Orbital (NBO) second order perturbative energies E n→π * obtained from NBO 7.0 [7] calculations coupled with ORCA 4.1.2 [8,9] at CCSD/cc-pVDZ level of theory were taken as the stabilization energies due to n → π * interactions. Such interactions, which play an important role in molecular reactivity and conformation (for instance, the Bürgi-Dunitz trajectory [10] preferred during nucleophilic attacks at a carbonyl carbon), comprise delocalization of lone-pair electrons (n) of an electronegative atom into an empty π * -antibonding orbital of an aromatic ring or a carbonyl group [11][12][13]. and PIMD simulations using a sGDML model for aspirin trained on CCSD/cc-pVDZ and CCSD/cc-pVTZ reference calculations, which display statisticaly the same behaviour. This means that intramolecular interactions in small molecules are already statistically well captured with cc-pVDZ basis set.

Supplementary Note 5. Molecular dynamics settings
The molecular dynamics simulations were done using the i-PI code [14,15] with the sGDML force field interface [16]. The integration time-step was set to 0.2 fs using the NVT ensemble and the total simulation time was 0.5 ns. All the simulations were performed at room temperature.
The path integral molecular dynamics (PIMD) simulations were performed using a baseline of 16 beads, value which in most cases gives converged thermodynamic properties of molecules at room temperature [17][18][19]. Such convergence was corroborated by running shorter PIMD trajectories (100 ps instead of 0.5 or 1 ns) using 32 beads. In the particular case of the methyl rotor dynamics in toluene, given the unexpected localization, the convergence of the results was validated running simulations using up to 64 beads.

Supplementary Note 6. Methyl rotor analysis: Different levels of theory
As mentioned in the main text, the rotor localization does not depend on the level of theory used to describe the PES of the toluene rotor localization. Supplementary Figure 2 shows direct evidence of this statement. The most evident case is that Hartree-Fock (HF/augcc-pVQZ) displays the same trends as CCSD(T)/cc-pVDZ even though their rotational barriers are 0.0008 kcal mol −1 and 0.028 kcal mol −1 , respectively. The same results are obtained using DFT with PBE0 and PBE0+MBD. The many body dispersion (MBD) [20,21] treatment of the van der Waals interaction [22,23] increases the PBE's rotational barrier given its tight relationship to the correlation energy. Despite the differences in the four levels of theory used for the simulation, the methyl rotor hindering was proven to be level-of-theory agnostic.
All the PIMD simulations in Supplementary Figure 2 were done using 16 beads, but they were validated with shorter simulations with upto 64 beads.

Supplementary Note 7. Origin of the rotor energetic barrier
From the electronic structure point of view, previous studies suggested that [24][25][26][27] the Me rotational energy barrier E MeRot originates from the difference of the π-bond order between the two ring bonds C 2 −C 3 and g. temperature, thermostat, and level of theory) for molecules 1 and 2, each trajectory is randomly sampled to collect a dataset of snapshots of molecular configurations (rectangles in orange and purple). Each one of the molecular datasets has N molecular configurations. B) From the two datasets, one for each molecule, a dataset of randomly selected molecular pairs is assembled separated at a distance R. Then, the intermolecular interaction energy is computed for the N molecular pairs, averaged and plotted as a function of R. The same procedure is repeated for several values of R.

Supplementary Note 8. Noncovalent interactions data preparation
In order to analyse the noncovalent interaction energy in a systematic way, a set of molecular pairs separated by regular intermolecular distance (R) intervals has to be created. An approximation to such curve can be obtained by sampling MD trajectories of individual molecules and then create molecular pairs from them. Supplementary Figure 4 shows the followed procedure. First, from a MD trajectory (classical or quantum) for molecule 1 simulated at a given temperature and level of theory, a set of N randomly sampled molecular configurations is created. This is repeated for molecule 2. The two molecules can be different or the same. This is illustrated in Supplementary Figure 4-A. From the two datasets, a new dataset of randomly selected molecular pairs is assembled enforcing a given intermolecular separation R and relative molecular orientation depending of which molecular arrangement wants to be analyzed, for example parallel benzene dimer, and separated at a distance R. Once the dataset of molecular pairs was created, the intermolecular interaction energy for each of the N molecular pairs is computed. Here, we have used symmetry-adapted perturbation theory (SAPT)/aug-cc-pVDZ [28][29][30] as implemented in Psi4 [3][4][5]. The resulting ensemble of interaction energies is then averaged and plotted as a function of R. This same procedure is repeated for several values of R to create an approximate interaction energy curve, as shown in Supplementary Figure 4-B. As mention before, the MD trajectories can be obtained from classical or PIMD simulations, allowing to approximate the change of intermolecular interactions due to the inclusion of nuclear quantum effects.
In order to generate the results in Supplementary Figure 5 of the main text, classical MD and PIMD (with 16 beads) simulations were performed at 300K and sGDML@CCSD(T) level of theory. Then following the method described above and in Supplementary Figure 4 the interaction energy curves were generated for the three configurations of benzene pairs. Each one of the ensembles of benzene pairs considered for a given separation R was formed by 100 molecular pairs. The resulting energy curves are the defined by the expectation energy value of each ensemble. Additionally, calculations estimating the increase in noncovalent interactions were performed for the methane dimer and the methane-benzene complex as shown in Supplementary Figure 5.