Importance Of Vibrational Zero-Point Energy Contribution To The Relative Polymorph Energies Of Hydrogen-Bonded Species

In press, in Crystal Growth and Design. A short paper with an important message.  In principle, solid-state quantum chemical methods provide the tools for both the prediction of crystal polymorphs and the calculation of the relative energies of characterized crystal forms to determine why one version forms preferentially. That said, modern solid-state quantum chemical methods are dominated by density functional theory which, although work is being performed to address this issue, are fundamentally incapable of accounting for all of the energy terms in the complete lattice energy calculation of solid-state materials because dispersion forces are not accounted for (there are methods around this, be they empirical or by way of solid-state Moeller-Plesset Perturbation Theory, which we may even have the computing power to handle someday…).

In molecular polymorphs, the energies of the crystal cells per molecule may be quite similar to one another because, being molecules with polar and non-polar regions, specific functional groups tend to bind preferentially to other specific functional groups, which all may involve similar interaction energies.  The point is that the lattice energy differences between different polymorphs may be quite small.  In such cases, the vibrational zero-point energy of the crystal cells may be very important contributors to the experimentally determined polymorph energy differences.  Such is found to be the case for the alpha- and gamma- polymorphs of glycine.

Specifically in the glycine example and other polymorphs for which proper thermochemistry is available, we even have a means to estimating what SHOULD be the lattice energy of the crystal without relying on theoretical models (and their inherent limitations).  For glycine, the experimental enthalpies for both crystal forms have been measured.  We have the experimental vibrational spectrum available against which to compare the theoretical work, from which we can determine the zero-point energy for the unit cell (simply 1/2 the sum of the vibrational mode energies).  With that information, we can determine (to a first approximation, there are many other factors to consider that play lesser roles in the final value) the experimental lattice energies.  These values then provide a benchmark for determining the abilities of theoretical models to reproduce this most fundamental of the solid-state quantum chemical properties.

Sharon A. Rivera1,2, Damian G. Allis1,3, Bruce S. Hudson1

1. Department of Chemistry, Syracuse University, Syracuse, NY 13224-4100
2. School of Chemistry, University of New South Wales, Sydney, NSW 2052, Australia
3. Nanorex, Inc., Bloomfield Hills, MI 48302-4100

Abstract: The relative stability of polymorphic crystal forms is a challenging conceptual problem of considerable technical interest.  Current estimates of relative polymorph energies concentrate on lattice energy.  In this work the contribution of differences in zero-point energy and vibrational enthalpy to the enthalpy difference for polymorphs is investigated.  The specific case investigated is that of alpha- and gamma-glycine, for which the experimental enthalpy difference is known.  Periodic lattice DFT computations are used to provide the vibrational spectrum at the Gamma-point.  It is confirmed that these methods provide reasonable descriptions of the inelastic neutron scattering spectra of these two crystals.  It is found that the difference in the zero-point energy is about 1.9 kJ/mol and that the vibrational thermal population difference is 0.9 kJ/mol in the opposite sense.  The overall vibrational contributions to the enthalpy difference are much larger than the observed value of ca. 0.3 kJ/mol.  The vibrational contribution must be largely compensated by the lattice energy difference.  The polymorphs of glycine differ in the pattern of their hydrogen bonds, a feature common to many polymorphs of interest.  The consequent difference in the N-H stretching frequencies is a contributor to the zero-point correction but the major effect stems from changes in the bending vibrations.

The Solid-State Terahertz Spectrum of MDMA (Ecstasy) – A Unique Test for Molecular Modeling Assignments

In press, in Chemical Physics Letters (CPL).  Yes, the blog has taken a bit of a turn, from high explosives to illicit drugs.  I expect my google rating to rise sharply with this post.  The protonated form of the 3,4-methylene-dioxymethamphetamine (Ecstasy, but I'll keep the post legit, so it is herein referred to as MDMA) molecule (herein referred to as MDMA:H+) and MDMA:H+ in its crystal cell with a chloride ion (Cl, the crystal form herein referred to as MDMA:HCl) is shown below in yet another fantastic NanoEngineer-1 rendering (if I do say so myself).

This CPL article is, to some extent, a response to those in the terahertz community who continue to attempt spectral assignments of crystalline and poly-crystalline samples using isolated-molecule quantum chemical calculations.  The MDMA:HCl sample and MDMA molecule, as a pair, are a very interesting case study of theory and experiment for reasons detailed below.  The spectra, shown below from a previous version of the paper (but the same spectra), show quite a bit of detail that will make sense shortly.

Panel A shows the isolated-molecule calculation for the neutral MDMA molecule at a B3LYP/6-31G(d) level of theory (in red).  You will note that this simulated spectrum is in very good agreement with experiment (in black), reproducing all of the major features and showing a number of smaller features that account for shoulders.  This agreement was the basis for the assignment of the MDMA:HCl spectrum reported in: G Wang, J Shen, Y Jia. "Vibrational spectra of ketamine hydrochloride and 3,4-methylenedioxymethamphetamine in terahertz range." Journal of Applied Physics 102 (2007) 013106/1-06/4.

The new theoretical analysis reported in the CPL article was instigated by this assignment in this previous publication.  Relevant to the measured sample and the previously reported assignment, two points arise that require address.

1. The previous calculation, as reported, was of the neutral MDMA molecule and is reasonably close to the MDMA spectrum shown above (in red.  This calculation was redone for the CPL article for comparative purposes).  As the experimental THz sample was of solid-state MDMA:HCl, the appropriate form of the molecule to run is not the neutral MDMA molecule, but the protonated form, MDMA:H+.  The protonated form has a different vibrational spectrum (shown in green in Panel A) than the isolated molecule form.  At the very least, the isolated-molecule to consider for the MDMA:HCl sample must be the protonated form.  Interestingly, the re-calculation at B3LYP/6-31G(d) reported for the CPL article predicts a fifth vibrational mode at 48.0 cm-1 that was not reported in the previous study.  We do not know if the previous group missed that peak in the write-up, decided that (since it is a low-intensity mode) it was not worth reporting, or if their starting molecular geometry was somehow different so that the other four modes were predicted to be in the same region and this mode was somehow turned off.

2. The solid-state spectrum shown in Panel B at a BP/DNP level of theory does not agree as well as the isolated-molecule MDMA B3LYP/6-31G(d) calculation. That being said, THAT IS NOT THE POINT.  The goal of a simulated spectrum IS NOT to obtain the closest spectral agreement with experiment.  The goal IS to explain the solid-state spectrum with the best theoretical model possible that, hopefully, is as close to the experimental result as possible.  In this case, the solid-state BP/DNP spectrum contains a finite number of vibrational modes that do group according to features in the THz spectrum, making the assignment reasonably straightforward.  Interestingly, the two most intense modes in the solid-state BP/DNP calculations involve the motions of the Cl…H+-N chains, which CANNOT be accounted for in an isolated-molecule calculation of either the neutral MDMA molecule or the protonated MDMA:H+.

In summary, as taken from the CPL paper:

With all of these considerations taken into account in this re-examination of the MDMA.HCl THz spectrum, it is found that this system serves as a fortuitous example of one whose THz spectrum is predicted quite precisely by two very different approaches, but is only described accurately by one that considers the crystal environment and the actual state of the molecule in its solid-state form.

Damian G. Allis1,2, Patrick M. Hakey1, and Timothy M. Korter1

1. Department of Chemistry, Syracuse University, Syracuse NY 13244-4100 USA
2. Nanorex, Inc., P.O. Box 7188, Bloomfield Hills, MI 48302-7188 USA

Abstract: The terahertz (THz, far-infrared) spectrum of 3,4-methylene-dioxymethamphetamine hydrochloride (Ecstasy) is simulated using solid-state density functional theory.  While a previously reported isolated-molecule calculation is noteworthy for the precision of its solid-state THz reproduction, the solid-state calculation predicts that the isolated-molecule modes account for only half of the spectral features in the THz region, with the remaining structure arising from lattice vibrations that cannot be predicted without solid-state molecular modeling.  The molecular origins of the internal mode contributions to the solid-state THz spectrum, as well as the proper consideration of the protonation state of the molecule, are also considered.

Vicinal Deuterium Perturbations on Hydrogen NMR Chemical Shifts in Cyclohexanes

In press, in the Journal of the American Chemical Society. To answer a student question (on the first day of class, no less) about the geometry/sp3 hybridization of carbon in my 2nd semester organic chemistry class, Prof. Baldwin began by writing the time-independent Shroedinger Equation on the chalkboard.

Needless to say, the wavefunction that was my future collapsed instantaneously.

To summarize considerable computational effort (and ignoring many, many equations you can find in the paper), the shifts in hydrogen (H) peak positions in the 1H NMR spectrum of cyclohexane that result from deuteration (D) can be explained by considering the contribution of the vibrational stretching mode differences of the C-H and C-D bonds. Interestingly, the trend determined from the computational studies is fairly insensitive to the use of electron correlation methods (B3LYP, MPW91, MP2 were all employed with 6-31++G(2d,p), 6-311++G(2d,2p), and Aug-cc-pVTZ basis sets), with Restricted Hartree-Fock (RHF) results very similar and much faster to calculate. The RHF results served as the basis for the analysis in the paper.

The geometry optimization of a molecule will provide you static bond lengths (C-H, C-D, C-C for cyclohexane). In the case of just the optimizations, H and D look identical to the electronic wavefunction (the basis of the Born-Oppenheimer Approximation). If the bond lengths stretched like harmonic oscillators (the approximation we almost always use when performing normal mode analyses in quantum chemical codes), the calculated bond lengths would be the correct (physical) average bond lengths. Of course, no chemical bond is a perfect harmonic oscillator, a fact that explains transition probabilities for overtone and combination bands that are forbidden within the harmonic approximation. For atom pairs that stretch with a significant anharmonic character (such as C-H/D), the actual description of the bond length is a very important factor to consider when trying to extract important trends based on the actual bond length behavior, and is as described in the figure below.

This figure shows the bottom of the potential wells for idealized harmonic and Morse oscillators. The figure, like the paper, only considers the contribution from the vibrational zero-point level (v0). The Heq/Deq position at the bottom of the potential wells is the calculated bond length after geometry optimization. If the stretching potential is harmonic (red), the average position of the H/D atom is always at this position regardless of the system energy (the H atom is sampling as much of the short-side of the potential as it is the long-side, which, you guessed it, averages to the middle).

Now consider the anharmonic oscillators (our Morse functions). The H atom, being lighter, has a higher vibrational zero-point energy (VZPE). Its average position (or mean displacement from the C atom) due to VZPE (Heq + v0) is slightly longer. The deuterium, with twice the mass, has a lower VZPE. As a result of lying lower along the Morse potential, its average VZPE-corrected position (Deq + v0) is closer to the calculated static position than for H. The difference in the contributed length to the average positions is shown in green. This position-based description is really just for show and neglects all of the quantum theory, but it does make the point. You can consider the green to be the effective additional distance within the Morse potential that the two atoms sample compared to the harmonic potential as part of their stretching mode, with the average position differing due to this larger spatial sampling in each oscillation.

So, "real C-D" is shorter than "real C-H". Accordingly, the 1H NMR peak position of C-H in the presence of C-H bonds (which are further away) is different than in the presence of C-D bonds (which are closer).

This is also the first paper I've been on that did not directly involve my making an image (hence the involved image above).

Daniel J. O'Leary1, Damian G. Allis2,3, Bruce S. Hudson2, Shelly James2, Katherine B. Morgera2, and John E. Baldwin2

1. Department of Chemistry, Pomona College, Claremont, California 9171-6338
2. Department of Chemistry, Syracuse University, Syracuse, New York 13244-4100
3. Nanorex, Inc., P.O. Box 7188, Bloomfield Hills, Michigan 48302-7188

Abstract: The substitution of a deuterium for a hydrogen is known to perturb the NMR chemical shift of a neighboring hydrogen atom. The magnitude of such a perturbation may depend on the specifics of bonding and stereochemical relationships within a molecule. For deuterium-labeled cyclohexanes held in a chair conformation at -80oC or lower, all four possible perturbations of H by D as H-C-C-H is changed to D-C-C-H have been determined experimentally, and the variations seen, ranging from 6.9 to 10.4 ppb, have been calculated from theory and computational methods. The predominant physical origins of the NMR chemical shift perturbations in deuterium-labeled cyclohexanes have been identified and quantified. The trends defined by the Δδ perturbation values obtained through spectroscopic experiments and by theory agree satisfactorily. They do not match the variations typically observed in vicinal JH-H coupling constants as a function of dihedral angles.