3.2.4. Virial Ratios, |V(**r***c*)|/G(**r***c*)

Analysis of the virial ratios at bond critical points serves as a more quantitative description of the nature of the interactions than the Laplacians of the electron density and the bond degree parameters. See the works of Grabowski [28] and of Rozas et al. [68] for a formal description of how the virial ratio is related to bonding. In short, local depletion of electron density (local dominance of the repulsive kinetic energy), which is indicative of ionic or long interactions have 0 < |V(**r***c*)|/G(**r***c*) < 1, local concentration of electron density (local dominance of the attractive potential energy), indicative of covalent interactions have |V(**r***c*)|/G(**r***c*) > 2, and the 1 < |V(**r***c*)|/G(**r***c*) < 2 interval describes interactions with mixed contributions.

Figure 7 shows the distribution of the virial ratios for all dimers found in this work, which cover the [0.61, 1.49] interval for primary hydrogen bonds and salt bridges (Figure 7A–C) and the [0.56, 0.92] interval for secondary HBs, exotic and dihydrogen contacts (Figure 7A–C). Notice that as in the analysis of the previous descriptors, the fitted distributions go a little beyond the actual limits. It is quite revealing that a large number of contacts, especially those involving the carbonyl group have virial ratios larger than 1, which confers them a high degree of covalency while not being formal bonds. Interestingly, these include the charged carboxylate which may naively be thought as being involved in highly ionic contacts. Virial ratios larger than 1 transcend the carbonyl group, which is indeed the case for the following HBs: N· · · H–O, C=O· · · H–O, S· · · H–O, N· · · H–S for *n* · · · *n* dimers, C=O<sup>−</sup> · · · H–O, N· · · H–N+, S· · · H–N+, C=O · · · H–N<sup>+</sup> for *<sup>n</sup>* · · · *<sup>z</sup>* and for all salt bridges in *z* · · · *z*. This high covalency of the *a priori* ionic contacts has been reported for other cases, including for example the microsolvation of charged species [49,53,72]. Most secondary HBs, H· · · H, and exotic contacts have virials smaller than the water dimer reference.

Surprisingly, the thiol group in a number of cases is involved in stronger interactions than the H–O· · · H–O and N· · · H–O contacts.

**Figure 6.** Bond degree parameters for intermolecular contacts for all dimers of cysteine found in this work. The top row shows only primary hydrogen bonds including salt bridges. The bottom row shows secondary hydrogen bonds, dihydrogen bonds, and all exotic interactions. The left column is reserved for the *n* · · · *n* dimers (subfigures (**A**,**D**)), the middle column for *n* · · · *z* (subfigures (**B**,**E**)) and the right column for *z* · · · *z* (subfigures (**C**,**F**)). All values taken from the B3LYP–D3/6–311++G(*d*, *p*) potential energy surfaces with water represented as a continuum solvent. Solid vertical lines mark the QTAIM boundaries separating locally stabilizing from the locally destabilizing interactions. Dashed vertical lines mark the reference value for the gas phase water dimer.

#### 3.2.5. NBO and NCI Picture of Intermolecular Interactions

Intermolecular interactions have been successfully studied under the NBO formalism in a wide range of problems [30,33]. In the particular case of the cysteine dimers, we proceeded to identify the localized donor Lewis orbitals from which charge is transferred to acceptor orbitals according to the *φ<sup>d</sup>* → *φ<sup>a</sup>* scheme. Tables 3 and 4 list the involved orbitals for each one of the 80 types of interactions found in this work, Figure 1 provides the corresponding surfaces for those dimers with populations larger than 5%. Once identified, we quantified the strength of the orbital interaction by second order perturbation theory on the Fock matrix as given by −*E* (2) *<sup>d</sup>*→*<sup>a</sup>* <sup>=</sup> *<sup>q</sup><sup>d</sup>* |h*φ<sup>d</sup>* |F|*φa*i|<sup>2</sup> /(*E<sup>a</sup>* − *E<sup>d</sup>* ). With this procedure, the strength of the interaction is directly related to the magnitude of *E* (2) *d*→*a* .

**Figure 7.** Virial ratios at bond critical points for intermolecular contacts for all dimers of cysteine found in this work. The top row shows only primary hydrogen bonds including salt bridges. The bottom row shows secondary hydrogen bonds, dihydrogen bonds, and all exotic interactions. The left column is reserved for the *n* · · · *n* dimers (subfigures (**A**,**D**)), the middle column for *n* · · · *z* (subfigures (**B**,**E**)) and the right column for *z* · · · *z* (subfigures (**C**,**F**)). All values taken from the B3LYP–D3/6– 311++G(*d*, *p*) potential energy surfaces with water represented as a continuum solvent. Vertical solid lines mark the QTAIM boundaries separating long range from intermediate character interactions. Dashed vertical lines mark the reference value for the gas phase water dimer.

Donor → acceptor orbital interactions resulting in primary hydrogen bonds and salt bridges in *n* · · · *n*, *n* · · · *z* and *z* · · · *z* cysteine dimers include everything from the very weak to the very strong, covering the wide [0.06, 50.30] kcal/mol interval (Figure 8A–C). Salt bridges exhibit an uncommonly complex distribution of energies with several shoulders. Interestingly, the strongest contacts are not salt bridges but rather N· · · H–O primary hydrogen bonds, listed as interaction **10** in Table 3 and shown in Figure 8A. This interaction type, which is also the strongest intermolecular contact found in alanine dimers [13], arises from *n*<sup>N</sup> → *σ* ∗ H−O charge transfer in *n* · · · *n* dimers. Next in the strength hierarchy are the highly ionic C=O<sup>−</sup> · · · H–O and N· · · H–N<sup>+</sup> contacts arising from *<sup>n</sup>*<sup>O</sup> <sup>→</sup> *<sup>σ</sup>* ∗ H−O and *n*<sup>N</sup> → *σ* ∗ H−N charge transfers in *n* · · · *z* dimers. These are listed as interactions **2, 12** with the corresponding distributions shown in Figure 8B. C=O<sup>−</sup> · · · H–N<sup>+</sup> salt bridges in *<sup>z</sup>* · · · *<sup>z</sup>* dimers come only third in the interaction energy scale. They arise from *n*<sup>O</sup> → *σ* ∗ H−N charge transfer, are listed as interaction **6** and the corresponding distributions are shown in Figure 8C. These sets of interactions are present on the structures with populations higher than 5%. In a manner consistent with the QTAIM descriptors analyzed above, a large number of primary HBs and salt bridges have interaction energies larger than 6.63 [30] kcal/mol, the orbital interaction energy for the reference water dimer, however, no secondary hydrogen bond, no exotic contact and no dihydrogen bond exceed the reference.

**Figure 8.** Donor · · · acceptor NBO energies for intermolecular contacts for all dimers of cysteine found in this work. The top row shows only primary hydrogen bonds including salt bridges. The bottom row shows secondary hydrogen bonds, dihydrogen bonds, and all exotic interactions. The left column is reserved for the *n* · · · *n* dimers (subfigures (**A**,**D**)), the middle column for *n* · · · *z* (subfigures (**B**,**E**)) and the right column for *z* · · · *z* (subfigures (**C**,**F**)). All values taken from the B3LYP–D3/6– 311++G(*d*, *p*) potential energy surfaces with water represented as a continuum solvent. The dashed vertical lines mark the reference value for the gas phase water dimer.

As stated above, primary HBs and salt bridges determine the structure of the dimers. According to Table 3, they always arise from orbital interactions of the *n*<sup>X</sup> → *σ* ∗ H−Y type with X,Y = O, N, S. As also stated above, secondary HBs, dihydrogen bonds and exotic contacts are usually a consequence of the structure. Notwithstanding, the orbitals involved in the weaker interactions offer a quite interesting and rather uncommon picture. First, notice that all exotic O· · · O contacts (**53–58** in Table 4) put the two negative ends of the fragments with various degrees of negative character in direct contact, with the most severe case being interaction **54** connecting two formal negative charges with no intermediaries. Second, notice that the 139 interactions grouped into **58, 61–63, 65, 71, 76, 78–80** may all be described by the general *n*<sup>X</sup> → *σ* ∗ Y−H charge transfer scheme with X, Y = O, N, S. These correspond to what David et al. [13] have called inverted hydrogen bonds because the lone pair on X donates electron charge to an antibonding *σ* ∗ Y−H orbital which is inverted from the usual *σ* ∗ H−Y , in other words, the charge donation occurs between two orbitals overlapping from one negative atom to another negative atom with no bridging proton. Anti electrostatic hydrogen bonds of the type described in this paragraph have been reported by Weinhold and Klein [73].

Figure 9 shows the obtained non covalent surfaces as well as the thoroughs of the reduced gradients for the largest binding energy dimers in the *n* · · · *n*, *n* · · · *z*, *z* · · · *z* potential energy surfaces, we include the corresponding interacting NBOs to help visualize the source of the NCI surfaces. See Figures S3–S5 in the supplementary material for the corresponding surfaces in all dimers with populations larger than 5%. The standard NCI color code [34,35] assigns green surfaces to weakly bonding contacts and blue surfaces to strong interactions, for the thoroughs, negative values of sign(*λ*2)*ρ* reveal bonding interactions [34] which are weaker when sign(*λ*2)*ρ* ≈ 0. NCI reveals that the two salt bridges in the *z* · · · *z* dimers are actually very similar (in fact, they cannot be told apart in the thoroughs) and that in most cases, large stabilizing surfaces arising from individual contacts transferring tiny amounts of charge to the interstitial region have significant contributions to the overall binding of the dimers. Unexpectedly, these charge transfer contributions are major contributors to the charged cases. Notice that charge transfer to the interstitial region appears to be the norm when several molecular units are stabilized via non covalent interactions: these fluxional surfaces of charge have been found to be a major player in the molecular interpretation of hydrophobicity [46], in the initial recognition and attachment of viruses to cell receptors [47,48,71], in the microsolvation and encapsulation of charged and neutral species, in the microscopic structure of ionic liquids [51], etc.

**Figure 9.** NCI surfaces and thoroughs describing the intermolecular contacts for the highest binding energy cysteine dimers. NBO pictures are also included to ease visualization of the interactions.

#### **4. Discussion and Context**

Accurate description and characterization of chemical bonding is a notoriously hard problem in chemistry, whose difficulty is magnified when dealing with weak intermolecular non covalent interactions. When studying molecules and their interactions, a large portion of the conceptual framework developed by experimentalists and theoreticians invokes a number of useful ideas that correspond to non observable quantities (partial atom charges, orbital interactions, virial ratios at BCPs, etc.); thus, there are no quantum mechanical operators whose expected values may be used to calculate them and therefore, approximate methods, however accurate, are used to determine these quantities. This approach has a fundamental problem: each quantity may be obtained by several methods and the results quite often vary among them. In this context, it is impressive and certainly reassuring that QTAIM, NBO and NCI, which are conceptually and methodologically substantially

different, afford a consistent, complementary picture of intermolecular bonding in the dimers of cysteine.

The large number of interactions, isomers, and types of binary contacts are intimately connected to the problem of molecular evolution and to the complexity of life observable on this planet. By virtue of the large number of accessible states, molecular systems where specific cysteine to cysteine contacts are observed are thermodynamically favored because of the ever increasing entropy of the universe, in other words, this type of systems will evolve towards equilibrium states with large structural diversity. Since the interactions dissected here are responsible for the molecular interactions between all aminoacid pairs, this argument of entropy driving molecular evolution readily applies to large proteins and biomolecules.

#### **5. Summary and Conclusions**

An intensive exploration of the potential energy surfaces for the interaction of neutral and charged cysteine monomers to form dimers in an aqueous environment represented by a continuum afforded a large number of isomers, amounting to 746 well characterized local minima. The isomers with the largest population are distributed within small energy differences of the putative global minima. Ten potential energy surfaces were explored in total for the *n* · · · *n*, *n* · · · *z*, and *z* · · · *z* combinations with two neutral (*n*) and two zwitterionic (*z*) forms. A number of strongly bound dimers were found, with interaction energies exceeding 20 kcal/mol in several cases and with interaction distances covering the very small to the very large in the [1.50, 4.19] Å interval. The nature of intermolecular bonding interactions was dissected using QTAIM, NCI, and NBO, three conceptually different methods, which for the present case afford consistent, complementary pictures. A total of 80 types of different intermolecular contacts were found in this complex and large universe of dimers. As a general rule, primary hydrogen bonds and salt bridges are the strongest of the interactions and determine the molecular geometry, conversely, secondary hydrogen bonds, exotic X· · · Y (X = C, N, O, S) and H· · · H dihydrogen contacts are weaker and most often a consequence of the structure. All interactions, even the highly ionic, may be described by the *φ<sup>d</sup>* → *φ<sup>a</sup>* orbital charge transfer scheme, leading to accumulation or depletion of electron density at the bond critical points as revealed by topological analysis of the electron densities. The large binding energies mentioned above are the result of unusually strong charge assisted hydrogen bonds and salt bridges. We found that the highly ionic salt bridges have large degrees of covalency; thus, a simplistic electrostatic attraction between positively and negatively charged fragments does not suffice for a proper account of bonding interactions whenever the zwitterions are involved. The weaker secondary hydrogen bonds, exotic X· · · Y and dihydrogen contacts in no few cases are stronger than, for example, the archetypal hydrogen bond in the water dimer. Moreover, independent of how strong or weak individual interactions are, their collective action cannot be ignored because they lead to the formation of large attractive non-covalent surfaces in the interstitial region between the fragments. A few antielectrostatic contacts [73] as well as a few inverted hydrogen bonds [13] in which the charge transfer occurs between the two negative ends of the fragments, were found; thus, they appear to be of common occurrence in nature.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/molecules27248665/s1, Figure S1: Electrostatic potential surfaces for the cysteine monomers.; Figure S2: Distributions of the Laplacians of the electron densities at bond critical points for all dimers.; Figures S3–S5: Additional NBO and NCI plots of descriptors of bonding interactions for the dimers with populations larger than 5%; Table S1: Binding energies and energy differences for neutral dimers.; Table S2: Binding energies and energy differences for mixed dimers. Table S3: Binding energies and energy differences for zwitterionic dimers. Cartesian coordinates for the entire set of 746 dimers are also provided.

**Author Contributions:** Conceptualization, S.G. (Sara Gómez), C.C. and A.R.; Data curation, S.G. (Santiago Gómez), S.G. (Sara Gómez), J.D. and D.G.; Formal analysis, S.G. (Santiago Gómez), S.G. (Sara Gómez), J.D., D.G., C.C. and A.R.; Funding acquisition, C.C. and A.R.; Investigation, S.G. (Santiago Gómez), S.G. (Sara Gómez), J.D. and D.G.; Methodology, S.G. (Sara Gómez) and A.R.; Project administration, A.R.; Resources, C.C. and A.R.; Supervision, S.G. (Sara Gómez) and A.R.; Validation, S.G. (Santiago Gómez), J.D. and D.G.; Visualization, S.G. (Santiago Gómez); Writing—original draft, S.G. (Santiago Gómez), S.G. (Sara Gómez), J.D., D.G. and A.R.; Writing—review and editing, S.G. (Santiago Gómez), S.G. (Sara Gómez), J.D., D.G., C.C. and A.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Universidad de Antioquia, grant "Estrategia para la sostenibilidad".

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data are available within the article or Supplementary Materials.

**Acknowledgments:** Internal support from Universidad de Antioquia via "Estrategia para la sostenibilidad" is acknowledged. We gratefully acknowledge the Center for High Performance Computing (CHPC) at SNS for providing the computational infrastructure.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

