Non-covalent bonding patterns are commonly harvested as a design principle in the field of catalysis, supramolecular chemistry, and functional materials to name a few. Yet, their computational description generally neglects finite temperature and environment effects, which promote competing interactions and alter their static gas-phase properties. Recently, neural network potentials (NNPs) trained on density functional theory (DFT) data have become increasingly popular to simulate molecular phenomena in condensed phase with an accuracy comparable to ab initio methods. To date, most applications have centered on solid-state materials or fairly simple molecules made of a limited number of elements. Herein, we focus on the persistence and strength of chalcogen bonds involving a benzotelluradiazole in condensed phase. While the tellurium-containing heteroaromatic molecules are known to exhibit pronounced interactions with anions and lone pairs of different atoms, the relevance of competing intermolecular interactions, notably with the solvent, is complicated to monitor experimentally but also challenging to model at an accurate electronic structure level. Here, we train direct and baselined NNPs to reproduce hybrid DFT energies and forces in order to identify what the most prevalent non-covalent interactions occurring in a solute-Cl–THF mixture are. The simulations in explicit solvent highlight the clear competition with chalcogen bonds formed with the solvent and the short-range directionality of the interaction with direct consequences for the molecular properties in the solution. The comparison with other potentials (e.g., AMOEBA, direct NNP, and continuum solvent model) also demonstrates that baselined NNPs offer a reliable picture of the non-covalent interaction interplay occurring in solution.

Non-covalent interactions (NCIs) impact all areas of chemistry. Among many examples, they govern the structural stability and activity of proteins and DNA,1,2 influence the regioselectivity and enantioselectivity in organocatalytic reactions,3–7 and are used as a design principle to build functional materials.8 Notable NCIs include hydrogen bonds, ππ interactions, anion/cation–π interactions, and σ-hole interactions. While textbook examples commonly depict non-covalently bonded patterns using a static dimer picture, realistic solvated chemical environments involve molecules undergoing a large number of competing NCIs. Favored conformations then arise from a subtle equilibrium between entropic and enthalpic contributions relative to all the possible NCIs. Yet, identifying which interaction plays the most prominent role is not always trivial and even less when solvent molecules intervene.

Achieving an accurate description of competing NCIs in a condensed environment and evaluating their relative importance have thus been a long-lasting goal of computational chemistry.9–12 Traditional static computations in implicit solvent provide an accurate estimation of the strength of individual NCIs, but at the cost of neglecting the role of thermal fluctuations, full entropic effects, and detailed descriptions of the solvation mechanism.13 

The joint experimental and computational investigations by Elmi and co-workers further stress the difficulties of reconciling theory and experiment when taking the solvent into consideration.14,15 They have, for instance, demonstrated that NCIs are strongly attenuated in a solvated environment.16 Despite the intense efforts associated with developing a posteriori dispersion corrections17–35 or with reducing the cost of wavefunction-based methods,36–42 the finite-temperature description of competing NCIs remains coarse if the solvent is represented as a continuum. In particular, interaction energies estimated by dispersion-corrected computations in implicitly modeled solvent can be one magnitude larger than experimental energies in solutions.16,43–45 While the gas-phase data are reproduced accurately by the underlying dispersion-corrected methods, the observed discrepancy suggests the inadequate description of the dispersion compensation and solute–solvent interactions by current implicit solvent models.43,44 As the influence of solvation on the dispersion forces is still being actively investigated,46 proper care should be placed on the unambiguous quantification and characterization of NCIs in solution when implicit solvent is used.47,48

To overcome the shortcomings associated with the simplified picture provided by implicit solvent models and achieve a thorough statistically converged dynamic representation of all the relevant driving forces behind the formation of minima stabilized by various NCIs, it is necessary to rely upon finite temperature molecular dynamic (MD) simulations explicitly, including the solvent molecules. Yet, the description of competing non-covalent interactions is a challenge for both classical and ab initio MD, although for different reasons. Classical force fields, for instance, suffer from the inaccurate description of the electrostatic contributions beyond the monopole approximation. While polarizable force fields (e.g., AMOEBA,49 SIBFA,50 NEMO,51 and GEM52) aim at addressing this issue, their main focus currently remains on the description of proteins and metalloproteins with little emphasis placed on organic molecules or on systems containing somewhat exotic elements. Alternatively, quantum mechanical/molecular mechanical (QM/MM) or full ab initio MD provide an accurate description of individual NCIs, but simulating the free-energy landscapes of non-covalent patterns in a solvated environment is limited by the insufficiently long timescales. The challenge associated with achieving statistically converged simulations at quantum chemical accuracy, thus, remains present.

Stemming from the advances in neural networks (NNs),53 high dimensional neural network-based potentials (NNPs)54,55 trained on density functional theory (DFT) or post-Hartree–Fock energies (and possibly forces), as pioneered by Behler and Parrinello, provide an efficient and accurate alternative to traditional ab initio MD. These approaches are increasingly used in the simulation of solid-state systems,56,57 molecules in the gas phase, liquids,58,59 or molecules on the interfaces.60 Their application in the simulations of molecular systems in solution, however, remains scarce.61,62 Recently, some of us demonstrated the advantages of applying NNPs to go beyond the semi-empirical description of reactions involving hydrogen peroxide, acid, and phenol in an amphiphilic environment.63 The neural network-based simulations showed the importance of relying upon accurate potentials and achieving long timescales for properly capturing the dynamic of the protonation states of the strong acid and the competition between different hydrogen bonds.

Here, we target the description of competing NCIs in solution. Specifically, we focus on chalcogen bonds as a relatively unexplored interaction64 with increasing applicability for the design of organocatalysts, anion transport through lipidic membranes, supramolecular architectures, and functional materials.65–68 The existence of chalcogen bonds is commonly understood as a result of the anisotropic charge density distribution around a covalently bound main group atom, which leads to electron-deficient σ-hole regions that interact with electron-rich moieties. Chalcogen bonds are non-negligible in strength and of marked importance, as common chalcogen-containing molecules are interacting through their σ-holes with lone pairs on N and O atoms and other prevalent Lewis bases.

A necessary prerequisite for warranting the relevance of chalcogen bonding in, for instance, molecular recognition or catalysis is to ensure strong self-association. A relevant heteroaromatic compound inclined to form strong chalcogen bonds in solution through N⋯ chalcogen interactions is 4,5,6,7-tetrafluorobenzo-2,1,3-telluradiazole (further referred only as benzotelluradiazole; see Fig. 1) and the Cl anion. The association constant for this complex in tetrahydrofuran (THF) was recently measured by means of UV–vis absorbance spectroscopy.69 The free energy corresponding to this interaction was estimated to be 7.0 kcal/mol, which aligns well with the chalcogen bond computed at the B97-D3/def2-TZVP level with THF represented by the Polarizable Continuum Model (PCM).69 Apart from the interaction between tellurium and the Cl anion, benzotelluradiazole can, in principle, promote alternative anion–π interactions as well as those involving the lone pair of the solvent molecules.

FIG. 1.

Chalcogen bond between the benzotelluradiazole and Cl. Left and middle: the chemical structure of the bonded complex, right: solvated system used in the simulations—H, C, N, O, F, Cl, and Te atoms are depicted in white, gray, blue, red, yellow, green, and magenta, respectively.

FIG. 1.

Chalcogen bond between the benzotelluradiazole and Cl. Left and middle: the chemical structure of the bonded complex, right: solvated system used in the simulations—H, C, N, O, F, Cl, and Te atoms are depicted in white, gray, blue, red, yellow, green, and magenta, respectively.

Close modal

To provide a realistic and atomistic dynamical description of chalcogen bonds in the condensed phase, we exploit baselined and direct NNPs similar to previous work.63 The choice of the reliable reference method for the baselined NNPs is essential to ensure the description of a broad range of NCIs. As seen in the recent benchmark studies by Řezáč et al.70–72 and Mehta and co-workers,73 the characterization of the interaction strength of NCIs and chalcogen bonding, in particular, is problematic for most popular DFT approximations owing to the possible role of the density errors and its worsening by dispersion corrected schemes. While the most reliable characterization is achieved by modern double hybrid DFT functionals, hybrid DFT functionals were generally shown to be a reasonable compromise.73 Here, we choose PBE074-D3BJ25 as the reference level. We also interface the i-PI dynamic driver with the xTB code75 to enhance the versatility of the baselined NNPs approach and enable the simulations of all the chemical elements present in the investigated system.

The objective of the simulation workflow discussed below is to address the following compelling questions regarding the nature and behavior of chalcogen bonds in solution: How strong and dominant is the interaction in tetrahydrofuran? Is the chalcogen bond and associated σ-hole responsible for the observed association between Cl and benzotelluradiazole or is it a consequence of more complicated interplay between the NCIs? Is the directionality of the chalcogen bond preserved in solution? How well are the relevant solvation effects captured by implicit solvent models?

These questions are addressed through careful comparison of molecular dynamics simulations obtained from two NNPs approaches, the AMOEBA polarizable force field, and static implicit solvent computations.

The workflow used herein (Fig. 2) is a variant of the baselined (BNNP) and direct (DNNP) NNPs training procedure published previously.63 The procedure begins with constructing a robust database of structures covering the relevant part of the conformational space (i.e., chalcogen bond and anion–π interactions), determined from umbrella sampling and temperature replica-exchange simulations at a baseline electronic structure level, which here corresponds to a GFN0-xTB75,76 Hamiltonian. Energies and forces are then recomputed at a more accurate level (PBE074 with D3BJ correction25) that serve as reference data. Behler–Parrinello NNPs54,77 are then trained to reproduce the reference data either directly or using a difference between the baseline and reference method.78 Molecular dynamics together with enhanced sampling techniques are run with each NNPs or in a combined multiple time step (MTS) scheme.79–81 

FIG. 2.

Scheme of a simulation workflow and software packages utilized in this work.

FIG. 2.

Scheme of a simulation workflow and software packages utilized in this work.

Close modal

At each time step of the dynamics, the forces and energies are computed using the LAMMPS software linked to n2p2 (DNNP) and, in addition, Atomistic Simulation Environment (ASE) linked to xTB (BNNP) using the current coordinates of the system. The results are then communicated to the i-PI driver that integrates the equations of motion and propagates the dynamics. The new coordinates of the system are then given back to the external software to update the energies and forces. The process is repeated as a loop for each time step of the dynamics. Similarly, PLUMED is connected to i-PI in order to compute the biased forces and energies whenever constraints are incorporated.

Details on the simulation protocol, preparation, and selection of the structures and symmetry functions chosen for the NNPs training and molecular dynamics trajectory are provided in Secs. II B–II G. A technical description of the electronic structure methods and force field parameterization employed as a comparison is provided in the supplementary material.

The superior accuracy and stability of kernel-based models and NNPs, trained from a baselined electronic structure level (i.e., Δ-ML78), has been largely demonstrated.63,82–88

This work uses GFN0–xTB as a baseline potential as the approach is well-suited for simulating solvated systems featuring a large diversity of chemical elements and environments. The GFNn–xTB family of Hamiltonians75 goes beyond the limits of pairwise parameterization imposed by several semiempirical approaches, such as the related Density Functional based Tight Binding (DFTB).89 In particular, the element-specific empirical fitting enables derivation of parameters that cover the majority of the Periodic Table (up to Z = 86), making GFNn–xTB suitable for computing the broad range of compounds and elements necessary here. GFN0–xTB, which is a non-self-consistent variant of GFNn–xTB, is also significantly faster than other semiempirical approaches.

To integrate GFN0-xTB into the working simulation machinery, we provide a calculator linking the xTB code and the ASE.90 The main advantage in calling xTB via ASE is the existence of a socket interface between ASE and i-PI, which share the information between xTB and ASE without direct modification of the codes. The ASE calculator for the xTB version 6.2 is available on Github (https://github.com/lcmd-epfl/xtb_ase_io_calculator).

The system is composed of one benzotelluradiazole solvated in 45 tetrahydrofuran (THF) molecules and one Cl anion, yielding 599 atoms in a 18.428 Å cubic box (see Fig. 1). The box was minimized and equilibrated to the experimental THF density using the General Amber Force Field in the sander module from AmberTools 16.91 The solute molecules were kept frozen in the geometry optimized at the M06-2X/def2-TZVP level during the equilibration run. Steered molecular dynamics92 were used to generate structures with a Te–Cl distance ranging from 2.5 to 7 Å. 19 structures from the steered MD runs were then used to perform umbrella sampling (US)93 simulations using the GFN0–xTB potential. Each window was simulated for 50 ps using the Abin code.94 The training set was constructed from 2000 structures selected from the US simulation with Farthest Point Sampling (FPS).95 Additional sets of 300 and 1500 structures, respectively, were generated by performing high temperature (1000 K) and two Replica Exchange MD (REMD)96 simulations within the 273.15–700 K temperature range. The first REMD simulation was initiated from structures featuring a chalcogen-bonded Cl. To assess the relevance of the anion–π interaction between Cl and benzotelluradiazole, the second REMD simulation was initiated with Cl placed above the benzotelluradiazole plane. Following a first round of NN training (vide infra), the training set was enlarged with 1200 additional structures taken from a trial US simulation with the trained BNNP. The complete set thus contained 5000 structures with their respective DFT energies and forces (2000 from GFN0–xTB US, 1800 from GFN0–xTB REMD, and 1200 from US using the BNNP trained with the previous 3800 structures). Such a set ensured configurational diversity and sufficient coverage of all the relevant regions of the potential energy landscape.

We used radial (G2) and angular (G4) Behler–Parrinello atomic symmetry functions (ASFs)54,97 as descriptors in the NNPs. ASFs were built for each elements using a cutoff equal to 16 a.u. for BNNP and 14 a.u. for DNNP. We initially generated 273 radial and 2064 angular ASFs per element for the BNNP and 784 and 2064 radial and angular ASFs for the DNNP. The radial functions for the DNNP were generated on a denser grid to provide a better resolution and to stabilize the potential. We used CUR decomposition as described in Ref. 95 to select the most relevant fingerprints for the system. Previous studies63,95 demonstrated that 64 symmetry functions per element yield sufficiently low training errors for systems containing 1 to 4 different elements. Regarding the larger number of species in the mixture, we monitored the root-mean-square error (RMSE) in the training of energies and forces for sets containing 64 to 256 symmetry functions per element. The results are summarized in the supplementary material (see Fig. S2). Based on the training performance, we selected 128 ASFs per element for the BNNP and 256 ASFs per element for the DNNP, yielding 896 (respectively, 1792) symmetry functions for the whole system. BNNPs require a lower number of ASFs than DNNP due to the use of baselined potential, which makes PES flatter and simplifies the training.

The DNNP and BNNP were trained using the Behler–Parrinello neural network potential architecture implemented in the n2p2 version 2.0.3.77 The input values to the DNNP are PBE0-D3BJ energies and forces computed for every structure in the training set. The BNNP is trained to reproduce the difference between the DFT and GFN0-xTB values and, therefore, served as a correction applied to the semiempirical computations. All NNs contained two hidden layers and 22 nodes per layer. The training set structures were normalized before the training, and the symmetry functions were centered and scaled by standard deviation.77 The complete training dataset of 5000 structures was split in 80% for training and 20% for testing. The weights were optimized using a Kalman filter. We use 100% of the energies and 0.5% of the force components per configuration in the training to balance the training accuracy and speed. An ensemble of three neural networks trained on different subsets of the training set was used to compute the BNNP energy and forces, and a single neural network was used for DNNP computations.

All the MD trajectories were propagated using the i-PI code, which was connected with the xTB–ASE interface, LAMMPS drivers for forces and energy computations,98,99 and Plumed 2.6.2,100,101 which served to apply restraints for US. To decrease the cost of BNNP, a multiple-time step (MTS) scheme combining DNNP and BNNP (MTS-NNP) was employed with the BAOAB integrator. In this approach, the dynamics were propagated by cheaper, though less reliable (vide infra), DNNP with a shorter time step. The deficiencies of DNNP were corrected by more computationally demanding BNNP, which was evaluated less frequently, i.e., with a larger simulation time step. In the MTS-NNP, the forces and energies were computed every 3 fs by the GNF0-xTB and BNNP (i.e., the outer time step) and every 0.5 fs by the DNNP (i.e., inner time step). The simulations using only DNNP were integrated with the time step 1 fs using 1 NNP. Stability of the trajectories was improved using both a velocity rescaling thermostat with a frequency of 10 fs and a thermostat based on the generalized Langevin equation102 (GLE) with the same parameters as in Ref. 63. The extrapolation regime of the NNPs during the MD was monitored by printing the extrapolation warnings at every step. The extrapolation warning was triggered when some of the symmetry functions reached value outside the interval on which it was trained.

Potential of Mean Force (PMF) simulations were used to evaluate the strength of the interaction between benzotelluradiazole and Cl. We applied umbrella sampling and selected the distance between the chloride ion and the Te atom as a reaction coordinate (RC). The RC was decomposed into 19 windows, ranging from 2.5 to 7 Å with a width of 0.25 Å. A harmonic restraint of 50 kcal/mol Å2 was applied in all windows. To improve the overlap among the windows in DNNP simulations, additional windows were added at 3.50, 5.00–5.20 Å with a restraint 100 kcal/mol Å2. We ran molecular dynamics with a total simulation time of 200 ps per window for MTS-NNP and 500 ps per window using DNNP alone, corresponding to a total simulation time of 4.8 and 9.5 ns, respectively. For the sake of comparison, we simultaneously ran the same simulation with the AMOEBA polarizable force field (10 ns/per window, total 190 ns). The first ∼20 ps of simulation for each window were used for equilibration and were thus removed from analysis. The free energy profile was reconstructed using the Weighted Histogram Analysis Method (WHAM) as implemented by Grossfield.103 The convergence of the PMF was verified by running the US in the reverse direction (see Figs. S9–S11). US with AMOEBA and DNNP exactly overlapped in both directions demonstrating sufficient convergence of the simulations. Reverse MTS-NNP PMF was ∼1 kcal/mol lower than the forward. Longer sampling might be therefore needed to achieve perfect overlap, but the difference between PMF was in an acceptable span. The confidence of the MTS-NNP results were assessed by running the simulations in two replicas, which demonstrated the reproducibility of the sampling (see Fig. S13).

The main purpose of this work is to build efficient and accurate NNPs to simulate complex mixtures encompassing a variety of atom type and competing NCIs, such as the chalcogen bond investigated herein. From the chemical perspective, our goal is to unravel the origin of the attractive interaction between benzotelluradiazole and Cl solvated in THF and to investigate the significance and structural properties of chalcogen bond in the presence of competing NCIs. We first evaluate the accuracy of the trained DNNP and BNNP and validate their robustness on a test set. We then analyze the outcome of umbrella sampling simulations with these potentials and estimate the strength of the interaction in THF prior to characterizing its nature. Finally, we describe the structural properties of the chalcogen bond in THF and assess the role played by the solvent and by other NCIs through comparisons with simulations performed with the AMOEBA polarizable force field and with static computations in PCM.

Figure 3 depicts the parity plots (i.e., reference vs predicted values) for BNNP and DNNP for the structures in the test set. Both NNPs achieved very low RMSE in the energy, i.e., ∼0.51 meV/atom for DNNP and ∼0.35 meV/atom for BNNP. The RMSE in the prediction of forces is equal to 118 and 80 meV/Å/atom for DNNP and BNNP, respectively.

To analyze the performance of the two NNPs variants and especially their ability to describe the variety of elements, the overall forces parity plots are split into individual plots for each element (see Figs. S4–S6 in the supplementary material). The reference and predicted forces show a good correlation, although larger errors are obtained for the elements in the solute, which are less represented than those in the solvent. The smallest error is obtained for the hydrogen atoms (48 meV/Å for BNNP and 85 meV/Å for DNNP) and the largest for tellurium (244 meV/Å for BNNP and 335 meV/Å for DNNP). Both the DNNP and BNNP forces display a better correlation with the reference DFT data than the original forces computed with GFN0-xTB, which are generally too low in absolute value. Similarly, both NNPs provide a very good description of the explicit solvent, as can be seen from very low errors in forces for the C, O, and H atoms.

FIG. 3.

Parity plot for energy and forces of DNNP (a) and (c) and BNNP (b) and (d). Results are reported after 100 training epochs. Data are obtained for the 1000 structures contained in the test set.

FIG. 3.

Parity plot for energy and forces of DNNP (a) and (c) and BNNP (b) and (d). Results are reported after 100 training epochs. Data are obtained for the 1000 structures contained in the test set.

Close modal

Both NNPs conserve the energy in the NVE simulations (see Figs. S7–S8 in the supplementary material), providing stable trajectories. During the short non-biased NVT runs (i.e., 50 ps), DNNP, however, triggers slightly more extrapolation warnings than the BNNP. Importantly, DNNP with the MTS setup does not extrapolate significantly while substantially decreasing the cost of the simulation. Following these promising results, we evaluate the stability and behavior of the DNNP and MTS-NNP in longer US simulations.

The US simulations are more challenging for the NNPs, as the restraining potential can easily push the structures toward regions underrepresented in the training set. While both potentials yield stable dynamics at long timescales (i.e., more than hundreds of ps), the extrapolation warnings triggered in the US simulations are indeed larger than in non-biased runs for the DNNP simulations (see Figs. S12–S14). Alternatively, the MTS-NNP ensures lower extrapolation risk due to the robustness of the underlying BNNP. Given the relatively weak strength of non-covalent interactions, even a mild difference in the PES may significantly alter the observed interplay. The reliable prediction of energies and forces in the regions outside the training set, i.e., the ability of NNPs to generalize, is therefore crucial. To evaluate the impact of the extrapolation on the predictive power of the NNPs, we build an additional validation set containing 500 structures selected by FPS from the DNNP US trajectories and compute the reference DFT energy and forces. The structures are extracted from the regions with significant extrapolation warnings, i.e., regions with the Te–Cl distance above 3.5 Å.

The predictions of forces [Figs. 4(a) and 4(b)] for both NNPs are in a very good agreement with the PBE0-D3BJ reference, showing a RMSE of 72 and 104 meV/Å/atom for BNNP and DNNP, respectively. The main difference between the NNPs is, however, visible in the prediction of the energy [Fig. 4(c)]. While the BNNP energy correlates well with the reference and RMSE is comparable with the test set (RMSE 0.52 meV/atom, 1.5 times higher), the DNNP strongly underestimates the energies of the selected out-of-sample structures (RMSE 4.82 meV/atom). The error in the energy predictions in these regions is therefore out of control. This behavior further highlights the need for a sufficiently large training set and careful validation of the DNNP, as the apparent stability of the simulations can be misleading and cause misinterpretation of the results. These errors are especially relevant in the context of non-covalent interactions owing to their sensitivity to the quality of the PES. While the issue could be resolved by retraining of the neural networks, including structures from the extrapolated regions, the perfect stabilization of the DNNP is out of the scope of this work. Importantly, MTS-NNP profits from the speed of the DNNP while preventing the extrapolation by corrections with more robust BNNP. Table I demonstrates the computational gain resulting from the application of the MTS-NNP approach compared to the original hybrid DFT, GFN0-xTB, BNNP, and DNNP computations. The MTS-NNP is more than 25 × 104 times faster than original DFT computations and four times faster than the BNNP alone. The sixfold speedup is not reachable in this case as a result of the non-negligible cost of the direct NNP. However, both the GFN0-xTB and NN computations can benefit from parallelization of the corresponding codes, which can further improve the speed of the actual simulations.

FIG. 4.

Energy (c) and forces (a) and (b) parity plots evaluated on an additional test set structures using both NNPs. RMSE for energy: 4.82 (DNNP) and 0.52 meV/atom (BNNP); RMSE for forces: 104 and 72 meV/Å/atom for DNNP and BNNP, respectively.

FIG. 4.

Energy (c) and forces (a) and (b) parity plots evaluated on an additional test set structures using both NNPs. RMSE for energy: 4.82 (DNNP) and 0.52 meV/atom (BNNP); RMSE for forces: 104 and 72 meV/Å/atom for DNNP and BNNP, respectively.

Close modal
TABLE I.

CPU time (core seconds) required to perform MD with a 0.5 fs time step. DFT timing is determined based on a single computing node with two 14 cores Intel Broadwell processors running at 2.6 GHz. GFN0-xTB and NN timings are computed from single-core execution using a node with an Intel Xeon processor running at 2.2 GHz. The BNNP comm. label indicates a prediction obtained with a committee of three BNNPs, and DNNP indicates results from a single direct NN. The MTS cost is computed from a ∼53 CPU s timing for a 1–6 MTS scheme with a 3 fs outer time step, consisting of 1 GFN0-xTB + 8 NN computations (3 BNNPs + 5 DNNPs).

GFN0-xTB MTS-NNP
PBE0-D3BJ+ BNNP comm.GFN0-xTBcomm.DNNP
227 × 103 37 29.7 8.8 2.5 
GFN0-xTB MTS-NNP
PBE0-D3BJ+ BNNP comm.GFN0-xTBcomm.DNNP
227 × 103 37 29.7 8.8 2.5 

The MTS-NNP is therefore perfectly suited for the long MD using enhanced sampling as it provides reliable generalization while significantly decreasing the computational cost. In Secs. III B and III C, we use MTS-NNP as a reference method to analyze the results of the US simulations, assess the impact of possible undesirable extrapolations in DNNP, and evaluate performance of AMOEBA, GFN0-xTB, and PCM in comparison with the NNPs.

The strength of the interaction between benzotelluradiazole and Cl is evaluated via the PMF [see Fig. 5(a)] reconstructed along the Te–Cl distance. The MTS-NNP predicts a favorable association of benzotelluradiazole and Cl with a free energy of 11 kcal/mol. The atomistic resolution of the free energy profiles reveals the directionality and nature of the bonding patterns associated with these energetics. The deep minimum located around the 2.75 Å Te–Cl distance coincides with an N–Te–Cl bond angle in the 160°–180° range [see Fig. 5(b)]. The short distance and nearly linear angle confirms the presence of a strong and directional interaction between benzotelluradiazole and Cl. The directionality does, however, persist only in the relatively short range, i.e., up to around 4.50 Å of Te–Cl distance. At the longer distances, the angle deviates from the liner arrangement and the mode of interaction changes. The evolving interactions occurring in this region are modulated by the explicit solvation responsible for the energy maximum associated with the transition state between the benzotelluradiazole–Cl complex and the solvent-separated structure (vide infra).

FIG. 5.

(a): Comparison of PMF along the Te–Cl coordinate for DNNP, MTS-NNP, AMOEBA, and GFN0-xTB umbrella sampling simulations and relaxed scan in PCM. (b): Evolution of the N–Te–Cl-angle along the collective variable. The average angle is computed for every window separately. Its standard deviation is represented with the shaded regions.

FIG. 5.

(a): Comparison of PMF along the Te–Cl coordinate for DNNP, MTS-NNP, AMOEBA, and GFN0-xTB umbrella sampling simulations and relaxed scan in PCM. (b): Evolution of the N–Te–Cl-angle along the collective variable. The average angle is computed for every window separately. Its standard deviation is represented with the shaded regions.

Close modal

The PMF profiles modeled with other potentials [see Fig. 5(a)] substantially differ from the MTS-NNP one. The profile at the bare GFN0-xTB level is shifted toward lower Te–Cl distances. This behavior is likely caused by the wrong description of the repulsive and electrostatic part of the potential, a limitation that is well known as GFN0-xTB does not use self-consistent evaluation of the charges.75 The difference between GFN0-xTB and MTS-NNP stresses the necessity of using a more accurate DFT level for achieving a proper description of the benzotelluradiazole–Cl interaction.

DNNP, AMOEBA, and the relaxed scan using a PCM all identify an energy minimum at a distance in agreement with MTS-NNP. Yet, the free energy obtained with AMOEBA is dramatically lower (about 1.8 kcal/mol) than the 11 kcal/mol of MTS-NNP, which is itself slightly higher than the experimental value (i.e., 1.5 times stronger). Akin to related halogen bonds, charge transfer contributions are expected to be part of the interaction energy.104–109 While the exact role of charge transfer in σ-hole interactions is still under debate,110–112 the original AMOEBA force field expression suffers from a lack of explicit charge transfer description,49 which could explain its collapse in the long-range. Similarly to MTS-NNP, the AMOEBA profile displays a maximum coinciding with the transition state, but at a shorter distance, as expected from the weaker interaction energy.

The 10 kcal/mol free energy values of the DNNP and the relaxed PCM scan are in better agreement with MTS-NNP. However, the smooth energy profile of the latter contrasts with the wiggly behavior of the DNNP profile, which is likely an outcome of the inconsistency discussed in Sec. III A. Finally, the relaxed scan in PCM reaches a plateau in the Te–Cl distance around 5 Å indicative of an absence of interaction above this range. The absence of maximum is directly related to the missing explicit solvent molecules.

Analysis of the N–Te–Cl angle [Fig. 5(b)] further confirms the presence of directional chalcogen bonding in these PMF profiles in the short range regions. In longer distances, however, the directionality with Cl is lost with the angle dropping below 160°similarly to MTS-NNP. In both the DNNP and AMOEBA simulations, the chalcogen bond becomes loose already around the Te–Cl distance of 3.5 Å. The implicit solvent energy profile using the PCM predicts a shift of behavior at around 4.00 Å.

Overall, the free energy profiles confirm the presence of a chalcogen bond between the benzotelluradiazole and Cl. The MTS-NNP result further highlights the strength and importance of this interaction, especially in the short range regions (below 4.50 Å). Yet, monitoring the PMF and angle is not sufficient to assess what is happening once the chalcogen bond is broken and to deeper understand the observed discrepancies between the MTS-NNP and other methods. The loss of directionality does not necessarily imply that Cl stop interacting with the Te atom and with benzotelluradiazole at larger distances. While not measurable by UV–vis spectroscopy, other interactions might influence the chalcogen bond and affect its properties, namely, anion–π interactions and competition between Te–Cl and Te–O chalcogen bonds with the anion and solvent molecules, respectively. The existence of competing interactions between benzotelluradiazole, THF, and the chloride anion is further analyzed by the density distributions of oxygen atoms from solvent molecules and Cl around the heteroatom discussed in Sec. III C.

Figure 6 summarizes the density distributions of O and Cl atoms in the windows showcasing the largest differences between the MTS-NNP and DNNP PMF curves as given in Fig. 5 (i.e., at Te–Cl distances of 3.00, 3.50, 4.75, and 5.25 Å).

FIG. 6.

3D density distribution of oxygen (red contour) and Cl (green contour) atoms around benzotelluradiazole. The density isocontours correspond to the density values of 0.005 (transparent) and 0.01 (opaque). A reference frame for density computation is defined by the N–Te–N atoms. C, N, F, and Te atoms are depicted in gray, blue, yellow, and magenta, respectively.

FIG. 6.

3D density distribution of oxygen (red contour) and Cl (green contour) atoms around benzotelluradiazole. The density isocontours correspond to the density values of 0.005 (transparent) and 0.01 (opaque). A reference frame for density computation is defined by the N–Te–N atoms. C, N, F, and Te atoms are depicted in gray, blue, yellow, and magenta, respectively.

Close modal

The MTS-NNP results confirms a directional chalcogen bond interaction between Te and Cl, which is apparent especially in the 3.00 Å distance window, which pictures the maximum of Cl density located close to the position of σ-hole. In the longer Te–Cl distances, interaction toward Cl weakens and the affinity toward the oxygen atoms of THF is amplified instead. The windows at Te–Cl distances of 4.75 and 5.25 Å finally demonstrate complete replacement of Te and Cl chalcogen bond by directional chalcogen bonds toward oxygen atoms, which occupy both σ-hole regions. The directionality of the newly formed Te–O bonds is further confirmed by the N–Te–O angle distributions (see Fig. S18 in the supplementary material). This behavior further stresses the importance of Te–Cl chalcogen bonds in the short range with its gradual replacement by Te–O chalcogen bonds at the Te–Cl distances above 4.00 Å, which explains the observed trends in the N–Te–Cl angle distributions. The density distribution at the 5.25 Å window depicts the transition from the Te–Cl complex to the fully solvated Cl corresponding to the maximum of the PMF profile located around 5.2 Å [see Fig. 5(a)].

At 3.00 Å, the density distributions with DNNP and AMOEBA already differ with those obtained with the MTS-NNP. They promote an interaction with Cl but also one with a THF oxygen atom, which is not present using MTS-NNP. At Te–Cl distances larger than 3.00 Å, the directional interaction with Cl is entirely replaced by the Te–O chalcogen bond interactions. With AMOEBA, the interaction between Te and Cl decreases progressively to ultimately vanish at distances greater than 5.25 Å. This behavior, which is aligned with the PMF, demonstrates the limitations of AMOEBA for describing chalcogen bonds involving tellurium. The density distributions with DNNP suggest an overly strong affinity for the solvent and, to some extent, too strong anion–pi interactions that are not observed in the MTS-NNP simulations. These artifacts arise from a poor extrapolation of the neural network potential in these regions and further emphasize the importance of its careful validation.

In general, the density distributions presented here show that the two σ-hole regions on the tellurium atom are good donors of chalcogen bond toward both the Cl anion and the THF oxygen atom. Benzotelluradiazole, therefore, interacts not only with the Cl but also with the solvent, which is not apparent from implicit solvent computations. The presence of the explicit solvent molecules is thus essential to properly describe the breaking of the chalcogen bond arising from competing interactions. More importantly, the solvent containing atoms in possession of lone pairs (e.g., O and N) can directly compete with the chalcogen bond acceptors and influence the interaction strength. The impact of the different solvents on the chalcogen bonded molecules was already observed in the formation of aggregates and crystallization.113 For instance, similar tellurium-containing compound crystallizes in a form of Te–O bonded macrocycles with a stoichiometry dictated directly by the solvent effects.114 Interestingly, THF molecules were even found inside the formed cavities further advocating for their interaction with tellurium. Chalcogen bonds between selenium containing species and oxygen and nitrogen atoms of solvent molecules were also confirmed by nuclear magnetic resonance (NMR) measurements.115 

Overall, this section provides key chemical and technical insights regarding both the nature and importance of chalcogen bonds and the reliability of their condensed phase simulations using distinct potentials, ranging from polarizable force fields to neural network surrogates for hybrid DFT. This comparison demonstrates the larger reliability and robustness of the MTS-NNP approach and the difficulties of the DNNP and polarizable force field in the description of challenging charge transfer complexes. From the chemical perspective, it is evident that chalcogen bonds will unavoidably form with the solvent molecules and will consequently enforce a more structured solvation shell near the solute, which cannot be captured by implicit solvent models.

This computational work investigates the persistence of chalcogen bonds between benzotelluradiazole and the Cl anion solvated in tetrahydrofuran. To ensure both quantum chemical accuracy and statistical convergence, we simulate the complex multielement systems made of 599 atoms relying upon Behler–Parrinello neural network potentials. We compare the performance of a MTS-NNP and a direct NNP, which are both trained to reproduce hybrid DFT energies and forces. The two variants show an impressive stability along the nanosecond timescale molecular dynamic simulations. The stability of the baseline-free DNNP is especially encouraging and, however, still requires significant effort to improve its generalization properties. From the technical perspective, we demonstrate the versatility of the MTS-NNP approach combining baselined and direct NNP in the simulation of systems incorporating several types of competing non-covalent interactions.

The application of this method brings important chemical insight into the behavior of chalcogen bonded molecules in the explicit solvent, which would not be possible to obtain with standard state-of-the-art approaches, i.e., ab initio molecular dynamics or polarizable force fields. Specifically, we confirm the persistence of the chalcogen bond in solution and observe an important competition between the chalcogen bond toward Cl and the solvent oxygen atoms. While the Te–Cl chalcogen bond is a dominant non-covalent interaction in the short range region up to Te–Cl distances of 4.00 Å, the interaction between the Te atom and oxygen of THF is prevailing in larger distances. This behavior coincides with the loss of directionality of the Te–Cl interaction and highlights the non-negligible impact of oxygen or nitrogen containing solvents on the topology of the chalcogen bond.

This work demonstrates that neural networks can be successfully applied to model the energy profiles of complex molecular systems interacting with the solvent. They provide important atomistic resolution of the solvation mechanism that helps rationalizing and harvesting the nature and behavior of chalcogen bonds in solution. Similar approaches would be useful to assess the role of Lewis basic solvents in the co-crystallization of chalcogen-bearing compounds and the potential reduction of catalysis activity or association constants induced by the solvent. More generally, it does open efficient and reliable routes to simulate homogeneous organocatalysts, host–guest systems, and other molecular systems exploiting a subtle balance of non-covalent interactions.

See the supplementary material for computational details, information on AMOEBA parameterization, NNP validation, convergence of PMF, and supplementary plots.

The authors are grateful to the EPFL for financial support and computational resources. V.J. acknowledges the funding from the Swiss National Science Foundation (SNSF, Grant No. 200020_175496), and F.C. and C.C. acknowledge the funding from the European Research Council (ERC, Grant No. 817977) within the framework of European Union’s Horizon 2020. The National Center of Competence in Research (NCCR) “Sustainable chemical process through catalysis (Catalysis)” of SNSF is acknowledged for financial support to R.L. The authors thank Raimon Fabregat for the assistance with graphical material, Daniel Hollas for the help with the Abin setup, and Kevin Rossi and Michele Ceriotti for scientific discussions.

The authors have no conflicts to disclose.

The xTB-ase calculator is available on Github at https://github.com/lcmd-epfl/xtb_ase_io_calculator. The data and the neural networks parameters are freely available on the Materials Cloud at https://archive.materialscloud.org/record/2022.42.

1.
E.
Frieden
, “
Non-covalent interactions: Key to biological flexibility and specificity
,”
J. Chem. Educ.
52
,
754
(
1975
).
2.
J.
Černý
and
P.
Hobza
, “
Non-covalent interactions in biomacromolecules
,”
Phys. Chem. Chem. Phys.
9
,
5291
5303
(
2007
).
3.
H. J.
Davis
and
R. J.
Phipps
, “
Harnessing non-covalent interactions to exert control over regioselectivity and site-selectivity in catalytic reactions
,”
Chem. Sci.
8
,
864
877
(
2017
).
4.
M.
Orlandi
,
J. A. S.
Coelho
,
M. J.
Hilton
,
F. D.
Toste
, and
M. S.
Sigman
, “
Parametrization of non-covalent interactions for transition state interrogation applied to asymmetric catalysis
,”
J. Am. Chem. Soc.
139
,
6803
6806
(
2017
).
5.
R. S. J.
Proctor
,
A. C.
Colgan
, and
R. J.
Phipps
, “
Exploiting attractive non-covalent interactions for the enantioselective catalysis of reactions involving radical intermediates
,”
Nat. Chem.
12
,
990
1004
(
2020
).
6.
S. J.
Zuend
and
E. N.
Jacobsen
, “
Mechanism of amido-thiourea catalyzed enantioselective imine hydrocyanation: Transition state stabilization via multiple non-covalent interactions
,”
J. Am. Chem. Soc.
131
,
15358
15374
(
2009
).
7.
C.
Uyeda
and
E. N.
Jacobsen
, “
Transition-state charge stabilization through multiple non-covalent interactions in the guanidinium-catalyzed enantioselective claisen rearrangement
,”
J. Am. Chem. Soc.
133
,
5062
5075
(
2011
).
8.
P. C.
Ho
,
J. Z.
Wang
,
F.
Meloni
, and
I.
Vargas-Baca
, “
Chalcogen bonding in materials chemistry
,”
Coord. Chem. Rev.
422
,
213464
(
2020
).
9.
K. E.
Riley
and
P.
Hobza
, “
Noncovalent interactions in biochemistry
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
1
,
3
17
(
2011
).
10.
B.
Meyer
,
S.
Barthel
,
A.
Mace
,
L.
Vannay
,
B.
Guillot
,
B.
Smit
, and
C.
Corminboeuf
, “
Dori reveals the influence of noncovalent interactions on covalent bonding patterns in molecular crystals under pressure
,”
J. Phys. Chem. Lett.
10
,
1482
1488
(
2019
).
11.
P.
Hobza
,
R.
Zahradník
, and
K.
Müller-Dethlefs
, “
The world of non-covalent interactions: 2006
,”
Collect. Czech. Chem. Commun.
71
,
443
531
(
2006
).
12.
A. M.
Maharramov
,
K. T.
Mahmudov
,
M. N.
Kopylovich
, and
A. J.
Pombeiro
,
Non-Covalent Interactions in the Synthesis and Design of New Compounds
(
John Wiley & Sons
,
2016
).
13.
H.-J.
Schneider
, “
Quantification of noncovalent interactions – Promises and problems
,”
New J. Chem.
43
,
15498
15512
(
2019
).
14.
A.
Elmi
and
S. L.
Cockroft
, “
Quantifying interactions and solvent effects using molecular balances and model complexes
,”
Acc. Chem. Res.
54
,
92
103
(
2021
).
15.
I. K.
Mati
,
C.
Adam
, and
S. L.
Cockroft
, “
Seeing through solvent effects using molecular balances
,”
Chem. Sci.
4
,
3965
3972
(
2013
).
16.
L.
Yang
,
C.
Adam
,
G. S.
Nichol
, and
S. L.
Cockroft
, “
How much do van der Waals dispersion forces contribute to molecular recognition in solution?
,”
Nat. Chem.
5
,
1006
1010
(
2013
).
17.
A. D.
Becke
and
E. R.
Johnson
, “
Exchange-hole dipole moment and the dispersion interaction
,”
J. Chem. Phys.
122
,
154104
(
2005
).
18.
A. D.
Becke
and
E. R.
Johnson
, “
A density-functional model of the dispersion interaction
,”
J. Chem. Phys.
123
,
154101
(
2005
).
19.
E. R.
Johnson
and
A. D.
Becke
, “
A post-Hartree–Fock model of intermolecular interactions
,”
J. Chem. Phys.
123
,
024101
(
2005
).
20.
A. D.
Becke
and
E. R.
Johnson
, “
Exchange-hole dipole moment and the dispersion interaction: High-order dispersion coefficients
,”
J. Chem. Phys.
124
,
014104
(
2006
).
21.
A. D.
Becke
and
E. R.
Johnson
, “
Exchange-hole dipole moment and the dispersion interaction revisited
,”
J. Chem. Phys.
127
,
154108
(
2007
).
22.
T.
Sato
and
H.
Nakai
, “
Density functional method including weak interactions: Dispersion coefficients based on the local response approximation
,”
J. Chem. Phys.
131
,
224104
(
2009
).
23.
S.
Grimme
,
J.
Antony
,
S.
Ehrlich
, and
H.
Krieg
, “
A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu
,”
J. Chem. Phys.
132
,
154104
(
2010
).
24.
S.
Grimme
, “
Density functional theory with London dispersion corrections
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
1
,
211
228
(
2011
).
25.
S.
Grimme
,
S.
Ehrlich
, and
L.
Goerigk
, “
Effect of the damping function in dispersion corrected density functional theory
,”
J. Comput. Chem.
32
,
1456
1465
(
2011
).
26.
S.
Grimme
,
A.
Hansen
,
J. G.
Brandenburg
, and
C.
Bannwarth
, “
Dispersion-corrected mean-field electronic structure methods
,”
Chem. Rev.
116
,
5105
5154
(
2016
).
27.
A.
Tkatchenko
and
M.
Scheffler
, “
Accurate molecular van der Waals interactions from ground-state electron density and free-atom reference data
,”
Phys. Rev. Lett.
102
,
073005
(
2009
).
28.
A.
Tkatchenko
,
R. A.
Distasio
,
R.
Car
, and
M.
Scheffler
, “
Accurate and efficient method for many-body van der Waals interactions
,”
Phys. Rev. Lett.
108
,
236402
(
2012
).
29.
S. N.
Steinmann
and
C.
Corminboeuf
, “
A system-dependent density-based dispersion correction
,”
J. Chem. Theory Comput.
6
,
1990
2001
(
2010
).
30.
S. N.
Steinmann
and
C.
Corminboeuf
, “
Comprehensive benchmarking of a density-dependent dispersion correction
,”
J. Chem. Theory Comput.
7
,
3567
3577
(
2011
).
31.
S. N.
Steinmann
and
C.
Corminboeuf
, “
A generalized-gradient approximation exchange hole model for dispersion coefficients
,”
J. Chem. Phys.
134
,
044117
(
2011
).
32.
C.
Corminboeuf
, “
Minimizing density functional failures for non-covalent interactions beyond van der Waals complexes
,”
Acc. Chem. Res.
47
,
3217
3224
(
2014
).
33.
S.
Ehrlich
,
J.
Moellmann
,
W.
Reckien
,
T.
Bredow
, and
S.
Grimme
, “
System-dependent dispersion coefficients for the DFT-D3 treatment of adsorption processes on ionic surfaces
,”
ChemPhysChem
12
,
3414
3420
(
2011
).
34.
E.
Caldeweyher
,
S.
Ehlert
,
A.
Hansen
,
H.
Neugebauer
,
S.
Spicher
,
C.
Bannwarth
, and
S.
Grimme
, “
A generally applicable atomic-charge dependent London dispersion correction
,”
J. Chem. Phys.
150
,
154122
(
2019
).
35.
E.
Caldeweyher
,
J.-M.
Mewes
,
S.
Ehlert
, and
S.
Grimme
, “
Extension and evaluation of the D4 London-dispersion model for periodic systems
,”
Phys. Chem. Chem. Phys.
22
,
8499
8512
(
2020
).
36.
F.
Neese
,
F.
Wennmohs
, and
A.
Hansen
, “
Efficient and accurate local approximations to coupled-electron pair approaches: An attempt to revive the pair natural orbital method
,”
J. Chem. Phys.
130
,
114108
(
2009
).
37.
A. E.
DePrince
and
C. D.
Sherrill
, “
Accurate noncovalent interaction energies using truncated basis sets based on frozen natural orbitals
,”
J. Chem. Theory Comput.
9
,
293
299
(
2013
).
38.
A. E.
DePrince
and
C. D.
Sherrill
, “
Accuracy and efficiency of coupled-cluster theory using density fitting/Cholesky decomposition, frozen natural orbitals, and a t1-transformed Hamiltonian
,”
J. Chem. Theory Comput.
9
,
2687
2696
(
2013
).
39.
C.
Riplinger
and
F.
Neese
, “
An efficient and near linear scaling pair natural orbital based local coupled cluster method
,”
J. Chem. Phys.
138
,
034106
(
2013
).
40.
R. M.
Parrish
,
C. D.
Sherrill
,
E. G.
Hohenstein
,
S. I. L.
Kokkila
, and
T. J.
Martínez
, “
Communication: Acceleration of coupled cluster singles and doubles via orbital-weighted least-squares tensor hypercontraction
,”
J. Chem. Phys.
140
,
181102
(
2014
).
41.
M.
Sparta
and
F.
Neese
, “
Chemical applications carried out by local pair natural orbital based coupled-cluster methods
,”
Chem. Soc. Rev.
43
,
5032
5041
(
2014
).
42.
B. S.
Fales
,
E. R.
Curtis
,
K. G.
Johnson
,
D.
Lahana
,
S.
Seritan
,
Y.
Wang
,
H.
Weir
,
T. J.
Martínez
, and
E. G.
Hohenstein
, “
Performance of coupled-cluster singles and doubles on modern stream processing architectures
,”
J. Chem. Theory Comput.
16
,
4021
4028
(
2020
).
43.
R.
Pollice
,
M.
Bot
,
I. J.
Kobylianskii
,
I.
Shenderovich
, and
P.
Chen
, “
Attenuation of London dispersion in dichloromethane solutions
,”
J. Am. Chem. Soc.
139
,
13126
13140
(
2017
).
44.
R.
Pollice
,
F.
Fleckenstein
,
I.
Shenderovich
, and
P.
Chen
, “
Compensation of London dispersion in the gas phase and in aprotic solvents
,”
Angew. Chem., Int. Ed.
58
,
14281
14288
(
2019
).
45.
A.
Hansen
,
C.
Bannwarth
,
S.
Grimme
,
P.
Petrović
,
C.
Werlé
, and
J.-P.
Djukic
, “
The thermochemistry of London dispersion-driven transition metal reactions: Getting the ‘right answer for the right reason
,”
ChemistryOpen
3
,
177
189
(
2014
).
46.
J. M.
Schümann
,
J. P.
Wagner
,
A. K.
Eckhardt
,
H.
Quanz
, and
P. R.
Schreiner
, “
Intramolecular London dispersion interactions do not cancel in solution
,”
J. Am. Chem. Soc.
143
,
41
45
(
2021
).
47.
S.
Grimme
, “
Comment on: “On the accuracy of DFT methods in reproducing ligand substitution energies for transition metal complexes in solution: The role of dispersive interactions” by H. Jacobsen and L. Cavallo
,”
ChemPhysChem
13
,
1407
1409
(
2012
).
48.
K. E.
Riley
,
J.
Vondrasek
, and
P.
Hobza
, “
Performance of the DFT-D method, paired with the PCM implicit solvation model, for the computation of interaction energies of solvated complexes of biological interest
,”
Phys. Chem. Chem. Phys.
9
,
5555
5560
(
2007
).
49.
J. W.
Ponder
,
C.
Wu
,
P.
Ren
,
V. S.
Pande
,
J. D.
Chodera
,
M. J.
Schnieders
,
I.
Haque
,
D. L.
Mobley
,
D. S.
Lambrecht
,
R. A.
DiStasio
, Jr.
 et al, “
Current status of the AMOEBA polarizable force field
,”
J. Phys. Chem. B
114
,
2549
2564
(
2010
).
50.
M.
Devereux
,
N.
Gresh
,
J.-P.
Piquemal
, and
M.
Meuwly
, “
A supervised fitting approach to force field parametrization with application to the SIBFA polarizable force field
,”
J. Comput. Chem.
35
,
1577
1591
(
2014
).
51.
A.
Holt
,
J.
Boström
,
G.
Karlström
, and
R.
Lindh
, “
A NEMO potential that includes the dipole–quadrupole and quadrupole–quadrupole polarizability
,”
J. Comput. Chem.
31
,
1583
1591
(
2010
).
52.
H.
Torabifard
,
O. N.
Starovoytov
,
P.
Ren
, and
G. A.
Cisneros
, “
Development of an AMOEBA water model using GEM distributed multipoles
,”
Theor. Chem. Acc.
134
,
1
10
(
2015
).
53.
R. M.
Balabin
and
E. I.
Lomakina
, “
Neural network approach to quantum-chemistry data: Accurate prediction of density functional theory energies
,”
J. Chem. Phys.
131
,
074104
(
2009
).
54.
J.
Behler
and
M.
Parrinello
, “
Generalized neural-network representation of high-dimensional potential-energy surfaces
,”
Phys. Rev. Lett.
98
,
146401
(
2007
).
55.
J.
Behler
, “
First principles neural network potentials for reactive simulations of large molecular and condensed systems
,”
Angew. Chem., Int. Ed.
56
,
12828
12840
(
2017
).
56.
N.
Artrith
,
T.
Morawietz
, and
J.
Behler
, “
High-dimensional neural-network potentials for multicomponent systems: Applications to zinc oxide
,”
Phys. Rev. B
83
,
153101
(
2011
).
57.
G.
Imbalzano
and
M.
Ceriotti
, “
Modeling the Ga/As binary system across temperatures and compositions from first principles
,”
Phys. Rev. Mater.
5
,
063804
(
2021
).
58.
T.
Morawietz
,
A.
Singraber
,
C.
Dellago
, and
J.
Behler
, “
How van der Waals interactions determine the unique properties of water
,”
Proc. Natl. Acad. Sci. U. S. A.
113
,
8368
8373
(
2016
).
59.
B.
Cheng
,
G.
Mazzola
,
C. J.
Pickard
, and
M.
Ceriotti
, “
Evidence for supercritical behaviour of high-pressure liquid hydrogen
,”
Nature
585
,
217
220
(
2020
).
60.
S. K.
Natarajan
and
J.
Behler
, “
Neural network molecular dynamics simulations of solid–liquid interfaces: Water at low-index copper surfaces
,”
Phys. Chem. Chem. Phys.
18
,
28704
28725
(
2016
).
61.
M.
Yang
,
L.
Bonati
,
D.
Polino
, and
M.
Parrinello
, “
Using metadynamics to build neural network potentials for reactive events: The case of urea decomposition in water
,”
Catal. Today
387
,
143
149
(
2022
).
62.
C.
Schran
,
F. L.
Thiemann
,
P.
Rowe
,
E. A.
Müller
,
O.
Marsalek
, and
A.
Michaelides
, “
Machine learning potentials for complex aqueous systems made simple
,”
Proc. Natl. Acad. Sci. U. S. A.
118
,
e2110077118
(
2021
).
63.
K.
Rossi
,
V.
Jurásková
,
R.
Wischert
,
L.
Garel
,
C.
Corminbœuf
, and
M.
Ceriotti
, “
Simulating solvation and acidity in complex mixtures with first-principles accuracy: The case of CH3SO3H and H2O2 in phenol
,”
J. Chem. Theory Comput.
16
,
5139
5149
(
2020
).
64.
C. B.
Aakeroy
,
D. L.
Bryce
,
G. R.
Desiraju
,
A.
Frontera
,
A. C.
Legon
,
F.
Nicotra
,
K.
Rissanen
,
S.
Scheiner
,
G.
Terraneo
,
P.
Metrangolo
 et al, “
Definition of the chalcogen bond (IUPAC Recommendations 2019)
,”
Pure Appl. Chem.
91
,
1889
1892
(
2019
).
65.
P.
Wonner
,
A.
Dreger
,
L.
Vogel
,
E.
Engelage
, and
S. M.
Huber
, “
Chalcogen bonding catalysis of a nitro-michael reaction
,”
Angew. Chem., Int. Ed.
58
,
16923
16927
(
2019
).
66.
S.
Benz
,
M.
Macchione
,
Q.
Verolet
,
J.
Mareda
,
N.
Sakai
, and
S.
Matile
, “
Anion transport with chalcogen bonds
,”
J. Am. Chem. Soc.
138
,
9093
9096
(
2016
).
67.
M.
Macchione
,
M.
Tsemperouli
,
A.
Goujon
,
A. R.
Mallia
,
N.
Sakai
,
K.
Sugihara
, and
S.
Matile
, “
Mechanosensitive oligodithienothiophenes: Transmembrane anion transport along chalcogen-bonding cascades
,”
Helv. Chim. Acta
101
,
e1800014
(
2018
).
68.
K. T.
Mahmudov
,
M. N.
Kopylovich
,
M. F. C.
Guedes da Silva
, and
A. J. L.
Pombeiro
, “
Chalcogen bonding in synthesis, catalysis and design of materials
,”
Dalton Trans.
46
,
10121
10138
(
2017
).
69.
G. E.
Garrett
,
G. L.
Gibson
,
R. N.
Straus
,
D. S.
Seferos
, and
M. S.
Taylor
, “
Chalcogen bonding in solution: Interactions of benzotelluradiazoles with anionic and uncharged Lewis bases
,”
J. Am. Chem. Soc.
137
,
4126
4133
(
2015
).
70.
J.
Řezáč
, “
Non-covalent interactions atlas benchmark data sets: Hydrogen bonding
,”
J. Chem. Theory Comput.
16
,
2355
2368
(
2020
).
71.
J.
Řezáč
, “
Non-covalent interactions atlas benchmark data sets 2: Hydrogen bonding in an extended chemical space
,”
J. Chem. Theory Comput.
16
,
6305
6316
(
2020
).
72.
K.
Kříž
,
M.
Nováček
, and
J.
Řezáč
, “
Non-covalent interactions atlas benchmark data sets 3: Repulsive contacts
,”
J. Chem. Theory Comput.
17
,
1548
1561
(
2021
).
73.
N.
Mehta
,
T.
Fellowes
,
J. M.
White
, and
L.
Goerigk
, “
CHAL336 benchmark set: How well do quantum-chemical methods describe chalcogen-bonding interactions?
,”
J. Chem. Theory Comput.
17
,
2783
2806
(
2021
).
74.
C.
Adamo
and
V.
Barone
, “
Toward reliable density functional methods without adjustable parameters: The PBE0 model
,”
J. Chem. Phys.
110
,
6158
6170
(
1999
).
75.
C.
Bannwarth
,
E.
Caldeweyher
,
S.
Ehlert
,
A.
Hansen
,
P.
Pracht
,
J.
Seibert
,
S.
Spicher
, and
S.
Grimme
, “
Extended tight-binding quantum chemistry methods
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
11
,
e1493
(
2021
).
76.
P.
Pracht
,
E.
Caldeweyher
,
S.
Ehlert
, and
S.
Grimme
, “
A robust non-self-consistent tight-binding quantum chemistry method for large molecules
,” ChemRxiv:8326202.v1 (
2019
).
77.
A.
Singraber
,
T.
Morawietz
,
J.
Behler
, and
C.
Dellago
, “
Parallel multistream training of high-dimensional neural network potentials
,”
J. Chem. Theory Comput.
15
,
3075
3092
(
2019
).
78.
R.
Ramakrishnan
,
P. O.
Dral
,
M.
Rupp
, and
O. A.
von Lilienfeld
, “
Big data meets quantum chemistry approximations: The Δ-machine learning approach
,”
J. Chem. Theory Comput.
11
,
2087
2096
(
2015
).
79.
M.
Tuckerman
,
B. J.
Berne
, and
G. J.
Martyna
, “
Reversible multiple time scale molecular dynamics
,”
J. Chem. Phys.
97
,
1990
(
1992
).
80.
N.
Luehr
,
T. E.
Markland
, and
T. J.
Martínez
, “
Multiple time step integrators in ab initio molecular dynamics
,”
J. Chem. Phys.
140
,
084116
(
2014
).
81.
V.
Kapil
,
J.
VandeVondele
, and
M.
Ceriotti
, “
Accurate molecular dynamics and nuclear quantum effects at low cost by multiple steps in real and imaginary time: Using density functional theory to accelerate wavefunction methods
,”
J. Chem. Phys.
144
,
054111
(
2016
); arXiv:1512.00176.
82.
L.
Shen
,
J.
Wu
, and
W.
Yang
, “
Multiscale quantum mechanics/molecular mechanics simulations with neural networks
,”
J. Chem. Theory Comput.
12
,
4934
4946
(
2016
).
83.
A. P.
Bartók
,
S.
De
,
C.
Poelking
,
N.
Bernstein
,
J. R.
Kermode
,
G.
Csányi
, and
M.
Ceriotti
, “
Machine learning unifies the modeling of materials and molecules
,”
Sci. Adv.
3
,
e1701816
(
2017
).
84.
B.
Meyer
,
B.
Sawatlon
,
S.
Heinen
,
O. A.
Von Lilienfeld
, and
C.
Corminboeuf
, “
Machine learning meets volcano plots: Computational discovery of cross-coupling catalysts
,”
Chem. Sci.
9
,
7069
7077
(
2018
).
85.
L.
Shen
and
W.
Yang
, “
Molecular dynamics simulations with quantum mechanics/molecular mechanics and adaptive neural networks
,”
J. Chem. Theory Comput.
14
,
1442
1455
(
2018
).
86.
G.
Sun
and
P.
Sautet
, “
Toward fast and reliable potential energy surfaces for metallic Pt clusters by hierarchical delta neural networks
,”
J. Chem. Theory Comput.
15
,
5614
5627
(
2019
).
87.
R.
Fabregat
,
A.
Fabrizio
,
B.
Meyer
,
D.
Hollas
, and
C.
Corminboeuf
, “
Hamiltonian-reservoir replica exchange and machine learning potentials for computational organic chemistry
,”
J. Chem. Theory Comput.
16
,
3084
3094
(
2020
).
88.
M.
Bogojeski
,
L.
Vogt-Maranto
,
M. E.
Tuckerman
,
K.-R.
Müller
, and
K.
Burke
, “
Quantum chemical accuracy from density functional approximations via machine learning
,”
Nat. Commun.
11
,
5223
(
2020
).
89.
M.
Gaus
,
A.
Goez
, and
M.
Elstner
, “
Parametrization and benchmark of DFTB3 for organic molecules
,”
J. Chem. Theory Comput.
9
,
338
354
(
2013
).
90.
A.
Hjorth Larsen
,
J.
Jørgen Mortensen
,
J.
Blomqvist
,
I. E.
Castelli
,
R.
Christensen
,
M.
Dułak
,
J.
Friis
,
M. N.
Groves
,
B.
Hammer
,
C.
Hargus
 et al, “
The atomic simulation environment—A python library for working with atoms
,”
J. Phys.: Condens. Matter
29
,
273002
(
2017
).
91.
D. A.
Case
,
T. E.
Cheatham
 III
,
T.
Darden
,
H.
Gohlke
,
R.
Luo
,
K. M.
Merz
Jr
.,
A.
Onufriev
,
C.
Simmerling
,
B.
Wang
, and
R. J.
Woods
, “
The Amber biomolecular simulation programs
,”
J. Comput. Chem.
26
,
1668
1688
(
2005
).
92.
S.
Izrailev
,
S.
Stepaniants
,
B.
Isralewitz
,
D.
Kosztin
,
H.
Lu
,
F.
Molnar
,
W.
Wriggers
, and
K.
Schulten
, “
Steered molecular dynamics
,” in
Computational Molecular Dynamics: Challenges, Methods, Ideas
(
Springer
,
1999
), pp.
39
65
.
93.
G. M.
Torrie
and
J. P.
Valleau
, “
Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling
,”
J. Comput. Phys.
23
,
187
199
(
1977
).
94.
D.
Hollas
,
J.
Suchan
,
M.
Ončák
, and
P.
Slavíček
, Photox/Abin: Pre-release of Version 1.1,
2018
.
95.
G.
Imbalzano
,
A.
Anelli
,
D.
Giofré
,
S.
Klees
,
J.
Behler
, and
M.
Ceriotti
, “
Automatic selection of atomic fingerprints and reference configurations for machine-learning potentials
,”
J. Chem. Phys.
148
,
241730
(
2018
).
96.
Y.
Sugita
and
Y.
Okamoto
, “
Replica-exchange molecular dynamics method for protein folding
,”
Chem. Phys. Lett.
314
,
141
151
(
1999
).
97.
J.
Behler
, “
Atom-centered symmetry functions for constructing high-dimensional neural network potentials
,”
J. Chem. Phys.
134
,
074106
(
2011
).
98.
V.
Kapil
,
M.
Rossi
,
O.
Marsalek
,
R.
Petraglia
,
Y.
Litman
,
T.
Spura
,
B.
Cheng
,
A.
Cuzzocrea
,
R. H.
Meißner
,
D. M.
Wilkins
 et al, “
i-PI 2.0: A universal force engine for advanced molecular simulations
,”
Comput. Phys. Commun.
236
,
214
223
(
2019
).
99.
A.
Singraber
,
J.
Behler
, and
C.
Dellago
, “
Library-based lammps implementation of high-dimensional neural network potentials
,”
J. Chem. Theory Comput.
15
,
1827
1840
(
2019
).
100.
The PLUMED consortium
, “
Promoting transparency and reproducibility in enhanced molecular simulations
,”
Nat. Methods
16
,
670
673
(
2019
).
101.
G. A.
Tribello
,
M.
Bonomi
,
D.
Branduardi
,
C.
Camilloni
, and
G.
Bussi
, “
PLUMED 2: New feathers for an old bird
,”
Comput. Phys. Commun.
185
,
604
613
(
2014
).
102.
M.
Ceriotti
,
G.
Bussi
, and
M.
Parrinello
, “
Colored-noise thermostats à la carte
,”
J. Chem. Theory Comput.
6
,
1170
1180
(
2010
).
103.
A.
Grossfield
, WHAM: the weighted histogram analysis method, http://membrane.urmc.rochester.edu/wordpress/content/wham.
104.
W.
Wang
,
B.
Ji
, and
Y.
Zhang
, “
Chalcogen bond: A sister noncovalent bond to halogen bond
,”
J. Phys. Chem. A
113
,
8132
8135
(
2009
).
105.
L. P.
Wolters
,
P.
Schyman
,
M. J.
Pavan
,
W. L.
Jorgensen
,
F. M.
Bickelhaupt
, and
S.
Kozuch
, “
The many faces of halogen bonding: A review of theoretical models and methods
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
4
,
523
540
(
2014
).
106.
C.
Wang
,
D.
Danovich
,
Y.
Mo
, and
S.
Shaik
, “
On the nature of the halogen bond
,”
J. Chem. Theory Comput.
10
,
3726
3737
(
2014
).
107.
J.
Řezáč
and
A.
de la Lande
, “
On the role of charge transfer in halogen bonding
,”
Phys. Chem. Chem. Phys.
19
,
791
803
(
2017
).
108.
B.
Inscoe
,
H.
Rathnayake
, and
Y.
Mo
, “
Role of charge transfer in halogen bonding
,”
J. Phys. Chem. A
125
,
2944
2953
(
2021
).
109.
F.
Nunzi
,
D.
Cesario
,
F.
Tarantelli
, and
L.
Belpassi
, “
Halogen bond interaction: Role of hybridization and induction
,”
Chem. Phys. Lett.
771
,
138522
(
2021
).
110.
J.
Thirman
,
E.
Engelage
,
S. M.
Huber
, and
M.
Head-Gordon
, “
Characterizing the interplay of Pauli repulsion, electrostatics, dispersion and charge transfer in halogen bonding with energy decomposition analysis
,”
Phys. Chem. Chem. Phys.
20
,
905
915
(
2018
).
111.
S. M.
Huber
,
E.
Jimenez-Izal
,
J. M.
Ugalde
, and
I.
Infante
, “
Unexpected trends in halogen-bond based noncovalent adducts
,”
Chem. Commun.
48
,
7708
7710
(
2012
).
112.
T.
Clark
,
M.
Hennemann
,
J. S.
Murray
, and
P.
Politzer
, “
Halogen bonding: The σ-hole
,”
J. Mol. Model.
13
,
291
296
(
2007
).
113.
N.
Biot
and
D.
Bonifazi
, “
Chalcogen-bond driven molecular recognition at work
,”
Coord. Chem. Rev.
413
,
213243
(
2020
).
114.
P. C.
Ho
,
P.
Szydlowski
,
J.
Sinclair
,
P. J. W.
Elder
,
J.
Kübel
,
C.
Gendy
,
L. M.
Lee
,
H.
Jenkins
,
J. F.
Britten
,
D. R.
Morim
, and
I.
Vargas-Baca
, “
Supramolecular macrocycles reversibly assembled by Te…O chalcogen bonding
,”
Nat. Commun.
7
,
11299
(
2016
).
115.
A.
Daolio
,
P.
Scilabra
,
M. E.
Di Pietro
,
C.
Resnati
,
K.
Rissanen
, and
G.
Resnati
, “
Binding motif of ebselen in solution: Chalcogen and hydrogen bonds team up
,”
New J. Chem.
44
,
20697
20703
(
2020
).

Supplementary Material