Computational Study on the Conformational Preferences of Neutral, Protonated and Deprotonated Glycine Dimers

Abstract

A conformational analysis has been carried out for monoprotonated, unprotonated and deprotonated glycine dimers in the gas phase and an aqueous solution. MP2/6-311++(d,p), B3LYP/6-311++(d,p) and M06/6-311++(d,p) optimizations were performed for more than 200 initial conformations comprising nonionic (COOH–CH₂–NH₂) (N) and zwitterionic (COO⁻–CH₂–NH₃⁺) (Z) structures for neutral monomers. All the methods indicate that Z monomers are preferred over N ones for the neutral and deprotonated dimers in aqueous solutions, whereas the reverse trend is observed in the gas phase (including also protonated dimers). NC and ZC structures coexist in aqueous solutions for the protonated glycine dimer. The preferred geometries are significantly different depending on the media and total dimer charge. Moreover, several minima display close energies in each series (media and total dimer charge). New conformers, not previously reported, are found to be significantly populated in those conformational mixtures. Dimers containing Z monomers are associated with larger absolute solvation energies and are more prone than N-containing ones to experience protonation and deprotonation in the gas phase, whereas the reverse trend is observed in the aqueous solution. The Quantum Theory of Atoms in Molecules (QTAIM) analysis reveals that uncharged dimers display trifling electron density transfer between monomers, whereas it is significant in anionic and cationic dimers.

1. Introduction

Dimers of amino acids and other clusters of biomolecules have received large attention as they are simple models for studying intermolecular hydrogen bonds (IHB), which play a major role for secondary and ternary protein structures [1,2,3], molecular recognition [4], and other important biochemical processes [5]. Moreover, the ability of diverse dimers of amino acids to act as molecular linkers prompted technical interest on their molecular structure and how pH alters it [6].

Glycine, the smallest amino acid and a model for more complex compounds, has received considerable attention in structural chemistry. Studies on glycine cover a wide range of properties and features, interconversion between zwitterions and non-ionic forms being one of them. Thus, zwitterions have not been found in isolated glycine monomers or dimers, whereas they were detected in glycine–water gas phase complexes. In fact, Ramaekers has reported the existence of zwitterions in n:1 gas phase water–glycine complexes, according to IR spectroscopy [7]. Moreover, the stabilization of glycine zwitterions in gas phase has been analysed in several studies, either microsolvanting it with water molecules [8] or adding metal ions [9]. This is not a straightforward problem because of the large number of possible structures for glycine–water complexes. Thus, a wide variety of structures was created by Hong [10] for a (1:1) glycine–water complex. Complexity increases substantially when five water molecules are included to induce the presence of zwitterions [11].

Could zwitterions be present in glycine homodimers? This question was considered in previous theoretical and/or experimental papers dealing with neutral and protonated glycine dimers in a gas phase. In contrast, we have not found any theoretical or experimental structural studies of anions.

A first detailed study on the molecular structure of the protonated cationic glycine homodimer, Gly₂H⁺, was carried out by Price et al. [12]. Using blackbody infrared radiative dissociation (BIRD), they found that its most stable structure contains a cation bonded to a non-zwitterionic monomer. According to their study, Gly₂H⁺ displays N–H···N IHB. In contrast, Raspopov and McMahon, combining high-pressure mass spectrometry and MPW1PW91/6-31G(d) calculations, reported that the most stable Gly₂H⁺ structure displays N–H···O IHB [13]. The energy of this structure is below the one reported by Price et al. Later, Wu and McMahon compared the calculated infrared spectra for the three most stable geometries of Gly₂H⁺ with the experimental one. They concluded that the dominant species is formed by a neutral non zwitterionic structure bonded to a cation bonded through a H₂N⁺–H···O IHB [14]. This structure was also proposed as the most stable one by Atkins et al., re-examining Gly₂H⁺ conformational preferences in the gas phase [15], and by Armentrout et al., using MP2(full)/6-311+G(2d,2p)//B3LYP/6-311+G(d,p) calculations [16]. Nevertheless, the latter authors also report a new low-energy geometry-displaying structure in N–H···N IHB.

The stability of the neutral glycine dimer, Gly₂, seems to be clear in the gas phase [17,18]. Nevertheless, Huang et al. concluded that monomers are preferred over dimers in supersaturated water solutions [19], in agreement with Hamad et al. in an ab initio molecular dynamics study [20]. Other experimental and theoretical studies [15,19] led to researchers conclude that Gly₂ dimers are stabilized by Na⁺ cations, which is feasible in biological systems. Moreover, several computational conformational analyses for Gly₂ were carried out at different theoretical levels. Most of them only considered the gas phase with non-zwitterionic structures [17,18], while the most recent ones also consider zwitterionic forms and aqueous solutions [20,21,22].

The present study expands previous ones, considering neutral, protonated, and deprotonated dimers, built with non-ionic or zwitterionic monomers. Calculations are carried out both for the gas phase and the aqueous solution. In all cases, a wide set of suitable geometries, with diverse kinds of IHBs, has been included in our conformational sampling. In fact, involving a large number of initial structures is crucial when dealing with very smooth molecular energy surfaces, as those of glycine dimers.

Describing the conformational preferences of these systems is a primary goal of this paper, but we also aim to analyse the effects provided by the non-ionic or zwitterionic structure of monomers on several properties, like glycine dimer acid/base features, the electron density transference between monomers, or the strength of the intermolecular hydrogen bonds. To these ends we used the quantum theory of atoms in molecules (QTAIM) [23].

2. Materials and Methods

In order to model different pH conditions, we have studied monoprotonated, unprotonated and deprotonated glycine dimers. They were built by joining different monomers: cation (COOH–CH₂–${\mathrm{NH}}_3^{+}$), C, nonionic neutral (COOH–CH₂–NH₂), N, zwitterion (COO⁻–CH₂–${\mathrm{NH}}_3^{+}$), Z, and anion (COO⁻–CH₂–NH₂), A. Therefore, we have considered the following dimers: NC and ZC for protonated species; NN, ZZ, and NZ for neutral ones; and NA and ZA for deprotonated dimers. For each of these dimers, we have built about 30 different geometries based on former papers by Chocholousova et al., who performed an extense set of MP2/6-31G(d,p) optimizations in the gas phase, [17], Carvalho et al. [18], and Friant-Michel et al. [22]. However, we extending them by including other types of IHB between monomers. We have also taken into account: (i) different possibilities to form N–H···O=C intramonomeric hydrogen bond (imHB); and (ii) internal rotations around the N–C–C=O dihedral angles in each monomer. These conformations can be roughly classified as: cyclic, with one or diverse kinds of IHBs, linear (with only one IHB), stacked (referring to roughly parallel main planes for the monomer skeletons) and gauche arrangement of monomers (denoted sometimes as T-shaped geometries).

The Gaussian 09 package [24] was used to perform complete geometry optimizations with the usual density functional (DFT) B3LYP, the metahybrid DFT M06, and MP2. The 6-311++G(d,p) basis set was used in all cases. Diffuse and polarization functions were included to allow better description of HBs [25,26]. Summing all in the seven groups of dimers, more than 200 geometries were optimized for each computational level and medium. This differs from what was done in a previous work, where all the geometries employed were only optimized at the B3LYP/6-31+G(d,p) in the gas phase summing all in the six groups of dimers [22]. When the geometries optimized with both methods were significantly different, we resorted to M06 to contrast results. Also, we performed single-point MP3 calculations and CCSD optimizations, when MP2 optimizations reversed the energy trends observed at other computational levels. All the energetic values shown in tables include ZPVE and energy thermal corrections. The vibrational analysis performed confirmed that optimized structures are minima with no imaginary frequency.

An aqueous solution was modelled, making use of the polarizable continuum model (PCM) to account for the solvation effects [27]. In all cases, PCM optimizations were performed using the same initial geometries employed in gas phase calculations. As the calculation with PCM frequencies is not accurate enough we have not estimated ZPVE and thermal corrections for PCM optimized structures [28]. Moreover, we have not assumed that gas phase corrections could be transferred into an aqueous solution because the geometries are, sometimes, very different.

The AIMPAC package [29] was used to perform the QTAIM electron density analysis. All atomic properties were computed with |L(Ω)| < 2 × 10⁻³ au, which usually guarantees their quality [30,31]. We also checked that electronic molecular population and energy were reproduced within 10⁻³ au and 2 kJ·mol⁻¹, respectively, by those obtained as summation of atomic properties, that is: ΣN(Ω) and ΣE(Ω). Our discussion also considers electron density at bond-critical points, ρ_b [23].

3. Results and Discussion

3.1. Conformational Analysis

The most stable structures obtained for each phase (gas and water solution) and group of dimers (NN, NZ, ZZ, NA, ZA, NC, ZC) at the MP2/6-311++G(d,p) level are described below (Table S1 reports Cartesian coordinates). There are significant differences between gas phase-optimized structures and those obtained in the aqueous solution. Both MP2 and B3LYP calculations provide more negative energies for zwitterions than for nonionic structures in the neutral monomers of the dimers in the aqueous solutions for all pH ranges, whereas the reverse trend is observed in the gas phase. Moreover, MP2 and B3LYP optimizations produce the same preferred conformers, excluding ZA series in the aqueous medium and NC series in the gas phase. As expected H₂N···H–O and H₂N···H–C hydrogen bonds are not observed, whereas H₂N···H–N bonds appear in several of the most stable structures.

3.1.1. Neutral Dimers in the Gas Phase

As previously found [17,18,19,20,21,22], the most stable dimer in the gas phase is formed by two non-ionic monomers (NN) and displays a cyclic structure with two O–H···O IHBs (Figure 1). The diverse arrangements of the two –CH₂–NH₂ residues give rise to different rotamers whose relative energies span in an 8.9 kJ·mol⁻¹ range according to MP2 calculations (10.2 kJ·mol⁻¹ with B3LYP). At both computational levels, the most stable rotamer is that where the N lone pair, Lp, is antiperiplanar to the C–C bond in both residues, while N–C and C=O bonds are synperiplanar. Diverse structures with a different arrangement of the main skeleton were also found as relative minima. They include dimers where the monomers are stacked (displaying N–H···O imHBs), arranged in T-shape (establishing O–H···O and N–H···O IHBs), or form an eight-membered cycle twisted with O–H···O and C–H···O IHBs, etc.

The energy sequence changes significantly from B3LYP to MP2 optimizations. Overall, our results confirm that NN cyclic dimers are more stable than stacked ones, contrary to what had been found by exploring the potential energy surface of this system [17] with the empirical potential of Cornell et al. [32].

The sequence of stabilities for the groups of dimers according to B3LYP is: NN > NZ > ZZ (

Table 1. Relative energies, ∆₁E, between the most stable dimer for the diverse zwitterion/non-ionic forms of cationic, neutral, and anionic species, ∆₁E. All values (in kJ·mol⁻¹) are relative to the most non-ionic form (NC, NN or NA, respectively), but CCSD ones.

Method	Gas			Aqueous Solution
	Cation	Neutral	Anion	Cation	Neutral	Anion
B3LYP	10.4	60.6 a (47.2) b	7.8	1.7	−34.2 a (−20.7) b	−15.7
MP2	13.8	47.5 a (55.2) b	12.7	−7.6	−34.1 a (−19.6) b	−23.6
M06	-	68.6 a(61.8) b,c	-	−1.7	−29.8 a (1.0) b,d	-
MP3	-	85.3 a (73.6) b	-	-	-	-
CCSD	-	9.9 e	-	-	-	-

^a ZZ-NN values. ^b NZ-NN values. ^c 61.8 kJ·mol⁻¹ corresponds to the M06 optimization of the same angular structure obtained as the most stable one according to MP2 and B3LYP optimizations. Nevertheless, the most stable NZ conformer according to the M06 study is a stacked dimer, whose NZ-NN relative energy is 59.0 kJ·mol⁻¹. ^d The positive value (1.0) corresponds to the same conformer obtained as the most stable one, with B3LYP and MP2 geometry optimizations (it displays N-H···N IHB). The most stable NZ dimer obtained with M06 for aqueous solution shows stacked structure, and its NZ-NN value is −11.5 kJ·mol⁻¹. ^e ZZ-NZ values.

). Although MP2 optimizations reverse the stability of NZ and ZZ, CCSD optimizations and single-point MP3 calculations on MP2 structures provide the same energy sequence obtained with B3LYP. We notice that whereas O–H···O=C IHBs appear in the most preferred geometries of NN dimers, N–H···O=C and C–H···O=C IHBs are displayed in the most stable NZ and ZZ structures (Figure 1).

Significant differences with previous studies, affecting to zwitterionic dimers (NZ and ZZ) are found:

• The most stable ZZ structures in the gas phase are stacked. The most favored one, ZZ-G1, (Figure 1) is different from those obtained in previous work [22]. The second and third most stables structures we obtain in this group, ZZ-G2 and ZZ-G3, (Figure 2) which are only 1.7 and 6.1 kJ·mol⁻¹ higher in energy than ZZ-G1, are closer in geometry to those reported previously as the most stable ones [22]. This sequence of conformational energies is found with B3LYP and MP2 levels.
• The geometry of the NZ dimer obtained as the most stable at the MP2 level, NZ-G1, differs significantly from the stacked geometry reported previously, NZ-GS, at the same computational level [22], which basically coincides with one of the B3LYP-optimized structures found here (Figure 3). Moreover, we notice that all our trials to perform MP2 optimizations for NZ stacked structures in the gas phase yielded NN ones. Nevertheless, at the B3LYP level NZ stacked structures were optimized, the most stable of them is only 10.7 kJ·mol⁻¹ higher.

3.1.2. Neutral Dimers in the Gas Phase

Aqueous solution computational studies for glycine dimers are scarcer. Friant-Michel and Ruiz-López [22] discussed three stacked ZZ structures, one linear NZ structure, and two cyclic NN ones. Hamad and Catlow [20] simulated three cyclic, doubly hydrogen-bonded ZZ structures in an aqueous solution using ab initio molecular dynamics. In the present work, the whole set of structures considered in the gas phase was optimized with PCM, both at B3LYP and MP2 levels. The following trends can be highlighted:

• The energy sequence changes to: ZZ < NZ < NN (Table 1).
• PCM optimized geometries are very different from those obtained in the gas phase. This is exemplified by the three lowest energy ZZ structures (ZZ-W2 and ZZ-W3 relative energies are, respectively, 2.5 and 3.0 kJ·mol⁻¹ at the MP2 level) (Figure 1 and Figure 2).
• Overall, ZZ aqueous solution preferred structures to be more opened than their gas phase counterparts. We notice they keep two N–H···O IHBs but lose one C–H···O and the imHB because of more favorable intermolecular interactions with water molecules.
• Preferred geometries keep more electron density on the most electronegative atoms.

In NZ dimers, the predominant structures in water solution are also different to those of the gas phase (Figure 1). They even display different kinds of IHBs. Thus, the most stable NZ geometry in an aqueous solution has one N–H···N IHB (which cannot be established between the two ${\mathrm{N}\mathrm{H}}_3^{+}$ units in ZZ). This conformer is about 4.4 kJ·mol⁻¹ more stable than that obtained by PCM-optimization of the gas phase preferred geometry. We also observe that whenever a solvated dimer includes one ⁺NH₃ unit, the most stable structure shows one N⁺–H···N IHB. This is even observed in ZC dimers in the gas phase. This allows solvation of carbonyl groups by water molecules, increasing dimer stability.

There are multiple NN dimers in an aqueous solution with similar energies, but very different in geometries (stacked, cyclic, linear and T-shaped) and showing different kinds of IHBs. We also observe that their most stable geometries differ from the gas phase (planar) to water solution (stacked) (Figure 1).

3.1.3. Anionic Dimers in the Gas Phase

Deprotonation of neutral forms gives rise to ZA or NA dimers. In the gas phase the NA dimers are preferred over the ZA ones (Table 1), and the reverse trend is observed in the aqueous solution. The most stable NA structure in gas phase displays O–H···O and N–H···O IHBs while the most stable ZA structure has two IHBs connecting ${\mathrm{N}\mathrm{H}}_3^{+}$ and OCO⁻ units (Figure 4).

The relative energy of the next stable NA form only exceeds 1 kJ·mol⁻¹ and shows a C–H···O IHB instead of the N–H···O one (Figure 5). Other conformers with N–H···O IHB and more lineal structures are higher in energy.

ZA-G1 and ZA-G2 conformers (nearly isoenergetic and with two N–H···O IHBs sharing the same N atom) get some kind of primacy among ZA dimers. In general, rotations around dihedral angles that keep the same IHBs do not provoke large energy differences. In contrast, two different types of IHB (N–H···O and C–H···O) are observed in other ZA dimers (ZA-G3 in Figure 5).

3.1.4. Anionic Dimers in an Aqueous Solution

We find the same situation described some lines above. When the ${\mathrm{N}\mathrm{H}}_3^{+}$ group exists (ZA dimer), the N⁺–H···N IHB is again present in the most stable geometry. When the ${\mathrm{N}\mathrm{H}}_3^{+}$ group disappears (NA dimer), the presence of only one O–H···O IHB gives rise to the most efficient solvation.

ZA in the aqueous media is one of the cases where B3LYP and MP2 do not lead to the same preferred geometries. According to MP2 the most stable structure, ZA-W1, bears one strong N⁺–H···N IHB and a weak C–H···C interaction, whereas ZA-W2, with N⁺–H···O and C–H···O IHBs is more than 6 kJ·mol⁻¹ less stable. The larger stability of ZA-W1 over ZA-W2 may be attributed to the availability of four oxygen atoms to establish IHBs with solvent molecules, while in ZA-W2 there are only two oxygens and one ${\mathrm{N}\mathrm{H}}_3^{+}$ available to establish them. In contrast, according to B3LYP ZA-W1 similar geometries (bearing N⁺–H···N IHBs) are around 4 kJ·mol⁻¹ less stable than a structure similar to ZA-W2, which is obtained as the preferred one. The conformational description provided by M06 calculations agrees with that of MP2. This leads us to question that maybe B3LYP does not perform a correct evaluation for N⁺–H···N IHBs.

The NA series of dimers (less stable than ZA in the aqueous solution) contains four conformers with very close energy (Figure 6). Specifically, NA-W1, NA-W2, and NA-W3, only differ in the dihedral angle containing the O–H···O IHB.

3.1.5. Cationic Dimers in the Gas Phase

Previous studies on protonated glycine dimers indicate the preference for NC forms in the gas phase [12,13,14,15,16]. Although Price et al. [12] reported that the most stable structure contains a N–H···N IHB, a consensus seems to be reached in the remaining studies [13,14,15,16] supporting the presence of N–H···O IHB in the most stable dimer. Finally, Armentrout et al., employing MP2(full)/6-311+G(2d,2p)//B3LYP/6-311+G(d,p), reported that a low-energy NC geometry (not fully converged in the frequency calculation) displays one N⁺–H···N IHB, although the most stable geometry displays N⁺–H···O IHB [16]. All the geometries considered in previous studies were included in our conformational sample for NC and ZC dimers. Our results, both B3LYP and MP2, agree with previous ones, indicating that NC structures are more stable than ZC in the gas phase. Moreover, the most stable of NC dimers displays N⁺–H···N IHB (Figure 7).

In contrast, the most stable structure for ZC contains N⁺–H···O and C–H···O IHBs. The second-most stable ZC geometry (nearly isoenergetic) rotates the N⁺–H···O IHB, resulting in a reinforcement of this bond while the C–H···O IHB is broken.

The formation of NC dimers with other IHB pattern costs at least 3.5 kJ·mol⁻¹, according to MP2 calculations. Thus, we observe a second NC geometry with a N⁺–H···O IHB, (Figure 8). Whereas the MP2 results for NC dimers in the gas phase contradict previous experimental and theoretical studies, our B3LYP study agrees with them. In fact, the N⁺–H···O structure is 8.2 kJ·mol⁻¹ below the N⁺–H···N one, and the geometry of the former coincides basically with ZC-G2 while that of the latter is substantially equal to ZC-G1 (Figure 7 and Figure 8).

3.1.6. Cationic Dimers in an Aqueous Solution

In the aqueous solution, MP2 results favor the ZC forms over NC ones. The contrary trend is obtained with B3LYP calculations. Although, M06 results agree with MP2 ones, the energy differences are certainly small. So, taking into account the limitations of PCM calculations, as well as the large series of NC and ZC possible conformers, we think that the simultaneous presence of ZC and NC dimers cannot be discarded in an aqueous solution.

The most stable geometry according to MP2 has ZC structure with N⁺–H···O IHB, similar to that obtained as the most stable within the ZC series in the gas phase (Figure 7), although the C–H···O bond path is replaced by a C···O one. The next conformer in stability (ZC-W2) results from a distortion of the Z monomer that allows the establishment of one imHB N⁺–H···O HB. In contrast, when the C–H···O and N⁺–H···O IHBs are kept in a water-solvated structure, the energy rises by 8 kJ·mol⁻¹ (ZC-W3), which is less stable compared to the most stable of NC dimers (NC-W1). This NC conformer displays N⁺–H···N IHB (Figure 7). Other NC conformers keeping these IHBs (e.g., NC-W2) show very similar relative energies, while ZC-W1-like geometries (as NC-W3) lead to higher energies (by 15 kJ·mol⁻¹ at least) (Figure 8).

3.2. Energy Trends

3.2.1. Zwitterionic versus Non Ionic Monomers

The energy difference between the most stable zwitterion and non-ionic structure for each global charge, ∆₁E, (Table 1), is positive in the gas phase, meaning that the most stable form is non-ionic, in agreement with experimental results. This trend is reversed in the aqueous solution, where the zwitterionic structures are preferred (negative Δ₁E values in Table 1). Cationic dimers are an exception with NC and ZC most stable solvated structures showing very similar energies at all the computational levels here employed. It is also noticeable that the energy difference between dimers containing N and Z monomers decreases in absolute value, both in the gas phase and in the aqueous solution, for anionic and cationic dimers with regard to the neutral ones (Table 1).

3.2.2. Solvation Energies

The difference between the energy obtained in the aqueous solution modelled by PCM and that of the gas phase, which is always negative, is represented by ∆₂E (

Table 2. Solvation energies, −∆₂E, (kJ·mol⁻¹) obtained subtracting the energy of the most stable aqueous solvated structure minus the most stable gas phase one ^a.

Dimers	B3LYP	MP2	M06
NC	295.3	288.9
ZC	304.5	311.6
NN	55.4	59.4	13.2
NZ	134.5	134.7	81.7
ZZ	159.9	142.0	110.1
NA	251.3	241.1
ZA	271.9	279.3

^a Thermal corrections are not included.

). We notice that the absolute value of ∆₂E is bigger in zwitterionic dimers than in non ionic ones. This means that Z-containing dimers are more stabilized by the aqueous solutions than their N-counterparts.

The lowest |∆₂E| values are observed for neutral dimers, whereas cations are more stabilized in water solution than anions. It is expected that hydration energies of cations are more negative than those of anions of similar size [33]. The reason for this trend is that solvating a cation in water implies breaking less IHBs, rather than solvating an anion of similar size. Nevertheless, the energy involved in breaking IHBs among solvent molecules is not taken into account by PCM (even in the non-electrostatic cavitation term [27]). Looking for an explanation, we have considered the polarization and distortion terms shown in

Table 3. Deformation and Polarization terms in MP2 solvation energies, and non-electrostatic terms for aqueous solution. All values in kJ mol⁻¹.

Dimers	Deformation	Polarization	Non Electrostatic a
NC	65.2	−354.1	30.3
ZC	75.0	−386.6	31.1
NN	57.8	−117.2	34.3
NZ	117.2	−251.8	33.8
ZZ	100.3	−242.2	32.9
NA	57.5	−301.8	38.3
ZA	147.0	−426.3	33.3

^a This term is not included in Δ₂E in Table 2.

. This partitioning reveals the leading role of polarization. Nevertheless, the largest polarization is shown by the anionic forms (ZA) which also experience the largest geometry deformation. As both effects counteract each other, the largest solvation energy corresponds to ZC, even when effects on solvent are not considered explicitly.

The sequence of polarization energies (Table 3) can be explained considering three factors: (i) the total charge of the molecule, which when other factors remain constant should lead to a sequence favoring anions over cations, and all types of them over neutral species; (ii) the charge separation within the system, favoring zwitterionic forms over non-ionic ones; and (iii) the geometry of the dimer, producing larger polarization energies when the most charged groups are external. Thus, conformers with N–H···N IHBs are more prone to increase their polarization energy than those including O–H···O IHBs. This explains why the polarization energy is larger for the most stable conformer of NC dimers than that of NA ones.

3.2.3. Protonation and Deprotonation Abilities

Neutral zwitterions are more prone to be protonated or deprotonated than non-ionic forms in the gas phase, whereas the reverse trend is observed in water solutions (

Table 4. Protonation and deprotonation energies (in kJ·mol⁻¹) for most stable conformers of diverse glycine dimers according to MP2 and B3LYP calculations.

Process	MP2	B3LYP	MP2	B3LYP
Protonation	Gas	Gas	Aq.sol.	Aq.sol.
NZ→CZ=ZC	991.7	987.8	1204.5	1195.2
NZ→NC	1005.4	998.2	1196.9	1196.8
ZZ→ZC	984.0	1001.3	1189.9	1181.7
NN→NC	950.2	951.0	1221.5	1232.5
Deprotonation	Gas	Gas	Aq.sol.	Aq.sol.
NZ→AZ=ZA	1320.4	1313.9	1211.4	1210.0
NZ→NA	1306.3	1306.1	1235.0	1225.7
ZZ→ZA	1328.1	1300.4	1226.0	1223.5
NN→NA	1362.9	1353.3	1210.3	1190.0

). This is observed both when we compare the protonation/deprotonation abilities of N and Z monomers in NZ dimers and when the comparison is established between NN and ZZ dimers. The only case where the computational level modifies the results is the protonation of NZ dimers in water solutions. This exception can be related to the very different geometries obtained between B3LYP and MP2 optimizations for the most stable structures of ZC forms (Figure 9).

3.3. Electron Density Analysis

3.3.1. Bond Properties

All the dimers obtained as most stable in each series are stabilized through IHBs. Depending on the series and phase, three kinds of strong IHB are observed: O–H···O, N–H···O, or N–H···N. Sometimes, they are combined with C–H···O weak IHBs. All of them are represented by bond-critical points, BCPs, in the QTAIM molecular graphs [23]. BCPs corresponding to C–H···C and C···O bond paths are also observed in some particular cases for aqueous solvated dimers (ZA and ZC, respectively). Furthermore, bond paths associated with N–H···O and O–H···N imHBs are observed in several monomers (Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8).

Electron densities at the BCP, ρ_b, are usually taken as relative indicators of bond strength within a series of similar bonds. Some examples preclude considering a direct relationship between preferred forms and ρ_b values. Thus, in the gas phase, the most stable ZZ conformer displays larger ρ_b values than the NN one (Figure 1), in spite of the NN form being the preferred one for neutral dimers in the gas phase (Table 1). By contrast, in water solution, the preferred neutral dimers (ZZ) display larger ρ_b values than NN.

We notice that O–H···O IHBs are only present in dimers that contain one neutral non-ionic monomer, N, and that their ρ_b values are larger when there is no other IHB in the dimer (NA-W1). Moreover, ρ_b values for O–H···O IHBs decrease progressively as the strength of the other IHB in the dimer rises. Thus, ρ_b for O–H···O is 0.080 au in NA-G2 (accompanied by one C–H···O IHB), 0.073 au in NZ-G1 (accompanied by N–H···O IHB), and 0.040 au in NN-G1 (accompanied by other O–H···O IHB). This effect can be related to geometry distortions of the O–H···O IHB. A good example is provided by the series of most stable NA dimers in the aqueous solution: NA-W1, NA-W2 and NA-W3, where the most stable form (with a 90° dihedral angle between both carboxyl groups) displays the largest ρ_b value and the smallest H···O bond distance (Figure 4 and Figure 6). Internal rotation between both carboxyl groups increases the H···O bond length and reduces ρ_b values.

N–H···O is the most commonly observed IHB in the diverse glycine dimer structures here obtained. Their ρ_b values span between 0.081 au (similar to O–H···O largest one) in ZC-G1 and 0.013 au in NN-G2. These values depend on the H···O bond length and N–H···O bond angle. Thus, the minimum values are observed in stacked geometries. A particular case is observed for the most stable structures of ZA dimer in the gas phase (ZA-G1 and ZA-G2), where the two N-H bonds of the zwitterionic monomer are involved in IHBs with the two oxygens of the anion (Figure 4 and Figure 5). Analyzing ρ_b values, we observe that the preferred structure has intermediate values (0.058 and 0.024 au) between those observed for ZA-G3 (0.085 and 0.009 au).

N–H···N IHBs display ρ_b values similar to those of most of N–H···O ones (around 0.06 au, see for instance ZC-W1 and NC-W1). This IHB is present in the most stable structures of glycine dimer whenever it contains only one ${\mathrm{N}\mathrm{H}}_3^{+}$ group. This IHB is more often observed in water solutions as it provides larger separation between COO⁻/COOH groups that allows solvation of both groups by water molecules, together with a relatively strong IHB.

ImHBs paths are related to the presence of certain coplanar arrangements. They appear between the carboxylate unit and one of the N–H bonds belonging to a ⁺NH₃ group, or between the nitrogen lone pair and the COOH group (only observed in NC-G2 and NN-W1).

3.3.2. Electron Density Transference between Monomers

The computation of QTAIM atomic properties allow us to analyse the electron density transfer between monomers that may take place in dimer formation (Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8). Atomic electron populations, N(Ω), for the atoms involved in intermolecular hydrogen bonds were analysed for the gas phase and the aqueous solvated dimers.

We notice that most of neutral dimers display zero or negligible electron density transference between monomers. In some cases null transference is due to symmetric structures, like NN-G1.

There are four exceptions where the transference is significant, all of them have been found in NZ dimers in gas or aqueous solutions. In three of them the electron density goes from Z to N monomer (NZ-G1, NZ-G2, and NZ-W2), but in other cases we observe the reverse transference (NZ-W1). We observe that dimers where the electron transference is, from Z to N, display: one N–H···O IHB with H donated by the zwitterionic ${\mathrm{N}\mathrm{H}}_3^{+}$ group, and one O–H···O IHB with neutral monomer donating the carboxylic H. By contrast, NZ-W1 displays one N–H···N interaction. Comparison of ρ_b values of N–H···O and O–H···O IHBs in NZ-G1 and NZ-G2, which display very similar electron transferences, indicates that the distortion of O–H···O IHB has little effect on the electron density transference.

Contrary to what has been found in neutral dimers, all the low-energy anionic and cationic dimers show significant electron density transferences. In all cases the anionic monomers donate electron density in an amount that is more affected by the kind of IHBs formed than for the non-ionic or zwitterionic character of the other monomer. In cationic dimers, the electron density is taken from the neutral monomer, and once more the amount of transference depends on the set of IHBs established.

Diverse studies previously carried out in a series of simple hydrogen bond complexes, as well as conformational analysis [34,35,36], concluded that relative energies could be explained on the basis of the summation of atomic electronic populations of hydrogens. This is not the case, as could expect taking into account the disparity of IHBs (and other factors) here involved. We notice that the most stable geometries in the gas phase display the largest summation of the atomic electron population for hydrogens. Nevertheless, this trend is not followed in the aqueous solution, where the most stable forms display the largest N(O) values. Summation of N(H) values in the neutral monomer display a clear dependence with the charge of the other attached monomer. They are reinforced by attached anions and depleted by cations.

Two noticeable effects related to the medium that affect all sort of dimers are: (a) The summation of N(H) values in each dimer is always bigger in the gas phase than in the aqueous solution; (b) On the contrary, the summation of N(O) values is lower in the gas phase. All of this is a consequence of solvation as it is favored by increased absolute values of atomic charges.

4. Conclusions

The large collection of structures optimized for different series of glycine dimers (neutral, anionic and cationic) has provided new most stable geometries, not previously described, for the ZZ series in the aqueous solution.

The stability sequence for stacked, cyclic, and linear dimers depends on the phase and total charge of the dimer. According to MP2 results, the most stable dimers that include one ${\mathrm{N}\mathrm{H}}_3^{+}$ group show N···H–N IHBs. This conformational preference is not confirmed at the B3LYP level in all the cases. Particularly, in the NC series, B3LYP optimizations, indicate that the most stable dimer displays O···H–N IHB, in agreement with all published interpretations of i.r. Gly₂H⁺ spectra [12,14,16].

Energy values show that: (i) zwitterions are not favored in the gas phase; (ii) water solvation is more exergonic when dimers contain Z monomers than when include N ones; and (iii) In the gas phase, Z monomers are more feasible to experience protonation and deprotonation than N ones, while the reverse trend is observed in the aqueous solution.

In general, O···H–O IHBs display the largest electron densities at the corresponding BCPs (up to 0.08 au), whereas this quantity is usually lower for O···H–N (spanning in a wide range depending of the presence of other IHBs in the dimer) and N···H–N IHBs (around 0.06 au).

N···H–O and N···H–C IHBs were not found in any optimized structure, but O···H–C bond paths are found accompanying other IHBs. Overall, ρ_c values are usually smaller in the aqueous solution than in the gas phase, indicating that solvation weakens the IHBs, as usually observed [37]. In contrast, the electron population of the hydrogen which belongs to this IHB increases in the aqueous solution with regard to that of the gas phase.

Most of the charged dimers display important electron density transference between monomers (around 0.1 au). The anionic monomer acts as electron density donor (ZA or NA), whereas the cation acts as acceptor (ZC or NC). In contrast, most of the neutral dimers display zero or negligible electron density transference between monomers.