Henry Rzepa's blog

Previously, I explored the unusual structure of a molecule with a hydrogen bonded interaction between a phenol and a pyridine. The crystal structure name was RAKQOJ and it had been reported as having almost symmetrical N…H…O hydrogen bonds. This feature had been determined using neutron diffraction crystallography, which is thought very reliable at determining proton positions. Another compound with these characteristics is 3-methyl-5-phenylpyrazole or MEPHPY01.[1] Here the neutron study showed it to apparently have the structure represented below, where the solid N-H lines indicate a proton equidistant between two nitrogens.

Inspection of the ORTEP plot shows a very odd feature; the thermal elipsoids (red arrow) for two of the N-H-N motifs are more or less spherical, indicating little thermal motion (the temperature of the determination is not noted, and is assumed as probably room temperature) but the other two (magenta arrows) are highly elongated in the direction of motion between the two nitrogen atoms. This feature was largely unexplained at the time of publication (1975) and indeed to this day. Here I offer a possible insight into this enigma.

The conventional structure is shown below showing four N-H bonds and four H…N hydrogen bonds.

So now for the results of some calculations. Computed at various B3LYP(±GD3BJ)/Def2-SVPP/Def2-TZVPP levels (Table, FAIR data DOI: 10.14469/hpc/10406), the located minimum in the total energy, saddle=0, corresponds to the conventional proton-localized structure shown above, where all four hydrogens are firmly attached to the four nitrogen atoms by a regular bond and the distances are 1.032 for the NH and 1.855Å for the hydrogen bond it forms. A zwitterionic isomer comprises the ion-pair shown below, examples of each component of which are known in the CSD (crystal structure database). 

There are three ways of distributing this motif, of which only 1 is stable to proton transfer. Structure 2 has a higher degree of charge separation whilst 3 superficially appears to reduce the degree of charge separation compared to 1. In fact, the three-dimensional structure of 1 allows the negative ion to stack above the positive ion, thus actually achieving minimal separation of charges.

The stacking also depends on the type of calculation. If dispersion correction is included, the aromatic faces stack directly above each other (as above). If omitted, the stacking actually corresponds more closely to that observed in the reported crystal structure, since the attraction between faces occurs not only within a structure but between adjacent structures in the solid state (something not modelled when the dispersion correction is applied only to a single unit).

The lesson learnt from the previous post is that the position of protons as determined by quantum-chemical geometry location using minimisation of the total computed energy might be misleading. Better perhaps to use the computed free energy? When this is done, as we saw in the previous post, the transition state for proton transfer as located in the total energy surface can actually have a free energy that is lower than that of the total energy minimum. So, for MEPHPY01, a stationary point in which all four hydrogens correspond to the apparently symmetrical experimental neutron diffraction structure emerges as saddle=3, corresponding to three force constants being negative. The bond lengths for this geometry occur in pairs, two with NH 1.25/N…H 1.30 and two with 1.28/1.28, revealing interesting asymmetry.

The normal vibrational modes for these three -ve force constants are shown below. The first (ν 1315i cm-1) shows all four hydrogens exchanging between nitrogens, a quadruple proton transfer. The second (ν 896i cm-1) shows a double proton transfer between one pair exchanging between two nitrogens and the last (ν 801i cm-1) is similar in form, but shows the other pair exchanging between a second different pair of nitrogens. These last two vibrational modes correspond to the very thermal ellipsoids seen in the crystal structure diagram at the top, where one pair of hydrogens show little motion and the other pair involves much greater motion between a pair of nitrogen atoms.This would correspond to formation of a species exhibiting two conventional NH…H hydrogen bonds and two symmetrical N…H…N units.

Two further stationary points corresponding to saddle=2 and saddle=1 can also be located (Table).

-1986.690265(14.1)

Now here is the wacky thing. At the gas phase SVPP basis set ± dispersion levels, these lower-order saddle points are actually HIGHER in free energy than the third order saddle point! Conventional wisdom is that the higher the order of the saddle point, the higher should its energy be! I am not aware of anyone reporting an inverse observation before. The effect however is solvation and also basis-set dependent, since adding dichloromethane as a continuum solvent changes the free energy minimum from the third to the second-order saddle point. It might well also be dependent on the density functional method.

What are we to conclude? The free energy barriers for all the proton transfer saddle points computed above are not that small, being ≥ 10 kcal/mol. But at room temperatures, these exchanges will in fact be fast kinetic processes and the measured neutral diffraction structure may well emerge as averaged in some way. The free energies of the higher order saddle points suggests the dynamics of this system may in fact be very complex and very different from any “normal” hydrogen bonded system. This is clearly not the final word yet, but it does hint that the proton transfer dynamics of 3-Methyl-5-phenylpyrazole may be a system very well worth looking at again! And indeed exploring how robust the effects noted above are to different density functionals.

This post has DOI: 10.14469/hpc/10512

References

1. F.H. Moore, A.H. White, and A.C. Willis, "3-Methyl-5-phenylpyrazole: a neutron diffraction study", Journal of the Chemical Society, Perkin Transactions 2, pp. 1068, 1975. https://doi.org/10.1039/p29750001068

The previous examples of four atom systems displaying two layers of aromaticity illustrated how 4 (B₄), 8 (C₄) and 12 (N₄) valence electrons were partitioned into 4n+2 manifolds (respectively 2+2, 6+2 and 6+6). The triplet state molecule B₂C₂ with 6 electrons partitioned into 2π and 4σ electrons, with the latter following Baird’s aromaticity rule.[1],[2]. Now for the final missing entry; as a triplet C₂N₂ has 10 electrons, which now partition into 4 + 6. But would that be 4π + 6σ or 4σ + 6π? Well, in a way neither! Read on.

Bonding MOs for C2N2. Click image to load 3D model
π3, 1 electron	π2, 1 electron
σ3 2 electron	σ2, 2 electron
π1 2 electron	σ1, 2 electron

The calculations (ωB97XD/Def2-TZVPP and CCSD(T)/Def2-TZVPP) are collected at FAIR DOI: 10.14469/hpc/10346. These show a partitioning into 5σ + 5π, a species that is not a minimum but undergoes a non-planar distortion.

However, the first excited state (the triplet) IS planar and is only 12.5 kcal/mol above the planar 5+5 precursor. It is now partitioned into 6σ and 4π, with the latter conforming to Baird’s rule for open shell triplets.[1],[2] So this is unlike C₂B₂, which showed 2π + 4σ partitioning with the σ series following Baird’s rule. Now we have two examples in which one of the σ or the π-manifolds follow Baird’s rule and the other follows Hückel’s rule. The systems themselves are somewhat contrived, but they show the simple fun and games that can be had with these aromaticity rules.

This post has DOI: 10.14469/hpc/10350

References

1. N.C. Baird, "Quantum organic photochemistry. II. Resonance and aromaticity in the lowest 3.pi..pi.* state of cyclic hydrocarbons", Journal of the American Chemical Society, vol. 94, pp. 4941-4948, 1972. https://doi.org/10.1021/ja00769a025
2. M. Rosenberg, C. Dahlstrand, K. Kilså, and H. Ottosson, "Excited State Aromaticity and Antiaromaticity: Opportunities for Photophysical and Photochemical Rationalizations", Chemical Reviews, vol. 114, pp. 5379-5425, 2014. https://doi.org/10.1021/cr300471v

Normally, aromaticity is qualitatively assessed using an electron counting rule for cyclic conjugated rings. The best known is the Hückel 4n+2 rule (n=0,1, etc) for inferring diatropic aromatic ring currents in singlet-state π-conjugated cyclic molecules^‡ and a counter 4n rule which infers an antiaromatic paratropic ring current for the system. Some complex rings can sustain both types of ring currents in concentric rings or regions within the molecule, i.e. both diatropic and paratropic regions. Open shell (triplet state) molecules have their own rule; this time the molecule has a diatropic ring current if it follows a 4n rule, often called Baird’s rule. But has a molecule which simultaneously follows both Hückel’s AND Baird’s rule ever been suggested? Well, here is one, as indeed I promised in the previous post.

The species shown above has two carbons and two borons in a ring. These have a total of 14 valence electrons, eight of which occupy the C-B bonds, leaving six contributing to circulating ring currents. These partition into two π-electrons which then form a Hückel 4n+2 aromatic (n=0) and four σ-electrons which then form a Baird 4n aromatic (n=1) as a triplet. The triplet for this molecule is indeed its lowest state, 38.9 kcal/mol or 45.4 kcal/mol in free energy lower than the two lowest energy singlet states. These arise by placing two electrons in either of the two orbitals σ2 or σ3 each singly occupied in the triplet state (FAIR Data collection: 10.14469/hpc/10267)

Bonding MOs for C2B2. Click image to load 3D model
σ3	σ2
σ1
π1

So here we see a different sort of doubly aromatic molecule, to add to C₄, B₄ and N₄. With two electrons less than C₄, it is now doubly aromatic as a triplet state, this time conforming to two different electron counting rules. It would be good to know if any other examples showing this pattern are known.

^‡Hückel’s rule originally applied to p-π electrons in a cycle, such as benzene. Nowadays it is also used for σ in-plane electrons in a cycle.

This post has DOI: 10.14469/hpc/10271.

I discussed in the previous post the small molecule C₄ and how of the sixteen valence electrons, eight were left over after forming C-C σ-bonds which partitioned into six σ and two π. So now to consider B₄. This has four electrons less, and now the partitioning is two σ and two π (CCSD(T)/Def2-TZVPPD calculation, FAIR DOI: 10.14469/hpc/10157). Again both these sets fit the Hückel 4n+2 rule (n=0).

Since B₄ has only two rather than six delocalized σ-orbitals, the contributions to the central B-B bond are weaker and so the B-B bond is much longer.

Bonding MOs for B4. Click image to load 3D model
σ1, -0.335 au
π1, -0.372 au

Next, N₄.

π-Bonding MOs for N4. Click image to load 3D model
π3	π2
π1
σ-Bonding MOs for N4. Click image to load 3D model
σ3	σ2
σ1

The pattern for N₄ is different in several aspects. Firstly the π-system has six bonding electrons distributed over only four atoms. This makes the electron repulsions too high and the species is no longer stable, having one large imaginary force constant corresponding to an out-of-plane distorsion. Secondly the lowest energy σ orbital is highly localised onto two nitrogens rather than being delocalised around the ring periphery. So all those electrons crammed into a small space have taken their toll.

Thus far we have identified three species, B₄, C₄ and N₄ with interesting sets of respectively 4,8 and 12 electrons, all partitioned into 4n+2 collections. But what happens if one cannot do that; lets say 6 and 10 electrons? Hang around to find out!

The title comes from the abstract of an article[1] analysing why Biotin (vitamin B7) is such a strong and effective binder to proteins, with a free energy of (non-covalent) binding approaching 21 kcal/mol. The author argues that an accumulation of both CH-π and CH-O together with more classical hydrogen bonds and augmented by a sulfur centered hydrogen bond, oxyanion holes and water solvation, accounts for this large binding energy.

Here, I thought I would present a visualisation of the surroundings of biotin using the method of NCI (non-covalent-interaction) analysis, which looks at the behaviour of the electron density in the “weak” (i.e. non-covalent) regions of the biotin. This provides a more objective measure of the important interactions, independent of what we might consider important by virtue of having labels attached (such as e.g. “hydrogen bond”).

1. I started by getting the coordinates of streptavidin (DOI: 10.2210/pdb3RY2/pdb) a protein where biotin has been co-crystallised.[2]
2. Loaded into the CCDC Mercury program, I selected the molecule biotin itself and then added to the selection its close contacts with various groups in the streptavidin protein. These additions were truncated and capped with a methyl group to allow a wavefunction for the assembly to be calculated.
3. Hydrogens were then added to this structure to complete atom valencies, using “idealised” positions and ensuring that when rotamers were possible, they were set up to form hydrogen bonds.
4. A calculation (DOI: 10.14469/hpc/9982 at the ωB97XD/Def2-TZVPP/SCRF=water level) was performed.
5. The heavy atom coordinates (i.e. not hydrogens) are unaltered from the X-ray structure. Since atom positions as measured by X-ray diffraction and as computed using a DFT procedure are slightly different, the original coordinates were also subjected to three cycles of DFT-based geometry optimisation (DOI: 10.14469/hpc/9983) to better reflect the electron density in the molecule.
6. The resulting wavefunctions in the form of an .fchk file (for both unoptimised and partially optimised geometries) were then used to compute a grid of total electron density points
7. The density, in the form of a cube of points, was fed to Jmol using the commands
   load biotin_den.cub; isosurface parameters [0.5 1 0.0005 0.05 0.95 1.00] NCI ""; color isosurface "bgyor" range -0.04 0.04;
   and the resulting NCI surface was written out using the command write biotin.jvxl for inclusion here.
8. This is the NCI plot obtained from the raw coordinates from the PDB file.
   
9. This is the NCI plot obtained from the coordinates from the PDB file after three geometry optimisation cycles. Can you spot any differences?
   
   
10. These models are now available for you to explore by clicking on the images above.
    • Blue regions represent “strong” or classical hydrogen bonds. There are four of these in the NCI diagrams above and they are all compact, another characteristic of strong hydrogen bonds.
    • The hydrogen bond to sulfur is somewhat weaker, and appears in the display as a compact, albeit now cyan-coloured surface.
    • The remaining regions are both diffuse and green and represent weaker “interactions”. They are less compact than the classical hydrogen bonds. They do not represent a bond so much as an attractive region in the molecule and hence the term non-classical. Most are CH groups close to the π-surface of an aromatic ring, but some are also CH…O interactions.

Do go ahead and load the 3D surface. You should particularly explore the CH-π regions and note that they are not necessarily associated with a particular CH bond, but with several of these combining to form an interaction with an aromatic π region.

What might emerge is the realisation that binding interactions are not always between specific atoms as in classical hydrogen “bonds”, but also constitute “stabilising regions” between the ligand and the protein. You will probably spot several of these regions that are not actually listed in the article itself.[1] I suggest that we do not refer to CH…π bonds such as in the quoted title of this post but instead as CH…π regions.

It would be great if the entire complex could be subjected to an NCI analysis. Wavefunctions for >2000 atoms can be obtained nowadays, but it would require a bit of work to ensure the density can be computed accurately enough and at high enough cubic resolution to be useful in the context of NCI analysis.

This blog has DOI: 10.14469/hpc/9984

References

1. D.B. McConnell, "Biotin’s Lessons in Drug Design", Journal of Medicinal Chemistry, vol. 64, pp. 16319-16327, 2021. https://doi.org/10.1021/acs.jmedchem.1c00975
2. I. Le Trong, Z. Wang, D.E. Hyre, T.P. Lybrand, P.S. Stayton, and R.E. Stenkamp, "Streptavidin and its biotin complex at atomic resolution", Acta Crystallographica Section D Biological Crystallography, vol. 67, pp. 813-821, 2011. https://doi.org/10.1107/s0907444911027806