On the effects of the basis set superposition error on the change of QTAIM charges in adduct formation. Application to complexes between morphine and cocaine and their main metabolites

Abstract

Eleven different complexes were found between heroin (or its main metabolite, morphine) and cocaine (or its main metabolite, ecgonine methyl ester) in a B3LYP computational study. Ten of these complexes display intermolecular hydrogen bonds, while only one is a stacked structure. All of them display positive complexation energies in the gas phase that turn into negative values in aqueous solution according to PCM calculations. Complexation energies become even more negative when explicit solvation water molecules are taken into account. Thus, complexes between these compounds could be present in biological media. The electron density of the 11 complexes was analysed within the QTAIM framework distinguishing three terms for every atomic property: (i) geometry distortion; (ii) BSSE estimated by extending the counterpoise (CP) method; and (iii) binding. We notice that: (i) geometry distortion effects are basically localised in the atoms involved in intermolecular bonds (interaction sites); (ii) counterpoise corrections for atomic properties are very small; and (iii) binding CP-corrected effects spread throughout the whole complex. The later are the most significant and indicate the important role played by hydrogens outside the interaction sites as electron density sources and sinks in complex formation.

---

1. Introduction

Cocaine and morphine were consumed millennia ago. In the late Bronze Age, leaves of Erythroxylum coca were consumed by the aboriginal people of South America. This tradition was based on its nutritional richness.¹ In Byzantium times, the medical tradition for the treatment of pain was through providing powder or drops of Papaver somniferum latex.² In modern times, the traditional chewing of coca leaves can be still found within the Andes Mountains area. In contrast, the consumption of the main phytochemical components of these plants is a worldwide affecting health problem. The consumption of cocaine, morphine and their derivatives develops drug dependence^3,4 through a mechanism which is not fully understood yet. It is known the cocaine side effects occur through the modification of the biochemical activity of the dopaminergic, serotonergic and noradrenergic pathways in the brain,^5,6 while those of morphine occur through the dopaminergic pathway, modulated by the μ-opiate and GABA receptors.⁷ The combined consumption of cocaine and heroin, (as in the mixture known as speedball) shows a synergic interaction which increases dramatically the dopamine release in the synaptic cleft,⁸ but not much more is known about the higher addiction observed for such kind of combinations.^4,6,9 Therefore, there are new different lines of research that are bringing more information to get insight on this matter.¹⁰ It is not possible to discard that complexes with high biochemical activity can be formed between both drugs or between crossed metabolites.

The first aim of this paper was to study the stability of the complexes formed between heroin (or morphine) with cocaine (or ecgonine methyl ester). To this end, we have performed a conformational analysis for all the complexes that can be formed between the four different pairs of molecules and evaluated their binding energy using the a posteriori counterpoise (CP) correction of the Basis Set Superposition Error (BSSE), proposed by Boys and Bernardi.¹¹ The electron density computed for the conformers was analysed by means of the Quantum Theory of Atoms in Molecules (QTAIM).¹² This work also explores how the computed electron density reorganization experienced by monomers upon complexation is influenced by BSSE. This analysis has been done introducing atomic contributions for the CP correction. Thus, atomic electron density reorganizations were split into terms due to binding CP-corrected interactions, geometry deformation, and BSSE.

2. Computational details

Heroin, H, and morphine, M, share (5a,6a)-7,8-didehydro-4,5-epoxy-17-methylmorphinan structure (Fig. 1) and only differ in the functional groups at positions 3 and 6. In both positions, heroin has two acetyl groups while those of morphine are hydroxyls. Cocaine, C, and ecgonine methyl ester, E, share methyl-(1R,2R,5S)-8-methyl-8-azabicyclo[3.2.1]octane-2-carboxylate structure (Fig. 1) and only differ in the functional group attached at position 3. It is a benzoyloxy group in cocaine and a hydroxyl group in ecgonine methyl ester. Formation energies for diverse complexes between heroin (or morphine) with cocaine (or ecgonine methyl ester) were calculated. The electron density reorganization involved in the complexation processes was analyzed using QTAIM and considering the effect of BSSE through the method proposed herein: a QTAIM partitioning method for CP correction.

Fig. 1 Molecular structure and atom numbering for the 8-azabicyclo[3.2.1]octane-2-methyl-ester structure, which is shared by cocaine and ecgonine methyl ester (left), and for the (5a,6a)-7,8-didehydro-4,5-epoxy-17-methylmorphinan structure (right), shared by heroin and morphine. In heroin, O1′ and C1′ correspond to the atoms in the substituent attached to O3, while O1′′ and C1′′ refer to atoms in the group attached to O6.

2.1. DFT calculations

All the calculations for the complexes studied here and their corresponding monomers were carried out using the GAUSSIAN 09 software revision D.01.¹³ The structures of these complexes were fully optimized with the CP correction using the density functional theory (DFT). The correlation functional used in these calculations was that proposed by Lee, Yang and Parr¹⁴ together with the Becke's three parameter exchange functional (B3LYP¹⁵) using the 6-31++G(d,p) 6d basis set. The vibrational frequencies and zero-point energy (ZPE) corrections were computed at the same DFT level, in order to ensure that the optimized structures were true minima.

In order to evaluate the stability of complexes formed between cocaine and heroin and/or their metabolites. The most stable complex, HE1, (Fig. 2) was modeled with PCM and with the insertion of three water molecules in the surroundings of the intermolecular hydrogen bond. To this end, solvent effects were considered by using the integral equation formalism polarizable continuum model (IEF-PCM).¹⁶ It treats solvation as a perturbation, defined in terms of the surface apparent charges, which, in turn, depend on the solute charge distribution estimated from the solute wave function.

Fig. 2 QTAIM molecular graphs of all the complexes study herein, organized by increasing CP-corrected complexation energies. The complex structures were calculated at B3LYP/6-31++G(d,p) 6d in gas phase.

Furthermore, the van der Waals interactions were also modeled. To this end, the 6-31++G** 6d basis set and the methods of the density functional theory B3LYP and M06 were used to perform the following calculations of the most stable complex formed between heroin and ecgonine methyl ester. Four different optimization procedures were done: (i) conventional optimization, (ii) optimization including the van der Waals interactions through the potential of Grimme,¹⁷ (iii) introducing the correction to the basis set superposition error by the counterpoise method and (iv) the sum of the inclusion of Grimme potential with counterpoise method.

2.2. QTAIM calculations

The QTAIM analysis of the electron density, ρ(r), for all of these systems was carried out with the AIMALL software.¹⁸ This analysis allows obtaining the atomic basins and their corresponding properties, such as the electron population, N(Ω), and energy, E(Ω). Integration errors were estimated through the differences between molecular properties and those obtained by summation of the atomic ones, N − ∑N(Ω) and E − ∑E(Ω). Their absolute values were always smaller than 2.5 × 10⁻³ au and 3.4 kJ mol⁻¹, respectively (ESI material, Table S1†). The absolute values obtained for L(Ω) were always smaller than 1.5 × 10⁻³ au.

2.3. A QTAIM partitioning method for the BSSE

Basis set superposition error (BSSE) has been long time recognised as a limitation affecting the calculation of interaction energies using a finite basis set. In these calculations, the energy of the complex is lowered because each of the interacting fragments is described with more basis functions in the complex than when it is isolated. In order to correct this error, Boys and Bernardi proposed the CP correction method.¹¹ If this (or other) BSSE correction is not employed (as could be done in the infinite basis set limit), the calculated complexation energy, Δ_cE, for two interacting systems, A and B, would be given by eqn (1), where each E^Z_Y(X) symbol denotes the energy of the X fragment in the Y geometry when computed with the Z basis set. The corresponding deformation, Δ_defE, and binding, Δ_bE, energies due, respectively, to: (i) geometrical distortion of the molecules, A and B, in order to become bonded in the formed complex; and (ii) intermolecular interaction once both fragments are in the geometry of the complex, are given by (2) and (3).

Δ_cE(AB) = E^AB_AB(AB) − E^A_A(A) − E^B_B(B) (1)

Δ_defE(A,B) = [E^A_AB(A) − E^A_A(A)] + [E^B_AB(B) − E^B_B(B)] (2)

Δ_bE(A,B) = E^AB_AB(AB) − E^A_AB(A) − E^B_AB(B) (3)

In some cases, two components can be considered for each monomer within Δ_defE. This happens when the geometry adopted by one (or both) of the monomers in the AB complex is closer to a conformer which is not the lowest energy one. Then Δ_defE defined by (2) can be split into four terms (eqn (4)), two are provided by the relative energies of the conformers in the isolated monomers, usually called conformational terms, Δ_confE. The others are due to the geometry distortions experienced by those conformers in the complex, Δ_gdE.

Δ_defE(AB) = Δ_confE(A) + Δ_gdE(A) + Δ_confE(B) + Δ_gdE(B) (4)

Once the complex AB has been calculated using a finite basis set, the BSSE comes up and, according to CP correction, the complexation and binding energies are overestimated to an extent given by eqn (5). The CP-corrected energy for the AB complex formation between fragment A and B is obtained by (6), and the corresponding CP-corrected binding energy by (7).

Δ_BSSEE(AB) = [E^A_AB(A) − E^AB_AB(A)] + [E^B_AB(B) − E^AB_AB(B)] (5)

Δ^CP_cE(AB) = [E^AB_AB(AB) − E^A_A(A) − E^B_B(B)] + Δ_BSSEE(AB) (6)

Δ^CP_bE(AB) = [E^AB_AB(AB) − E^A_AB(A) − E^B_AB(B)] + Δ_BSSEE(AB) (7)

QTAIM allows a disjoint partitioning of any molecular property, P, of a certain molecular system, S, into their atomic counterparts, P(Ω), (eqn (8)) where Ω denotes an atom in the system.

Thus, molecular energy, E, can be split into atomic ones, E(Ω). In the same vein, atomic contributions to CP-uncorrected or corrected complexation or binding energies, as well as deformation energies can be obtained through eqn (9)–(13). In these equations, either E^A_A(Ω) = 0 or E^B_B(Ω) = 0, depending on which molecule contains the Ω atom. Thus, any of these expressions can be reduced to the difference between the atomic energy within the complex, X, or in the isolated molecule, M, enclosing the atom Ω under study. Moreover, it is possible to define an atomic contribution to the CP correction, Δ_BSSEE(Ω_X) as the difference between those atomic energies obtained for Ω in the isolated molecule with the basis set employed for the monomer, E^M_X(Ω_M), and that used for the whole system, E^M_X(Ω_M), both of them computed in the geometry of the complex (eqn (14)).

Δ_BSSEE(Ω_X) = E^M_X(Ω_M) − E^X_X(Ω_M) (14)

BSSE may affect not only to atomic energies, but to any other atomic property, P, and the corresponding atomic CP correction, Δ_BSSEP(Ω_X), can be obtained through the generalization of eqn (14), yielding eqn (15).

Δ_BSSEP(Ω_X) = P^M_X(Ω_M) − P^X_X(Ω_M) (15)

Hereafter, specific values of QTAIM atomic properties will be denoted as P^β_α(Ω_χ), noticing that the density property P was integrated within the atomic basin Ω in the χ molecular system displaying α geometry when using the β basis set. We will pay particular attention to atomic electron populations, making extensive use of: N^X_X(Ω_X), N^M_M(Ω_M), N^X_X(Ω_M), and N^M_X(Ω_M) values. This allows estimating (to our knowledge for the first time) to which extent the depletion or enlargement experienced by the electron population of atoms involved (or not) in intermolecular bonds is affected by BSSE. Thus, the usual variations of atomic electron population computed for a complexation process, Δ_cN(Ω), can be corrected by estimating the effects of BSSE on it by considering Δ_BSSEN(Ω_X), giving rise to the CP-corrected term, Δ_c^CPN(Ω). Moreover, Δ_c^CPN(Ω) can be subsequently split into deformation, Δ_defN(Ω), and CP-corrected binding terms, as shown in eqn (16). These terms are given by eqn (17)–(19). While the two formers are physically meaningful, the latter only estimates the effect of a computational shortcoming connected to the basis set size.

Δ_defN(Ω) = N^M_X(Ω_X) − N^M_M(Ω_M) (17)

Δ_BSSEN(Ω_X) = N^M_X(Ω_M) − N^X_X(Ω_M) (19)

Moreover, when the geometry of a fragment in a certain complex corresponds to a distortion of a non most stable conformer, the deformation energy involves a conformational term, Δ_confE, given by the relative energy of that conformer.

3. Results and discussion

3.1. On the stability of complexes formed between cocaine and heroin and/or their metabolites

In order to understand complexation energy from its atomic contribution into the QTAIM framework, the conformational analysis was carried out to obtain the most stable complexes between heroin (or morphine) with cocaine (or ecgonine methyl ester) but its complete conformational distribution. Therefore, the conformational analysis was done by a manual procedure. The criterion to find complexes was defined on the basis that intermolecular hydrogen bonds and stacking interactions stabilize these monomers due to the presence of these functional groups: hydroxyl in M and E; ammonium in M, H, C and E; acetyl in H, E and C; aromatic rings in M, H and C; and vinylic hydrogens in M and H (Fig. 1). The former interaction was searched between hydrogen donor atoms, such as O–H, N–H and =C–H, with acceptor atoms, such as =O and the oxygen in –O–H. The latter interaction was searched between aromatic groups and all the previously functional groups described. Thus, we have found eleven different complexes between heroin and morphine with cocaine and ecgonine methyl ester (Fig. 2). These complexes are named making reference to their constituting monomers. Whereas no energy minimum was obtained for the adduct containing heroin and cocaine, two different minima (HE1 and HE2) are found for the complex between heroin and ecgonine methyl ester, morphine and cocaine were found attached in four different arrangements (MC1–MC4), and five complexes are formed by morphine and ecgonine methyl ester (ME1–ME5). These 11 structures display intramolecular (only within cocaine and ecgonine methyl ester monomers) and intermolecular hydrogen bond, IHB, but in one case, MC3, where the monomers display stacking interactions with no IHB.

The diverse structures found for each of the three different pairs of molecules (HE, MC and ME) differ in the atoms involved in the main intermolecular bond, as shown in Table 1.

Table 1 Distances, R (in Å), and corresponding BCP parameters (multiplied by 10³ and in au) electron density, ρ_BCP, its Laplacian, ∇²ρ_BCP, and total energy density, H_BCP for hydrogen bonds and stacking interactions displayed by the complexes found^a,b

Complex	Interaction	R	ρBCP	∇ρBCP	HBCP
a Most significant IHBs are in boldface.b Atomic numbering for IHBs and stacking interactions refer to the first monomer in the acronym before the interaction sign (⋯) and to the second one after it. These values were calculated at B3LYP/6-31++G(d,p) 6d in gas phase.
Intramolecular hydrogen bond
HE1	(O1′⋯H–N8)E	1.791	40.17	113.68	−1.35
HE2	(O1′⋯H–N8)E	1.771	42.04	119.10	−1.47
MC1	(O1′⋯H–N8)C	1.781	41.05	117.09	−1.36
MC2	(O1′⋯H–N8)C	2.058	22.95	67.97	−0.50
MC3	(O1′⋯H–N8)C	1.782	40.91	116.68	−1.35
MC4	(O1′⋯H–N8)C	1.931	29.78	85.64	−0.96
ME1	(O1′⋯H–N8)E	1.778	41.29	117.38	−1.40
ME2	(O1′⋯H–N8)E	1.783	40.86	116.03	−1.37
ME3	(O1′⋯H–N8)E	1.937	29.37	84.17	−0.97
ME4	(O1′⋯H–N8)E	1.893	32.30	91.52	−1.12
ME5	(O1′⋯H–N8)E	1.934	29.61	84.94	−0.99
Intermolecular hydrogen bond (IHB)
HE1	C1′′=O1′′⋯H–O3	1.960	24.34	67.71	−0.55
	C2′′–H⋯O3′	2.487	7.94	30.75	1.09
HE2	C1′=O1′⋯H–O3	1.887	28.29	81.26	−0.18
	C2′–H⋯O3′	2.610	6.33	25.07	1.06
MC1	O6–H⋯O1′′=C1′′	1.883	25.96	85.33	0.87
	O6⋯H6a–C6	2.563	7.42	25.80	0.81
	O6⋯H4e–C4	3.081	2.68	11.00	0.65
MC2	O3⋯H–N8	2.218	14.10	46.23	0.29
	O3–H3⋯O1′=C1′	2.210	13.81	58.85	1.55
	O⋯H8a–C8	2.501	8.15	28.76	0.83
	O⋯H–C5	3.113	2.66	10.54	0.67
	O6⋯H–C5	3.035	2.61	10.66	0.70
MC4	O3–H3⋯O1′=C1′	2.179	13.79	56.85	1.28
	O3⋯H–N8	2.426	8.31	32.56	0.94
ME1	O3⋯H–O3	2.039	19.20	56.12	−0.29
	O6⋯H6a–C6	4.573	0.11	0.41	0.04
ME2	O6⋯H–O3	1.976	22.31	64.51	−0.29
ME3	O3–H3⋯O1′=C1′	2.163	14.43	59.60	1.31
	O3⋯H8–N8	2.353	9.80	37.30	0.87
	O⋯Ha–C8	2.700	5.15	19.86	0.91
ME4	O3–H3⋯O1′=C1′	2.108	15.04	54.97	0.76
	O⋯Hg^−–C2′	3.472	1.16	4.98	0.37
	O6⋯Hg^+–C2′	3.608	0.98	3.84	0.29
ME5	O6–H6⋯O1′=C1′	2.098	16.41	63.68	1.05
	O6⋯H–N8	2.392	8.82	34.30	0.92
Stacking interactions
MC3	O3⋯C2r	3.397	4.59	14.79	0.60
	O3⋯O1′	3.597	2.46	11.20	0.63
	H5⋯H3r	2.810	1.99	6.56	0.44
	O⋯C3r	4.053	1.39	5.25	0.34
	O6⋯H3r–C3r	3.319	1.87	7.04	0.50
	H2⋯O1′=C1′	3.338	1.64	6.67	0.45

---

All the structures obtained for these adducts display significant positive values for complexation energies (even uncorrected ones). They become slightly larger when zero point vibrational corrections, ZPE, are considered (Table 2). The least positive values of complexation energies correspond to the two complexes established between heroin and ecgonine methyl ester, whose Δ_c^CPE values are, at least, 10 kJ mol⁻¹ less positive than any ME complex.

Table 2 Uncorrected, Δ_cE, and CP-corrected, Δ_c^CPE, complexation energies, basis set superposition error, BSSE, and CP-corrected binding energies for the complexes here studied. Geometry distortion, Δ_gdE, and conformational terms, Δ_confE, for their constituting monomers. Contributions to Δ_cE to due zero point energy corrections, ΔZPE, are also shown. All values are in kJ mol^−1a,b

Complex	ΔZPE	ΔcE^a	ΔBSSEE	Δc^CPE	ΔgdE^b		ΔconfE
					H/M	C/E	H/M	C/E
a ΔZPE corrections not included.b Total deformation energy, as defined by eqn (2), is obtained by adding the corresponding ΔconfE and ΔgdE terms. These values were calculated at B3LYP/6-31++G(d,p) 6d in gas phase.
HE1	4.50	74.94	2.33	77.28	3.94	0.87	0.00	1.84
HE2	3.39	81.36	2.42	83.77	12.79	1.00	0.00	1.84
MC1	3.42	80.95	4.30	85.24	0.45	2.32	5.38	0.00
MC2	3.33	88.21	4.57	92.78	2.47	6.67	13.12	0.00
MC3	0.46	103.38	3.30	106.68	0.53	0.60	0.00	0.00
MC4	1.01	104.86	3.15	108.01	1.18	2.48	0.00	0.00
ME1	3.09	91.47	2.67	94.15	2.29	0.62	12.05	0.00
ME2	3.24	93.91	2.60	96.51	1.31	1.38	0.00	1.84
ME3	2.80	93.69	3.63	97.32	3.14	2.60	12.05	0.00
ME4	2.90	98.86	3.01	101.87	1.09	1.13	5.38	0.00
ME5	2.78	102.97	3.09	106.06	0.87	2.34	5.38	0.00

---

3.1.1. The solvent effects on stability of the most stable complex, HE1. In contrast, both Δ_cE and Δ_c^CPE values become negative when the calculations are carried out in aqueous solution modelled with PCM. Thus, Δ_c^CPE is −19.2 kJ mol⁻¹ for HE1. As PCM does not consider IHB with solvent molecules, we thought that it was more reliable to perform PCM calculations for complex formation including explicitly some water molecules. This task was only carried out for the complex formed with HE1 with 3 explicit water molecules, HE1·3W. The water molecules were inserted in the surroundings of the IHB region, two around the donor and one on the acceptor atom (Fig. 3). These positions of the added water molecules to the atoms involved in IHB were chosen to stabilize the IHB through charge transferred from water molecules.¹⁹ The complex was optimized in gas phase and then was subsequently subjected to a PCM single point calculation. The same computational procedure was considered for the complexes formed by one heroin molecule attached to two water molecules, H·2W, and that formed by ecgonine methyl ester with three water molecules, E·3W (Fig. 3). Heroin and ecgonine methyl ester were arranged in the corresponding initial geometry of H·2W and E·3W complexes in the same conformation displayed by those monomers in the HE1·3W optimized structure. Thus, the BSSE-uncorrected energy variation for the process indicated by eqn (20) was obtained.

H·2W + E·3W → HE1·3W + 2W (20)

Fig. 3 B3LYP/6-31++G(d,p) 6d optimized structure for the HE1·3W complex (H monomer on the left). Dotted lines indicate hydrogen bonds. The complex structures with the inclusion of water molecules was calculated at B3LYP/6-31++G(d,p) 6d in aqueous solution through PCM.

BSSE corrections for this energy were obtained considering the three-step process, detailed in Scheme 1. The two fragments for CP corrections can be clearly identified in the products (steps I and II) or in the reagents (step III) for every step.

Scheme 1

Table 3 details Δ^CPE obtained for each step in gas phase and PCM modeled aqueous solution. The global complexation process (19) displays a positive value for Δ_c^CPE in gas phase (9.68 kJ mol⁻¹), but turns into a negative one (−73.30 kJ mol⁻¹) in aqueous solution.

Table 3 CP-Corrected complexation energies, Δ^CPE, for the different steps in Scheme 1 are shown. These energy values were obtained from calculations at B3LYP/6-31++G(d,p) 6d in gas phase and in aqueous solution. In both cases, the inclusion of explicit water molecules at the atoms involved in the intermolecular hydrogen bond were taken into account as is shown in Fig. 3

Step	Δ^CPE
	GAS	PCM
H·2W → H·W + W	17.45	14.10
E·3W → E·2W + W	−16.49	−21.46
H·W + E·2W → HE1·3W	9.68	−73.30

---

The analysis of charge transfer involved in the formation of hydrogen bonds in water clusters have been studied by several research groups.¹⁹ Albrecht et al. properly described the charge transfer during the hydrogen bond formation in water clusters,¹⁹ this is a consequence of interacting the lone pair on the accepting oxygen with an empty hydride antibond of the donating O–H atoms, . Occupancy of the antibonding orbital raises the energy of the O–H species while stabilizing the oxygen. The stabilization increases as the cluster size increases because of hydrogen bond cooperativity. This can be described through changing cluster geometry: decreasing O⋯O distances coupled with increasing O–H covalent bond lengths indicate a weakening of the covalent OH bond and strengthening of the H⋯O hydrogen bond interaction. Furthermore, there is an increasingly linear O–H⋯O angle as the size of the ring increases. We notice that the explicit inclusion of water molecules stabilizes the complex formed between heroin and one of cocaine metabolites, ecgonine methyl ester. This stabilization is well supported by the geometrical variations observed on the bond distances and angles of the atoms involved in the intermolecular hydrogen bond between heroin and ecgonine methyl ester. The obtained bond distances in gas phase are shown and in aqueous phase are in parenthesis. O⋯O is 2.88737 Å (2.85010 Å), O–H 0.97687 Å (0.97837 Å), H⋯O hydrogen bond interaction 1.95990 Å (1.90066 Å), and O–H⋯O angle 157.647° (162.831°). Thus, the calculations here reported indicate that the formation of this kind of complexes in human fluids cannot be discarded. Moreover, it could be explored as a reason behind the stronger addiction observed for combined abuse of heroin and cocaine.

3.1.2. The intramolecular van der Waals interactions on the most stable complex, HE1. The general characteristics of the structures of heroin and ecgonine methyl ester, which are part of the complex HE1, are described by polar and nonpolar groups. In the first case, heroin has two carboxylate groups and an ammonium group, and ecgonine methyl ester has a hydroxyl group, a carboxylate group and an ammonium group. On the other hand, the remaining structures that comprise these molecules are aliphatic, in both systems, and an aromatic, in heroin. These molecules bind by means of polar groups. However, the interactions with non-polar groups are not considered in a conventional optimization. Therefore, the study of interactions van der Waals is relevant.

For this purpose, the results of calculating the interaction energy for HE1 complex will be analyzed in this order: (i) conventional optimization with DFT methods: B3LYP and M06, (ii) with the inclusion of potential Grimme, (iii) with the inclusion of the correction to the basis set superposition error by the counterpoise method and (iv) with the inclusion of Grimme potential with the counterpoise method (Table 4).

Table 4 The interaction energy values obtained in gas phase for HE1 complex were calculated at 6-31++G** 6d with different methods. The energy value for complex with the basis set of complex in geometry of complex, E^X_X(X), monomers with the basis set of the complex in geometry of complex, E^X_X(M), monomers with the basis set of isolated monomer in geometry of isolated monomer, E^M_M(M) as well as the uncorrected, Δ_cE, and CP-corrected, Δ_c^CPE, complexation energies and Δ_BSSEE are shown. The terms E^X_X(X), E^M_M(M) and Δ_cE correspond to eqn (1) and Δ_BSSEE, E^X_X(M) and Δ_c^CPE correspond to eqn (5)–(7), respectively. The energy values of the systems are in au and energy derived values are in kJ mol⁻¹

Method	Complex	Monomers	E^XX(X)	E^XX(M)	E^MM(M)	ΔcE	Δ^CPcE	ΔBSSEE
B3LYP	HE1	H	−1917.48289	—	−1245.34854	72.46	—	—
		E		—	−672.16195
M06	HE1	H	−1916.25656	—	−1244.55290	71.76	—	—
		E		—	−671.73099
B3LYP-D3	HE1	H	−1917.58611	—	−1245.41173	63.45	—	—
		E		—	−672.19855
M06-D3	HE1	H	−1916.27612	—	−1244.56514	67.27	—	—
		E		—	−671.73660
B3LYP + CP	HE1	H	−1917.48289	−1245.34837	−1245.34797	69.01	70.99	1.98
		E		−672.16184	−672.16149
M06 + CP	HE1	H	−1916.25634	−1244.55177	−1244.55135	67.93	69.94	2.01
		E		−671.73121	−671.73086
B3LYP-D3 + CP	HE1	H	−1917.58619	−1245.41079	−1245.41035	59.01	61.24	2.23
		E		−672.19873	−672.19831
M06-D3 + CP	HE1	H	−1916.27587	−1244.56401	−1244.56358	63.48	65.60	2.12
		E		−671.73685	−671.73647

---

The value of the interaction energy obtained with B3LYP is 2.14 kJ mol⁻¹ greater than M06. This is because M06 includes in its definition 27% of potential Hartree–Fock exchange. Instead, B3LYP includes 20% of the difference between the potential of Hartree–Fock exchange regarding the potential of the local density. This difference in the weight of Hartree–Fock exchange potential shows that the interaction, through the polar groups, is stabilized by the electronic exchange.

When the dispersion terms of the nonpolar groups in HE1 complex are taken into account, those stabilize the interaction energy. In the case of the energy of interaction calculated with B3LYP with the inclusion of the potential Grimme, it was obtained that the latter is 10.45 kJ mol⁻¹ more stable than the one obtained with only B3LYP. In the case of M06, a similar trend was observed, which is that the interaction energy calculated with M06 and the van der Waals interactions is 4.48 kJ mol⁻¹ more stable than that obtained only with M06. From these results, the importance of including the terms of dispersion in the calculation of the interaction energy is observed. Although the potential of Grimme not describe the real origin of the dipole–dipole interaction between the molecules forming the complex, but rather its effect on the effective potential of midfield of the respective density functional.

Furthermore, by correcting the interaction energy overlaying the BSSE by the method of counterpoise, it was found that the optimized geometry with B3LYP with the counterpoise method is 2.91 kJ mol⁻¹ more stable than that obtained with only B3LYP. Analogously to the interaction energy values obtained with M06, the difference is 1.82 kJ mol⁻¹.

Finally, at the optimization of the complex geometry with either B3LYP or M06 with the potential of Grimme and the counterpoise method was obtained that the interaction energy obtained with B3LYP is 14.89 kJ mol⁻¹ more stable than that obtained with only B3LYP. Similarly, to M06, the difference is 8.28 kJ mol⁻¹.

3.2. Complexation energy terms in gas phase complexes

Retaking our methodology objective, assessing BSSE effects on electron density rearrangements involved in complex formation, and for the sake of simplicity, we have considered gas phase complexes established without including water molecules.

Complexation energies, not including ΔZPE values, as well as all their components mentioned and CP corrections indicated in Section 2.3 (Table 2) are successfully reproduced (within 2.0 kJ mol⁻¹) with atomic QTAIM energies, E(Ω).

The values shown in Table 2 and formula introduced in Section 2.3 allow obtaining CP-corrected binding, Δ_b^CPE, and deformation energies, Δ_defE, for the eleven complexes found. We highlight that all binding energies are positive, precluding negative Δ_c^CPE values. Clearly, Δ_b^CPE values are positive because the repulsion between the global positive charges exhibited by both monomers exceeds the stabilization due to IHBs formation or that obtained from stacking interactions in MC3. The least unfavourable binding corresponds to HE complexes (68.1 and 70.6 kJ mol⁻¹ for, respectively, HE2 and HE1). In contrast, the most unfavourable binding corresponds to the stacked complex MC3 (105.6 kJ mol⁻¹).

We also notice that the largest Δ_gdE values correspond to heroin and cocaine in HE2 and MC2, respectively. Both molecules display just one low energy conformer, while morphine displays four and the ecgonine methyl ester two.^20,21 Consequently, conformational terms, Δ_confE are important only for some of the complexes containing E or M monomers. In fact, the largest deformation energies, Δ_defE, correspond to MC2 (22.3 kJ mol⁻¹) and ME3 (17.8 kJ mol⁻¹) because of the conformational change experienced by the morphine monomer.

We believe it is worth to include here a short comment on how aqueous solution stabilizes these complexes, yielding negative Δ_c^CPE values in PCM calculations. The global charge of complexes doubles that of monomer and results in a more effective solvent polarisation for the former which compensates the still positive Δ_b^CPE values observed with PCM. Moreover, the inclusion of explicit water molecules gives rise to larger distances among the most positive atoms in the complex. We highlight these are not the nitrogen atoms, as incorrectly indicated by the usual Lewis structures employed to represent them. As shown in previous electron density analysis,^20,21 the positive charge of any of the four monomers is basically shared among the ammonium hydrogen and those hydrogens of the methylammonium group, with a contribution from the bridgehead carbons in H and M and of bicycle hydrogens in C and E.

3.3. Variation of atomic properties involved in complex formation

The method introduced in Section 2.3 allows distinguishing between binding and deformation (even splitting it into conformational and geometry distortion) effects on the electron density reorganization accompanying complexation. Moreover, we can make a CP correction for binding effects on electron density. In this section we report and discuss separately the effects of conformational and geometry distortion, BSSE, and CP corrected binding. We detail the amount of CP-corrections and where they become significant. Once more, for the sake of clarity, we restrict our discussion to gas phase complexes with no water molecules.

3.3.1. Effects due to deformations of monomers. In order to analyze deformations (including conformational reorganizations and geometry distortions) and their effects, we look at variations of bond lengths, R_ij, and some QTAIM atomic properties: electron population, N(Ω), energy, E(Ω), the module of the first moment of the electron density, M(Ω), and atomic volume, v(Ω).

Table 5 summarizes the bonds whose length is distorted by more than 5 × 10⁻³ Å. They are concentrated in those regions where IHB are established. Thus, the heroin monomer is mainly distorted around O6 in HE1, while its largest distortions concentrate around O3 in HE2, as the IHBs formed are, respectively, O6–C1′′–O1′′⋯H3–O3 and O3–C1′–O1′⋯H3–O3. They are reinforced in both cases by an additional (and weak) C–H⋯O3 IHB (Fig. 4). In the same vein, the E monomer is mainly distorted, in HE1 and HE2, around the hydroxyl group (O3). The map of distortions is completed by those of the intramolecular HB in C and E monomers, which is present in all the complexes here obtained. We also notice that the bonds including the furan oxygen of H and M monomers (C4–O and C5–O) display moderate distortions in most of the 11 complexes.

Table 5 Most significant variations of bond lengths experienced by the monomers (indicated by the letter preceding complex acronym) of the complexes here studied, Δ_defR_ij, (in Å multiplied by 10³). Variations whose absolute values are below 5 × 10⁻³ Å and those affecting C–H bonds are not shown. These values were calculated at B3LYP/6-31++G(d,p) 6d in gas phase

Bond in H/M	H, HE1	H, HE2	M, MC1	M, MC2	M, MC3	M, MC4	M, ME1	M, ME2	M, ME3	M, ME4	M, ME5
C3–O3	—	8	—	10	—	12	17	−1	10	8	−1
C3–C4	—	—	—	−5	—	—	—	—	—	—	—
C4–O	5	—	—	7	—	−5	−0	7	—	—	7
C5–O	−6	—	—	7	—	6	5	−7	8	—	−4
C6–O6	11	—	—	5	—	—	—	18	5	—	13
C6–C7	—	—	—	—	—	—	—	—	−5	—	−3
O6–H	—	—	8	—	—	—	—	—	—	—	—
O3–C1′	—	−31	—	—	—	—	—	—	—	—	—
C1′=O1′	—	14	—	—	—	—	—	—	—	—	—
O6–C1′′	−25	−1	—	—	—	—	—	—	—	—	—
C1′′=O1′′	13	0	—	—	—	—	—	—	—	—	—

---

Bond in E/C	E, HE1	E, HE2	C, MC1	C, MC2	C, MC3	C, MC4	E, ME1	E, ME2	E, ME3	E, ME4	E, ME5
C2–C3	—	—	—	—	—	—	—	5	—	—	—
C3–O3	−4	−6	7	—	—	—	—	—	—	—	—
O3–H	10	12	—	—	—	—	5	6	—	—	—
N8–H	—	—	—	−10	—	−8	—	—	−10	−7	−9
O1′⋯H(N8)	27	7	—	278	—	151	—	18	165	102	154
O3–C1′′	—	—	−17	—	—	—	—	—	—	—	—
C1′′=O1′′	—	—	10	—	—	—	—	—	—	—	—

---

Fig. 4 Significant Δ_defN(Ω) (in boldface and in au multiplied by 10³), Δ_defE(Ω) (in italics and expressed in kJ mol⁻¹), Δ_defM(Ω) (in au multiplied by 10³), and Δ_defv(Ω) (in parenthesis and in au) values in HE1. See ESI material† for these values in the remaining complexes studied in this work. Δ_defN(Ω), Δ_defE(Ω), Δ_defM(Ω) and Δ_defv(Ω) are not shown when their absolute values are below 10⁻³ au, 1 kJ mol⁻¹, 10⁻³ au, and 1 au respectively. All the values were calculated at B3LYP/6-31++G(d,p) 6d in gas phase.

The principal variations experienced by the main atomic properties due to the geometry deformation involved in the complexation process, Δ_defN(Ω), Δ_defE(Ω), Δ_defM(Ω), and Δ_defv(Ω), are shown for HE1, as an example, in Fig. 4, while the corresponding values for all the systems studied herein are reported as ESI material.† All these conformational variations are analyzed below for each monomer.

Heroin (only present in two complexes) distorts to interact through its acetyl groups (at position 6 in HE1 and the one at position 3 in HE2) as H-acceptor. As indicated above, the largest Δ_defR_ij values are shown by two bonds in the acetyl group: C=O (lengthens) and C–O (shrinks). In accordance, Δ_defN(Ω) values for HE1 (Fig. 4) and HE2 (ESI material†) indicate that electron density is transferred from the C=O group to reinforce the CO bond at position 3 or 6 of the heroin backbone. In fact, the methyl group in the acetyl unit only increases scarcely its electron population after this deformation (4 × 10⁻³ au in HE1 and 7 × 10⁻³ au in HE2). This distortion of the geometry also reduces the polarisation of both atoms of the carbonyl group and increases that of the attached sp³ oxygen. Overall, the C–O–C=O unit involved in IHB destabilizes more than that the whole H monomer (28 vs. 3.6 kJ mol⁻¹ in HE1, and 41 vs. 13.7 kJ mol⁻¹ in HE2). The difference between both quantities is due to stabilizations spread among the remaining atoms of the monomer.

Morphine is included by nine of the eleven complexes herein studied. When this molecule is isolated, it shows six stable conformers.²⁰ They arise from the two possible orientations of the 3-OH group with regard to the attached aromatic ring (synperiplanar, s, or antiperiplanar, a) and from the three different dispositions of the 6-OH group with regard to its attached aliphatic ring (named g, t, and g′ for approximate values of 60°, 180° or −60° in the C5–C6–O–H dihedral angle). Four of them are observed within the nine complexes. The most stable (ag′) corresponds basically to the M monomer in MC3, MC4, and ME2 complexes, while sg′ closely reminds the geometry displayed in MC1, ME4 and ME5, and those of higher energy are included in ME1 and ME3 (sg) and MC2 (st). In all the complexes, the M monomer displays small distortions from the corresponding conformer, as can be anticipated from Δ_gdE values (Table 2). MC1 and MC3 display the least deformed M monomers (Table 5). In the rest of the complexes, the largest |Δ_defR_ij| values correspond to the C–O bond with the hydroxyl group at positions 3 or 6, depending on which is the interaction site for M in the complex. The electron density of the corresponding oxygen (O3 or O6) depletes, and it is mostly transferred to its attached carbon, whose Δ_defN(Ω) values span from 5 × 10⁻³ to 31 × 10⁻³, and rarely to the hydroxyl hydrogen, whose Δ_defN(Ω) values are below 6 × 10⁻³ au. The most deformed M monomers are those bonded through their both hydroxyl groups (MC2, MC3, and ME1). In these, cases Δ_defN(C3) and Δ_defN(C6) are both positive.

Cocaine is in four complexes. Its distortions are oriented toward forming two types of intermolecular interactions: stacking or IHB. Moreover, C can work as H donor or H acceptor or both things depending on the complex considered. The stacked complex MC3 displays the lowest Δ_defE(C) value, 0.60 kJ mol⁻¹, with no significant bond distortion (all are smaller than 3 × 10⁻³ Å), either in C or in M (Table 5). The largest |Δ_defR_ij| values are placed within the benzoyloxy group attached to C3 in MC1, where this monomer is H acceptor exclusively. When cocaine also acts as H donor (MC2 and MC4) N8–H8 shrinks significantly and the N8–H8⋯O1′=C1′ intramolecular hydrogen bond lengthens. Δ_defN(Ω) and Δ_defE(Ω) values within the benzoyloxy group follow the same trends detailed above for the acetyl groups of heroin. The electron density lost by the C1′′=O1′′ unit is mainly transferred to O3, and not to the phenyl group. In contrast, the electron density lost by N8 (and attached methyl) in MC2 and MC4 is mainly transferred to the ammonium hydrogen atom.

Ecgonine methyl ester is in seven complexes. It distorts from two of its most populated conformers.²¹ One observed in ME1, ME3, ME4, and ME5, while the other is in HE1, HE2, and ME2. Its complexes can be also classified considering the role E plays in IHBs: H donor in HE1, HE2, ME1 and ME2; H acceptor in ME4; and both roles in ME3 and ME5. When E acts exclusively as H-donor its interaction site is the 3-OH group, but when it also acts as H-acceptor the H-donor site is the N8–H8 bond. ME1 shows the lowest deformation energy for E (Table 2) with no bond distorted by more than 5 × 10⁻³ Å (Table 5). Slightly larger Δ_defE(E) values are displayed in the other complexes where the E monomer is exclusively H-donor, and the most distorted bonds are concentrated around the 3-OH group, with the intramolecular hydrogen bond scarcely affected. In contrast, when E is H-acceptor, the largest distortions are experienced by the N8–H bond and the N8–H8⋯O1′=C1′ intramolecular hydrogen bond lengthens by more than 0.1 Å. Δ_defN(N8) is positive in the seven E complexes, but there is a quantitative difference between the cases where the E interaction site is 3-OH (1 × 10⁻³ au < Δ_defN(N8) < 3 × 10⁻³ au) and those where N8–H is involved in IHBs (Δ_defN(N8) > 8 × 10⁻³ au).

When the distortions of the four monomers are compared we notice two general trends: (i) heroin acetyl groups produce larger geometry distortion of the aliphatic part of their common backbone than the hydroxyl groups of morphine, as shown by larger Δ_defE values; (ii) cocaine and ecgonine methyl ester monomers that interact affecting their intramolecular hydrogen bond display larger deformation energies than those that interact through 3-OH hydroxyl.

We highlight that, in all cases, if we exclude the stacked complex MC3, distortions are clearly localized in specific areas of each monomer, which coincide with the interaction sites for complex formation. Thus, in morphine, the hydroxyl group and the attached C atom at positions 3 or 6 amount for 35 ± 12% of Σ|Δ_defN(Ω)|, with almost all monomers within the 1σ interval and only those of MC2 (17%), ME2 (52%) and ME3 (21%) out of it. The corresponding contribution for Δ_defE(Ω) is of 44 ± 10%, and only MC2, ME2 and ME3 are more than 1σ from the average. When the same study is carried out along H monomers, just four atoms, the acetyl group at C6 represent 50% of Σ|Δ_defN(Ω)|, 66% of Σ|Δ_defE(Ω)|, 57% of Σ|Δ_defM(Ω)|, and 41% of Σ|Δ_defv(Ω)| in HE1. In the same vein, the corresponding four atoms of the acetyl group at C3 represent 41% of Σ|Δ_defN(Ω)|, 60% of Σ|Δ_defE(Ω)|, 38% of Σ|Δ_defM(Ω)| and 25% of Σ|Δ_defv(Ω)| in HE2. The smaller localization of the distortions observed for H monomer in HE2 can be explained considering that the interaction site is attached to an aromatic ring. In cocaine, we distinguish three cases: (i) in MC1 the interaction site is part of the benzoyloxy group and C3, O3, C1′′, and O1′′ give rise to 50% and 63% of Σ|Δ_defN(Ω)| and Σ|Δ_defE(Ω)|, respectively; (ii) in MC2 and MC4, cocaine mainly interacts through N8 and H8 and their respective contributions are 29%, 22% and 24%, 35%; and (iii) the atoms of the two C–OH bonds of MC3 display the largest Δ_defN(Ω) and Δ_defE(Ω) values. In fact, they account for almost one half (51%) of Σ|Δ_defN(Ω)| and Σ|Δ_defE(Ω)|. Nevertheless, this percentage diminish up to 32% in Δ_defM(Ω) and even to 12% in Δ_defv(Ω), indicating the less localized character of stacking interactions.²² Two different behaviours have been distinguished above for ecgonine methyl ester: interaction as O–H hydrogen donor (HE1, HE2, ME1, and ME2), where C3 and O3 atoms display the largest Δ_defP(Ω) values and simultaneous interaction as H-acceptor and N–H hydrogen donor (ME3, ME4, and ME5), where the most significant Δ_defP(Ω) values correspond to N8, H8, O1′, and C1′. Once again the participation of |Δ_defP(Ω)| for these atoms accounts for significant percentages of the corresponding Σ|Δ_defP(Ω)| quantity.

3.3.2. BSSE effect on atomic properties. As previously detailed, it can be expected that BSSE affects to any atomic property. In this section, we analyze, making use of Δ_BSSEP(Ω) quantities defined in eqn (15), the effect of BSSE on Δ_cN(Ω), Δ_cE(Ω), Δ_cM(Ω) and Δ_cv(Ω) for the complexes here studied. Overall, we notice that the Δ_BSSEP(Ω) in all these complexes were always small. Their highest absolute values for each monomer are confined between 0.5 × 10⁻³ and 2.2 × 10⁻³ au for Δ_BSSEN(Ω), 1.8 and 14.4 kJ mol⁻¹ for Δ_BSSEE(Ω), 1.0 × 10⁻³ and 13.5 × 10⁻³ au for Δ_BSSEM(Ω), and 0.1 to 0.4 au for Δ_BSSEv(Ω). Thus, even for the most affected properties in relative values, Δ_cE(Ω) and Δ_cM(Ω), the influence of BSSE partitioned into its atomic contribution for complex formation by IHB does not seem to be of high significance.

Even though, as depicted by Fig. 5 for HE1, there are scarce atoms displaying differences between Δ_cP(Ω) and Δ_c^CPP(Ω) (or between Δ_bP(Ω) and Δ_b^CPP(Ω)) within the significant figures. As could be expected, part of the largest |Δ_BSSEP(Ω)| values correspond to atoms placed in the region where the number of basis functions centred close enough enlarges in the complex. That is, atoms placed where the two monomers approach each other, no matter if these atoms are involved or not in intermolecular bonds. Thus, we notice that the most positive Δ_BSSEN(Ω) values in HE1 correspond to H(O3) in E (the H involved in IHB) and to H(C2_r) in H, which is not involved in any intermolecular bond, but placed in the contact region between the monomers of this complex. Nevertheless, other Δ_BSSEN(Ω) values do not have a so simple explanation. Thus, we notice that Δ_BSSEN(C2) in H is negative in spite it is close to the contact region between monomers and does not participate in any intermolecular bond. It is the unique atom bonded to H(C2_r), but we think that Δ_BSSEN(Ω) cannot be considered the result of electron density displacement through chemical bonds, as other N(Ω) variations accompanying protonation,²⁴ or other chemical processes like complexation itself. In fact, BSSE effects are consequences of the spatial reinforcement or depletion of the number of basis functions whose numerical values are high enough.

Fig. 5 Significant Δ_BSSEN(Ω) (in boldface and in au multiplied by 10³), Δ_BSSEE(Ω) (in italics and expressed in kJ mol⁻¹) Δ_BSSEM(Ω) (in au multiplied by 10³), and Δ_BSSEv(Ω) (in parenthesis and in au) values in HE1. Reddish colour atoms indicate atoms where Δ_c^CPE(Ω) > Δ_cE(Ω), while greenish ones are the atoms extra-stabilized according to CP-uncorrected atomic energies, that is: Δ_c^CPE(Ω) < Δ_cE(Ω). See ESI material† for these values in the remaining complexes studied in this work. Δ_BSSEN(Ω), Δ_BSSEE(Ω), Δ_BSSEM(Ω) and Δ_BSSEv(Ω) are not shown when their absolute values are below 5 × 10⁻⁴ au, 1 kJ mol⁻¹, 5 × 10⁻⁴ au, and 0.1 au, respectively. All the values were calculated at B3LYP/6-31++G(d,p) 6d in gas phase.

3.3.3. CP corrected atomic properties related to intermolecular interactions. CP corrections to BSSE for atomic properties detailed above were employed, as indicated in eqn (18) to obtain CP-corrected atomic binding energies, Δ_b^CPE(Ω), and the variations due to binding for atomic electron population, Δ_b^CPN(Ω), atomic module of the first moment of the electron density, Δ_b^CPM(Ω), and atomic volume, Δ_b^CPv(Ω). The results obtained are shown in Fig. 6 for HE1, as an example, and presented in ESI material† for the remaining complexes. As commented above, adding geometry distortion and binding contributions provides complexation ones.

Fig. 6 Most significant Δ_b^CPN(Ω) (in boldface and in au and multiplied by 10³), Δ_b^CPE(Ω) values (in italics and expressed in kJ mol⁻¹), Δ_b^CPM(Ω) (in au multiplied by 10³), and Δ_b^CPv(Ω) (in parenthesis and in au) in HE1. See ESI material† for corresponding values for the remaining complexes studied here in. Summation of these atomic variations (excluding Δ_b^CPM(Ω)) are shown for whole monomers (M), atoms and groups whose variation is not detailed (M′), and hydrogens not detailed or included in the methyl group attached to C3′′ (H(M′)). The arrow indicates the total electron density transference between both monomers (in au multiplied by 10³). All the values were calculated at B3LYP/6-31++G(d,p) 6d in gas phase.

Although the figure details the variations of atomic properties experienced by those atoms involved in intermolecular interactions, we notice that, in contrast with variations due to geometry distortion, which are localized, binding introduces significant modifications in the atomic properties of the whole complex with regard to those in isolated monomers. This is shown in Fig. 6 by comparing the summations of Δ_b^CPN(Ω) and Δ_b^CPE(Ω) for each monomer and the corresponding sum extended only to those atoms not involved in intermolecular bonds. As previously shown for much more simple complexes,²⁵ the electron density transferred between monomers in HE1 (0.028 au) (or in any other complex here studied) is exceeded by the loss of electron density shown by the atoms not involved in IHBs in the electron density donor (heroin, in this case) (0.037 au). We highlight the relevant role played by the hydrogens of this part of the molecule, which all together lose 0.040 au. Looking at the acceptor monomer, ecgonine methyl ester in this case, we observe that most of the electron density transferred results in increasing the electron density of that part of the molecule which is again not involved in IHBs (0.020 au). Moreover, the hydrogens of this part of E are the atoms that, in total, increase their electron density the most (0.016 au in HE1). That is, the hydrogens on the periphery of the monomers play a very important role for understanding electron density reorganization in complex formation. This role reminds that of sink or source of electron density, shown initially by Stutchbury and Cooper three decades ago,²⁶ for the hydrogens of methyl groups in protonations and deprotonations, and confirmed by our group for the protonation of a wide series of compounds.^24,27

Among the complexes here obtained, ME2 is the simplest case for analyzing how electron density reorganizes because of intermolecular binding. In fact, this is the only complex where we only find one intermolecular bond path (the IHB connecting O6 in M to H3 in E) in QTAIM analysis (Table 1). For this analysis we consider the following four disjoint fragments in ME2: (i) the atomic basins of the H-acceptor C–O–H unit of morphine (C6, O6 and H6) in the IHB; (ii) the remaining atoms of morphine; (iii) the atomic basins of the H-donor unit H–O3–C3 of ecgonine methyl ester; and (iv) the remaining atomic basins of ecgonine methyl ester. They are denoted as fragments I, II, III, and IV in what follows. For each of these fragments we have performed atomic integrations of the proper density function for the zero-flux atomic basins of these fragments, defined from electron densities computed for the complex and the monomer with the geometry and basis functions of the complex (respectively [ρ^X_X(r)]_X and [ρ^X_X(r)]_M), as indicated by eqn (18). The results indicate that binding is the responsible of transferring 0.023 au from M to E. The electron density of both fragments involved in the IHB (I and III) is reinforced in, 0.050 and 0.009 au, respectively. The electron density of fragment IV is also reinforced by binding (0.013 au), while the atoms not involved in IHB in M monomer (fragment II) are the ones losing electron density (0.072 au). Once again, hydrogens play a leading role in the electron density reorganization as they provide most of the electron density (0.066 au) lost by fragment II and receive almost one half of that received by fragment IV (0.012 au). These results also lead us to conclude that representing the IHB formation as a certain kind of nucleophilic attack from the H-acceptor to the H of the donor unit is unreliable.

It is also worth to discuss the electron density reorganization accompanying MC3 formation, as it is the unique complex formed by stacking interactions. In this case, the electron density flows from C to M. The transference is certainly small (0.004 au), as previously found in stacking complexes,^23,28 but exceeds that obtained for other IHB complexes in this study. In contrast with IHB complexes, both molecular fragments defined by collecting the atoms not involved in stacking bond paths, display negative values ΣΔ_b^CPN(Ω) values (−0.028 au in C and −0.027 au in M). Again, both fragments formed by the atoms connected by bond paths reinforce significantly their electron population (0.022 au in C and 0.031 au in M). In this case, the set of hydrogens not included in intermolecular bond paths plays the role of electron density donors (0.021 au in C and 0.035 au in M). The electron density taken from these hydrogens reinforces that of the atoms involved in stacking interactions and also that of heavy atoms in the rest of each monomer.

As shown in the three cases analyzed in detail, the electron density can be transferred from the H/M monomer to the E/C one (HE1, HE2, MC2, ME1, and ME2), or it can follow the opposite way (MC1, MC3, MC4, ME4, and ME5). Moreover, this electron density transference can be negligible (ME3).

4. Conclusions

B3LYP optimizations carried out for complexes formed between heroin or morphine with cocaine or ecgonine methyl ester lead to obtain eleven different minima described in this work. While all the structures obtained display positive values for complexation energy in gas phase, these values become negative in aqueous solution simulated with PCM. Explicit inclusion of water molecules stabilizes the complexes more. Thus, the calculations here reported indicate that the formation of this kind of complexes in human fluids cannot be discarded, and could be explored as a reason behind the stronger addiction observed for combined abuse of heroin and cocaine.

Moreover, the calculated interaction energy for the most stable complex, HE1, with M06/6-31++G** 6d is more stable than B3LYP/6-31++G** 6d due to the definition of the former density functional, which gives more weight to the Hartree–Fock exchange potential. Furthermore, the modeling of the dipole–dipole interactions between the molecules of the complex allows a better description of the interaction, which is reflected in a more stable interaction energy of the complex. On the other hand, when the interaction energy is calculated with the inclusion of the dispersion terms with the BSSE correction by counterpoise method is observed that the difference between the latter value and that obtained only with the density functional is equal to the sum of its parts, this is to say that the difference of the interaction energy value obtained with the Grimme potential from the obtained only with the respective density functional plus the difference of the interaction energy value obtained with the counterpoise method from the obtained only with the respective density functional. Whose errors relative to the values obtained with B3LYP and M06 are of 7.22% and 0.36%, respectively.

The reorganization of electron density in the formation of these eleven complexes in gas phase was subjected to a QTAIM analysis distinguishing between deformation and binding terms for atomic properties. While deformation terms are computed for structures obtained with the same basis set, binding terms are affected by basis set superposition error. In this work, we have extended the counterpoise method to obtain CP corrections for binding atomic properties, as the difference between a certain atomic property computed with the electron density of the optimized structure of the complex and with that computed for isolated monomer with the geometry and basis set of the complex. These corrections were found to be very small in the eleven complexes here reported.

The values of the four atomic properties here considered (electron population, total energy, module of the first moment of the electron density and volume) indicate that: (i) geometry distortions affect basically the atoms involved in intermolecular bondpaths and their closest neighbours, allowing the definition of interactions sites for each complex; and (ii) in contrast, binding modifies the atomic properties of the whole molecule in a significant extent. Thus adding the variation of atomic properties due to CP corrected binding out of the interaction sites of one monomer usually exceeds the corresponding value for the interaction sites.

CP corrected binding variations of atomic electron population for complexes established by intermolecular hydrogen bonds, indicate that, overall, those atoms (especially hydrogens) placed outside of the interaction sites of each monomer provide (or receive) more electron density than that transferred from one monomer to another. Thus the electron density reorganization accompanying the formation of intermolecular hydrogen bonds can be described by distinguishing the monomers as electron density donor and acceptor. The non-interacting part of the donor provides electron density to such an extent that allows reinforcing the electron density of its interaction sites and transferring electron density to the acceptor. The acceptor distributes the electron density gained reinforcing that of the interaction sites, but also increasing significantly the atomic populations of several atoms in the rest of the molecule. Therefore, representing the IHB formation as a certain kind of nucleophilic attack from the H-acceptor to the H of the donor unit is unreliable.

Moreover, in the unique stacking complex found for these compounds, the electron density of atoms involved in intermolecular bond paths and heavy atoms in the rest of the monomers is reinforced by depleting that of hydrogens.

Acknowledgements

We thank Centro de Supercomputación de Galicia, CESGA, for free access to its computational facilities.

Notes and references

1. (a) T. D. Dillehay, J. Rossen, D. Ugent, A. Karathanasis, V. Vásquez and P. J. Netherly, Antiquity, 2010, 84, 939; (b) T. D. Dillehay, J. Rossen and P. J. Netherly, Am. Sci., 1997, 85, 46; (c) M. A. Rivera, A. C. Aufderheide, L. W. Cartmell, C. M. Torres and O. Langsjoen, J. Psychoact. Drugs, 2005, 37, 455; (d) M. E. Penny, A. Zavaleta, M. Lemay, M. R. Liria, M. L. Huaylinas, M. Alminger, J. McChesney, F. Alcaraz and M. B. Reddy, Food Nutr. Bull., 2009, 30, 205.
2. I. A. Ramoutsaki, H. Askitopoulou and E. Konsolaki, Int. Congr. Ser., 2002, 1242, 43.
3. (a) J. W. Cornish and C. P. O'Brien, Annu. Rev. Public Health, 1996, 17, 259; (b) N. W. Withers, L. Pulvirenti, G. F. Koob and J. C. Gillin, J. Clin. Psychopharmacol., 1995, 15, 63; (c) G. F. Koob and M. Le Moal, Science, 1997, 278, 52; (d) D. E. Selley, C. C. Cao, T. Sexton, J. A. Schwegel, T. J. Martin and S. R. Childers, Biochem. Pharmacol., 2001, 62, 447; (e) C. E. Inturrisi, M. Schultz, S. Shin, J. G. Umans, L. Angel and E. J. Simon, Life Sci., 1983, 33, 773; (f) J. M. Andersen, A. Ripel, F. Boix, P. T. Normann and J. Mørland, J. Pharmacol. Exp. Ther., 2009, 331, 153.
4. G. F. Koob and M. Le Moal, Neuropsychopharmacology, 2001, 24, 97.
5. (a) R. D. Blakely, L. J. De Felice and H. C. Hartzell, J. Exp. Biol., 1994, 196, 263; (b) I. Sora, F. S. Hall, A. M. Andrews, M. Itokawa, X. F. Li, H. B. Wei, C. Wichems, K. P. Lesch, D. L. Murphy and G. R. Uhl, Proc. Natl. Acad. Sci. U. S. A., 2001, 98, 5300; (c) H. E. Melikian, R. D. Blakely, J. K. McDonald, H. Gu, G. Rudnick and K. R. Moore, J. Biol. Chem., 1994, 269, 12290; (d) F. S. Hall, X. F. Li, J. R. Thompson, I. Sora, D. L. Murphy, K. P. Lesch, M. Caron and G. R. Uhl, Neuroscience, 2009, 162, 870; (e) I. Sora, C. Wichems, N. Takahashi, X. F. Li, Z. Zeng, R. Revay, K. P. Lesch, D. L. Murphy and G. R. Uhl, Proc. Natl. Acad. Sci. U. S. A., 1998, 95, 7699; (f) M. C. Ritz, R. J. Lamb, S. R. Goldberg and M. J. Kuhar, Science, 1987, 237, 1219; (g) B. J. Hoffman, E. Mezey and M. J. Brownstein, Science, 1991, 254, 579; (h) J. E. Kilty, D. Lorang and S. G. Amara, Science, 1991, 254, 578; (i) S. Shimada, S. Kitayama, C. L. Lin, A. Patel, E. Nanthakumar, P. Gregor, M. Kuhar and G. Uhl, Science, 1991, 254, 576.
6. (a) J. Camí and M. Farré, N. Engl. J. Med., 2003, 349, 975; (b) F. H. Gawin, Science, 1991, 251, 1580.
7. S. W. Johnson and R. A. North, J. Neurosci., 1992, 12, 483.
8. (a) G. E. Torres, R. R. Gainetdinov and M. G. Caron, Nat. Rev. Neurosci., 2003, 4, 13; (b) S. S. Negus, M. B. Gatch and N. K. Mello, J. Pharmacol. Exp. Ther., 1998, 285, 1123; (c) F. Leri, J. Bruneau and J. Stewart, Addiction, 2003, 98, 7.
9. (a) G. N. Stowe, L. F. Vendruscolo, S. Edwards, J. E. Schlosburg, K. K. Misra, G. Schulteis, A. V. Mayorov, J. S. Zakhari, G. F. Koob and K. D. Janda, J. Med. Chem., 2011, 54, 5195; (b) S. Wee, M. J. Hicks, B. P. De, J. B. Rosenberg, A. Y. Moreno, S. M. Kaminsky, K. D. Janda, R. G. Crystal and G. F. Koob, Neuropsychopharmacology, 2012, 37, 1083.
10. (a) T. Cunha-Oliveira, A. C. Rego, J. Garrido, F. Borges, T. Macedo and C. R. Oliveira, Toxicology, 2010, 276, 11; (b) J. M. P. J. Garrido, M. P. M. Marques, A. M. S. Silva, T. R. A. Macedo, A. M. Oliveira-Brett and F. Borges, Anal. Bioanal. Chem., 2007, 388, 1799.
11. (a) H. B. Jansen and P. Ross, Chem. Phys. Lett., 1969, 3, 140; (b) S. F. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553.
12. (a) R. F. W. Bader, Monatsh. Chem., 2005, 136, 819; (b) R. F. W. Bader, Chem. Rev., 1991, 91, 893; (c) R. F. W. Bader, Atoms in Molecules: A Quantum Theory, Oxford University Press, New York, 1990; (d) P. L. A. Popelier, Atoms in Molecules: An Introduction, Prentice Hall, 1999.
13. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09, Revision A.02, Gaussian, Inc., Wallingford CT, 2009.
14. (a) C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785; (b) B. Miehlich, A. Savin, H. Stoll and H. Preuss, Chem. Phys. Lett., 1989, 157, 200.
15. (a) A. D. Becke, Phys. Rev. A, 1988, 38, 3098; (b) A. D. Becke, Chem. Phys., 1993, 98, 5648.
16. (a) B. Mennuccia and J. Tomasi, J. Chem. Phys., 1997, 106, 5151; (b) M. T. Cances, B. Mennucci and J. Tomasi, J. Chem. Phys., 1997, 107, 3032; (c) M. Cossi, V. Barone, B. Mennucci and J. Tomasi, Chem. Phys. Lett., 1998, 286, 253.
17. S. Grimme, J. Antony, S. Ehrlich and S. Krieg, J. Chem. Phys., 2010, 132, 154104.
18. T. A. Keith, AIMAll, Version 13.02.26, TK Gristmill Software, Overland Park, KS, USA, 2012, http://www.aim.tkgristmill.com.
19. (a) L. Albrecht and R. J. Boyd, J. Phys. Chem. A, 2012, 116, 3946; (b) S. S. Xantheas, Chem. Phys., 2000, 258, 225; (c) R. A. Klein, in NIC Symp, NIC Series 32, ed. G. Munster, D. Wolf and D. Kremer, J. v. N. Inst. Comp., Julich, 2006, p. 65; (d) B. F. King and F. J. Weinhold, Chem. Phys., 1995, 103, 333; (e) F. J. Weinhold, J. Mol. Struct.: THEOCHEM, 1997, 181, 398.
20. D. A. Rincon, M. N. D. S. Cordeiro and R. A. Mosquera, Int. J. Quantum Chem., 2010, 110, 2472.
21. D. A. Rincon, M. N. D. S. Cordeiro, R. A. Mosquera and F. Borges, Chem. Phys. Lett., 2009, 467, 249.
22. D. A. Rincon, M. N. D. S. Cordeiro and R. A. Mosquera, J. Phys. Chem. A, 2009, 113, 13937.
23. R. A. Mosquera, M. J. González Moa, L. Estévez, M. Mandado and A. M. Graña, An Electron Density-Based Approach to the Origin of Stacking Interactions, in Quantum Biochemistry, ed. C. F. Matta, Wiley VCH, Weinheim, 2010, vol. 1, pp. 365–387.
24. (a) N. Otero, M. J. González Moz, M. Mandado and R. A. Mosquera, Chem. Phys. Lett., 2006, 428, 249; (b) D. Ferro-Costas and R. A. Mosquera, J. Phys. Chem. A, 2013, 117, 257.
25. A. Vila and R. A. Mosquera, Chem. Phys., 2003, 291, 73.
26. N. C. J. Stutchbury and D. L. Cooper, J. Chem. Phys., 1983, 79, 4967.
27. (a) M. J. González Moa and R. A. Mosquera, J. Phys. Chem. A, 2005, 109, 3682; (b) J. L. López, A. M. Graña and R. A. Mosquera, J. Phys. Chem. A, 2009, 113, 2652; (c) A. Vila and R. A. Mosquera, Tetrahedron, 2001, 57, 9415.
28. (a) M. J. González Moa, M. Mandado and R. A. Mosquera, J. Phys. Chem. A, 1998, 2007, 111; (b) L. Estévez, N. Otero and R. A. Mosquera, J. Phys. Chem. A, 2009, 113, 11051; (c) L. Estévez, M. Sánchez-Lozano and R. A. Mosquera, RSC Adv., 2014, 4, 25018.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6ra22736h

---