In this study, we have developed a python script code aimed at automatizing the procedure of computing the interaction energy by quantum mechanical calculations followed by the analysis of the different interaction energy contributions for a large set of complexes. Specifically, the developed code is applied to the investigation of 250 complexes formed by the glycine amino acid, whose geometry has been taken from Ropo et al. [44], and 250 drugs selected from the ZINC database [43]. The general workflow of the procedure driven by the code, represented in Fig. 1, consists of different steps: (i) choice of the initial intermolecular distance between the two molecules of the complex; (ii) determination of the orientation between the molecules of the complex for which the interaction is the most favourable one; (iii) calculation of the potential energy curve for the most favourable orientation at a low level of theory; (iv) high-level energy calculations for five geometries selected around the low-level energy minimum; and (v) EDA employing the high-level electron densities for those five geometries. A detailed explanation of these steps will be given for an specific complex formed by the glycine amino acid and the 3-(2-methoxyphenoxy)propane-1,2-diol drug (see Fig. 1A). From now on, they will also be called molecule 1 and molecule 2, respectively, for simplicity.
In order to avoid a very large repulsive or a very low attractive interaction between the molecules of the complex, an appropriate initial intermolecular distance between the molecules has to be chosen in step (i). This distance is defined between the geometric centres of the two molecules, and it is calculated in the following manner. First, the interatomic distances between all the atoms within each of the molecules is calculated. Then, the largest interatomic distance for each molecule is selected (\({d}_{1}\) and \({d}_{2}\)). These two distances are then employed as the diameters of two generated spheres centred on the geometric centres of the molecules (see Fig. 2A). In order to avoid overlaps between atoms of different molecules, each sphere is enlarged by adding to its radius twice the van der Waals radius of the largest atom of the corresponding molecule. In the case of the example complex shown in Fig. 2A, for molecule 1 the largest distance (or diameter d1) corresponds to the distance between the oxygen atom forming the C = O double bond of the carboxyl group and the hydrogen atom of the amino group, which is 4.33 Å. Then, twice the van der Waals radius of the carbon atom (1.70 Å), which is the largest atom of the molecule, is added to sphere radius. In the case of molecule 2, the maximum distance is found between two hydrogen atoms (see Fig. 2A), and the atom with the largest van der Waals radius is again the carbon atom. Finally, the distance \(R\) between the geometrical centres of the two molecules is computed as \(R={r}_{1}+{2r}_{1,vdW}+{r}_{2}+{2r}_{2,vdW}\).
Then, for the previously determined intermolecular distance, both molecules are randomly rotated until obtaining 15 different orientations. For each orientation, a single point energy at PM6 level of theory (low level) using the Gaussian09 software[45] has been performed, and the lowest-energy orientation is chosen (step ii). Therefore, by the combination of the steps (i) and (ii), a geometry of the complex that is likely very close to the global minimum of the low-level method is chosen. For this favourable geometry, the potential energy curve is computed at PM6 level in the step (iii) by performing energy calculations at different intermolecular distances, from which the energy of the isolated molecules is subtracted to obtain the interaction energy. Specifically, the intermolecular distance is shortened and elongated in steps of 0.1 Å until the interaction energy is higher than 20 kcal/mol, in the repulsive region, and smaller than 10% of the absolute value of the interaction energy at the energy minimum, in the region of large intermolecular distances. In this way, it is assured that the PM6 potential energy curve considers the intermolecular distances around the energy minimum, which could be thermally populated in a dynamic scenario because they present a favourable, or at least not highly repulsive, interaction energy.
Once the potential energy curve is computed at low level of theory, the interaction energy for five geometries around the low-level energy minimum is refined by higher-level computations. Specifically, the low-level energy minimum geometry, two consecutive geometries at shorter distances, and two geometries taken every two steps at larger distances are selected from the PM6 curve. For those geometries, the input files for computing the interaction energy, step (iv), with the M06-2X functional and 6-31G* basis set using the Gaussian09 software [45] are created with the MoBioTools kit [46, 47]. In the interaction energy calculation, the basis set superposition error is corrected by applying the counterpoise correction [48]. Finally, in the step (v), the electron densities of the complex and the isolated molecules are read by the EDA-NCI program [41, 42, 49] to obtain the different interaction energy contributions, namely, Pauli, electrostatic, dispersion and induction interactions. The sum of the last two terms constitutes the polarisation energy.
The energy contributions are represented in Fig. 2B for the five intermolecular distances of the example complex. As can be seen, for this particular case, the electrostatic, polarisation, dispersion and induction contributions go to more negative values, while the Pauli repulsion term becomes more positive, when the molecules get closer. Overall, the DFT total interaction energy presents a minimum of -7.36 kcal/mol at 3.49 Å. In this case, the position of the DFT minimum coincides with that of the PM6 one, although this is not the case for all the complexes investigated. The EDA calculation also provides the electron density polarisation induced by the intermolecular interaction. Thus, the functional groups that modify its electron density because of the interaction can be easily identified. Figure 2C displays the polarisation density for the glycine/3-(2-methoxyphenoxy)propane-1,2-diol complex, where it can be seen that there exist an electron density migration from the amino acid to the alcohol group of the drug. All the steps described above for this sample case have been performed for 250 complexes and the results are shown in the Supporting Information.