Assessing the persistence of chalcogen bonds in solution with neural network potentials

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.⁸ 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.¹³ 

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.¹⁶ Despite the intense efforts associated with developing a posteriori dispersion corrections^17–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,⁴⁶ 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,⁴⁹ SIBFA,⁵⁰ NEMO,⁵¹ and GEM⁵²) 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),⁵³ 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.⁶⁰ 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.⁶³ 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 interaction⁶⁴ 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.⁶⁹ 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).⁶⁹ 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.

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.⁶³ 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,⁷³ 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.⁷³ Here, we choose PBE0⁷⁴-D3BJ²⁵ as the reference level. We also interface the i-PI dynamic driver with the xTB code⁷⁵ 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.⁶³ 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-xTB^75,76 Hamiltonian. Energies and forces are then recomputed at a more accurate level (PBE0⁷⁴ with D3BJ correction²⁵) that serve as reference data. Behler–Parrinello NNPs^54,77 are then trained to reproduce the reference data either directly or using a difference between the baseline and reference method.⁷⁸ Molecular dynamics together with enhanced sampling techniques are run with each NNPs or in a combined multiple time step (MTS) scheme.^79–81 

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., Δ-ML⁷⁸), 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 Hamiltonians⁷⁵ goes beyond the limits of pairwise parameterization imposed by several semiempirical approaches, such as the related Density Functional based Tight Binding (DFTB).⁸⁹ 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.⁹⁰ 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.⁹¹ The solute molecules were kept frozen in the geometry optimized at the M06-2X/def2-TZVP level during the equilibration run. Steered molecular dynamics⁹² 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)⁹³ simulations using the GFN0–xTB potential. Each window was simulated for 50 ps using the Abin code.⁹⁴ The training set was constructed from 2000 structures selected from the US simulation with Farthest Point Sampling (FPS).⁹⁵ Additional sets of 300 and 1500 structures, respectively, were generated by performing high temperature (1000 K) and two Replica Exchange MD (REMD)⁹⁶ 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 (G₂) and angular (G₄) 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 studies^63,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.⁷⁷ 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.⁷⁷ 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 equation¹⁰² (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 Ų 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 Ų. 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.¹⁰³ 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.

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 × 10⁴ 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.

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).

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.⁷⁵ 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,⁴⁹ 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 Å).

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.¹¹³ For instance, similar tellurium-containing compound crystallizes in a form of Te–O⁻ bonded macrocycles with a stoichiometry dictated directly by the solvent effects.¹¹⁴ 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.¹¹⁵ 

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.