Author: Henry Rzepa

Personal Impressions from WATOC 2020 – Dispersion and non Born-Oppenheimer models.
WATOC 2020 was just held in 2022 in Vancouver Canada, over one week. With many lectures held in parallel, it is not possible for one person to cover anything like the topics presented, so this is a personal view of some of those talks that I attended. As happens with many such events, common themes gradually emerge and here I highlight just two that struck me as important for the future of computational chemistry.
1. Dispersion. This goes back to Fritz London and his formula: E_disp = -(C6/R⁶), where where coefficient C6 depends on the expectation values of the instantaneous dipole moments and average atomic excitation energies. The nature of this formula suggests that it decays rapidly with the distance between any pair of nuclei, R. But an increasing body of evidence is suggesting that such simple approaches (implemented as a correction in many e.g. DFT methods and known as e.g. D3+BJ, or D4 etc) may be underestimating the long range dispersion attractions. One nice example is what is known as the exfoliation of layers of graphite, where the forces holding the layers together can be measured quite accurately and which emerge as a great deal greater than the simple formulae suggest. It appears we now have a renaissance in developing new more accurate dispersion energy methods which include various higher order terms and are being applied to a variety of discrete molecule and solid state systems. One space to look out for!
2. .Non Born-Oppenheimer behaviour. It is a mainstay of most solutions of the Schroedinger equation where the nuclei are treated as classical point charge objects with fixed positions in an electronic field described by a wavefunction. But there is now considerable activity in developing methods that generate an extended Hessian (2nd derivative matrix) describing the forces that depends on both the classical nuclear coordinates of non-hydrogen atoms and the expectation values of quantum proton coordinates. This matrix is diagonalised to obtain the coupled vibrational frequencies which now naturally include the anharmonicity of the now quantum-treated protons and recovers the electron-proton correlation. It impacts most directly on so-called proton tunnelling and isotope effects, which can slice off 2-4 kcal/mol from barriers, but is now seen as a manifestation of electron-proton correlations in non-Born Oppenheimer potentials. The classical approach is to shave these energies off using eg Eckart potentials, but is now being replaced by e.g. a nuclear-electronic orbital method (NEO) which calculate the barriers from first principles. Typical types of reactions that are affected by non-BO behaviour are proton coupled electron transfers (PCET, see here for an example) which are increasingly seen as important in many biological processes.
I have tried to highlight just two themes that emerged from WATOC of personal interest to me; of course there was a great deal of new and exciting stuff that I have not mentioned. The next WATOC will be in Oslo in 2025, and no doubt new and exciting themes will emerge there as well!
July 11, 2022
Dioxane tetraketone – an ACS molecule of the week with a mystery.

I have long been fascinated by polymers of either carbon dioxide,^† or carbon monoxide, or combinations of both. One such molecule, referred to as dioxane tetraketone when it was featured on the ACS molecule-of-the-week site and also known as the anhydride of oxalic acid, or more formally 1,4-dioxane-2,3,5,6-tetraone, has been speculated upon for more than a century.[cite]10.1002/cber.19080410335[/cite]

The history of chemistry has many molecules whose existence has been speculated upon, but where attempted syntheses have failed and for which sound theoretical reasons often only emerged many years later.[cite]10.1023/A%3A1009270411806[/cite]

The synthesis of dioxane tetraketone was finally achieved in 1998[cite]10.1039/A803430C[/cite] at low temperatures (243K), although it was noted that in CDCl₃/Et₂O solutions at 273K it quickly decomposed to give equal quantities of carbon monoxide and dioxide. The characterisation was by ¹³C NMR, for which a single signal at 144.9 ppm was observed. The predicted value using the ACD/CNMR Predictor 2.0 Program (a so-called additive rule-based method) was 154 ppm (the value obtained using a similar tool available in Chemdraw is 150.9 ppm). The monomer oxirane-2,3-dione was also eliminated because of its predicted ¹³C shift using the same method of 167 ppm (155.3 using Chemdraw). Here I thought I would check these chemical shifts using a DFT-based method and also look at the barrier to the decomposition to see if it corresponds to a facile reaction at 273K (FAIR Data DOI: 10.14469/hpc/10619).

Firstly the NMR, using eg ωB97XD/aug-cc-pvdz/SCRF=chloroform. The calculated value of 148.2 ppm compares well with the observed value of 144.9 ppm. The value calculated for oxirane-2,3-dione was 156.6 ppm, rather lower than the ACD/Predictor method but in agreement with the Chemdraw implementation. The predicted IR spectrum (not reported) is shown below, should it ever be measured for this species.

Next, the reaction energy profile, this time calculated using ωB97XD/Def2-TZVPP for the reaction mechanism shown below.

The IRC reveals that the mechanism (black arrows) is followed, in a concerted process that reveals absolutely no sign of any ionic intermediate (red) which could then lead to oxirane-2,3-dione (blue). The barrier ΔG^‡ is 36.9 kcal/mol (it is lower than the total energy inferred below because the entropy is very positive, one molecule being converted to four during the reaction) which is far to high to correspond to a reaction that easily occurs at 273K. The value in water as solvent is very similar, again indicating that the ionic route is not enhanced by a polar solvent. The transition state has another feature of interest. It has C₂ chiral symmetry, typical of a pericyclic reaction with Möbius topology, as indeed would be appropriate for an eight electron process.

So what about that mystery then? Well, experimentally dioxane tetraketone decomposes at 273K, which would correspond to a free energy barrier of around 14-15 kcal/mol. The calculated value is far higher, too high to be simply an error in the DFT method. So here is a suggestion. CDCl₃, unless very carefully purified, contains HCl, which could very easily catalyse the reaction. So if another solvent were to be tried, lets say acetonitrile in which any trace of acid has been removed, would solutions of dioxane tetraketone then persist at room temperatures for far longer? An experiment perhaps to be tried!

^†Perhaps the most fascinating is the cyclic trimer of carbon dioxide, which arguably has pretensions to be aromatic. It has very recently been synthesized.[cite]10.1021/acs.joc.6b00647[/cite],[cite]10.1021/jp802872p[/cite]

June 22, 2022

Checking a conclusion we made in 1987: Tetrahedral intermediates formed by nitrogen and oxygen attack of aromatic hydroxylamines on acetyl cyanide

Minds (and memories) can work in wonderful ways. In 1987[cite]10.1021/jo00389a050[/cite] we were looking at the properties of “stable” tetrahedral intermediates formed in carbonyl group reactions. The reaction involved adding phenylhydroxylamine to acetyl cyanide. NMR signals for two new species were detected, and we surmised one was due to N-attack on the carbonyl and one was due to O-attack, in each case to form a stable tetrahedral intermediate. To try to identify which was which, ¹⁵N labelled hydroxylamine was used and then the ¹⁵N-¹³C coupling constants were measured, which could either be ^1-bondJ (for N-attack) or ^2-bondJ (for O-attack).

Well, 35 years later, literally in a dream on the morning of 7th June, 2022, these results came back to me and the dream involved wondering whether we had gotten the assignments of the N- and O-species the correct way around. You see we had assigned the larger of the ¹⁵N-¹³C couplings to the two bond (O-attack, species 3 below) rather than one-bond (N-attack, species 4 below) coupling. In 1987, the art of accurately computing such couplings was still in its infancy, but now in 2022 it is quick and easy to do. So here I report the results, which 35 years on allows a check of those assignments.

The necessary calculations are assembled at FAIR DOI: 10.14469/hpc/10593 conducted at the ωB97XD/aug-cc-pvdz/scrf=acetonitrile level. Firstly, it is important that the conformational space of these molecules is explored, since they contain a plethora of interesting anomeric effects. I will not discuss this process, simply quoting what I believe to be the lowest energy conformation for both isomers.

#	Property	Species 3	Species 4
1	ΔG₂₉₈	-608.600542	-608.598472
2	ΔG₂₁₅	-608.586956	-608.585163
3	NBO E(2)	14.3,19.4,10.9,8.1	10.0,11.2,9.9
4	δ_C obs	94.3 ppm	85.0 ppm
5	δ_C calc	97.2 (Δδ 2.9)	88.1 (Δδ 3.1)
6	J_N-C obs	²J ±2.5 Hz	¹J ±1.3 Hz
7	J_N-C calc	²J +1.7	¹J +0.8

The relative free energies ΔΔG₂₉₈ favour 3 over 4 by 1.3 kcal/mol at 298K (9:1). The article notes that 3 is significantly favoured over 4 at higher temperatures (i.e. ~298K) but that the concentration of 4 increases at lower temperatures.
At 215K, ΔΔG₂₁₅ reduces to 1.1 kcal/mol, but this equates to 13:1 at this temperature. ΔΔG₂₁₅ would need to be about 0.8 kcal/mol for 4 to increase (6.5:1), but these are small errors in energy and a more accurate calculation would have to be done to get this aspect correct.
The NBO E(2) terms indicating overlap between a lone pair and an acceptor orbital (the anomeric effect), show a dazzling variety of interactions for such a small molecule. Species 3 shows four significant interactions, species 4 one less.
The chemical shifts measured for 3 and 4 –
– are matched by the calculation, the error being similar for both species.
The ¹⁵N-¹³C coupling constants –
– are again matched, with the ¹J coupling being about half the value of the ²J coupling for both obs and calculated values.

The nature of modern scientific research, and the funding available for it, means that old work is rarely re-investigated using more recent techniques. In this case, the reinvestigation does not require the molecules to be re-synthesized again, merely that a retrospective computational layer be applied. As a result of my dream of four days ago, this process has produced an interesting new layer which thankfully confirms the original conclusions.

June 11, 2022

3-Methyl-5-phenylpyrazole: a crystallographic enigma?

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.[cite]10.1039/p29750001068[/cite] 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).

stationary points	B3LYP+GD3BJ/gas	B3LYP+GD3BJ/DCM	B3LYP/gas
SVPP, saddle=0, neutral	-1984.448365(0.0)	-1984.460070(0.0)	-1984.241838(0.0)
SVPP, saddle=0, ion-pair	-1984.426827(13.5)	-1984.439305(13.0)	-1984.218379(14.7)
TZVPP, saddle=0, neutral	-1986.712743(0.0)	-1986.725277 (0.0)	-1986.509333(0.0)
TZVPP, saddle=0, ion-pair	-1986.688467(15.2)	-1986.703562 (13.6)	-1986.481384(17.5)

SVPP, saddle=1	-1984.428898(12.1)	-1984.442004(11.3)	-1984.220238(13.6)
SVPP, saddle=2	-1984.430085(11.5)	-1984.442986(10.7)	-1984.220968(13.1)
SVPP, saddle=3	-1984.431457(10.6) 1315, 896, 801	-1984.441075(11.9) 970, 603, 294	-1984.223176 (11.7) 1335, 782, 760
TZVPP, saddle=1	-1986.691518(13.2)	–
TZVPP, saddle=2	-1986.689979(14.3)	–
TZVPP, saddle=3	-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

May 19, 2022

Geometries of proton transfers: modelled using total energy or free energy?

Proton transfers are amongst the most common of all chemical reactions. They are often thought of as “trivial” and even may not feature in many mechanistic schemes, other than perhaps the notation “PT”. The types with the lowest energy barriers for transfer often involve heteroatoms such as N and O, and the conventional transition state might be supposed to be when the proton is located at about the half way distance between the two heteroatoms. This should be the energy high point between the two positions for the proton. But what if a crystal structure is determined with the proton in exactly this position? Well, the first hypothesis is that using X-rays as the diffracting radiation is unreliable, because protons scatter x-rays very poorly. Then a more arduous neutron diffraction study is sometimes undertaken, which is generally assumed to be more reliable in determining the position of the proton. Just such a study was undertaken for the structure shown below (RAKQOJ)[cite]10/c3zxh2[/cite], dataDOI: 10.5517/cc57db3 for the 80K determination. The substituents had been selected to try to maximise the symmetry of the O…H…N motif via pKa tuning (for another tuning attempt, see this blog). The more general landscape this molecule fits into[cite]10.1039/C1RA00219H[/cite] is shown below:

The results obtained for the position of the proton for RAKQOJ were fascinating. They were very dependent on the temperature of the crystal! At room temperatures (using X-rays), the proton was measured as 1.09Å from the oxygen and 1.47Å from the nitrogen (neutral form above). At 20K, the OH distance was 1.309Å and the HN 1.206Å (~ionic form above). Indeed, the very title of this article is First O-H-N Hydrogen Bond with a Centered Proton Obtained by Thermally Induced Proton Migration. The authors give a number of reasons for this behaviour (their ref 17[cite]10/c3zxh2[/cite] and also[cite]10.1039/C1RA00219H[/cite]), but one they do not mention is thermally induced changes in the dielectric constant of the crystal with temperature, given that in one position for the proton the molecule is ionic and in the other neutral. So I decided to model the system as a function of solvent. In this model, the solvent dielectric is used to approximate the crystal dielectric. My first choice of energy function is to compute geometries using the B3LYP+GD3BJ/Def2=TZVPP/SCRF=solvent method to see what might emerge and as a possible prelude to trying other functionals. FAIR data for these calculations are collected at DOI: 10.14469/hpc/10368.

Solvent	ε	ΔG₂₉₈ for O…HN	r_O…H	r_HN	ΔG₂₉₈ for OH…N	r_OH	r_H…N	ΔG₂₉₈ TS (PT)	r_OH	r_HN
Water	78.4	-2893.387188 -2893.334325^♠	1.4913	1.0827	-2893.386705 -2893.334333^♠	1.0364	1.5696	-2893.387668 -2893.336183^♠	1.1852	1.2899
Dichloro methane	8.9	-2893.385173	1.4566	1.0945	-2893.385662	1.0309	1.5878	-2893.386022	1.2072	1.2642
Chloroform	4.7	-2893.382254	1.4227	1.1082	-2893.384514	1.0261	1.6049	-2893.384773	1.2321	1.2388
Dibutyl ether	3.1	-2893.380813	1.3778	1.1302	-2893.383511	1.0213	1.6235	-2893.382918	1.2667	1.2078
Toluene	2.4	-2893.379752	1.3248	1.1635	-2893.382915	1.0178	1.6385	-2893.379773	1.2851	1.1934
Gas phase	0	n/a			-2893.377949	1.0009	1.7387	n/a
Expt (RT) [cite]10/c3zxh2[/cite]	?	n/a				1.09	1.47	n/a
Expt (20K) [cite]10/c3zxh2[/cite]	?	n/a	1.309	1.206	n/a

^♠ At 20K

Results:

The geometries for each model are obtained by minimising the total energy of the system as a function of the 3N-6 geometric variables (coordinates).
The geometries show that for all solvents, TWO minima in the total energy are obtained, one for the ionic and one for the neutral form. This is called a double-well energy potential. Even a non-polar solvent such as toluene produces a solvation energy of ~3.1 kcal/mol compared to the gas phase, which is sufficient to induce a double-well potential.
Without solvent (gas phase), only the neutral geometry is obtained.
In the most polar solvent water, the double well potential looks like this:

The ionic well is about 0.4 kcal/mol lower in total energy (and ~0.3 kcal/mol in free energy, see table above) than the neutral form, with a barrier connecting neutral to ionic only 1.0 kcal/mol. A transition state + intrinsic reaction coordinate (IRC) can be easily located on this total energy potential, confirming the double-well form.
When free energies ΔG are computed, which include thermal effects such as entropy and zero-point energy, the transition state emerges as 0.3 kcal/mol less than the total energy of the ionic form (red entries, Table). In effect, the free energy potential surface is INVERTED compared to the total energy surface and the “transition state” becomes the lowest point on the energy surface. So this point is a minimum in the free energy but a maximum in the total energy, the result of adding thermal effects to the total energy.
In dichloromethane, the free energy of the neutral form is now lower by 0.3 kcal/mol than the ionic form. The OH bond is starting to get shorter and the NH one longer. The transition state is now 0.22 kcal/mol lower than the neutral form. With chloroform, the OH and HN bonds have become ~equal in length, the proton is symmetrically disposed.
By the time dibutyl ether as solvent is reached, the transition state is no longer lower in ΔG than the neutral form, moving on to being 2.0 kcal/mol higher for toluene. So as the solvent polarity decreases, we see a change in the potential from a single well in ΔG, in which the proton is centred, to a very asymmetric well in which the proton is attached to the oxygen.
Can we match the observed neutron diffraction results to the calculations? As the temperature decreases, the neutron diffraction shows the start of proton transfer from oxygen to nitrogen to form an ionic species. The calculations show that this can be modelled by an increase in the effective dielectric constant of the medium. The computed “transition state” for proton transfer somewhere between dibutyl ether and toluene (as a dielectric media) emerges as approximately the best model for the structure of this species. At this dielectric, the calculated ΔG is no longer quite the lowest free energy point in the potential. This might be due to the many approximations used in this model such as minimisation of total energy, the partition function method used to calculate entropy, the nature of the DFT functional, the continuum solvation model, the basis set, etc.

Conclusions:

These results were obtained with the approximation that minimising the total molecular energy produces a computed geometry that can be compared to the experimental neutron diffraction structures. But can one do better? Obtaining molecular geometries by minimising the computed free energies would be non-trivial. Firstly, minimisation would depend on availability of first derivatives of the energy function with respect to coordinates, in this case ΔG. These are not available for any DFT codes. The result would itself be temperature dependent (as indeed are the experimental results shown above). Furthermore, ΔG is computed from normal vibrational modes and these are only appropriate when the first derivatives of the function are zero, at which point the so-called six rotations and translations of the molecule in free space also have zero energy. So we need vibrations to compute derivatives, but we need derivatives to compute vibrations in this classical approach.

It would be great for example if the approximate model of the potential for a hydrogen transfer used above as based on minimising total energies for derivatives could be checked against a model based on geometries optimised using free energies instead. Such procedures do exist,[cite]10.1063/1.2715941[/cite] using molecular dynamics trajectory methods.

This post has DOI: 10.14469/hpc/10382 [cite]10/hqsm[/cite]

April 18, 2022

C2N2: a 10-electron four-atom molecule displaying both Hückel 4n+2 and Baird 4n selection rules for ring aromaticity.

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.[cite]10.1021/ja00769a025[/cite],[cite]10.1021/cr300471v[/cite]. 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 C₂N₂.
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.[cite]10.1021/ja00769a025[/cite],[cite]10.1021/cr300471v[/cite] 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

April 7, 2022

Bonding MOs for C₂N₂. Click image to load 3D model
π3, 1 electron	π2, 1 electron

σ3 2 electron	σ2, 2 electron

π1 2 electron	σ1, 2 electron

Raw data and the evolution of crystallographic FAIR data. Journals, processed and raw structure data.

In my previous post on the topic, I introduced the concept that data can come in several forms, most commonly as “raw” or primary data and as a “processed” version of this data that has added value. In crystallography, the chemist is interested in this processed version, carried by a CIF file. However on rare occasions when a query arises about the processed component, this can in principle at least be resolved by taking a look at the original raw data, expressed as diffraction images. I established with much appreciated help from CCDC that since 2016, around 65 datasets in the CSD (Cambridge structural database) have appeared with such associated raw data. The problem is easily reconciling the two sets of data (the raw data is not stored on CSD) and one way of doing this is via the metadata associated with the datasets. In turn, if this metadata is suitably registered, one can query the metadata store for such associations, as was illustrated in the previous post on the topic. Here I explore the metadata records for five of these 65 sets to find out their properties, selected to illustrate the five data repositories thus far that host such data for compounds in the CSD database.

Raw data repository	Raw Data DOI	Raw data →CSD?	CSD→ Raw data?	⇐Journal⇒
Zenodo	10.5281/zenodo.4271549	No	No	10.1039/C6RA28567H
Imperial College research data repository	10.14469/hpc/2298	Yes	Yes	10.1021/acsomega.7b00482
RepoD, a Harvard Dataverse instance	10.18150/repod.6628285	No	No	10.1021/acs.cgd.0c01252
Cambridge university repository	10.17863/CAM.21968	No	No	10.1016/j.inoche.2018.08.024
Isis neutron and muon source data journal	10.5286/ISIS.E.RB1620465	No	No	10.1039/D0CC02418J

Ideally, one is looking for bidirectional links between the data as expressed in the metadata and in both directions. As you can see from the above, these links are present in only one of the five sets. More common is that both the raw and the processed data will contain links to the journal article where the data is discussed. Very much less commonly are there links from the journal article to the raw data, although such links are slightly more likely to exist from the journal to the processed data. If you click on the link in any of the last three columns, a copy of the metadata will download for you to inspect. There you can verify if the assertions made above are correct.

What the metadata records demonstrate above is a very small scale so-called PID graph (DOI: [cite]10.5438/jwvf-8a66[/cite] 10.5438/jwvf-8a66) where each DOI is a node in that graph and if a connection exists, it is shown by a line connecting the nodes. The PID graph can be extended to include a third type of node, the journal article and then it starts to get interesting! I will investigate if I can generate the PID graph for the above, although be prepared, it will not (yet) contain very many lines between nodes!

March 28, 2022

Sir Geoffrey Wilkinson: An anniversary celebration. 23 March, 2022, Burlington House, London.
The meeting covered the scientific life of Professor Sir Geoffrey Wilkinson from the perspective of collaborators, friends and family and celebrated three anniversaries, the centenary of his birth (2021), the half-century anniversary of the Nobel prize (2023) and 70 years almost to the day (1 April) since the publication of the seminal article on Ferrocene (2022).[cite]10.1021/ja01128a527[/cite]

The meeting was organised as “inverse hybrid” (to use the new terminology), with a maximum capacity in-person audience attending along with fourteen speakers, three of whom were remote and one who could not attend on the day but whose presentation was given on their behalf. I will not give abstracts for the talks here, but note two common themes that I thought emerged during the day.
1. All the speakers found themes in either their memories of Wilkinson and their time in his laboratories or their current research work that show how he continues to influence, along with the famous text book that he co-wrote, the modern world of chemistry. He truly left a remarkable legacy.
2. This is a personal observation, but in his day, Wilkinson was famously sceptical of the ability of molecular modelling to cast profound insights into the molecules his group were studying. Yesterday I think with only one or two exceptions, the talks were accompanied by “DFT modelling” helping to provide such insights, either into the reaction mechanisms via energy profiles or into the properties of the molecules themselves, including their spectroscopy.
A small exhibition of artefacts included his famous portrait, all the editions of the text book and other items from his desk.

Finally, I thought I might explore the famous controversy surrounding the model of ferrocene which is shown in the photos below. It is shown with the two cyclopentadienyl rings in a so-called “eclipsed” conformation. To cast light on this, I show a search of the Cambridge crystal database of all molecules with this sub-structure. There are 24,868 of them.

The histogram plot of the dihedral angle is shown below. The staggered geometry has a dihedral of 36° and you can see a small maximum at this point in the distribution below. But this is dwarfed by 0°, the value for the eclipsed orientation. The barrier to rotation is known to be very small, and this is reflected in the almost continuous distribution amongst those 24,868 molecules.
March 24, 2022
A four-atom molecule exhibiting simultaneous compliance with Hückel 4n+2 and Baird 4n selection rules for ring aromaticity.

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 C₂B₂.
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.

March 22, 2022
More aromatic species with four atoms. B4 and N4.

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 B₄.
Click image to load 3D model

σ1, -0.335 au

π1, -0.372 au

Next, N₄.

π-Bonding MOs for N₄.
Click image to load 3D model

π3 π2

π1

σ-Bonding MOs for N₄.
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!

March 19, 2022

π-Bonding MOs for N₄. Click image to load 3D model
π3	π2

π1

σ-Bonding MOs for N₄. Click image to load 3D model
σ3	σ2

σ1