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

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 are the most prevalent non-covalent interactions occurring in a solute-Cl$^-$-THF mixture. 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.


I. INTRODUCTION
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][4][5][6][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 noncovalently 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 has thus been a long-lasting goal of computational chemistry. [9][10][11][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 Cockroft and coworkers further stress the difficulties of rec-onciling 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 corrections [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34] or with reducing the cost of wavefunction-based methods, [35][36][37][38][39][40][41] 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,[42][43][44] 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. 42,43 As the influence of solvation on the dispersion forces is still being actively investigated, 45 proper care should be placed on the unambiguous quantification and characterization of NCIs in solution when implicit solvent is used. 46,47 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, 48 SIBFA, 49 NEMO, 50 GEM 51 ) 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, 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 time scales. The challenge associated with achieving statistically converged simulations at quantum chemical accuracy thus remains present.
Stemming from the advances in Neural Networks, 52 high dimensional neural network-based potentials (NNPs) 53,54 trained on 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, 55,56 molecules in the gas phase, liquids 57,58 or molecules on the interfaces. 59 Their application in the simulations of molecular systems in solution however remains scarce. 60,61 Recently, some of us demonstrated the advantages of applying NNPs to go beyond the semiempirical description of reactions involving hydrogen peroxide, acid and phenol in an amphiphilic environment. 62 The neural network-based simulations showed the importance of relying upon accurate potentials and achieving long timescale 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 63 with increasing applicability for the design of organocatalyts, anion transport through lipidic membranes, supramolecular architectures and functional materials. [64][65][66][67] 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 nonnegligible 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 Clanion. The association constant for this complex in tetrahydrofuran (THF) was recently measured by means of UV-vis absorbance spectroscopy. 68 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 Polarizable Continuum Model (PCM). 68 Apart from the interaction between tellurium and the Clanion, 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. 62 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 study byŘezáč et al. [69][70][71] and Goerigk and coworkers, 72 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 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. 72 Here we choose PBE0 73 -D3BJ 24 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 Cland 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.

II. METHODS AND COMPUTATIONAL DETAILS
a. Overview The workflow used herein (Fig. 2) is a variant of the baselined (BNNP) and direct (DNNP) NNPs training procedure published previously. 62 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 74 Hamiltonian. Energies and forces are then recomputed at a more accurate level (PBE0 73 with D3BJ correction 24 ) that serve as reference data. Behler-Parrinello NNPs 53,75 are then trained to reproduce the reference data either directly or using a difference between the baseline and reference method. 76 Molecular dynamics together with enhanced sampling techniques are run with each NNPs or in a combined multiple time step (MTS) scheme. 77 Details on the simulation protocol, preparation and selec- b. Baselined NNPs using the GFN0-xTB semiempirical potential The superior accuracy and stability of kernel-based models and NNPs, trained from a baselined electronic structure level (i.e., ∆-ML 76 ) has been largely demonstrated. 62,[78][79][80][81][82][83] 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 74 goes beyond the limits of pairwise parametrization imposed by sev-eral semiempirical approaches, such as the related DFTB. 84 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 Atomistic simulation environment (ASE) 85 . 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 modifiction 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).
c. Training set construction The system is composed of one benzotelluradiazole solvated in 45 tetrahydrofuran (THF) molecules and one Clanion, 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. 86 The solute molecules were kept frozen in the geometry optimized at the M06-2X/def2-TZVP level during the equilibration run. Steered Molecular Dynamics 87 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) 88 simulations using the GFN0-xTB potential. Each window was simulated for 50 ps using the Abin code. 89 The training set was constructed from 2 000 structures selected from the US simulation with Farthest Point Sampling (FPS). 90 Additional sets of 300 and 1500 structures respectively were generated by performing high temperature (1000 K) and two Replica Exchange MD (REMD) 91 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 Cland benzotelluradiazole, the second REMD simulation was initiated with Clplaced 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 simulations with the trained BNNP. The complete set thus contained 5 000 structures with their respective DFT energies and forces (2 000 from GFN0-xTB US, 1 800 from GFN0-xTB REMD, 1 200 from US using the BNNP trained with the previous 3 800 structures). Such a set ensured configurational diversity and sufficient coverage of all the relevant regions of the potential energy landscape.
d. Atomic Symmetry Functions (ASFs) definition We used radial (G 2 ) and angular (G 4 ) Behler-Parrinello atomic symmetry functions (ASFs) 53,92 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. 90 to select the most relevant fingerprints for the system. Previous studies 62,90 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 SM (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 (resp. 1792) symmetry functions for the whole system.
e. NNPs training The DNNP and BNNP were trained using the Behler-Parrinello neural network potential architecture implemented in the n2p2 version 2.0.3. 75 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 2 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. 75 The complete training data set of 5 000 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. An ensemble of 3 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.
f. Molecular dynamics 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 93,94 and Plumed 2.6.2 95,96 , 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. 97 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 98 (GLE) with the same parameters as in the Ref. 62. 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.
g. Umbrella sampling Potential of Mean Force (PMF) was 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 ns and 9.5 ns respectively. For the sake of comparison, we simultaneously ran the same simulation with 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 99 . The convergence of the PMF was verified by running the US in the reverse direction (see Fig. S7-S10). 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.

III. RESULTS AND DISCUSSION
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 Clsolvated 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. Fig. 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 meV/Å/atom and 80 meV/Å/atom for DNNP and BNNP, respectively.

Validation of the NNPs
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 Fig. S4 -S6 in SM). The reference and predicted forces show a good correlation, although larger er-rors 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.
All methods yield stable trajectories during short nonbiased MD simulations (i.e., 50 ps). DNNP however triggers slightly more extrapolation warnings than BNNP. Importantly, DNNP with the MTS set up 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-Cldistance above 3.5 Å.
The predictions of forces (Fig. 4, A and B) for both NNPs are in a very good agreement with the PBE0-D3BJ reference, showing a RMSE of 72 meV/Å/atom 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 consistent with the test set (RMSE 0.52 meV/atom), 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 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. 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 the next section, 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.

Role of NCIs in benzotelluradiazole-Classociation
The strength of the interaction between benzotelluradiazole and Clis evaluated via the PMF (see Fig. 5A) reconstructed along the Te-Cldistance. The MTS-NNP predicts a favorable association of benzotelluradiazole and Clwith 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-Cldistance coincides with a N-Te-Clbond 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-Cldistance. At the longer distances, the angle deviates from the liner arrangement and the mode of interaction changes. chalcogen bond is replaced by a different 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-Clcomplex and the solvent-separated structure (vide infra).
The PMF profiles modeled with other potentials (see Fig.  5A) substantially differ from the MTS-NNP one. The profile at the bare GFN0-xTB level is shifted towards lower Te-Cldistances. This behavior is likely caused by the wrong description of the repulsive and electrostatic part of the potential, a limitation which is well known as GFN0-xTB does not use self-consistent evaluation of the charges. 74 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-Clinteraction.
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 then 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. [100][101][102][103][104][105] While the exact role of charge transfer in σ -hole interactions is still under debate, 106-108 the original AMOEBA force field expression suffers from a lack of explicit charge transfer description, 48 which could explain its collapse in the long-range. An adjustment of the multipoles (including higher order terms) and Van der Waals terms 109 (by softening the VdW radii 110,111 ) could be a route to retrieve the missing charge transfer term. 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 the previous section. Finally, the relaxed scan in PCM reaches a plateau in the Te-Cldistance 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-Clangle (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 Clis lost with the angle dropping below 160 • similarly to MTS-NNP. In both the DNNP and AMOEBA simulations, the chalocogen bond becomes loose already around the Te-Cldistance 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 Clstop 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 af- 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 Cldensity located close to the position of σhole. In the longer Te-Cldistances, interaction towards Clweakens and the affinity towards the oxygen atoms of THF is amplified instead. The windows at Te-Cldistance of 4.75 and 5.25 Å finally demonstrate complete replacement of Te and Clchalcogen bond by directional chalcogen bonds towards 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. S15 in the SM). This behavior further stresses the importance of Te-Clchalcogen bonds in the short range with its gradual replacement by Te-O chalcogen bonds at the Te-Cldistances above 4.00 Å, which explains the observed trends in the N-Te-Clangle distributions. The density distribution at the 5.25 Å window depicts the transition from the Te-Cl complex to the fully solvated Clcorresponding to the maximum of the PMF profile located around 5.2 Å (see Fig. 5A).
At 3.00 Å, the density distributions with DNNP and AMOEBA already differ with those obtained with the MTS-NNP. They promote an interaction with Clbut also one with a THF oxygen atom, which is not present using MTS-NNP. At Te-Cldistances larger than 3.00 Å, the directional interaction with Clis entirely replaced by the Te-O chalcogen bond interactions. With AMOEBA, the interaction between Te and Cldecreases 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 are rather nonphysical exhibiting a too strong affinity for the solvent and to some extent too strong anion-π 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 towards both the Clanion and the THF oxygen atom. Benzotelluradiazole therefore interacts not only with the Clbut 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, 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. 112 For instance, similar tellurium-containing compound crystallizes in a form of Te-Obonded macrocycles with a stoichiometry dictated directly by the solvent effects. 113 Interestingly, THF molecules were even found inside the formed cavities further advocating for their interaction with tellurium. Chalcogen bonds be-FIG. 6. 3D Density distribution of oxygen (red contour) and Cl -(green contour) atoms around the 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. tween selenium containing species and oxygen and nitrogen atoms of solvent molecules were also confirmed by NMR measurements. 114 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.

IV. CONCLUSIONS
This computational work investigates the persistence of chalcogen bonds between benzotelluradiazole and the Clanion 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, 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 system incorporating several types of competing noncovalent interactions.
The application of this method brings important chemical insight into 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 towards Cland the solvent oxygen atoms. While the Te-Clchalcogen bond is a dominant non-covalent interaction in the short range region up to Te-Cldistances of 4.00 Å, the interaction between Te atom and oxygen of THF is prevailing in larger distances. This behavior coincides with the loss of directionality of the Te-Clinteraction and highlights the nonnegligible 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 which helps rationalizing and harvesting the nature and behavior of chalcogen bonds in solution. Similar ap-proaches would be useful for assessing 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.

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

ACKNOWLEDGMENTS
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, No. 514149), and F. C. and C. C. from the European Research Council (ERC, Grant Agreement No. 817977) within the framework of European Union's H2020. The National Center of Competence in Research (NCCR) "Sustainable chemical process through catalysis (Catalysis)" of SNSF is acknowledged for financial support of R. L. The authors thank Raimon Fabregat for the assistance with graphical material, Daniel Hollas for the help with Abin set up and Kevin Rossi for scientific discussions.

DATA AVAILABILITY STATEMENT
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 will be freely available on the Materials Cloud. II. AMOEBA parametrization S3 A. Parameter validation S10

SUPPLEMENTARY MATERIAL
III. Symmetry functions selection S11 IV. NNP validation S12 A. Training errors of the test set S12 B. Forces parity plots per element S13 V. Potential of Mean Force: Convergence verification S14 A. AMOEBA S14 B. DNNP S15 C. MTS-NNP S17 VI. Supplementary/Supplementary plots S18 A. Comparison of implicit solvent models S18 B. Extrapolation warnings S19 C. N-Te-O angle evolution S22 S2 I. COMPUTATIONAL DETAILS A. Electronic structure description GFN0-xTB S1 with periodic boundary conditions as implemented in xTB 6.2 S2 is used both for preliminary sampling and as a baseline potential for NNP simulations. GFN0-xTB is a non-iterative member of the extended tight-binding family of Hamiltonians developed by Grimme et al., which makes the approach significantly faster compared to self-consistent semiempirical methods. Reference energies and forces are computed at the PBE0-D3BJ S3,S4 level of theory as implemented in CP2K 6.1. S5 We use the TZV2P-MOLOPT basis set for H, C, N, O, F and Cl, S6 and a DZVP-MOLOPT-SR for Te with a Goedecker-Teter-Hutter pseudopotential (GTH PBE) S7 for all elements. The plane-wave cutoff is set to 700 Ry with a relative cutoff of 70 Ry. All computations employ a Coulomb operator truncated at R = 9 Å and the auxiliary density matrix method with a cpFIT3 fitting basis set for H, C, N, O, F and Cl and FIT9 for Te. S8 PBE-D3BJ wave functions were used as the initial guess for the PBE0-D3BJ computations to speed up the convergence. S4,S9 PBC are taken into account with a cubic box of 18.428 Å and a total charge of -1.0e for the simulation cell. The free energy curves in implicit solvent are computed using the PCM, S10 SMD, S11 and COSMO-RS S12 models. Structures for the computations are generated by a relaxed scan along the Te-Cldistance in a range from 2.0 to 7.75 Å with a step of 0.25 Å at the PBE0-D3BJ/aug-cc-pVTZ level with the solvent approximated by PCM and SMD, respectively, as implemented in Gaussian 16. S13 Free energy curves for the PCM and SMD models are constructed from 26 structures computed at PBE0-D3BJ/aug-cc-pVTZ single points with the respective solvent model using the quasi-harmonic correction to entropy S14 and enthalpy S15 as implemented in the GoodVibes program. S16 Energy estimates and COSMO-RS solvation energies are computed at the PBE0-dDsC/TZ2P level, as implemented in the ADF package , S17 using the same geometries as in the PCM computations. Free energy corrections are obtained using the PBE0-D3BJ/aug-cc-pVTZ level along with the quasi-harmonic correction to entropy and enthalpy as in previous computations.

B. Classical force fields simulations
To compare the NNPs performance with state-of-the-art force fields, we perform the umbrella sampling simulations with polarizable (AMOEBA) force field. Since the suitable AMOEBA parameters for THF and benzotelluradiazole do not exist, both systems were parametrized before proceeding with the simulations using the following protocol and mixture: one benzotelluradiazole molecule, 45 THF molecules, one Clanion and one Na + cation. Na + , which can influence the PMF, is constrained at a distance of 10 Å with a force constant of 100 kcal/mol cdotÅ 2 applied in all simulations. The molecular dynamics simulations are carried out using the Tinker 8.8 software. S18 The equations of motion are integrated by the RESPA multiple-time-step integrator with a time step of 2 fs. The Smooth Particle Mesh Ewald algorithm is employed with a grid dimension of 24 Å × 24 Å × 24 Å to calculate forces in the reciprocal space. The Ewald cutoff is set to 7 Å and the Van der Waals cutoff to 9 Å. The NVT thermodynamic ensemble is used with the temperature maintained to 298.15 K by the Bussi thermostat. For the AMOEBA force field, the convergence criteria for multipoles are set to 10 −5 . We define the same windows and spring constants as for the NNPs simulations, but the simulation time per window is increased to 10 ns. The first two nanoseconds are kept for equilibration and not considered in the analysis. Finally, the PMF is reconstructed and the statistical error estimated in both cases using the WHAM implementation of Alan Grossfield. S19 S3 II. AMOEBA PARAMETRIZATION AMOEBA parameters are generated for the benzotelluradiazole and THF molecules according to the procedure described by Ren et al. S20 with the following main steps: 1. The geometry of an initial structures are optimized at the MP2/6-311G** level using the Gaussian 09 software.
2. The distributed Multipole Analysis is performed to generate an initial set of atomic multipoles.
3. The multipoles are refined by fitting to the molecular electrostatic potential obtained from a single point calculation at the MP2/aug-cc-pvtz level of theory.
4. Van der Waals, bond, angle, dihedral, and out of plane bending parameters are assigned using the Tinker valence program.