2 Computational Details

All theoretical calculations were performed via the Gaussian 16 suite of programs[31]. The dispersion-corrected density functional theory (DFT-D3) methods were used to study the systems we chose. The structures of all the monomers and the complexes included in this study were fully optimized at the B3LYP-D3/6-311+G(d,p) level of theory. The vibrational frequencies of the optimized structures were carried out at the same level and used to obtain the zero-point vibrational energy (ZPE). All stationary points were characterized as a local energy minimum on their potential energy surfaces by verifying that all the corresponding frequencies were real. The interaction and the binding energy were obtained using the double-hybrid functional method B2PLYP-D3 with jun-cc-pVTZ as the basis set. Both B3LYP-D3 and B2PLYP-D3 methods used in present study are all DFT methods including a version of Grimme’s D3 dispersion model (DFT-D3)with BJ-damping[32, 33]. Due to the inclusion of both long and medium-range dispersion effects, these DFT methods have been widely used and show generally satisfactory performance in the calculations of non-covalent interactions in similar systems[34-36]. Jun-cc-pVTZ is a modified version of aug-cc-pVTZ by removing itsf -type diffusion basis function. This modification significantly facilitates the convergence of the geometry optimization process while the loss of accuracy is trivial[37].
The interaction energy, denoted by ΔE int, is defined as the difference between the complex and the sum of energies of monomer whose geometries come from the optimized structure of the complex. As for the binding energy, denoted by ΔE bin, the energy of the monomer used as the reference point is derived from the energetic minima of the isolated monomers. The deformation energy, denoted byE def, is calculated as the difference between ΔE bin and ΔE int of the complexes, which is the energy difference between the monomers in their equilibrium geometries and at their relaxed geometries in the complexes, and its value is positive since the complexation results in changes of the structure of monomers[38, 39]. The basis-set superposition error (BSSE) was obtained using the counterpoise correction method proposed by Boys and Bernardi[40] to correct both ΔE int and ΔE bin. The total interaction energy in the ternary system, denoted by ΔE (ABC) or ΔE total, is calculated as the difference between the energy of ternary complex and the energy sum of the monomers which is frozen in the geometry of the ternary complexes. The cooperative effects can be explained by the many-body interaction analysis [41, 42]. The two-body terms, denoted by ΔE (AB), ΔE (BC), and ΔE (A-C)far, are defined as the difference between the energy of each binary system and the energy sum of the monomers, which come from the geometry of ternary complexes. The interplay between the two interactions in the ternary can be estimated with the cooperative energy, denoted byE coop, which are obtained by the following formulas:E coopE (ABC)-ΔE (AB)-ΔE (BC)-ΔE (A-C)far. This methodology has been widely used to study how the different interactions affect each other in the complexes[43, 44].
Molecular electrostatic potentials (MEPs) maps of the isolated monomer were calculated with the Wave Function Analysis-Surface Analysis Suite (WFA-SAS) program[45] at the B3LYP-D3/6-311+G(d,p) level on the 0.001 a.u. contour of electronic density to locate the position and the value of the minima and maxima. To unveil the nature of the intermolecular interactions, the atoms-in-molecules (AIM) topological analysis were carried out with the use of the Multiwfn package[46] based on the wavefunctions generated at the B3LYP-D3/6-311+G(d,p) level. The following characteristics of bond critical points (BCPs) corresponding to the intermolecular interaction were analyzed including the electron density (ρBCP), the Laplacian of the electron density (∇2ρBCP), and the electronic energy density (HBCP), which is decomposed into electronic kinetic energy density (GBCP) and the electronic potential energy density (VBCP). The Visual Molecular Dynamics (VMD) program[47]was applied to visualize the weak interaction. The natural bond orbital (NBO) analysis was performed via the NBO 3.1 program[48]implemented in Gaussian 16 to estimate the orbital-orbital interactions as well as charge transfer (CT).